New Approach for the Inversion of Structured Matrices via Newton ’ s Iteration

Newton’s iteration is a fundamental tool for numerical solutions of systems of equations. The well-known iteration ( ) 1 2 , 0 i i i X X I MX i + = − ≥ rapidly refines a crude initial approximation 0 X to the inverse of a general nonsingular matrix. In this paper, we will extend and apply this method to n n × structured matrices M , in which matrix multiplication has a lower computational cost. These matrices can be represented by their short generators which allow faster computations based on the displacement operators tool. However, the length of the generators is tend to grow and the iterations do not preserve matrix structure. So, the main goal is to control the growth of the length of the short displacement generators so that we can operate with matrices of low rank and carry out the computations much faster. In order to achieve our goal, we will compress the computed approximations to the inverse to yield a superfast algorithm. We will describe two different compression techniques based on the SVD and substitution and we will analyze these approaches. Our main algorithm can be applied to more general classes of structured matrices.

× structured matrices M , in which matrix multiplication has a lower computational cost.These matrices can be represented by their short generators which allow faster computations based on the displacement operators tool.However, the length of the generators is tend to grow and the iterations do not preserve matrix structure.So, the main goal is to control the growth of the length of the short displacement generators so that we can operate with matrices of low rank and carry out the computations much faster.In order to achieve our goal, we will compress the computed approximations to the inverse to yield a superfast algorithm.We will describe two different compression techniques based on the SVD and substitution and we will analyze these approaches.Our main algorithm can be applied to more general classes of structured matrices.

Introduction
Frequently, matrices encountered in practical computations that have some special structures which can be exploited to simplify the computations.In particular, computations with dense structured matrices are ubiquitous in sciences, communications and engineering.Exploitation of structure enables dramatic acceleration of the computations and a major decrease in memory space, but sometimes it also leads to numerical stability problems.
The best-known classes of structured matrices are Toeplitz and Hankel matrices, but Cauchy and Vandermonde types are also quite popular.The computations with such matrices are widely applied in the areas of algebraic coding, control, signal processing, solution of partial differential equations and algebraic computing.For example, Toeplitz matrices arise in some major signal processing computations and Cauchy matrices appear in the study of integral equations and conformal mappings.The complexity of computations with n n × dense structured matrices dramatically decreases in comparison with the general n n × matrices, that is, from the order of 2  n words of storage space and n ω arithmetic operations (ops) with 2.37 3 ω < ≤ in the best algorithms, to ( ) O n words of storage space and to ( ) n ops (and sometimes to ( ) log O n n ops) (Table 1).The displacement rank r for an m n × Toeplitz and Hankel matrix is at most 2. A matrix is Toeplitz-like or Hankel-like if r is small or bounded by a small constant independent of m and n .Those matrices can be represented by ( ) m n r + entries of its short displacement generators instead of mn entries which leads to more efficient storage of the entries in computer memory and much faster computations [1].In this paper, we will focus our study on Newton's iteration for computing the inverse of a structured matrix and we will analyze the resulting algorithms.This iteration is well-known for its numerically stability and faster convergence.Newton's iteration ( ) for matrix inversion was initially proposed by Schultz in 1933 and studied by Pan and others.The authors of [2] described quadratically convergent algorithms for the refinement of rough initial approximations to the inverse of Toeplitz and Toeplitz-like matrices and to the solutions of Toeplitz and Toeplitz linear systems of equations.
The algorithms are based on the inversion formulas for Toeplitz matrices.The Cauchy-like case was studied in [3].Since those matrices can be represented by their short generators, which allow faster com-putations based on the displacement operators tool, we can employ Newton's iteration to compute the inverse of the input structured matrices.However, Newton's iteration destroys the structure of the structured matrices when used.Therefore, in Section 5 we will study two compression techniques in details, namely, truncation of the smallest singular values of the displacement and the substitution of approximate inverses for the inverse matrix into its displacement expression to preserve the structure.Finally, each iteration step is reduced to two multiplications of structured matrices and each iteration is performed by using nearly linear time and optimal memory space which leads to superfast algorithm, especially if the input matrix is a well-conditioned structured matrix.Our main algorithm of Section 6 can be applied to more general classes of structured matrices.We will support our analysis with numerical experiments with Toeplitz matrices in Section 7.

