# A Lanczos approach to the inverse square root of a large and sparse matrix

## 1 Introduction

The computation of the inverse square root of a matrix is a special problem in scientific computing. It is related to the matrix sign and polar decomposition .

One may define the matrix sign by:

$$\text{sign}(A)=A(A^2)^{\frac{1}{2}}$$
(1)   
where $`A`$ is a complex matrix with no pure imaginary eigenvalues.

In polar coordinates, a complex number $`z=x+iy`$, is represented by

$$z=|z|e^{i\varphi },\varphi =\mathrm{arctan}\frac{y}{x}$$
(2)   
In analogy, the polar decomposition of a matrix $`A`$ is defined by:

$$A=V(A^{}A)^{\frac{1}{2}},V^1=V^{}$$
(3)   
where $`V`$ is the polar part and the second factor corresponds to the absolute value of $`A`$.

The mathematical literature invloving the matrix sign function traces back to 1971 when it was used to solve the Lyapunov and algebraic Riccati equations .

In computational physics one may face a similar problem when dealing with Monte Carlo simulations of fermion systems, the so-called sign problem . In this case the integration measure is proportional to the determinant of a matrix and the polar decomposition may be helpful to monitor the sign of the determinant.

The example brought in this paper comes from the recent progress in formulating Quantum Chromodynamics (QCD) on a lattice with exact chiral symmetry .

In continuum, the massless Dirac propagator $`D_{cont}`$ is chirally symmetric, i.e.

$$\gamma _5D_{cont}+D_{cont}\gamma _5=0$$
(4)   
On a regular lattice with spacing $`a`$ the symmetry is suppressed according to the Ginsparg-Wilson relation: :

$$\gamma _5D+D\gamma _5=aD\gamma _5D$$
(5)   
where $`D`$ is the lattice Dirac operator.

An explicit example of a Dirac operator obeying this relation is the so-called overlap operator :

$$aD=1A(A^{}A)^{1/2},A=MaD_W$$
(6)   
where $`M`$ is a shift parameter in the range $`(0,2)`$, which I have fixed at one.

$`D_W`$ is the Wilson operator,

$$D_W=\underset{\mu =1}{\overset{4}{}}\gamma _\mu _\mu \frac{a}{2}\underset{\mu =1}{\overset{4}{}}\mathrm{\Delta }_\mu $$
(7)   
which is a nearest-neighbors discretization of the continuum Dirac operator (it violates the Ginsparg-Wilson relation). $`_\mu `$ and $`\mathrm{\Delta }_\mu `$ are first and second order covariant differences given by:

$`(_\mu \psi )_i=\frac{1}{2a}(U_{\mu ,i}\psi _{i+\widehat{\mu }}U_{\mu ,i\widehat{\mu }}^{}\psi _{i\widehat{\mu }})`$

$`(\mathrm{\Delta }_\mu \psi )_i=\frac{1}{a^2}(U_{\mu ,i}\psi _{i+\widehat{\mu }}+U_{\mu ,i\widehat{\mu }}^{}\psi _{i\widehat{\mu }}2\psi _i)`$

where $`\psi _i`$ is a fermion field at the lattice site $`i`$ and $`U_{\mu ,i}`$ an $`SU(3)`$ lattice gauge filed associated with the oriented link $`(i,i+\widehat{\mu })`$. These are unitary 3 by 3 complex matrices with determinant one. A set of such matrices forms a lattice gauge “configuration”.

$`\gamma _\mu ,\mu =1,\mathrm{},5`$ are $`4\times 4`$ Dirac matrices which anticommute with each-other.

Therefore, if there are $`N`$ lattice points in four dimensions, the matrix $`A`$ is of order $`12N`$. A restive symmetry of the matrix $`A`$ that comes from the continuum is the so called $`\gamma _5symmetry`$ which is the Hermiticity of the $`\gamma _5A`$ operator.

By definition, computation of $`D`$ involves the inverse square root of a matrix. This is a non-trivial task for large matrices. Therefore several algorithms have been proposed by lattice QCD physicists .

All these methods rely on matrix-vector multiplications with the sparse Wilson matrix $`D_W`$, being therefore feasible for large simulations.

In fact, methods that approximate the inverse square root by Legendre  and Chebyshev polynomials  need to know apriori the extreme eigenvalues of $`A^{}A`$ to be effective. This requires computational resources of at least one Conjugate Gradients (CG) inversion.

In  the inverse square root is approximated by a rational approximation, which allows an efficient computation via a multi-shift CG iteration. Storage here may be an obstacle which is remedied by a second CG step .

The Pade approximation used by  needs the knowledge of the smallest eigenvalue of $`A^{}A`$. Therefore the method becomes effective only in connection with the $`D`$ inversion .

