Asymptotically fast solution of toeplitz and related systems of linear equations

We present an inversion algorithm for the solution of a generic N X N Toeplitz system of linear equations with computational complexity O(Nlog²N) and storage requirements O(N). The algorithm relies upon the known structure of Toeplitz matrices and their inverses and achieves speed through a doubling method. All the results are derived and stated in terms of the recent concept of displacement rank, and this is used to extend the scope of the algorithm to include a wider class of matrices than just Toeplitz and also to include block Toeplitz matrices.     

Fast Wavelet Transform for Toeplitz Matrices and Property Analysis

Fast wavelet transform algorithms for Toeplitz matrices are proposed in this paper. Distinctive from the well known discrete trigonometric transforms, such as the discrete cosine transform （DCT） and the discrete Fourier transform （DFT） for Toeplitz matrices, the new algorithms are achieved by compactly supported wavelet that preserve the character of a Toeplitz matrix after transform, which is quite useful in many applications involving a Toeplitz matrix. Results of numerical experiments show that the proposed method has good compression performance similar to using wavelet in the digital image coding. Since the proposed algorithms turn a dense Toeplitz matrix into a band-limited form, the arithmetic operations required by the new algorithms are O（N） that are reduced greatly compared with O（N log N） by the classical trigonometric transforms.     

Fast numerical solution of fredholm integral equations with stationary kernels

A fast recursive matrix method for the numerical solution of Fredholm integral equations with stationary kernels is derived. IfN denotes the number of nodal points, the complexity of the algorithm isO(N ²), which should be compared toO(N ³) for conventional algorithms for solving such problems. The method is related to fast algorithms for inverting Toeplitz matrices.Applications to equations of the first and second kind as well as miscellaneous problems are discussed and illustrated with numerical examples. These show that the theoretical improvement in efficiency is indeed obtained, and that no problems with numerical stability or accuracy are encountered.     

A direct method to solve block banded block Toeplitz systems with non-banded Toeplitz blocks

A fast solution algorithm is proposed for solving block banded block Toeplitz systems with non-banded Toeplitz blocks. The algorithm constructs the circulant transformation of a given Toeplitz system and then by means of the Sherman-Morrison-Woodbury formula transforms its inverse to an inverse of the original matrix. The block circulant matrix with Toeplitz blocks is converted to a block diagonal matrix with Toeplitz blocks, and the resulting Toeplitz systems are solved by means of a fast Toeplitz solver.The computational complexity in the case one uses fast Toeplitz solvers is equal to ξ(m,n,k)=O(mn³)+O(k³n³) flops, there are m block rows and m block columns in the matrix, n is the order of blocks, 2k+1 is the bandwidth. The validity of the approach is illustrated by numerical experiments.     

Inversion of certain extensions of Toeplitz matrices

It is shown how an algorithm for inverting Toeplitz matrices using O(n²) operations can be modified to deal with a certain type of extension, named conjugate-Toeplitz matrices. The block partitioned case is also included.     

Remarks on a displacement-rank inversion method for Toeplitz systems

Comments are made regarding the implementation of a Toeplitz-matrix inversion algorithm described by Bitmead and Anderson in [1]. We show that although the algorithm is asymptotically efficient with O(N(logN)²) operations, it requires a 10⁶×10⁶ matrix to break even with the class of algorithms whose operation count is of the order of O(N²) (as found in [4]).     

Fast Structured Matrix Computations: Tensor Rank and Cohn–Umans Method

We discuss a generalization of the Cohn–Umans method, a potent technique developed for studying the bilinear complexity of matrix multiplication by embedding matrices into an appropriate group algebra. We investigate how the Cohn–Umans method may be used for bilinear operations other than matrix multiplication, with algebras other than group algebras, and we relate it to Strassen’s tensor rank approach, the traditional framework for investigating bilinear complexity. To demonstrate the utility of the generalized method, we apply it to find the fastest algorithms for forming structured matrix–vector product, the basic operation underlying iterative algorithms for structured matrices. The structures we study include Toeplitz, Hankel, circulant, symmetric, skew-symmetric, f-circulant, block Toeplitz–Toeplitz block, triangular Toeplitz matrices, Toeplitz-plus-Hankel, sparse/banded/triangular. Except for the case of skew-symmetric matrices, for which we have only upper bounds, the algorithms derived using the generalized Cohn–Umans method in all other instances are the fastest possible in the sense of having minimum bilinear complexity. We also apply this framework to a few other bilinear operations including matrix–matrix, commutator, simultaneous matrix products, and briefly discuss the relation between tensor nuclear norm and numerical stability.     