Definition of Structured Matrices
for every pair of its entries , i j t and for every pair of its entries , i j h and 1, 1 i j h − + .For a given vector ( ) , the matrix ( ) is called a Vandermonde matrix.Given two vectors s and t such that i j s t ≠ for all i and j , the n n × matrix ( ) is a Cauchy (generalized Hilbert) matrix where ( ) , ( ) , and ( ) respectively.We refer the reader to [4] and [5] for more details on the basic properties and features of structured matrices.

Displacement Representation of Structured Matrices
The concept of displacement operators and displacement rank, initially introduced by Kailath, Kung, and Morf in 1979, and studied by other authors such as Bini and Pan, is a powerful tool for dealing with matrices with structure.The displacement rank approach was intended for more restricted use when it was initially introduced in [6] (also [7]), namely, to measure how "close" to Toeplitz a given matrix is.Then the idea turned out to be even deeper and more powerful, thus it was developed, generalized and extended to other structured matrices.In this paper, we consider the most general and modern interpretation of the displacement of a matrix M (see Definition 2).The main idea is, for a given structured matrix M , we need to find an operator L that transforms the matrix M into a low rank matrix ( ) L M , such that one can easily recover M from its image ( ) L M and operate with low rank matrices instead.Such operators that shift and scale the entries of the structured matrices turn out to be appropriate tools for introducing and defining the matrices of Toeplitz-like, Hankel-like, Vandermonde-like, and Cauchy-like types which we will introduce in Definition 4. Definition 2 For any fixed field F (say the complex field C ) and a fixed pair { } , A B of operator matrices, we define the linear displacement operators and Stein type The image ( ) L M of the operator L is called the displacement of the matrix M .The operators of Sylvester and Stein types can be transformed easily into one another if at least one of the two associated operator matrices is nonsingular.This leads to the following theorem.
if the operator matrix A is nonsingular, and , and The operator matrices that we will use are the unit f-circulant matrix ( ) , and the diagonal matrix . Note that ( ) is unit lower triangular Toeplitz matrix.
. Then if we apply the displacement operator of Sylvester type (2) to Toeplitz matrix T using the shift operator matrices 1 Z and 0 Z , we will get: This operator displaces the entries of the matrix T so that the image can be represented by the first row and the last column.Furthermore, the image is a rank 2 matrix Notice that the operator matrices vary depending on the structure of the matrix M .We use shift operator matrices (such as f Z and T f Z for any scalar f ) for the structures of Toeplitz and Hankel type, scaling operator matrices (such as ( ) D x ) for the Cauchy structure, and shift and scaling operator matrices for the structure of Vandermonde type.One can restrict the choices of f to 0 and 1.However, other operator matrices may be used with other matrix structures (see Table 2).
Definition 3 For an m n × structured matrix M and an associated operator is called the L -rank or the displacement rank of the matrix M .The constant r in Definition 3 usually remains relatively small as the size of M grows large, that is, r is small relative to ( ) min , m n , say ( ) as m and n grow large.In this case, we call the matrix L -like and having structure of type L .Now let us recall the linear operators L of ( 2) and ( 3) that map a matrix M into its displacements: and where A and B are operator matrices of Table 2, ( ) trices.Then the matrix pair ( ) and l is small relative to m and n , then M is said to have L -structure or to be an L -structured matrix.
The displacement ( ) (as a matrix of small rank) can be represented with a short generator defined by only a small number of parameters.For instance, the singular value decomposition (SVD) of the matrix ( ) L M of size m n × gives us its orthogonal generator that consists of two generator matrices G and H having orthogonal columns and of sizes m r × and n r × respectively, that is a total of ( ) , G H is called an L -generator (or displacement generator) of length r for M .It is also a generator of length r for ( ) L M .Note that the pair ( ) , G H is nonunique for a fixed ( ) L M .In particular, for Toeplitz, Hankel, Vandermonde, and Cauchy matrices, the length r of the associated generators (4) and ( 5) is as small as 1 or 2 for some appropriate choices of the operator matrices A and B of Table 2.The four classes of structured matrices are naturally extended to Toeplitz-like, Hankel-like, Vandermonde-like, and Cauchylike matrices, for which r is bounded by a fixed (not too large) constant.
Definition 4 An m n ) where for a fixed pair of scalars e and f .A matrix M is Hankel-like if MJ or JM is Toeplitz- like where ( ) is the reflection matrix.
Therefore, the problems of solving Toeplitz-like and Hankel-like linear systems are reduced to each other.

