{'fraction_non_alphanumeric': 0.029754959159859977, 'fraction_numerical': 0.009626604434072346, 'mean_word_length': 4.176985804892782, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A fast minimal residual solver for overlap fermions

Computing quark propagators with overlap fermions requires the solution of a shifted unitary linear system. Jagels and Reichel have shown that for such systems it is possible to construct a minimal residual algorithm by short recurrences. The J\"ulich-Wuppertal group have found this algorithm to be the fastest among overlap solvers. In this paper we present a three-term recurrence for the Arnoldi unitary process. Using the new recurrence we construct a minimal residual solver which is the fastest among all Krylov subspace algorithms considered so far for the overlap inversion.


Introduction
Lattice theories with chiral quarks provide an accurate tool for studying the physics of strong interactions. The physical information of these theories is contained in quark propagators, which are then combined to construct meson, nucleon and even more complicated elementary particle propagators. Quark propagator computations amount to solving large linear systems of the type: where D ∈ C N ×N is a dense matrix operator representing the chiral Dirac operator on a regular four dimensional space-time lattice, x, b ∈ C N are the quark propagator and its source.
More than a decade ago two formulations of chiral fermion theories were discovered: the domain wall fermions [1,2] and overlap fermions [3]. These apparently different formulations are closely related to each other [4]. In particular, the truncated overlap variant of domain wall fermions [5] can be shown to be equivalent to overlap fermions in four dimensions at any lattice spacing [6]. Therefore, for computational purposes it is sufficient to consider the Neuberger Dirac operator, which is a shifted unitary matrix of the form [7]: Here D W is the Wilson-Dirac operator, which is a non-Hermitian sparse matrix. It is easy to see why V is a unitary matrix: starting from the singular value decomposition of D W = XΣY * , one gets V = XY * .
Another way to look into the Neuberger operator is to use operators Γ 5 =diag(I N/2 , −I N/2 ) 1 andΓ 5 =sign(H W ), where H W = Γ 5 D W is the Hermitian Wilson-Dirac operator. In this case it is easy to see that the overlap operator takes the form: This representation can be used to show that D * D commutes with Γ 5 . As a result, D * D is block diagonal in the chiral basis of Γ 5 . This observation has important applications: if we would like to solve linear systems with coefficient matrix D * D we can use Conjugate Gradients (CG) algorithm, which we know is optimal. Indeed, CG is very well suited for the molecular dynamics algorithm which generates the ensemble of gauge fields in lattice QCD with overlap fermions.
However, another important task is the computation of quark propagators which we consider here. This requires the solution of the linear systems of the type 1.1. One can still use CG on normal equations (CGNE) but this method may not be optimal in this case. In this paper we seek optimal solutions of the above system in the Krylov subspace: It is well known that the generalized minimal residual method (GMRES) is the optimal method that can be constructed in this subspace. However, it is well known as well that the underlying Arnoldi process produces long recurrences, a feature which limits the benefits of the GMRES algorithm for large problems. Nonetheless, Jagels and Reichel [8] found that for shifted unitary systems it is possible to construct an Arnoldi process with short recurrences. The method, the shifted unitary minimal residual algorithm (SUMR) was shown by Jülich-Wuppertal group [9] to be the fastest linear solver for quark propagator computations.
Here we propose an approach, which is based on a three-term recurrence of the unitary Arnoldi process. This process is used to construct two optimal Krylov subspace methods: the shifted unitary orthogonal method (SUOM) and its geometric optimal counterpart SHUMR.
We show that this algorithm converges faster then SUMR for all lattices and quark masses considered in this paper.
The paper is organised as follows: section 2 describes the new recurrence of the unitary Arnoldi process. Section 3 deals with the construction of the SUOM and SHUMR algorithms.
In section 4 we compare directly these algorithms to SUMR and other Krylov subspace methods. Conclusions are drawn in section 5.

A three-term recurrence for the Arnoldi unitary process
Our aim in this section is to construct an orthogonal basis of the Krylov subspace 1.7 with as few operations as possible which will then be useful for the solution of the linear system 1.1. As a starting point we use the well-known Arnoldi process [10], which is essentially a modified Gram-Schmidt process, given in Algorithm 1. After k steps of this algorithm, the next unnormalized Arnoldi vector is given by: The right hand side contains a linear combination over all computed Arnoldi vectors. This makes the overall complexity of the algorithm to grow quadratically with k. Moreover, the computer memory grows linearly with k, a requirement which is often prohibitive and thus undesirable for large problems like ours.
Fortunately, for unitary matrices one can do better. After k steps of Algorithm 1 we can write: where Q k = [q 1 , . . . , q k ] is the matrix of orthonormal Arnoldi vectors, H k is an k × k upper Hessenberg matrix andH k is a k + 1 × k matrix obtained by appending the row h k+1,k e T k to the matrix H k .

Proposition 2.1H k has orthonormal columns.
Proof. Multiplying both sides of 2.2 from left by Q * k we get: Multiplying both sides of 2.2 from left again but now by (V Q k ) * we get: From 2.3 and 2.2 one gets Q * k V * Q k = H * k and q * k+1 V Q k = h k+1 e T k . Thus, we find that: which proves the proposition. ✷

Corollary 2.2 From 2.5 it is clear that H k is a unitary matrix if its last column is normalised.
This property was used be Jagels and Reichel to write H k as a product of k elementary Givens rotations, which are then exploited to construct a coupled two-term recurrence Arnoldi algorithm (see Algorithm 3.1 of [8]). They note that solving the coupled recurrences would still lead to an algorithm where all vectors are required.
In contrast to their work, we give an algorithm which is defined by short recurrences. We start by noting that the LU-decomposition of the upper Hessenberg matrix H k can be written in the form: where L k is a lower bidiagonal matrix and U k is upper triangular. For our convenience, we take the diagonal elements of U k to be one. Substituting this decomposition into 2.5, multiplying the result from the right by U k and from the left by U * k we get: Since the right hand side is a tridiagonal Hermitian matrix so must be the left hand side. Hence, U k should be upper bidiagonal. This decomposition, eq. 2.6 was used for the first time by Rutishauser to compute the eigenvalues of orthogonal Hessenberg matrices [11]. We can easily extend this decomposition for the matrixH k : Multiplying both sides of 2.2 by U k we get: This way, the next Arnoldi vector can be computed using: These expressions allow us to construct the three-term recurrence Arnoldi unitary process as given in Algorithm 2.
Note that the algorithm needs an additional inner product than the usual (eg. Lanczos) threeterm recurrence algorithms. The reason is the appearance of the matrix V both in the diagonal and subdiagonal terms of 2.10. The algorithm is equivalent to the normal Arnoldi algorithm in exact arithmetic as the following proposition shows: k is a matrix with orthonormal columns.
Proof. The algorithm produces orthonormal vectors Q k and matrices L k and U k such that equations 2.9 is satisfied. Multiplying both sides of these equations by U −1 k we get new equations where L k U −1 k is upper Hessenberg. From Proposition 2.1 follows thatH k and thusL k+1 U −1 k is a matrix with orthonormal columns. ✷

Remark 2.4 Both Arnoldi processes produce the same basis for the Krylov subspace but Algorithm 2 is more efficient: its complexity is linear in k and its computer memory requirement is
constant at each step k.

Optimal Krylov subspace methods for the overlap inversion
The overlap operator is a non-Hermitian matrix operator. For such matrices GMRES is knwon to be the optimal algorithm. We call this class of algorithms geometrically optimal. Another algorithm which can be used for such problems is the full orthogonalisation method (FOM). In this case the k th residual vector, lies orthogonal to the Krylov subspace K k . Because of this algebraic property we call this class of algorithms algebraically optimal. It can be shown that when the norm-minimising process of GMRES is converging rapidly, the residual norms in the corresponding Galerkin process of FOM exhibit similar behaviour [12].
Both methods GMRES and FOM use the Arnoldi process to generate the iterates. In this section we will use the Arnoldi unitary process to construct two algorithms: the first one is the specialisation of FOM to shifted unitary systems, whereas the second is the specialisation of GMRES to these systems.

A full orthogonalisation strategy
We seek the solution in the Krylov subspace spanned by the Arnoldi vectors, i.e.
Assuming for simplicity that b = 1, and using eq. 3.1 the residual vector will be, Using the definition of D, eq. 1.2 and the relation 2.9 we get: where1l k is obtained by appending a row of k zeros to the unit matrix 1l. Projecting this equation onto the Krylov subspace K k one has: or equivalently: Note that the matrix on the left hand side is tridiagonal and we can write: The LU decomposition of the matrix T k is denoted by T k =L kŨk . Hence the solution can be written as: Then from and from the recurrence for z k , eq. 3.8 we get: Finally, denoting, we have: Using this result and the definition of the residual vector, eq. 3.1 we get: Using matricesL k ,Ũ k it is easy to show (see Appendix A) that the following recurrences hold: In this way we have specified the iterations for this linear solver, which is called the shifted Note that the code can be trivially changed in order to get the residual error norm instead of the error norm, or to use x 0 as a starting guess by defining the starting residual error as Another advantage of our MATLAB code is its easy inclusion into a C++ code which uses uBLAS libraries for inner products and Euclidean norms 3 . We have used these libraries to construct such a C++ function for use with overlap fermions. In this case, the matrix-vector multiplication is an external routine which applies the inverse square root of D * W D W to the current Arnoldi vector q followed by a D W multiplication.

A minimal residual strategy
As in the previous case our starting point is the Krylov subspace K k . But instead of projecting the residual vector we seek the minimum of its 2-norm on this subspace. The solution in this case is denoted byx k and is written formally as a solution of a Least Squares Problem (LSP): Requiring x = Q k y and using eq. 3.4 we get: Since Arnoldi vectors are orthonormal, the matrix Q k+1 can be ignored and we end up with a much smaller LSP. Note that the matrixH k =L k U −1 k has orthonormal columns (see Proposition 2.1), a property which can be used to get a short recurrence algorithm as in the case of the SUMR algorithm [8]. However, as in the case of SUOM, we follow a strategy that involves a tridiagonal matrix in the LSP. This way the solution of the smaller problem is given by: where In order to compute the iterates from those of the SUOM algorithm we split the solution vector as follows:z Using equations 3.7-3.9 we get: We solve this LSP by QR factorization of the matrixT k : with O k being a k + 1 × k + 1 unitary matrix and R k an upper tridiagonal matrix. Therefore we have: As it is usual in this case, the unitary matrix can be constructed using Givens rotations. At step k one can express it in the form: From this it is clear that: O k e k+1 = s k e k + c k e k+1 (3.26) which gives: The right hand side is minimal if its first term is minimal. If we assume that R k has full rank then this term must vanish. In this case the solution is given by: Using 3.18-3.21, the solution to the original problem can be written as: wherex In order to simplify the expression we define the matrix P k : and denote ω k = να k s k . This way we have: To complete the algorithm one should write the recurrence on p k . Multiplying both sides of 3.31 by R k from the right we get P k R k = Q k U k , which gives: where by µ k , ε k and θ k we denote the only non-zero entries of the last column of R k . From this equation we get: This completes the description of the method, which we call SHUMR in order to distinguish it from the SUMR algorithm. The details of the QR decomposition are given in Appendix B.
The MATLAB code of the algorithm is listed below. It differs from the SUOM code with lines between "start added lines" and "end added lines". One can make here the same remarks as made in the case of the SUOM.m function. In particular, the code can be easily modified into a C++ code using uBLAS libraries.

A numerical example
We note that SUMR and SHUMR algorithms differ in the underlying Arnoldi process: the SUMR algorithm uses two coupled two-term recurrences which, when are solved, yield the usual long recurrence of the Arnoldi algorithm; the SHUMR algorithm uses a three-term recurrence. In principle, it is possible to compare theoretically the effect of such a difference, but this goes beyond the purpose of this paper. We have chosen instead a direct numerical comparison in case of overlap fermions, since this is of great importance for practical lattice computations.
Before we do so we give a first example in the case of a small, i.e. 200 × 200 unitary matrix.
The example is similar to Example 2 of the paper of Jagels and Reichel [8].

Comparison of algorithms for the overlap inversion
In this section we compare the convergence of SUMR, SUOM and SHUMR algorithms in the case of overlap fermions. SUMR and SHUMR are geometrically optimal algorithms for shifted linear systems, whereas SUOM is algebraically optimal in the sense defined above. For completeness we display the convergence of Conjugate Residuals (CR), Conjugate Gradients on Normal Equations (CGNE) and a special variation of CGNE which we call CG-CHI. The latter exploits the fact that D * D is block diagonal and solves simultaneously the two decoupled chiral systems.
Note that the computation of D as applied to a vector is a numerical problem, which is by now well researched. A good review of these methods can be found in [13]. We use the Lanczos method [14,15], in the double pass version and without H W -eigenvalue projection.
In the figures below we show the convergence of algorithms as a function of Wilson matrixvector multiplication number on 8 3 16 quenched lattices at various couplings and quark masses. The first observation is that SHUMR, SUMR, SUOM and CR are more efficient than CGNE and CG-CHI algorithms. This is observed by the other groups as well [9]. Hence, we decided not run these algorithms further for smaller quark masses. The second observation is that SUMR, SUOM and CR converge neck-to-neck with CR being slightly worse at β = 6. The third observation is that SHUMR converges 10 − 15% faster than SUMR and SUOM. This is to be contrasted to the linear convergence of CR. A second observation at this mass is the emergence of the pattern that SHUMR converges faster than SUMR and the latter converges faster than SUOM. This pattern is confirmed in Figure 4 where the quark mass is lowered to Hence, the best algorithms are the optimal algorithms SHUMR and SUMR with SHUMR converging 10 − 15% faster in all cases. Although both algorithms minimise the residual vector norm they differ in the underlying Arnoldi process. The nature of short recurrences employed by SHUMR may impact numerical properties of the algorithm. Hence, generated Krylov subspaces are different and we may conclude that SHUMR explores it more efficiently. However, at present we have no theoretical tool to characterise the difference.

Conclusions
In    shown in this paper. In particular, we presented a new Arnoldi unitary process, Algorithm 2, which is intrinsically a short recurrence process, a result which appears for the first time.
On the application level we conclude that SHUMR algorithm outperforms all other knwon iterative solvers for quark propagator computations with overlap fermions. The MATLAB codes of SUOM and SHUMR algorithms given here allow an easy conversion to other object oriented code like C++ which uses uBLAS libraries.

Appendix B
In order to complete the derivation of the SHUMR algorithm one has to specify the Givens matrix G k and the last column of R k . Lett k be the last column ofT k andμ k be the last column ofR k . Then, we have:μ k = O ktk ,t k = νe k+1 + t kk e k − γ k−1 e k−1 where t kk = c 1 + c 2 l kk . From the definition of O k , eq. 3.25 it is clear that: It is easy to see that µ k = c k µ + s k ν ands k = c k ν/µ. Using c 2 k + |s k | 2 = 1 one has: For µ = 0 one has c k = 0 and s k = 1.