Factorizing symmetric indefinite matrices

The LDL^T factorization of a symmetric indefinite matrix, although efficient computationally, may not exist and can be unstable in the presence of round off error. The use of block diagonal 2×2 pivots is attractive, but there are some difficulties in determining an efficient and stable pivot strategy. Previous suggestions have required O(n>³) operations (either multiplications or comparisons) just to implement the pivot strategy. A new strategy is described which in practice only requires O(n²) operations. Indeed, the effort required by this pivot strategy is less than that required when using partial pivoting with an unsymmetric LU factorization, which is the usual way of factorizing indefinite matrices.     

On computing bounds for the least singular value of a triangular matrix

An algorithm for rapid computation of a lower bound for the least singular value of a triangular matrix is presented. It requiresO(N) operations whereN is the order of the matrix, and is based on the Perron-Frobenius theory of non-negative matrices. The input data are the diagonal elements and the off-diagonal elements of maximum modulus in each row.     

Inverses of tridiagonal Toeplitz and periodic matrices with applications to mechanics

The linear algebraic equation Ax = b with tridiagonal coefficient matrix of A is solved by the analytical matrix inversion. An explicit formula is known if A is a Toeplitz matrix. New formulas are presented for the following cases: (1) A is of Toeplitz type except that A(1, 1) and A(n, n) are different from the remaining diagonal elements. (2) A is p-periodic (p > 1), by which is meant that in each of the three bands of A a group of p elements is periodically repeated. (3) The tridiagonal matrix A is composed of periodic submatrices of different periods. In cases (2) and(3) the problem of matrix inversion is reduced to a second-order difference equation with periodic coefficients. The solution is based on Floquet's theorem. It is shown that for p = 1 the formulae found for periodic matrices reduce to special forms valid for Toeplitz matrices. The results are applied to problems of elastostatics and of vibration theory.     

An O(N) algorithm for constructing the solution operator to 2D elliptic boundary value problems in the absence of body loads

The large sparse linear systems arising from the finite element or finite difference discretization of elliptic PDEs can be solved directly via, e.g., nested dissection or multifrontal methods. Such techniques reorder the nodes in the grid to reduce the asymptotic complexity of Gaussian elimination from O(N ²) to O(N ^1.5) for typical problems in two dimensions. It has recently been demonstrated that the complexity can be further reduced to O(N) by exploiting structure in the dense matrices that arise in such computations (using, e.g., \(\mathcal {H}\) -matrix arithmetic). This paper demonstrates that such accelerated nested dissection techniques become particularly effective for boundary value problems without body loads when the solution is sought for several different sets of boundary data, and the solution is required only near the boundary (as happens, e.g., in the computational modeling of scattering problems, or in engineering design of linearly elastic solids). In this case, a modified version of the accelerated nested dissection scheme can execute any solve beyond the first in O(N _boundary) operations, where N _boundary denotes the number of points on the boundary. Typically, N _boundary ～ N ^0.5. Numerical examples demonstrate the effectiveness of the procedure for a broad range of elliptic PDEs that includes both the Laplace and Helmholtz equations.     

New inversion formulas for matrices classified in terms of their distance from Toeplitz matrices

The problem of solving linear equations, or equivalently of inverting matrices, arises in many fields. Efficient recursive algorithms for finding the inverses of Toeplitz or displacement-type matrices have been known for some time. By introducting a way of characterizing matrices in terms of their “distance” from being Toeplitz, a natural extension of these algorithms is obtained. Several new inversion formulas for the representation of the inverse of non-Toeplitz matrices are also presented.     

An implicit QR algorithm for symmetric semiseparable matrices

The QR algorithm is one of the classical methods to compute the eigendecomposition of a matrix. If it is applied on a dense n × n matrix, this algorithm requires O(n³) operations per iteration step. To reduce this complexity for a symmetric matrix to O(n), the original matrix is first reduced to tridiagonal form using orthogonal similarity transformations. In the report (Report TW360, May 2003) a reduction from a symmetric matrix into a similar semiseparable one is described. In this paper a QR algorithm to compute the eigenvalues of semiseparable matrices is designed where each iteration step requires O(n) operations. Hence, combined with the reduction to semiseparable form, the eigenvalues of symmetric matrices can be computed via intermediate semiseparable matrices, instead of tridiagonal ones. The eigenvectors of the intermediate semiseparable matrix will be computed by applying inverse iteration to this matrix. This will be achieved by using an O(n) system solver, for semiseparable matrices. A combination of the previous steps leads to an algorithm for computing the eigenvalue decompositions of semiseparable matrices. Combined with the reduction of a symmetric matrix towards semiseparable form, this algorithm can also be used to calculate the eigenvalue decomposition of symmetric matrices. The presented algorithm has the same order of complexity as the tridiagonal approach, but has larger lower order terms. Numerical experiments illustrate the complexity and the numerical accuracy of the proposed method. Copyright © 2005 John Wiley & Sons, Ltd.     