Remark 1
In the case where the operator L is associated with Toeplitz, Hankel, Vandermonde, or Cauchy matrices M , we call an L -like matrix Toeplitz-like, Hankel-like, Vandermonde-like, or Cauchy-like matrix, respectively, and we say that the matrix has the structure of Toeplitz, Hankel, Vandermonde, or Cauchy type respectively.To take advantage of the matrix structure, we are going to COMPRESS the structured input matrix M via its short L-generators based on Theorem 2 below, then we will OPERATE with L-generators rather than with the matrix itself, and finally RECOVER the output from the computed short L-generators.
Definition 6 We define the norm of a nonsingular linear operator L and its inverse 1 L − as follows: ∞ , or F , and where the supremum is taken over all matrices M having a positive L -rank of at most r .In this case the condition number κ of the operator L is defined as Theorem 2 [8] Let G and H be a pair of n l × matrices, ( ) , where i g and i h denote the i -th columns of G and H respectively.Let ( ) . Then we have one of the following Toeplitz type matrices: ( ) ( ) ( ) ( ) Theorem 3 Let M be a nonsingular matrix, then we have the following: 1) 2) , if A is a nonsingular matrix. 3)

Newton's Iteration for General Matrix Inversion
Given an n n × nonsingular matrix A and an initial approximation 0 X to its inverse 1 A − that satisfies the bound 0 1 I X A ρ = − < for some positive ρ and some fixed matrix norm, Newton's iteration (1) can be written in the form and rapidly improves the initial approximation 0 X to the inverse.0 X can be supplied by either some computational methods such as preconditioning conjugate gradient method or by an exact solution algorithm that requires refinement due to rounding errors.
Lemma 1 Newton's iteration (6) converges quadratically to the inverse of an n n × nonsingular matrix A .
Proof.From (6) , since We define the error matrices and hence the norm of the residual matrices i R is bounded above by 2 i ρ .This proves that Newton's iteration converges quadratically.
is guaranteed in (6), where x     stands for the smallest integer not exceeded by x .(7) also implies that so that the matrix i X approximates 1 A − with a relative error norm less than  .
Remark 2 Each step of Newton's iteration (6) involves the multiplication of A by i X , the addition of 2 to all of the diagonal entries of the product, which takes n additions, and the pre-multiplication of the resulting matrix by i X .The most costly steps are the two matrix multiplications; each step takes 3 n multiplications and .Then by quadratic convergence of Newton's iteration, we only need 5 iterative steps of (6) to compute 5 X that satisfies the bound ( ) and hence If we consider any linearly convergent scheme such as Jacobi's or Gauss-Seidels's iteration ([9], p. 611) for the same input matrix, the norm of the error matrix of the computed approximation to

Newton-Structured Iteration
In Section 4, we showed that for a general input matrix A , the computational cost of iteration ( 6) is quite high because matrix products must be computed at every step (see Remark 2).However, for structured matrices M , matrix multiplication has a lower computational cost of ( ) .In this case, Newton's iteration can be efficiently applied and rapidly improves a rough initial approximation to the inverse, but also destroys the structure of the four classes of structured matrices.Therefore, Newton's iteration has to be modified to preserve the initial displacement structure of the input structure matrix during the iteration and maintain matrix structure throughout the computation.The main idea here is to control the growth of the length of the short displacement generators so that we can operate with matrices of low rank and carry out the computations much faster.This can be done based on one of the following two techniques.The first technique is to periodically chop off the components corresponding to the smallest singular values in the SVDs of the displacement matrices defined by their generators.This technique can be applied to a Toeplitz-like input matrix.For a Cauchy-like input matrix C , there is no need to involve the SVD due to the availability of the inverse formula for 1 C − [10].The second technique is to control the growth of the length of the generators.This can be done by substituting the approximate inverses for the inverse matrix into its displacement.Here we assume that the matrices involved can be expressed via their displacement generators.Now, let us assume that the matrices M and 0 X have short generators ∇ -generators are used.In this paper, we restrict our presentation to ∇ -generators (see Theorem 1).

