Search

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]

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]

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 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]

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]

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]

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 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]

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]

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]

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 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]

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]

We present a multilevel method for discrete ill-posed problems arising from the discretization of Fredholm integral equations of the first kind. In this method, we use the Haarwavelet transform to define restriction and prolongation operators within a multigrid-type iteration. The choice of the Haar wavelet operator has the advantage of preserving matrix structure, such as Toeplitz, between grids, which can be exploited to obtain faster solvers on each level where an edge preserving Tikhonov regularization is applied. Finally, we present results that indicate the promise of this approach for restoration of signals and images with edges. [ABSTRACT FROM AUTHOR]

The convergence of multigrid methods can be analyzed based on a Fourier analysis of the method or by proving certain inequalities that have to be fulfilled by the smoother and by the coarse grid correction separately. Here, we analyze the multigrid method for the constant coefficient Poisson equation with a compact Fourier analysis using the formalism of multilevel Toeplitz matrices and their generating functions or symbols. The Fourier analysis is applied for determining the smoothing factor and the overall error of the combined smoothing and coarse grid correction error reduction of a twogrid step (TGS) by representing the twogrid step explicitly by a symbol. If the effects of the smoothing correction and the coarse grid correction are orthogonal to each other, then, in a twogrid step, the error is removed in one step, and the twogrid method can be considered as a direct solver. If the coarse linear system is identical to the original matrix, then the same projection and smoother again make the twogrid step a direct solver. Hence, the multigrid cycle has to be applied only once with one smoothing step on each level, and therefore the whole multigrid method can be considered as a direct solver. In this paper we want to identify multigrid as a direct solver in one-dimension with a coarse system derived by linear and constant interpolation, and in two-dimensions for a modified projection related to [U. Trottenberg, C.W. Osterlee, and A. Schüller, Multigrid, Academic Press, San Diego, 2001, paragraph A.2.3]. By studying not only smoothing and coarse grid correction separately but the symbol of the full twogrid step as well, we are furthermore able to derive better convergence estimates and information on how to find efficient combinations of projection and smoother. As smoothers we consider also approximate inverse smoothers, colored smoothers that take into account the different character of grid points, and rank-reducing smoothers that coincide with the given matrix on a large number of entries. Numerical examples show the effectiveness of the new approach. [ABSTRACT FROM AUTHOR]

We study the Newton method and a fixed-point method for solving the system of nonlinear equations arising from the Sinc--Galerkin discretization of the Burgers equations. In each step of the Newton method or the fixed-point method, a structured subsystem of linear equations is obtained and needs to be solved numerically. In this paper, preconditioning techniques are applied to solve such linear subsystems. The bounds for eigenvalues of the preconditioned matrices are derived and numerical examples are given to illustrate the effectiveness of the proposed methods. We also find that a combination of the Newton/fixed-point iteration with the preconditioned GMRES method is quite efficient for the Sinc--Galerkin discretization of the Burgers equations. [ABSTRACT FROM AUTHOR]

Large 2-level Toeplitz systems arise in a variety of applications (see, e.g., [R. H. Chan and M. Ng, SIAM Rev., 38 (1996), pp. 427-482]) for which efficient numerical methods for their solution are required. Some successful numerical techniques need the explicit knowledge of the generating function f of the considered system Tn(f)x = b, an assumption that usually is not fulfilled in real applications. In this paper we analyze and complete the procedure proposed in [D. Noutsas, S. Serra Capizzano, and P. Vassalos, Numer. Linear Algebra Appl., 12 (2005), pp. 231-239] for the 2-level case. In such a way, from the knowledge of the coefficients of Tn(f), we determine optimal preconditioning strategies for the solution of our systems. Finally, some numerical experiments are performed and discussed in connection with our theoretical analysis. [ABSTRACT FROM AUTHOR]

In this paper, we study factorized banded inverse preconditioners for matrices with Toeplitz structure. We show that if a Toeplitz matrix T has certain off-diagonal decay property, then the factorized banded inverse preconditioner approximates T-1 accurately, and the spectra of these preconditioned matrices are clustered around 1. In nonlinear image restoration applications, Toeplitz-related systems of the form I + T* D T are required to solve, where D is a positive nonconstant diagonal matrix. We construct inverse preconditioners for such matrices. Numerical results show that the performance of our proposed preconditioners are superior to that of circulant preconditioners. A two-dimensional nonlinear image restoration example is also presented to demonstrate the effectiveness of the proposed preconditioner. [ABSTRACT FROM AUTHOR]

We develop a new algorithm for the fast computation of discrete sums f(yj) := ∑Nk=1 αkK(yj −xk) (j = 1, . . . ,M) based on the recently developed fast Fourier transform (FFT) at nonequispaced knots. Our algorithm, in particular our regularization procedure, is simply structured and can be easily adapted to different kernels K. Our method utilizes the widely known FFT and can consequently incorporate advanced FFT implementations. In summary, it requires O(N log N + M) arithmetic operations. We prove error estimates to obtain clues about the choice of the involved parameters and present numerical examples in one and two dimensions. [ABSTRACT FROM AUTHOR]

Studies the solutions of symmetric positive definite Toeplitz systems. Derivation of solutions by recursive-based preconditioned conjugate gradient method; Numerical results illustrating the effectiveness of the approach.