Essentially optimal computation of the inverse of generic polynomial matrices

We present an inversion algorithm for nonsingular n×n matrices whose entries are degree d polynomials over a field. The algorithm is deterministic and, when n is a power of two, requires O^∼(n³d) field operations for a generic input; the soft-O notation O^∼ indicates some missing log(nd) factors. Up to such logarithmic factors, this asymptotic complexity is of the same order as the number of distinct field elements necessary to represent the inverse matrix.     

Contractions with rank one defect operators and truncated CMV matrices

The main issue we address in the present paper are the new models for completely nonunitary contractions with rank one defect operators acting on some Hilbert space of dimension N?∞. These models complement nicely the well-known models of Livšic and Sz.-Nagy-Foias. We show that each such operator acting on some finite-dimensional (respectively, separable infinite-dimensional Hilbert space) is unitarily equivalent to some finite (respectively semi-infinite) truncated CMV matrix obtained from the “full” CMV matrix by deleting the first row and the first column, and acting in CN (respectively ?²(N)). This result can be viewed as a nonunitary version of the famous characterization of unitary operators with a simple spectrum due to Cantero, Moral and Velázquez, as well as an analog for contraction operators of the result from [Yu. Arlinski?, E. Tsekanovski?, Non-self-adjoint Jacobi matrices with a rank-one imaginary part, J. Funct. Anal. 241 (2006) 383-438] concerning dissipative non-self-adjoint operators with a rank one imaginary part. It is shown that another functional model for contractions with rank one defect operators takes the form of the compression f(ζ)→PK(ζf(ζ)) on the Hilbert space L²(T,dμ) with a probability measure μ onto the subspace K=L²(T,dμ)?C. The relationship between characteristic functions of sub-matrices of the truncated CMV matrix with rank one defect operators and the corresponding Schur iterates is established. We develop direct and inverse spectral analysis for finite and semi-infinite truncated CMV matrices. In particular, we study the problem of reconstruction of such matrices from their spectrum or the mixed spectral data involving Schur parameters. It is pointed out that if the mixed spectral data contains zero eigenvalue, then no solution, unique solution or infinitely many solutions may occur in the inverse problem for truncated CMV matrices. The uniqueness theorem for recovered truncated CMV matrix from the given mixed spectral data is established. In this part the paper is closely related to the results of Hochstadt and Gesztesy-Simon obtained for finite self-adjoint Jacobi matrices.     

A fast method to block-diagonalize a Hankel matrix

In this paper, we consider an approximate block diagonalization algorithm of an n×n real Hankel matrix in which the successive transformation matrices are upper triangular Toeplitz matrices, and propose a new fast approach to compute the factorization in O(n ²) operations. This method consists on using the revised Bini method (Lin et al., Theor Comp Sci 315: 511–523, 2004). To motivate our approach, we also propose an approximate factorization variant of the customary fast method based on Schur complementation adapted to the n×n real Hankel matrix. All algorithms have been implemented in Matlab and numerical results are included to illustrate the effectiveness of our approach.     

On neighbouring matrices with quadratic elementary divisors

Summary Algorithms are presented which compute theQR factorization of an order-n Toeplitz matrix inO(n ²) operations. The first algorithm computes onlyR explicitly, and the second computes bothQ andR. The algorithms are derived from a well-known procedure for performing the rank-1 update ofQR factors, using the shift-invariance property of the Toeplitz matrix. The algorithms can be used to solve the Toeplitz least-squares problem, and can be modified to solve Toeplitz systems inO(n) space.     

Generating polynomials for matrices with Toeplitz or conjugate-Toeplitz inverses

A complex matrix A = [a_ij] has been called conjugate-Toeplitz if a_ij = c^i?1(a_i?j), where c( ) denotes conjugation. A necessary and sufficient condition is derived for a matrix H to have a conjugate-Toeplitz inverse. The elements of H are generated from the coefficients of certain polynomials. The result is simplified either when H is to have a Toeplitz inverse, or when H is banded and is to have a conjugate-Toeplitz inverse. Each of these consequences of the main theorem is a different generalization of a previous result on banded matrices having a Toeplitz inverse.