Search

This paper shows both analytically and numerically that the skin effect in systems of non-Hermitian subwavelength resonators is robust with respect to random imperfections in the system. The subwavelength resonators are highly contrasting material inclusions that resonate in a low-frequency regime. The non-Hermiticity is due to the introduction of a directional damping term (motivated by an imaginary gauge potential), which leads to a skin effect that is manifested by the system's eigenmodes accumulating at one edge of the structure. We elucidate the topological protection of the associated (real) eigenfrequencies and illustrate numerically the competition between the two different localization effects present when the system is randomly perturbed: the non-Hermitian skin effect and the disorder-induced Anderson localization. We show numerically that, as the strength of the disorder increases, more and more eigenmodes become localized in the bulk. Our results are based on an asymptotic matrix model for subwavelength physics and can be generalized also to tight-binding models in condensed matter theory. [ABSTRACT FROM AUTHOR]

The present work is devoted to the eigenvalue asymptotic expansion of the Toeplitz matrix Tn(a), whose generating function a is complex-valued and has a power singularity at one point. As a consequence, Tn(a) is non-Hermitian and we know that in this setting, the eigenvalue computation is a nontrivial task for large sizes. First we follow the work of Bogoya, Böttcher, Grudsky, and Maximenko and deduce a complete asymptotic expansion for the eigenvalues. In a second step, we apply matrixless algorithms, in the spirit of the work by Ekström, Furci, Garoni, Serra-Capizzano et al., for computing those eigenvalues. Since the inner and extreme eigenvalues have different asymptotic behaviors, we worked on them independently and combined the results to produce a high precision global numerical and matrixless algorithm. The numerical results are very precise, and the computational cost of the proposed algorithms is independent of the size of the considered matrices for each eigenvalue, which implies a linear cost when the entire spectrum is computed. From the viewpoint of real-world applications, we emphasize that the class under consideration includes the matrices stemming from the numerical approximation of fractional diffusion equations. In the final section a concise discussion on the matter and a few open problems are presented. [ABSTRACT FROM AUTHOR]

In this work, we propose a novel preconditioned Krylov subspace method for solving an optimal control problem of wave equations, after explicitly identifying the asymptotic spectral distribution of the involved sequence of linear coefficient matrices from the optimal control problem. Namely, we first show that the all-at-once system stemming from the wave control problem is associated to a structured coefficient matrix-sequence possessing an eigenvalue distribution. Then, based on such a spectral distribution of which the symbol is explicitly identified, we develop an ideal preconditioner and two parallel-in-time preconditioners for the saddle point system composed of two block Toeplitz matrices. For the ideal preconditioner, we show that the eigenvalues of the preconditioned matrix-sequence all belong to the set (-3/2, -1/2) ∪ (1/2,3/2) well separated from zero, leading to mesh-independent convergence when the minimal residual method is employed. The proposed parallel-in-time preconditioners can be implemented efficiently using fast Fourier transforms or discrete sine transforms, and their effectiveness is theoretically shown in the sense that the eigenvalues of the preconditioned matrix-sequences are clustered around ±1, which leads to rapid convergence. When these parallel-in-time preconditioners are not fastly diagonalizable, we further propose modified versions which can be efficiently inverted. Several numerical examples are reported to verify our derived localization and spectral distribution result and to support the effectiveness of our proposed preconditioners and the related advantages with respect to the relevant literature. [ABSTRACT FROM AUTHOR]