The method presented earlier by the author  relies on taking exactly the inverse square root from the Ritz values. These are the roots of the Lanczos polynomial approximating the inverse of $`A^{}A`$.

In that work the Lanczos polynomial was constructed by applying the Hermitian operator $`\gamma _5A`$. The latter is indefinite, thereby responsible for observed oscillations in the residual vector norm .

Here I use a Lanczos polynomial on the positive definite matrix $`A^{}A`$. In this case the residual vector norm decreases monotonically and leads to a stable method. This is a crucial property that allows a reliable stopping criterion that I will present here.

The paper is selfcontained: in the next section I will briefly present the Lanczos algorithm and set the notations. In section 3, I use the algorithm to solve linear systems, and in section 4, the computation of the inverse square root is given. The method is tested in section 5 and conclusions are drawn in the end.

## 2 The Lanczos Algorithm

The Lanczos iteration is known to approximate the spectrum of the underlying matrix in an optimal way and, in particular, it can be used to solve linear systems .

Let $`Q_n=[q_1,\mathrm{},q_n]`$ be the set of orthonormal vectors, such that

$$A^{}AQ_n=Q_nT_n+\beta _nq_{n+1}(e_n^{(n)})^T,q_1=\rho _1b,\rho _1=1/b_2$$
(8)   
where $`T_n`$ is a tridiagonal and symmetric matrix, $`b`$ is an arbitrary vector, and $`\beta _n`$ a real and positive constant. $`e_m^{(n)}`$ denotes the unit vector with $`n`$ elements in the direction $`m`$.

By writing down the above decomposition in terms of the vectors $`q_i,i=1,\mathrm{},n`$ and the matrix elements of $`T_n`$, I arrive at a three term recurrence that allows to compute these vectors in increasing order, starting from the vector $`q_1`$. This is the $`LanczosAlgorithm`$:

$$\begin{array}{c}\beta _0=0,\rho _1=1/b_2,q_0=o,q_1=\rho _1b\hfill \\ fori=1,\mathrm{}\hfill \\ v=A^{}Aq_i\hfill \\ \alpha _i=q_i^{}v\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ \beta _i=v_2\hfill \\ if\beta _i<tol,n=i,endfor\hfill \\ q_{i+1}=v/\beta _i\hfill \end{array}$$
(9)   
where $`tol`$ is a tolerance which serves as a stopping condition.

The Lanczos Algorithm constructs a basis for the Krylov subspace :

$$\text{span}\{b,A^{}Ab,\mathrm{},(A^{}A)^{n1}b\}$$
(10)   
If the Algorithm stops after $`n`$ steps, one says that the associated Krylov subspace is invariant.

In the floating point arithmetic, there is a danger that once the Lanczos Algorithm (polynomial) has approximated well some part of the spectrum, the iteration reproduces vectors which are rich in that direction . As a consequence, the orthogonality of the Lanczos vectors is spoiled with an immediate impact on the history of the iteration: if the algorithm would stop after $`n`$ steps in exact arithmetic, in the presence of round off errors the loss of orthogonality would keep the algorithm going on.

## 3 The Lanczos Algorithm for solving $`A^{}Ax=b`$

Here I will use this algorithm to solve linear systems, where the loss of orthogonality will not play a role in the sense that I will use a different stopping condition.

I ask the solution in the form

$$x=Q_ny_n$$
(11)   
By projecting the original system on to the Krylov subspace I get:

$$Q_n^{}A^{}Ax=Q_n^{}b$$
(12)   
By construction, I have

$$b=Q_ne_1^{(n)}/\rho _1,$$
(13)   
Substituting $`x=Q_ny_n`$ and using (8), my task is now to solve the system

$$T_ny_n=e_1^{(n)}/\rho _1$$
(14)   
Therefore the solution is given by

$$x=Q_nT_n^1e_1^{(n)}/\rho _1$$
(15)   
This way using the Lanczos iteration one reduces the size of the matrix to be inverted. Moreover, since $`T_n`$ is tridiagonal, one can compute $`y_n`$ by short recurences.

If I define:

$$r_i=bA^{}Ax_i,q_i=\rho _ir_i,\stackrel{~}{x}_i=\rho _ix_i$$
(16)   
where $`i=1,\mathrm{}`$, it is easy to show that

$$\begin{array}{c}\rho _{i+1}\beta _i+\rho _i\alpha _i+\rho _{i1}\beta _{i1}=0\hfill \\ q_i+\stackrel{~}{x}_{i+1}\beta _i+\stackrel{~}{x}_i\alpha _i+\stackrel{~}{x}_{i1}\beta _{i1}=0\hfill \end{array}$$
(17)   
Therefore the solution can be updated recursively and I have the following Algorithm1 for solving the system $`A^{}Ax=b`$:

$$\begin{array}{c}\beta _0=0,\rho _1=1/b_2,q_0=o,q_1=\rho _1b\hfill \\ fori=1,\mathrm{}\hfill \\ v=A^{}Aq_i\hfill \\ \alpha _i=q_i^{}v\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ \beta _i=v_2\hfill \\ q_{i+1}=v/\beta _i\hfill \\ \stackrel{~}{x}_{i+1}=\frac{q_i+\stackrel{~}{x}_i\alpha _i+\stackrel{~}{x}_{i1}\beta _{i1}}{\beta _i}\hfill \\ \rho _{i+1}=\frac{\rho _i\alpha _i+\rho _{i1}\beta _{i1}}{\beta _i}\hfill \\ r_{i+1}:=q_{i+1}/\rho _{i+1}\hfill \\ x_{i+1}:=y_{i+1}/\rho _{i+1}\hfill \\ if\frac{1}{|\rho _{i+1}|}<tol,n=i,endfor\hfill \end{array}$$
(18)   
## 4 The Lanczos Algorithm for solving $`(A^{}A)^{1/2}x=b`$

Now I would like to compute $`x=(A^{}A)^{1/2}b`$ and still use the Lanczos Algorithm. In order to do so I make the following observations:

Let $`(A^{}A)^{1/2}`$ be expressed by a matrix-valued function, for example the integral formula :

$$(A^{}A)^{1/2}=\frac{2}{\pi }_0^{\mathrm{}}𝑑t(t^2+A^{}A)^1$$
(19)   
From the previous section, I use the Lanczos Algorithm to compute

$$(A^{}A)^1b=Q_nT_n^1e_1^{(n)}/\rho _1$$
(20)   
It is easy to show that the Lanczos Algorithm is shift-invariant. i.e. if the matrix $`A^{}A`$ is shifted by a constant say, $`t^2`$, the Lanczos vectors remain invariant. Moreover, the corresponding Lanczos matrix is shifted by the same amount.

This property allows one to solve the system $`(t^2+A^{}A)x=b`$ by using the same Lanczos iteration as before. Since the matrix $`(t^2+A^{}A)`$ is better conditioned than $`A^{}A`$, it can be concluded that once the original system is solved, the shifted one is solved too. Therefore I have:

$$(t^2+A^{}A)^1b=Q_n(t^2+T_n)^1e_1^{(n)}/\rho _1$$
(21)   
Using the above integral formula and putting everything together, I get:

$$x=(A^{}A)^{1/2}b=Q_nT_n^{1/2}e_1^{(n)}/\rho _1$$
(22)   
There are some remarks to be made here:

a) As before, by applying the Lanczos iteration on $`A^{}A`$, the problem of computing $`(A^{}A)^{1/2}b`$ reduces to the problem of computing $`y_n=T_n^{1/2}e_1^{(n)}/\rho _1`$ which is typically a much smaller problem than the original one. But since $`T_n^{1/2}`$ is full, $`y_n`$ cannot be computed by short recurences. It can be computed for example by using the full decomposition of $`T_n`$ in its eigenvalues and eigenvectors; in fact this is the method I have employed too, for its compactness and the small overhead for moderate $`n`$.

b) The method is not optimal, as it would have been, if one would have applied it directly for the matrix $`(A^{}A)^{1/2}`$. By using $`A^{}A`$ the condition is squared, and one looses a factor of two compared to the theoretical case!

c) From the derivation above, it can be concluded that the system $`(A^{}A)^{1/2}x=b`$ is solved at the same time as the system $`A^{}Ax=b`$.

d) To implement the result (22), I first construct the Lanczos matrix and then compute

$$y_n=T_n^{1/2}e_1^{(n)}/\rho _1$$
(23)   
To compute $`x=Q_ny_n`$, I repeat the Lanczos iteration. I save the scalar products, though it is not necessary.

Therefore I have the following Algorithm2 for solving the system $`(A^{}A)^{1/2}x=b`$:

$$\begin{array}{c}\beta _0=0,\rho _1=1/b_2,q_0=o,q_1=\rho _1b\hfill \\ fori=1,\mathrm{}\hfill \\ v=A^{}Aq_i\hfill \\ \alpha _i=q_i^{}v\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ \beta _i=v_2\hfill \\ q_{i+1}=v/\beta _i\hfill \\ \rho _{i+1}=\frac{\rho _i\alpha _i+\rho _{i1}\beta _{i1}}{\beta _i}\hfill \\ if\frac{1}{|\rho _{i+1}|}<tol,n=i,endfor\hfill \\ \\ Set(T_n)_{i,i}=\alpha _i,(T_n)_{i+1,i}=(T_n)_{i,i+1}=\beta _i,otherwise(T_n)_{i,j}=0\hfill \\ y_n=T_n^{1/2}e_1^{(n)}/\rho _1=U_n\mathrm{\Lambda }_n^{1/2}U_n^Te_1^{(n)}/\rho _1\hfill \\ \\ q_0=o,q_1=\rho _1b,x_0=o\hfill \\ fori=1,\mathrm{},n\hfill \\ x_i=x_{i1}+q_iy_n^{(i)}\hfill \\ v=A^{}Aq_i\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ q_{i+1}=v/\beta _i\hfill \end{array}$$
(24)   
where by $`o`$ I denote a vector with zero entries and $`U_n,\mathrm{\Lambda }_n`$ the matrices of the eigenvectors and eigenvalues of $`T_n`$. Note that there are only four large vectors necessary to store: $`q_{i1},q_i,v,x_i`$.

## 5 Testing the method

I propose a simple test: I solve the system $`A^{}Ax=b`$ by applying twice the $`Algorithm2`$, i.e. I solve the linear systems

$$(A^{}A)^{1/2}z=b,(A^{}A)^{1/2}x=z$$
(25)   
in the above order. For each approximation $`x_i`$, I compute the residual vector

$$r_i=bA^{}Ax_i$$
(26)   
The method is tested for a SU(3) configuration at $`\beta =6.0`$ on a $`8^316`$ lattice, corresponding to an order $`98304`$ complex matrix $`A`$.

In Fig.1 I show the norm of the residual vector decreasing monotonically. The stagnation of $`r_i_2`$ for small values of $`tol`$ may come from the accumulation of round off error in the $`64`$-bit precision arithmetic used here.

This example shows that the tolerance line is above the residual norm line, which confirms the expectation that $`tol`$ is a good stopping condition of the $`Algorithm2`$.

## 6 Conclusions

I have presented a Lanczos method to compute the inverse square root of a large and sparse positive definite matrix.

The method is characterized by a residual vector norm that decreases monotonically and a consistent stopping condition. This stability should be compared with a similar method presented earlier by the author , where the underlying Hermitian but indefinite matrix $`\gamma _5A`$ led to appreciable instabilities in the norm of the residual vector.

In terms of complexity this algorithm requires less operations for the same accuracy than its indefinite matrix counterpart. This property is guaranteed by the monotonicity of the residual vector norm. Nontheless, the bulk of the work remains the same.

With the improvement in store the method is complete.

It shares with methods presented in  the same underlying Lanczos polynomial. As it is wellknown  CG and Lanczos methods for solving a linear system produce the same results in exact arithmetic. In fact, CG derives from the Lanczos algorithm by solving the coupled two-term recurences of CG for a single three-term recurence of Lanczos. However, the coupled two-term recurences of CG accumulate less round off. This makes CG preferable for ill-conditioned problems.

There are two main differences between the method presented here and those in :

a) Since CG and Lanczos are equivalent, they produce the same Lanczos matrix. Therefore, any function of $`A^{}A`$ translates for both algorithms into a function of $`T_n`$ (given the basis of Lanczos vectors). The latter function translates into a function of the Ritz values, the eigenvalues of $`T_n`$. That is, whenever the methods of papers  try to approximate the inverse square root of $`A^{}A`$, the underlying CG algorithm shifts this function to the Ritz values. It is clear now that if I take the inverse square root from the Ritz values exactly, I don’t have any approximation error. This is done in $`Algorithm2`$.

b) $`Algorithm2`$ sets no limits in the amount of memory required, whereas the multi-shift CG needs to store as many vectors as the number of shifts. For high accuracy approximations the multi-shift CG is not practical. However, one may lift this limit in expense of a second CG iteration (two-step CG) . Therefore $`Algorithm2`$ and the two-step CG have the same iteration workload, with $`Algorithm2`$ computing exactly the inverse square root.

Additionally, $`Algorithm2`$ requires the calculation of Ritz eigenpairs of $`T_n`$, which makes for an overhead proportional to $`n^2`$ when using the QR algorithm for the eigenvalues and the inverse iteration for the eigenvectors . Since the complexity of the Lanczos algorithm is $`nN`$, the relative overhead is proportional to $`n/N`$. For moderate gauge couplings and lattice sizes this is a small percentage.

I conclude that the algorithms of  may be used in situations where a high accuracy is not required and/or $`A`$ is well-conditioned.

Experience with overlap fermions shows that high accuracy is often essential . In such situations the $`Algorithm2`$ is best suited.