SVD-Based Compression of the Generators
For a fixed operator L and a matrix M , one can choose the orthogonal (SVD-based) L -generator matrices to obtain better numerical stability [11].
For real matrices, U is orthogonal iff is orthogonal, then i u form an orthogonal basis for m C .

Definition 8 Every m n
× matrix W has its SVD (singular value decomposition): ( ) .Based on the SVD of a matrix W, its orthogonal generator of the minimum length is defined by However, one can use an alternative diagonal scaling, in this case (11) can be written as Note that if ( ) for a matrix M , then the pair ( ) " produces the above results.Clearly, 1 3 σ = is the largest singular value and 2 1 σ = is the smallest singular value.

Outline of the Techniques
By Theorem 6, we assume that ∇ -displacement rank r for large i .This fact motivates the following approach to decrease the computational cost of the iteration in (6): We will replace i X in (6) with a nearby matrix i Y having , B A ∇ -displacement rank of at most r and then restart the iteration with i Y instead of i X .The advantage of this modification is the decrease of the computational cost at the i -th Newton step and at all subsequent steps.However, the disadvantage is a possible deterioration of the approximation, since we do not know 1 M − , and the transition from i X to i Y may occur in a wrong direction.But since i X and i Y are supposed to lie close to 1 M − , hence they both lie close to each other.This observation enables us to bound the deterioration of the approximation relative to the current approximation error, so that the quadratic improvement in the error bound at the next step of Newton's iteration will immediately compensate us for such a deterioration.Furthermore, the transition from i X to a nearby matrix i Y of a small displacement rank can be done at a low computational cost.In the next main algorithm we describe this approach for Sylvester type operators only.The same results can be extended to Stein operators L using ∇ -generators of length at most 3r for the matrices ( ) and then apply the compression subalgorithm 1 C or 2 C (transition from ∇ -generators of length at most r for the matrices , where 2,i e as in (13).
Proof.Using the estimate from ([9], p. 442), that is, for any two matrices A and E , ( and conclude that ( . Now, we use (2) of Lemma 2 and deduce , for j r > .If we combine the latest inequality with equation (1) of Lemma 2 Now let us use Lemma 2 and Lemma 3 to prove the bound in the next theorem.

Theorem 7 If
, B A ∇ is a nonsingular linear operator and 2,i e is a positive scalar defined in (13), then the Proof.By Definition 6 replacing L by , and a , B A ∇ -generator ( ) (the computation is not costly since it is performed with and since i k is small, see [11]).Set to zero the diagonal entries 1 , ,  Next, we will present the second compression technique based on substitution which does not depend on the SVD.By Theorem 4 we conclude that ( ) ( ) , and by Theorem 3 we also have ( ) ( ) ( ) if A is nonsingular.Clearly this expresses the displacement of 1 M − via the displacement of M , the operator matrices A and B , and the product of the inverse with 2r specified vectors, where ( ) ( ) , let us substitute i X for 1 M − on the right hand side and define a short ( ) , of length at most 3r for a matrix 1 i X + of Newton's iteration (6). , of length at most r for a matrix 1 i Y + satisfying the equation ( ) where ( ) 1 ˆi e + of ( 14), and for i e of (13).

COMPUTATIONS:
Compute the matrix products ( ) and ( ) Now the pair ( ) , ∇ -generator of length at most r for a matrix Here we assume the operator , B A L = ∇ , then the matrix 1 i Y + is uniquely defined by its L- W + and thus (18) and (19) completely define Newton's process without involving the matrix . This subalgorithm achieves compression because ( ) ( ) and preserves approximation.Since , we deduce from Theorem 3 that ( ) . The computation of the n r × matrices of ( 16) is reduced to multi- plication of the matrix , and for .Then we can write ( ) ( ) Now, take the norm of j E : Simplifying the latest inequality yield the following bound: ( ) ( ) If we use (20) for 1 j i = + , we get ( ) arithmetic operations.However, it is actually applied to structured matrices with smaller rank, so it is not really expensive.In addition, CLAPACK is one of the faster SVD algorithms that can be applied on real and complex matrices with less time complexity (see [11]).The tests results are summarized in the following example.
Example 6 N refers to number of Newton's steps, n is the matrix size, K is number of times we ran the test, ε is a selected bound, and ( ) ( ) ( ) ( ) (the first table of Table 3) or 2 − (the second table of Table 3) (see Figure 1).
2) The matrices (third table represents symmetric positive definite and the fourth table represents symmetric positive indefinite of Table 3) (see Figure 1).
3) Randomly generated Toeplitz matrices, whose entries are chosen randomly from the interval [ ) 0,1 and are uniformly distributed with mean 0 and standard deviation of 1 (see Table 4 and Figure 2).