An algorithm is implemented for the computation of inertial manifolds to arbitrary accuracy. It is applied to the Kuramoto­Sivashinsky equation to recover the high wave number components of elements on the global attractor. The computational results demonstrate that this algorithm is significantly more accurate than previous methods, when compared after just a few iterations. An essential feature of the algorithm is that it does not require that the entire manifold be constructed. Hence the algorithm is particularly suited for computing trajectories on the manifold, and its computational complexity is independent of the dimension. [ABSTRACT FROM AUTHOR]

In this paper, we are interested in the iterative solution of ill-conditioned Toeplitz systems generated by continuous nonnegative real-valued functions f with a finite number of zeros. We construct new w-circulant preconditioners without explicit knowledge of the generating function f by approximating f by its convolution f * KN with a suitable positive reproducing kernel KN. By the restriction to positive kernels we obtain positive definite preconditioners. Moreover, if f has only zeros of even order ≤ 2s, then we can prove that the property [This symbol cannot be presented in ASCII format]π-π t2k KN(t)dt ≤ CN-2k (k = 0, … ,s) of the kernel is necessary and sufficient to ensure the convergence of the PCG method in a number of iteration steps independent of the dimension N of the system. Our theoretical results were confirmed by numerical tests. [ABSTRACT FROM AUTHOR]

We give an overview of fast algorithms for solving least squares problems with Toeplitz structure, based on generalization of the classical Schur algorithm, and discuss their stability properties. In order to obtain more accurate triangular factors of a Toeplitz matrix as well as accurate solutions for the least squares problems, methods based on corrected seminormal equations (CSNE) can be used. We show that the applicability of the generalized Schur algorithm is considerably enhanced when the algorithm is used in conjunction with CSNE. Several numerical tests are reported, where different variants of the generalized Schur algorithm and CSNE are compared for their accuracy and speed. [ABSTRACT FROM AUTHOR]

Block-boundary value methods are considered and applied to solve initial value problems. Results on their convergence and stability properties are presented and their intermediate position between the class of multistep and Runge{Kutta methods is established. [ABSTRACT FROM AUTHOR]

Many ill-posed problems are solved using a discretization that results in a least squares problem or a linear system involving a Toeplitz matrix. The exact solution to such problems is often hopelessly contaminated by noise, since the discretized problem is quite ill conditioned, and noise components in the approximate null-space dominate the solution vector. Therefore we seek an approximate solution that does not have large components in these directions. We use a preconditioned conjugate gradient algorithm to compute such a regularized solution. A unitary change of coordinates transforms the Toeplitz matrix to a Cauchy-like matrix, and we choose our preconditioner to be a low rank Cauchy-like matrix determined in the course of Gu's fast modified complete pivoting algorithm. We show that if the kernel of the ill-posed problem is smooth, then this preconditioner has desirable properties: the largest singular values of the preconditioned matrix are clustered around one, the smallest singular values, corresponding to the lower subspace, remain small, and the upper and lower spaces are relatively unmixed. The preconditioned algorithm costs only O(n lg n) operations per iteration for a problem with n variables. The effectiveness of the preconditioner for filtering noise is demonstrated on three examples. [ABSTRACT FROM AUTHOR]

A method for computing the smallest eigenvalue of a symmetric positive definite Toeplitz matrix by application of Newton's method to the characteristic polynomial has been recently introduced by Mastronardi and Boley [SIAM J. Sci. Comput., 20 (1999), pp. 1921–1927]. Though considerably slower than methods developed by the authors [W. Mackens and H. Voss, SIAM J. Matrix. Anal. Appl., 18 (1997), pp. 521–534], [W. Mackens and H. Voss, Linear Algebra Appl., 275⁄276 (1998), pp. 401–415], [H. Voss, Linear Algebra Appl., 287 (1999), pp. 359–371] the new approach is conceptually much simpler. In this paper we improve the performance of the new method substantially while keeping its simplicity. [ABSTRACT FROM AUTHOR]

Blur removal is an important problem in signal and image processing. The blurring matrices obtained by using the zero boundary condition (corresponding to assuming dark background outside the scene) are Toeplitz matrices for one-dimensional problems and block-Toeplitz­Toeplitz- block matrices for two-dimensional cases. They are computationally intensive to invert especially in the block case. If the periodic boundary condition is used, the matrices become (block) circulant and can be diagonalized by discrete Fourier transform matrices. In this paper, we consider the use of the Neumann boundary condition (corresponding to a reflection of the original scene at the boundary). The resulting matrices are (block) Toeplitz-plus-Hankel matrices. We show that for symmetric blurring functions, these blurring matrices can always be diagonalized by discrete cosine transform matrices. Thus the cost of inversion is significantly lower than that of using the zero or periodic boundary conditions. We also show that the use of the Neumann boundary condition provides an easy way of estimating the regularization parameter when the generalized cross-validation is used. When the blurring function is nonsymmetric, we show that the optimal cosine transform preconditioner of the blurring matrix is equal to the blurring matrix generated by the symmetric part of the blurring function. Numerical results are given to illustrate the efficiency of using the Neumann boundary condition. [ABSTRACT FROM AUTHOR]