Solving the linear system (KK+)-1b is often the major computational burden when a forward-backward evolutionary equation must be solved in a problem, where K is the so-called all-at-once matrix of the forward subproblem after space-time discretization. An efficient solver requires a good preconditioner for KKT. Inspired by the structure of K, we precondition KKT by PαPαT with Pα being a block α-circulant matrix constructed by replacing the Toeplitz matrices in K by the α-circulant matrices. By a block Fourier diagonalization of Pα, the computation of the preconditioning step (PαPαT)-1r is parallelizable for all the time steps. We give a spectral analysis for the preconditioned matrix (PαPαT)-1(KKT) and prove that for any one-step stable time-integrator the eigenvalues of (PαPαT)-1 (KKT) spread in a mesh-independent interval [(1 + √2δ-1, (1-√2δ)-1] if the parameter α weakly scales in terms of the number of time steps Nt as α = δ/√Nt, where δ ∈ (0, 1/√2) is a free constant. Two applications of the proposed preconditioner are illustrated: PDE-constrained optimal control problems and parabolic source identification problems. Numerical results for both problems indicate that spectral analysis predicts the convergence rate of the preconditioned conjugate gradient method very well. [ABSTRACT FROM AUTHOR]

We propose recursive algorithms to update an orthogonal numerical basis matrix of the null space of the block row, (banded) block Toeplitz, and block Macaulay matrix, which is the multivariate generalization of the (banded) block Toeplitz matrix. These structured matrices are often constructed in an iterative way, and, for some applications, a basis matrix of the null space is required in every iteration. Consequently, recursively updating a numerical basis matrix of the null space, while exploiting the inherent structure of the matrices involved, induces large savings in the computation time. Moreover, we also develop a sparse adaptation of one of the recursive algorithms that avoids the explicit construction of the block Macaulay matrix and results in a considerable reduction of the required memory. We provide several numerical experiments to illustrate the proposed algorithms: for example, we solve four multiparameter eigenvalue problems via the null space of the block Macaulay matrix and notice that the recursive and sparse approach are, on average, 450 and 1300 times faster than the standard approach, respectively. [ABSTRACT FROM AUTHOR]

Log-concavity and almost log-convexity of the cycle index polynomials were proved by Bender and Canfield [J. Combin. Theory Ser. A, 74 (1996), pp. 57--70]. Schirmacher [J. Combin. Theory Ser. A, 85 (1999), pp. 127--134] extended them to q-log-concavity and almost q-log-convexity. Motivated by these, we consider the stronger properties total positivity from the Toeplitz matrix and Hankel matrix. By using exponential Riordan array methods, we give some criteria for total positivity of the triangular matrix of coefficients of the generalized cycle index polynomials, the Toeplitz matrix and Hankel matrix of the polynomial sequence in terms of the exponential formula, the logarithmic formula, and the fractional formula, respectively. Finally, we apply our criteria to some triangular arrays satisfying some recurrence relations, including Bessel triangles of two kinds and their generalizations, the Lah triangle and its generalization, the idempotent triangle, and some triangles related to binomial coefficients, rook polynomials, and Laguerre polynomials. We not only get total positivity of these lower-triangles, and q-Stieltjes moment properties and 3-q-log-convexity of their row-generating functions, but also prove that their triangular convolutions preserve the Stieltjes moment property. In particular, we solve a conjecture of Sokal on the q-Stieltjes moment property of rook polynomials. [ABSTRACT FROM AUTHOR]

We study pseudospectra of matrices arising in models of ion channel dynamics, bimolecular reactions, and Michaelis--Menten enzyme kinetics. We show that matrices for ion channel models and for bimolecular reactions are examples of a special class known as twisted Toeplitz matrices. Using Michaelis--Menten enzyme kinetics as an example, we suggest more complicated reactions cannot be modeled simply by the twisted Toeplitz ideas currently in the literature, and thus we propose a generalization to what we call block twisted Toeplitz matrices. Furthermore, we show that certain simple cases of ion channel models always exhibit a phenomenon known as localization. This result comes by studying the symbol of the operator and an application of a theorem of Trefethen and Chapman. We further observe that a larger magnitude of the so-called twist-ratio, a number that comes from the symbol, is associated with stronger localization. For certain bimolecular reactions, specialized to a dimerization, we show by the same theorem that there are large regions of parameter space where localization is guaranteed and that this is always guaranteed in the limit as the matrix size becomes very large. However, we also find examples for these reactions where the hypotheses of the theorem of Trefethen and Chapman do not hold, and yet in numerical experiments, we still observe localization. Finally we return to the first of our matrix models for a more detailed study to give further insight into important applications to cardiac ion channels and in particular sodium channels. We demonstrate behavior similar to the so-called cutoff phenomenon of Markov processes, and we show that the stationary distribution, which is so crucial to applications of these models to ion channels and especially sodium channels, can be extremely localized. [ABSTRACT FROM AUTHOR]

In [McDonald, Pestana, and Wathen, SIAM J. Sci. Comput., 40 (2018), pp. A1012--A1033], a block circulant preconditioner is proposed for all-at-once linear systems arising from evolutionary partial differential equations, in which the preconditioned matrix is proven to be diagonalizable and to have identity-plus-low-rank decomposition in the case of the heat equation. In this paper, we generalize the block circulant preconditioner by introducing a small parameter ϵ>0 into the top-right block of the block circulant preconditioner. The implementation of the generalized preconditioner requires the same computational complexity as that of the block circulant one. Theoretically, we prove that (i) the generalization preserves the diagonalizability and the identity-plus-low-rank decomposition; (ii) all eigenvalues of the new preconditioned matrix are clustered at 1 for sufficiently small\epsilon; (iii) GMRES method for the preconditioned system has a linear convergence rate independent of size of the linear system when\epsilon is taken to be smaller than or comparable to square root of time-step size. Numerical results are reported to confirm the efficiency of the proposed preconditioner and to show that the generalization improves the performance of block circulant preconditioner. [ABSTRACT FROM AUTHOR]

In this work, we perform a spectral analysis of flipped multilevel Toeplitz sequences, i.e., we study the asymptotic spectral behavior of {YnTn (f)} n, where Tn (f) is a real, square multilevel Toeplitz matrix generated by a function f ∈ L¹([-π, π]d) and Yn is the exchange matrix, which has 1's on the main antidiagonal. In line with what we have shown for unilevel flipped Toeplitz matrix sequences, the asymptotic spectrum is determined by a 2 x 2 matrix-valued function whose eigenvalues are ±|f|. Furthermore, we characterize the eigenvalue distribution of certain preconditioned flipped multilevel Toeplitz sequences with an analysis that covers both multilevel Toeplitz and circulant preconditioners. Finally, all our findings are illustrated by several numerical experiments. [ABSTRACT FROM AUTHOR]

The aim of this paper is to propose a robust n umerical solver, which is capable of efficiently solving a three-dimensional elliptic problem in a data-sparse quantized tensor format. In particular, we use the combined Tucker and quantized tensor train format (TQTT), which allows us to use astronomically large grid sizes. However, due to the ill-conditioning of discretized differential operators, such fine grids lead to numerical instabilities. To obtain a robust solver, we utilize the well-known alternating direction implicit method and modify it to avoid multiplication by differential operators. So as to make the method efficient, we derive an explicit TQTT representation of the iteration matrix and quantized tensor train representations of the inverses of symmetric tridiagonal Toeplitz matrices as an auxiliary result. As an application, we consider accurate solution of elliptic problems with singular potentials arising in electronic Schrodinger's equation. [ABSTRACT FROM AUTHOR]

The problem of finding the distance from a given n×n matrix polynomial of degree k to the set of matrix polynomials having the elementary divisor (λ-λ0)j,j⩾r, for a fixed scalar λ0 and 2⩽r⩽kn is considered. It is established that polynomials that are not regular are arbitrarily close to a regular matrix polynomial with the desired elementary divisor. For regular matrix polynomials the problem is shown to be equivalent to finding minimal structure preserving perturbations such that a certain block Toeplitz matrix becomes suitably rank deficient. This is then used to characterize the distance via two different optimizations. The first one shows that if λ0 is not already an eigenvalue of the matrix polynomial, then the problem is equivalent to computing a generalized notion of a structured singular value. The distance is computed via algorithms like BFGS and Matlab's globalsearch algorithm from the second optimization. Upper and lower bounds of the distance are also derived and numerical experiments are performed to compare them with the computed values of the distance. [ABSTRACT FROM AUTHOR]

The treatment of two-dimensional random walks in the quarter plane leads to Markov processes which involve semi-infinite matrices having Toeplitz or block Toeplitz structure plus a lowrank correction. We propose an extension of the framework introduced in [D. A. Bini, S. Massei, and B. Meini, Math. Comp., 87 (2018), pp. 2811–2830] which allows us to deal with more general situations such as processes involving restart events. This is motivated by the need for modeling processes that can incur in unexpected failures like computer system reboots. We present a theoretical analysis of an enriched Banach algebra that, combined with appropriate algorithms, enables the numerical treatment of these problems. The results are applied to the solution of bidimensional quasi-birth-death processes with infinitely many phases which model random walks in the quarter plane, relying on the matrix analytic approach. The reliability of our approach is confirmed by extensive numerical experimentation on several case studies. [ABSTRACT FROM AUTHOR]

Quadratic matrix equations of the kind A1X2 +A0X +A 1 = X are encountered in the analysis of Quasi--Birth-Death stochastic processes where the solution of interest is the minimal nonnegative solution G. In many queueing models, described by random walks in the quarter plane, the coefficients A1,A0,A 1 are infinite tridiagonal matrices with an almost Toeplitz structure. Here, we analyze some fixed point iterations, including Newton's iteration, for the computation of G and introduce effective algorithms and acceleration strategies which fully exploit the Toeplitz structure of the matrix coefficients and of the current approximation. Moreover, we provide a structured perturbation analysis for the solution G. The results of some numerical experiments which demonstrate the effectiveness of our approach are reported. [ABSTRACT FROM AUTHOR]

We present new algorithms for randomized construction of hierarchically semiseparable (HSS) matrices, addressing several practical issues. The HSS construction algorithms use a partially matrix-free, adaptive randomized projection scheme to determine the maximum off-diagonal block rank. We develop both relative and absolute stopping criteria to determine the minimum dimension of the random projection matrix that is sufficient for desired accuracy. Two strategies are discussed to adaptively enlarge the random sample matrix: repeated doubling of the number of random vectors and iteratively incrementing the number of random vectors by a fixed number. The relative and absolute stopping criteria are based on probabilistic bounds for the Frobenius norm of the random projection of the Hankel blocks of the input matrix. We discuss parallel implementation and computation and communication cost of both variants. Parallel numerical results for a range of applications, including boundary element method matrices and quantum chemistry Toeplitz matrices, show the effectiveness, scalability, and numerical robustness of the proposed algorithms. [ABSTRACT FROM AUTHOR]

Given a Lebesgue integrable function f over [-π; π], we consider the sequence of matrices {YnTn[f]}n, where Tn[f] is the n-by-n Toeplitz matrix generated by f and Yn is the anti-identity matrix. Because of the unitary nature of Yn, the singular values of Tn[f] and YnTn[f] coincide. However, the eigenvalues are affected substantially by the action of Yn. Under the assumption that the Fourier coefficients of f are real, we prove that {YnTn[f]}n is distributed in the eigenvalue sense as ±|f|. A generalization of this result to the block Toeplitz case is also shown. We also consider the preconditioning introduced by [J. Pestana and A. Wathen, SIAM J. Matrix Anal. Appl., 36 (2015), pp. 273{288] and prove that the preconditioned matrix-sequence is distributed in the eigenvalue sense as π1 under the mild assumption that f is sparsely vanishing. We emphasize that the mathematical tools introduced in this setting have a general character and can be potentially used in different contexts. A number of numerical experiments are provided and critically discussed. [ABSTRACT FROM AUTHOR]

When solving linear systems with nonsymmetric Toeplitz or multilevel Toeplitz matrices using Krylov subspace methods, the coefficient matrix may be symmetrized. The preconditioned MINRES method can then be applied to this symmetrized system, which allows rigorous upper bounds on the number of MINRES iterations to be obtained. However, effective preconditioners for symmetrized (multilevel) Toeplitz matrices are lacking. Here, we propose novel ideal preconditioners and investigate the spectra of the preconditioned matrices. We show how these preconditioners can be approximated and demonstrate their effectiveness via numerical experiments. [ABSTRACT FROM AUTHOR]

In cryo-electron microscopy, the three-dimensional (3D) electric potentials of an ensemble of molecules are projected along arbitrary viewing directions to yield noisy two-dimensional images. The volume maps representing these potentials typically exhibit a great deal of structural variability, which is described by their 3D covariance matrix. Typically, this covariance matrix is approximately low rank and can be used to cluster the volumes or estimate the intrinsic geometry of the conformation space. We formulate the estimation of this covariance matrix as a linear inverse problem, yielding a consistent least-squares estimator. For n images of size N-by-N pixels, we propose an algorithm for calculating this covariance estimator with computational complexity O(nN4 + p √κN6 logN), where the condition number κ is empirically in the range 10{200. Its efficiency relies on the observation that the normal equations are equivalent to a deconvolution problem in six dimensions. This is then solved by the conjugate gradient method with an appropriate circulant preconditioner. The result is the first computationally efficient algorithm for consistent estimation of the 3D covariance from noisy projections. It also compares favorably in runtime with respect to previously proposed nonconsistent estimators. Motivated by the recent success of eigenvalue shrinkage procedures for high-dimensional covariance matrix estimation, we incorporate a shrinkage procedure that improves accuracy at lower signal-to-noise ratios. We evaluate our methods on simulated datasets and achieve classification results comparable to state-of-the-art methods in shorter running time. We also present results on clustering volumes in an experimental dataset, illustrating the power of the proposed algorithm for practical determination of structural variability. [ABSTRACT FROM AUTHOR]

Standard Krylov subspace solvers for self-adjoint problems have rigorous convergence bounds based solely on eigenvalues. However, for non-self-adjoint problems, eigenvalues do not determine behavior even for widely used iterative methods. In this paper, we discuss time-dependent PDE problems, which are always non-self-adjoint. We propose a block circulant preconditioner for the all-at-once evolutionary PDE system which has block Toeplitz structure. Through reordering of variables to obtain a symmetric system, we are able to rigorously establish convergence bounds for MINRES which guarantee a number of iterations independent of the number of time-steps for the all-at-once system. If the spatial differential operators are simultaneously diagonalizable, we are able to quickly apply the preconditioner through use of a sine transform; and for those that are not, we are able to use an algebraic multigrid process to provide a good approximation. Results are presented for solution to both the heat and convection diffusion equations. [ABSTRACT FROM AUTHOR]

We propose a simple proof of a generalized Mirsky inequality comparing the differences of singular values of two matrices with the singular values of their difference. We then discuss the implication of this generalized inequality for the recovery of low-rank matrices via concave minimization. [ABSTRACT FROM AUTHOR]

The computation of the matrix exponential is a ubiquitous operation in numerical mathematics, and for a general, unstructured n x n matrix it can be computed in O(n3) operations. An interesting problem arises if the input matrix is a Toeplitz matrix, for example as the result of discretizing integral equations with a time invariant kernel. In this case it is not obvious how to take advantage of the Toeplitz structure, as the exponential of a Toeplitz matrix is, in general, not a Toeplitz matrix itself. The main contribution of this work are fast algorithms for the computation of the Toeplitz matrix exponential. The algorithms have provable quadratic complexity if the spectrum is real, or sectorial, or, more generally, if the imaginary parts of the rightmost eigenvalues do not vary too much. They may be eficient even outside these spectral constraints. They are based on the scaling and squaring framework, and their analysis connects classical results from rational approximation theory to matrices of low displacement rank. As an example, the developed methods are applied to Merton's jump-diffusion model for option pricing. [ABSTRACT FROM AUTHOR]

Recently, nonlocal models attract the wide interest of scientists. They mainly come from two applied scientific fields: peridynamics and anomalous diffusion. Even though the matrices of the algebraic equation corresponding to the nonlocal models are usually Toeplitz (denote a0 as the principal diagonal element, a1 as the trailing diagonal element, etc). There are still some differences for the models in these two fields. For the model of anomalous diffusion, a0/a1 is uniformly bounded; most of the time, a0/a1 of the model for peridynamics is unbounded as the step size h tends to zero. Based on the uniform boundedness of a0/a1, the convergence of the two-grid method is well established [R. H. Chan, Q.-S. Chang, and H.-W. Sun, SIAM J. Sci. Comput., 19 (1998), pp. 516-529; H. Pang and H. Sun, J. Comput. Phys., 231 (2012), pp. 693-703; M. H. Chen, Y. T.Wang, X. Cheng, and W. H. Deng, BIT, 54 (2014), pp. 623-647]. This paper provides the detailed proof of the convergence of the two-grid method for the nonlocal model of peridynamics. Some special cases of the full multigrid and the V-cycle multigrid methods are also discussed. The numerical experiments are performed to verify the convergence. [ABSTRACT FROM AUTHOR]

In this paper we propose an efficient method for a convex optimization problem which involves a large nonsymmetric and non-Toeplitz matrix. The proposed method is an instantiation of the alternating direction method of multipliers applied in Krylov subspace. Our method offers significant advantages in computational speed for the convex optimization problems involved with general matrices of large size. We apply the proposed method to the restoration of spatially variant blur. The matrix representing spatially variant blur is not block circulant with circulant blocks (BCCB). Efficient implementation based on diagonalization of BCCB matrices by the discrete Fourier transform is not applicable for spatially variant blur. Since the proposed method can efficiently work with general matrices, the restoration of spatially variant blur is a good application of our method. Experimental results for total variation restoration of spatially variant blur show that the proposed method provides meaningful solutions in a short time. [ABSTRACT FROM AUTHOR]

We consider the stiffness matrices arising from the Galerkin B-spline isogeometric analysis discretization of classical elliptic problems. By exploiting their specific spectral properties, compactly described by a symbol, we design an efficient multigrid method for the fast solution of the related linear systems. The convergence rate of general-purpose multigrid methods, based on classical stationary smoothers, is optimal (i.e., bounded independently of the matrix size), but it also worsens exponentially with respect to the spline degree. The symbol all ows us to give a detailed theoretical explanation of this exponential worsening in the case of the two-grid scheme. In addition, thanks to a specific factorization of the symbol, we are able to design an ad hoc multigrid method with an effective preconditioned CG or GMRES smoother at the finest level, in the spirit of the multiiterative idea. The convergence rate of this multi-iterative multigrid method is not only optimal but also robust (i.e., bounded substantially independently of the spline degree). This can again be explained by the symbol, in combination with the theory of generalized locally Toeplitz sequences. A selected set of numerical experiments confirms our symbol-based analysis, as well as the effectiveness of the proposed multi-iterative multigrid method, also for larger spline degree. [ABSTRACT FROM AUTHOR]

We consider the Poisson equation (I -- P)u = g, where P is the transition matrix of a quasi-birth-and-death process with infinitely many levels, g is a given infinite dimensional vector, and u is the unknown. Our main result is to provide the general solution of this equation. To this purpose we use the block tridiagonal and block Toeplitz structure of the matrix P to obtain a set of matrix difference equations, which are solved by constructing suitable resolvent triples. [ABSTRACT FROM AUTHOR]

When one computes an average of positive definite (PD) matrices, the preservation of additional matrix structure is desirable for interpretations in applications. An interesting and widely present structure is that of PD Toeplitz matrices, which we endow with a geometry originating in signal processing theory. As an averaging operation, we consider the barycenter, or minimizer of the sum of squared intrinsic distances. The resulting barycenter, the Kähler mean, is discussed along with its origin. Also, a generalization of the mean towards PD (Toeplitz-block) block-Toeplitz matrices is discussed. For PD Toeplitz-block block-Toeplitz matrices, we derive the generalized barycenter, or generalized Kaahler mean, and a greedy approximation. This ap proximation is shown to be close to the generalized mean with a significantly lower computational cost. [ABSTRACT FROM AUTHOR]

In this paper, we study the linear systems arising from the discretization of timedependent space-fractional diffusion equations. By using a finite difference discretization scheme for the time derivative and a finite volume discretization scheme for the space-fractional derivative, Toeplitz-like linear systems are obtained. We propose using the approximate inverse-circulant preconditioner to deal with such Toeplitz-like matrices, and we show that the spectra of the corresponding preconditioned matrices are clustered around 1. Experimental results on time-dependent and space-fractional diffusion equations are presented to demonstrate that the preconditioned Krylov subspace methods converge very quickly. [ABSTRACT FROM AUTHOR]

We investigate the spectrum of one-dimensional periodic lattices with geometric interaction in the asymptotic regimes of vanishing attenuation and long lattice lengths. Our mathe- matical analysis is motivated by a number of real-world problems where chains of radiative elements are embedded in waveguides in order to enhance the interlattice interactions and collective resonances. Under certain assumptions the problem reduces to determining the eigenvalues of a complex-valued, symmetric Toeplitz matrix whose elements follow a geometric progression. A recurrence formula is derived for the characteristic polynomial and becomes the basis for developing spectral asymptotics and resolving the eigenvalue patterns on the complex plane. Analytical asymptotic formulas are derived for the eigenvalues and their accuracy is assessed numerically. [ABSTRACT FROM AUTHOR]

We present a superfast divide-and-conquer method for finding all the eigenvalues as well as all the eigenvectors (in a structured form) of a class of symmetric matrices with off-diagonal ranks or numerical ranks bounded by r, as well as the approximation accuracy of the eigenvalues due to off-diagonal compression. More specifically, the complexity is O(r²n log n) + O(rn log² n), where n is the order of the matrix. Such matrices are often encountered in practical computations with banded matrices, Toeplitz matrices (in Fourier space), and certain discretized problems. They can be represented or approximated by hierarchically semiseparable (HSS) matrices. We show how to preserve the HSS structure throughout the dividing process that involves recursive updates and how to quickly perform stable eigendecompositions of the structured forms. Various other numerical issues are discussed, such as computation reuse and deation. The structure of the eigenvector matrix is also shown. We further analyze the structured perturbation, i.e., how compression of the off-diagonal blocks impacts the accuracy of the eigenvalues. They show that rank structured methods can serve as an effective and efficient tool for approximate eigenvalue solutions with controllable accuracy. The algorithm and analysis are very useful for finding the eigendecomposition of matrices arising from some important applications and can be modified to find SVDs of nonsymmetric matrices. The efficiency and accuracy are illustrated in terms of Toeplitz and discretized matrices. Our method requires significantly fewer operations than a recent structured eigensolver, by nearly an order of magnitude. [ABSTRACT FROM AUTHOR]

The pseudo-polar Fourier transform is a specialized nonequally spaced Fourier transform, which evaluates the Fourier transform on a near-polar grid known as the pseudo-polar grid. The advantage of the pseudo-polar grid over other nonuniform sampling geometries is that the trans-formation, which samples the Fourier transform on the pseudo-polar grid, can be inverted using a fast and stable algorithm. For other sampling geometries, even if the nonequally spaced Fourier transform can be inverted, the only known algorithms are iterative. The convergence speed of these algorithms and their accuracy are difficult to control, as they depend both on the sampling ge- ometry and on the unknown reconstructed object. In this paper, a direct inversion algorithm for the three-dimensional pseudo-polar Fourier transform is presented. The algorithm is based only on one-dimensional resampling operations and is shown to be significantly faster than existing iterative inversion algorithms. [ABSTRACT FROM AUTHOR]

For a family of near banded Toeplitz matrices, generalized characteristic polynomials are shown to be orthogonal polynomials of two variables, which include the Chebyshev polynomials of the second kind on the deltoid as a special case. These orthogonal polynomials possess a maximal number of real common zeros, which generate a family of Gaussian cubature rules in two variables. [ABSTRACT FROM AUTHOR]

We study the spectral properties of the stiffness matrices coming from the approximation of a d-dimensional second order elliptic differential problem by the ℚp Lagrangian finite element method (FEM); here, p = (p1,...,pd) and pj represents the polynomial approximation degree in the jth direction. After presenting a construction of these matrices, we investigate the conditioning and the spectral distribution in the Weyl sense, and we determine the so-called (spectral) symbol describing the asymptotic spectrum. We also study the properties of the symbol, which turns out to be a d-variate function taking values in the space of N(p) X N(p) Hermitian matrices, where N(p) = Πj=1dpj. Unlike the stiffness matrices coming from the p-degree B-spline isogeo-metric analysis approximation of the same differential problem, where a unique d-variate real-valued function describes all the spectrum, here the spectrum is described by N(p) different functions, i.e., the N(p) eigenvalues of the symbol, which are well-separated, far away, and exponentially diverging with respect to p and d. This very involved picture provides an explanation of: (a) the difficulties encountered in designing robust solvers, with convergence speed independent of the matrix size, of the approximation parameters p, and of the dimensionality d; (b) the convergence deterioration of known iterative methods, already observed in practice for moderate p and d. [ABSTRACT FROM AUTHOR]

We develop a fast and accurate finite element method for space-fractional dispersion equations in two space dimensions, which are expressed in terms of fractional directional derivatives in all the directions that are integrated with respect to a probability measure on the unit circle. The fast method significantly reduces the computational work of solving the discrete linear algebraic systems from O(N³) by a direct solver to O(N log N) per iteration and a memory requirement from O(N²) to O(N). Furthermore, the fast method reduces the evaluation of the stiffness matrix, which often constitutes a large portion of the CPU time, from O(N²) to O(N). The developed preconditioned fast Krylov subspace iterative solver significantly reduces the number of iterations in a Krylov subspace iterative method and may improve the convergence behavior of the solver. Numerical results show the utility of the method. [ABSTRACT FROM AUTHOR]

The equation with the time fractional substantial derivative and space fractional derivative describes the distribution of the functionals of the Lévy flights, and the equation is derived as the macroscopic limit of the continuous time random walk in unbounded domain and the Lévy flights have divergent second order moments. However, in more practical problems, the physical domain is bounded and the involved observables have finite moments. Then the modified equation can be derived by tempering the Lévy measure of the Lévy flights and the corresponding tempered space fractional derivative is introduced. This paper focuses on providing the high order algorithms for the modified equation, i.e., the equation with the time fractional substantial derivative and space tempered fractional derivative. More concretely, the contributions of this paper are as follows: (1) Detailed numerical stability analysis and error estimates of the schemes with first order accuracy in time and second order in space are given in complex space, which is necessary since the inverse Fourier transform needs to be made for getting the distribution of the functionals after solving the equation. (2) We further propose the schemes with high order accuracy in both time and space, and the techniques of treating the issue of keeping the high order accuracy of the schemes for nonhomogeneous boundary/initial conditions are introduced. (3) Multigrid methods are effectively used to solve the obtained algebraic equations which still have the Toeplitz structure. (4) We perform extensive numerical experiments, including verifying the high convergence orders and simulating the physical system which needs to numerically make the inverse Fourier transform to the numerical solutions of the equation. [ABSTRACT FROM AUTHOR]

Circulant preconditioning for symmetric Toeplitz linear systems is well established; theoretical guarantees of fast convergence for the conjugate gradient method are descriptive of the convergence seen in computations. This has led to robust and highly efficient solvers based on use of the fast Fourier transform exactly as originally envisaged in [G. Strang, Stud. Appl. Math., 74 (1986), pp. 171-176]. For nonsymmetric systems, the lack of generally descriptive convergence theory for most iterative methods of Krylov type has provided a barrier to such a comprehensive guarantee, though several methods have been proposed and some analysis of performance with the normal equations is available. In this paper, by the simple device of reordering, we rigorously establish a circulant preconditioned short recurrence Krylov subspace iterative method of minimum residual type for nonsymmetric (and possibly highly nonnormal) Toeplitz systems. Convergence estimates similar to those in the symmetric case are established. [ABSTRACT FROM AUTHOR]

The fractional diffusion equation is discretized by an implicit finite difference scheme with the shifted Grünwald formula, which is unconditionally stable. The coefficient matrix of the discretized linear system is equal to the sum of a scaled identity matrix and two diagonal-times-Toeplitz matrices. Standard circulant preconditioners may not work for such Toeplitz-like linear systems. The main aim of this paper is to propose and develop approximate inverse preconditioners for such Toeplitz-like matrices. An approximate inverse preconditioner is constructed to approximate the inverses of weighted Toeplitz matrices by circulant matrices, and then combine them together row-by-row. Because of Toeplitz structure, both the discretized coefficient matrix and the preconditioner can be implemented very efficiently by using fast Fourier transforms. Theoretically, we show that the spectra of the resulting preconditioned matrices are clustered around one. Thus Krylov subspace methods with the proposed preconditioner converge very fast. Numerical examples are given to demonstrate the effectiveness of the proposed preconditioner and show that its performance is better than the other testing preconditioners. [ABSTRACT FROM AUTHOR]

This paper shows that the matrix fraction description, given by a pair {A(z) B(z)} of matrix polynomials of z, for a linear time-invariant system may not be unique even if A(z) is monic, A(z) and B(z) have no common left factor, and the matrix coefficients corresponding to the highestorder terms of A(z) and B(z) are full row rank. The orders of all possible matrix fraction descriptions (MFDs) of a given system are completely characterized. Testing criteria for determining whether a matrix pair is an MFD of the system are derived, which involve rank tests of certain Toeplitz matrices derived from either the impulse response or output correlation functions of the system. A decision procedure is devised that generates sequentially all MFDs for a given system. Identification algorithms are introduced that estimate all MFDs of a given system from its input-output data or output data only. The results are then extended to cover ARMAX systems. [ABSTRACT FROM AUTHOR]

We consider half-infinite triangular Toeplitz matrices with slow decay of the elements and prove under a monotonicity condition that the elements of the inverse matrix, as well as the elements of the fundamental matrix, decay to zero. We provide a quantitative description of the decay of the fundamental matrix in terms of p-norms. The results add to the classical results of Jaffard and Vecchio and are illustrated by numerical examples. [ABSTRACT FROM AUTHOR]

When the heat equation and wave equation are approximated by ut=-Ku and utt = -Ku (discrete in space), the solution operators involve e-Kt, √K, cos(√Kt), and sinc(√Kt). We compute these four matrices and find accurate approximations with a variety of boundary conditions. The second difference matrix K is Toeplitz (shift-invariant) for Dirichlet boundary conditions, but we show why e-Kt also has a Hankel (anti-shift-invariant) part. Any symmetric choice of the four corner entries of K leads to Toeplitz plus Hankel in all functions f(K). Overall, this article is based on diagonalizing symmetric matrices, replacing sums by integrals, and computing Fourier coefficients. [ABSTRACT FROM AUTHOR]

We study the conditioning of nonlocal integral operators with singular and integrable kernels in fractional Sobolev spaces. These operators are used, for instance, in peridynamics formulation and nonlocal diffusion. In one dimension (1D), we present sharp quantification of the extremal eigenvalues in all three parameters: size of nonlocality, mesh size, and regularity of the fractional Sobolev space. We accomplish sharpness both rigorously and numerically. For the minimal eigenvalue, we obtain sharpness analytically by using a nonlocal characterization of Sobolev spaces. We verify this estimate by exploiting the Toeplitz property of the stiffness matrix. However, the analytical approach fails to give sharp quantification of the maximal eigenvalue. Hence, in 1D, we take an algebraic approach by directly working with the stiffness matrix entries, which have complicated expressions due to all three parameters. We systematically characterize the nonzero entries and dramatically simplify their expressions by using convenient algebra. We establish the zero row sum property of the stiffness matrix and negativity of the off-diagonal entries. Eventually, we arrive at sharpness through the use of the Gerschgorin circle theorem. [ABSTRACT FROM AUTHOR]

The main aim of this paper is to develop a fast algorithm for solving weighted Toeplitz regularized least squares problems arising from image restoration. Based on augmented system formulation, we develop new Hermitian and skew-Hermitian splitting (HSS) preconditioners for solving such linear systems. The advantage of the proposed preconditioner is that the blurring matrix, weighting matrix, and regularization matrix can be decoupled such that the resulting preconditioner is not expensive to use. We show that for a preconditioned system that is derived from a saddle point structure of size (m+n) x (m+n), the preconditioned matrix has an eigenvalue at 1 with multiplicity n and the other m eigenvalues of the form 1 -- λ with |λ| < 1. We also study how to choose the HSS parameter to minimize the magnitude of λ, and therefore the Krylov subspace method applied to solving the preconditioned system converges very quickly. Experimental results for image restoration problems are reported to demonstrate that the performance of the proposed preconditioner is better than the other testing preconditioners. [ABSTRACT FROM AUTHOR]

We present some superfast (O((m + n) log²(m + n)) complexity) and stable structured direct solvers for m× n Toeplitz least squares problems. Based on the displacement equation, a Toeplitz matrix T is first transformed into a Cauchy-like matrix C, which can be shown to have small off-diagonal numerical ranks when the diagonal blocks are rectangular. We generalize standard hierarchically semiseparable (HSS) matrix representations to rectangular ones, and construct a rectangular HSS approximation to C in nearly linear complexity with randomized sampling and fast multiplications of C with vectors. A new URV HSS factorization and a URV HSS solution are designed for the least squares solution. We also present two structured normal equation methods. Systematic error and stability analysis for our HSS methods is given, which is also useful for studying other HSS and rank structured methods. We derive the growth factors and the backward error bounds in the HSS factorizations, and show that the stability results are generally much better than those in dense LU factorizations with partial pivoting. Such analysis has not been done before for HSS matrices. The solvers are tested on various classical Toeplitz examples ranging from well-conditioned to highly ill-conditioned ones. Comparisons with some recent fast and superfast solvers are given. Our new methods are generally much faster, and give better (or at least comparable) accuracies, especially for ill-conditioned problems. [ABSTRACT FROM AUTHOR]

We study the tensor structure of two operations: the transformation of a given multidimensional vector into a multilevel Toeplitz matrix and the convolution of two given multidi-mensional vectors. We show that the low-rank tensor structure of the input is preserved in the output and propose efficient algorithms for these operations in the newly introduced quantized tensor train (QTT) format. Consider a d-dimensional 2n x ⋯ x 2n-vector x. If it is represented elementwise, the number of parameters is (2n)d. However, if we assume that x is given in a QTT representation with ranks bounded by p, the number of parameters is reduced to O (dp² log ra). Under this assumption we show how the multilevel Toeplitz matrix generated by x can be obtained in the QTT format with ranks bounded by 2p in O (dp² log ra) operations. We also describe how the convolution x * y of x and a d-dimensional n x ⋯ x n-vector y can be computed in the QTT format with ranks bounded by 2t in O (dt² logra) operations, provided that the matrix xy' is given in the same format with ranks bounded by t. We exploit an algorithm for the inexact matrix-vector multiplication in the QTT format to accelerate the convolution algorithm dramatically. The performance of our approach to convolution is demonstrated with numerical examples, including the computation of the Newton potential of a strong cusp on fine grids with up to 220 x 220 x 220 points in three dimensions. [ABSTRACT FROM AUTHOR]

We introduce a flexible optimization framework for nuclear norm minimization of matrices with linear structure, including Hankel, Toeplitz, and moment structures and catalog applications from diverse fields under this framework. We discuss various first-order methods for solving the resulting optimization problem, including alternating direction methods of multipliers, proximal point algorithms, and gradient projection methods. We perform computational experiments to compare these methods on system identification problems and system realization problems. For the system identification problem, the gradient projection method (accelerated by Nesterov's extrapolation techniques) and the proximal point algorithm usually outperform other first-order methods in terms of CPU time on both real and simulated data, for small and large regularization parameters, respectively, while for the system realization problem, the alternating direction method of multipliers, as applied to a certain primal reformulation, usually outperforms other first-order methods in terms of CPU time. We also study the convergence of the proximal alternating direction methods of multipliers used in this paper. [ABSTRACT FROM AUTHOR]

In this paper, we show large deviations for random step functions of type Zn(t) = 1/nΣ[nt] k=1 X2 k, where {Xk}k is a stationary Gaussian process. We deal with the associated random measures vn = 1/n Σn k. The proofs require a Szegö theorem for generalized Toeplitz matrices which is analogous to a result of Kac, Murdoch, and Szegö [J. Rational Mech. Anal., 2 (1953), pp. 767-800]. We also study the polygonal line built on Z2 kбk/n. The proofs require a Szegö theorem for generalized Toeplitz matrices which is analogous to a result of Kac, Murdoch, and Szegö [J. Rational Mech. Anal., 2 (1953), pp. 767-800]. We also study the polygonal line built on Zn(t) and show moderate deviations for both random families. [ABSTRACT FROM AUTHOR]

This paper presents a preconditioning strategy applied to certain types of kernel matrices that are increasingly ill-conditioned. The ill-conditioning of these matrices is tied to the unbounded variation of the Fourier transform of the kernel function. Hence, the basic idea is to differentiate the kernel in order to suppress the variation. The idea resembles some existing preconditioning methods for Toeplitz matrices, where the theory heavily relies on the underlying fixed generating function. The theory does not apply to the case of a fixed domain with increasingly fine discretizations because the generating function depends on the grid size. For this case, we prove equal distribution results on the spectrum of the resulting matrices. Furthermore, the proposed preconditioning technique also applies to non-Toeplitz matrices, thus eliminating the reliance on a regular grid structure of the points. The preconditioning strategy can be used to accelerate an iterative solver for solving linear systems with respect to kernel matrices. [ABSTRACT FROM AUTHOR]

In this paper a new O(N log³ N) solver for N x N Toeplitz-like systems, based on a divide and conquer technique, is presented. Similarly to the superfast algorithm MBA for the inversion of a Toeplitz-like matrix [R. R. Bitmead and B. D. O. Anderson, Linear Algebra Appl., 34 (1980), pp. 103-116; M. Morf, Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing, 1980, pp. 954-959], it exploits the displacement properties. In order to avoid the well-known numerical instability of the explicit inversion, the new algorithm relies on the triangular factorization and back-substitution formula for the system seen as a 2x2 block system with blocks of half size. This idea is the one used in [M. Stewart, SIAM J. Matrix Anal. Appl., 25 (2003), pp. 669-693] to improve the numerical stability of superfast methods based on the generalized Schur algorithm for positive definite Toeplitz matrices, but the algorithm we propose can be applied also to nonsymmetric Toeplitz-like systems. The stability of the algorithm is examined through numerical experiments. [ABSTRACT FROM AUTHOR]

The computation of the structured pseudospectral abscissa and radius (with respect to the Frobenius norm) of a Toeplitz matrix is discussed and two algorithms based on a low-rank property to construct extremal perturbations are presented. The algorithms are inspired by those considered in [N. Guglielmi and M. Overton, SIAM J. Matrix Anal. Appl., 32 (2011), pp. 1166-1192] for the unstructured case, but their extension to structured pseudospectra and analysis presents several difficulties. Natural generalizations of the algorithms, allowing us to draw significant sections of the structured pseudospectra in proximity of extremal points, are also discussed. Since no algorithms are available in the literature to draw such structured pseudospectra, the approach we present seems promising to extend existing software tools (Eigtool, Seigtool) to structured pseudospectra representation for Toeplitz matrices. We discuss local convergence properties of the algorithms and show some applications to a few illustrative examples. [ABSTRACT FROM AUTHOR]

We propose a superfast solver for Toeplitz linear systems based on rank structured matrix methods and randomized sampling. The solver uses displacement equations to transform a Toeplitz matrix T into a Cauchy-like matrix C, which is known to have low-numerical-rank offdiagonal blocks. Thus, we design a fast scheme for constructing a hierarchically semiseparable (HSS) matrix approximation to C, where the HSS generators have internal structures. Unlike classical HSS methods, our solver employs randomized sampling techniques together with fast Toeplitz matrixvector multiplications, and thus converts the direct compression of the off-diagonal blocks of C into the compression of much smaller blocks. A strong rank-revealing QR factorization method is used to generate/preserve certain special structures, and also to ensure stability. A fast ULV HSS factorization scheme is provided to take advantage of the special structures. We also propose a precomputation procedure for the HSS construction so as to further improve the efficiency. The complexity of these methods is significantly lower than some similar Toeplitz solvers for large matrix size n. Detailed flop counts are given, with the aid of a rank relaxation technique. The total cost of our methods includes O(n) flops for HSS operations and O(n log² n) flops for matrix multiplications via FFTs, where n is the order of T. Various numerical tests on classical examples, including ill-conditioned ones, demonstrate the efficiency, and also indicate that the methods are stable in practice. This work shows a practical way of using randomized sampling in the development of fast rank structured methods. [ABSTRACT FROM AUTHOR]

The paper establishes error orders for integral limit approximations to the traces of products of truncated Toeplitz operators generated by integrable real symmetric functions defined on the real line. These approximations and the corresponding error bounds are of importance in the statistical analysis of continuous-time stationary processes (asymptotic distributions and large deviations of Toeplitz-type quadratic functionals, estimation of the spectral parameters and functionals, etc.). An explicit second-order asymptotic expansion is found for the trace of a product of two truncated Toeplitz operators generated by the spectral densities of continuous-time stationary fractional Riesz-Bessel motions. The order of magnitude of the second term in this expansion is shown to depend on the long-memory parameters of the processes. Also, it is shown that the pole in the first-order approximation is removed by the second-order term, which provides a substantially improved approximation to the original functional. [ABSTRACT FROM AUTHOR]

The notion of compact Fourier analysis (CFA) is discussed. CFA allows description of multigrid (MG) in a nutshell and offers a clear overview on all MG components. The principal idea of CFA is to model the MG mechanisms by means of scalar generating functions and matrix functions (block symbols). The formalism of the CFA approach is presented by describing the symbols of the fine and coarse grid problems, the prolongation and restriction, the smoother, and the coarse grid correction, resp., smoothing corrections. CFA uses matrix functions and their features (e.g., product, inverse, adjugate, norm, spectral radius, eigenvectors, eigenvalues of multilevel ω-circulant matrices), and scalar functions and their roots. This leads to an elementary description and allows for an easy analysis of MG algorithms. A first application is to utilize CFA for deriving MG as a direct solver, i.e., an MG cycle that will converge in just one iteration step. Necessary and sufficient conditions that have to be fulfilled by the MG components are given for obtaining MG functioning as a direct solver. Furthermore, new general and practical smoothers and transfer operators that lead to efficient MG methods are introduced. In addition, we study sparse approximations of the Galerkin coarse grid operator yielding efficient and practicable MG algorithms (approximately direct solvers). Numerical experiments demonstrate the theoretical results. [ABSTRACT FROM AUTHOR]