rapidly refines a crude initial approximation 0 X
to the inverse of a general nonsingular matrix.In this paper, we will extend and apply this method to n n Hankel, Vandermonde, and Cauchy matrices respectively that can be represented by the 5-dim vector ( ) pair of general n n × input matrices A and 0 X .Of course, this cost will be reduced if the input matrix is structured (see Section 5).Example 3 (Newton's iteration versus linear convergent schemes.)Let 1

1 A
− decreases by factor 2 in each iteration step.The error norm decreases to below1 2 i − −in i iteration steps, so this will take 19 iteration steps to ensure the bound(7) for the same

r and 2 r,
for the two input matrices M and 0 X respectively.So if M and 0 X are structured matrices, say Toeplitz or Toeplitz-like, the computations required for the Newton's iteration can be carried out more efficiently.In general, for the four classes of structured matrices of Definition 1, we usually associate various linear displacement operators L of Sylvester type ( ) by operating with short ∆ -generators of the matrices M , i X , and i MX (or i X M ), we will have to modify Newton's iteration and use one of the above two techniques to control the length of a , B A ∆ -generator of i X , which seems to triple at every iterative step of (6).The same technique applies in the case


are singular values of the matrix W, U * and V * denote the Hermitian (conjugate) transpose of the orthogonal matrices U and V , respectively.That is, U and V are real matrices, then T and

Remark 3
In numerical computation of the SVD with a fixed positive tolerance ε to the output errors caused by rounding ( ε is the output error bound which is also an upper bound on all the output errors of the

Theorem 1 Y to the inverse 1 M 1 C or 2 C∇∇
with slight modifications to the technique and the theory.Algorithm 1 (Newton's Iteration Using Sylverster Type Operators) INPUT: An integer 0 r > , two matrices A and B defining the operator , length at most r , a sufficiently close initial approximation 0 − given with its , B A ∇ -generator of length at most r , an upper bound λ on the number of Newton's iteration steps, and a compression subalgorithm for the transition from a , B A∇ -generator of length at most 3r for an n n -generator of length at most r for a nearby matrix.-generator of length at most r for a matrix Y λ approximating 1

of Definition 8 .
For all j r > we have the proof. Subalgorithm 1 ( 1 C : Compression of Displacement by Truncation of their Smallest Singular Values.)INPUT: An integer 0 r > , two matrices A and B defining the operator , A B ∇ for an n n × nonsingular structured matrices M , computational cost of performing Subalgorithm 1 is dominated by the computation of the SVD of the displacement which requires ( )

2 C:∇
Compression of the displacement by substitution.)INPUT: An integer 0 r > , a pair of n n × operator matrices A and B defining the operator , -generator of length r for a nonsingular n n × matrix M where

 Theorem 9
For

Algorithm 1 balgoritm 1 .
was tested numerically for n n × Toeplitz input matrices M by applying the compression Su- We use Matlab 8.3.0 on a Window 7 machine to run our tests.The Matlab SVD subroutine was used to obtain the singular value decomposition.In general, the complexity of the SVD on an n n

1 )
-th condition number of the matrix M .We restricted our numerical experiments to the following classes of n n × Toeplitz matrices: Real symmetric tridiagonal Toeplitz matrices: ( )

For a symmetric positive
definite matrix M, the initial matrix used, whereas for unsymmetric and symmetric indefinite Toeplitz matrices M, Newton's iteration was initialized with

Table 1 .
Parameters and flops (floating point arithmetic operations) count.
generator of minimum length for the matrix M .
generator of length at most r for a matrix i Y such that

Table 3 .
First top two tables represent Tridiagonal Matrices of Class 1 and the last two tables represent matrices of Class 2.

Table 4 .
Random Matrices of Class 3.