The numerical and computational aspects of the overlap formalism in lattice quantum chromodynamics are extremely demanding due to a matrix-vector product that involves the sign function of the hermitian Wilson matrix. In this paper we investigate several methods to compute the product of the matrix sign-function with a vector, in particular Lanczos based methods and partial fraction expansion methods. Our goal is two-fold: we give realistic comparisons between known methods together with novel approaches and we present error bounds which allow to guarantee a given accuracy when terminating the Lanczos method and the multishift-CG solver, applied within the partial fraction expansion methods.

This paper provides two results on the numerical behavior of the classical Gram-Schmidt algorithm. The first result states that, provided the normal equations associated with the initial vectors are numerically nonsingular, the loss of orthogonality of the vectors computed by the classical Gram-Schmidt algorithm depends quadratically on the condition number of the initial vectors. The second result states that, provided the initial set of vectors has numerical full rank, two

The numerical and computational aspects of chiral fermions in lattice quantum chromodynamics are extremely demanding. In the overlap framework, the computation of the fermion propagator leads to a nested iteration where the matrix vector multiplications in each step of an outer iteration have to be accomplished by an inner iteration; the latter approximates the product of the sign function of the hermitian Wilson fermion matrix with a vector.In this paper we investigate aspects of this nested paradigm. We examine several Krylov subspace methods to be used as an outer iteration for both propagator computations and the Hybrid Monte-Carlo scheme. We establish criteria on the accuracy of the inner iteration which allow to preserve an a priori given precision for the overall computation. It will turn out that the accuracy of the sign function can be relaxed as the outer iteration proceeds. Furthermore, we consider preconditioning strategies, where the preconditioner is built upon an inaccurate approximation to the sign function. Relaxation combined with preconditioning allows for considerable savings in computational efforts up to a factor of 4 as our numerical experiments illustrate. We also discuss the possibility of projecting the squared overlap operator into one chiral sector.

We investigate optimal choices for the (outer) iteration method to use when solving linear systems with Neuberger's overlap operator in QCD. Different formulations for this operator give rise to different iterative solvers, which are optimal for the respective formulation. We compare these methods in theory and practice to find the overall optimal one. For the first time, we apply the so-called SUMR method of Jagels and Reichel to the shifted unitary version of Neuberger's operator, and show that this method is in a sense the optimal choice for propagator computations. When solving the "squared" equations in a dynamical simulation with two degenerate flavours, it turns out that the CG method should be used.

We consider the solution of the linear systemfor various real values of σ. This family of shifted systems arises, for example, in Tikhonov regularization and computations in lattice quantum chromodynamics. For each single shift σ this system can be solved using the conjugate gradient method for least squares problems (CGLS). In literature various implementations of the, so-called, multishift CGLS methods have been proposed. These methods are mathematically equivalent to applying the CGLS method to each shifted system separately but they solve all systems simultaneously and require only two matrix-vector products (one by A and one by A T ) and two inner products per iteration step. Unfortunately, numerical experiments show that, due to roundoff errors, in some cases these implementations of the multishift CGLS method can only attain an accuracy that depends on the square of condition number of the matrix A. In this paper we will argue that, in the multishift CGLS method, the impact on the attainable accuracy of rounding errors in the Lanczos part of the method is independent of the effect of roundoff errors made in the construction of the iterates. By making suitable design choices for both parts, we derive a new (and efficient) implementation that tries to remove the limitation of previous proposals. A partial roundoff error analysis and various numerical experiments show promising results.