An iterative method to compute the sign function of a non-Hermitian matrix and its application to the overlap Dirac operator at nonzero chemical potential

The overlap Dirac operator in lattice QCD requires the computation of the sign function of a matrix. While this matrix is usually Hermitian, it becomes non-Hermitian in the presence of a quark chemical potential. We show how the action of the sign function of a non-Hermitian matrix on an arbitrary vector can be computed efficiently on large lattices by an iterative method. A Krylov subspace approximation based on the Arnoldi algorithm is described for the evaluation of a generic matrix function. The efficiency of the method is spoiled when the matrix has eigenvalues close to a function discontinuity. This is cured by adding a small number of critical eigenvectors to the Krylov subspace, for which we propose two different deflation schemes. The ensuing modified Arnoldi method is then applied to the sign function, which has a discontinuity along the imaginary axis. The numerical results clearly show the improved efficiency of the method. Our modification is particularly effective when the action of the sign function of the same matrix has to be computed many times on different vectors, e.g., if the overlap Dirac operator is inverted using an iterative method.     

Numerical methods for the QCDd overlap operator. I. Sign-function and error bounds

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.     

A nested Krylov subspace method to compute the sign function of large complex matrices

We present an acceleration of the well-established Krylov–Ritz methods to compute the sign function of large complex matrices, as needed in lattice QCD simulations involving the overlap Dirac operator at both zero and nonzero baryon density. Krylov–Ritz methods approximate the sign function using a projection on a Krylov subspace. To achieve a high accuracy this subspace must be taken quite large, which makes the method too costly. The new idea is to make a further projection on an even smaller, nested Krylov subspace. If additionally an intermediate preconditioning step is applied, this projection can be performed without affecting the accuracy of the approximation, and a substantial gain in efficiency is achieved for both Hermitian and non-Hermitian matrices. The numerical efficiency of the method is demonstrated on lattice configurations of sizes ranging from 4⁴ to 10⁴, and the new results are compared with those obtained with rational approximation methods.     

Computing charge densities with partially reorthogonalized Lanczos

This paper considers the problem of computing charge densities in a density functional theory (DFT) framework. In contrast to traditional, diagonalization-based, methods, we utilize a technique which exploits a Lanczos basis, without explicit reference to individual eigenvectors. The key ingredient of this new approach is a partial reorthogonalization strategy whose goal is to ensure a good level of orthogonality of the basis vectors. The experiments reveal that the method can be a few times faster than ARPACK, the implicit restart Lanczos method. This is achievable by exploiting more memory and BLAS3 (dense) computations while avoiding the frequent updates of eigenvectors inherent to all restarted Lanczos methods.     

基于重启Lanczos过程的模型降阶方法

针对大规模的线性时不变系统,提出了基于重启Lanczos过程的模型降阶方法。首先,通过重启Lanczos过程分别得到原始系统的可控Gram矩阵的近似矩阵及可观Gram矩阵的近似矩阵。然后,根据原始系统的可控Gram矩阵及可观Gram矩阵所满足的Lyapunov方程构造映射Sylvester方程并求解,对解进行双正交化,得到降阶所需的变换矩阵,从而得到降阶系统。运用此方法对大规模线性时不变系统进行降阶,能够得到具有较高近似精度的稳定的降阶系统。最后,数值算例验证了此方法是行之有效的。     

Model order reduction methods for coupled systems in the time domain using Laguerre polynomials

In this paper, based on Laguerre polynomials, we present new methods for model reduction of coupled systems in the time domain. By appropriately selected projection matrices, a reduced order system is produced to retain the topology structure of the original system. Meanwhile, it preserves a desired number of Laguerre coefficients of the system’s output, thereby providing good approximation accuracy. We also study the two-sided projection method in the time domain, as well as the stability of reduced order systems. Two numerical examples are used to illustrate the efficiency of the proposed methods.     

Numerical methods for the QCD overlap operator: III. Nested iterations

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.     

On the application of an element-by-element lanczos solver to large offshore structural engineering problems

Iterative methods for solving large sets of linear equations have been used as an alternative to direct methods of solution since the early beginning of numerical analysis. The conjugate gradient method (CCM), one of the most widely used, seeks a solution that minimizes the potential energy of the finite element assemblage. Recently, the use of Lanczos algorithm for the solution of large sets of linear equations has been re-examined. Lanczos biorthogonalization procedure is an oblique projection method that provides a solution approximation whose residual is orthogonal to a Krylov subspace. It has been shown that Lanczos and CGM share several properties but the former has the advantage of not being necessary to compute the approximated solution at each iteration. Jacobi preconditioning can also be employed in order to accelerate convergence. The Lanczos procedure was implemented using an element-by-element (EBE) scheme. The applications spread over typical offshore engineering problems encompassing regular and irregular meshes. These problems are normally ill-conditioned when compared with continuum problems. For all the analyses addressed the element-by-element Lanczos procedures presented outstanding computational efficiency.     

Model-order reductions for MIMO systems using global Krylov subspace methods

This paper presents theoretical foundations of global Krylov subspace methods for model order reductions. This method is an extension of the standard Krylov subspace method for multiple-inputs multiple-outputs (MIMO) systems. By employing the congruence transformation with global Krylov subspaces, both one-sided Arnoldi and two-sided Lanczos oblique projection methods are explored for both single expansion point and multiple expansion points. In order to further reduce the computational complexity for multiple expansion points, adaptive-order multiple points moment matching algorithms, or the so-called rational Krylov space method, are also studied. Two algorithms, including the adaptive-order rational global Arnoldi (AORGA) algorithm and the adaptive-order global Lanczos (AOGL) algorithm, are developed in detail. Simulations of practical dynamical systems will be conducted to illustrate the feasibility and the efficiency of proposed methods.     

Bilinear Lanczos components for fast dimensionality reduction and feature extraction

Matrix-based methods such as generalized low rank approximations of matrices (GLRAM) have gained wide attention from researchers in pattern recognition and machine learning communities. In this paper, a novel concept of bilinear Lanczos components (BLC) is introduced to approximate the projection vectors obtained from eigen-based methods without explicit computing eigenvectors of the matrix. This new method sequentially reduces the reconstruction error for a Frobenius-norm based optimization criterion, and the resulting approximation performance is thus improved during successive iterations. In addition, a theoretical clue for selecting suitable dimensionality parameters without losing classification information is presented in this paper. The BLC approach realizes dimensionality reduction and feature extraction by using a small number of Lanczos components. Extensive experiments on face recognition and image classification are conducted to evaluate the efficiency and effectiveness of the proposed algorithm. Results show that the new approach is competitive with the state-of-the-art methods, while it has a much lower training cost.     

Zur iterativen Einschließung von positiven Inversen

In a partially ordered ring we give a method for the two-sided approximation of the inverse of an element which needs the same number of multiplications as known methods of convergence order two. The new method has the convergence order three.     

Arnoldi-based model reduction for fractional order linear systems

In this paper, the Arnoldi-based model reduction methods are employed to fractional order linear time-invariant systems. The resulting model has a smaller dimension, while its fractional order is the same as that of the original system. The error and stability of the reduced model are discussed. And to overcome the local convergence of Padé approximation, the multi-point Arnoldi algorithm, which can recursively generate a reduced-order orthonormal basis from the corresponding Krylov subspace, is used. Numerical examples are given to illustrate the accuracy and efficiency of the proposed methods.     

一类特殊类型子空间上Ritz对的性质及其应用

§1.引言设A∈R~(M×N),定义增广矩阵 其中上标T表示转置。不失一般性,假设M≥N,设σ_i,i=1,2,…,N是A的奇异值,u_i和v_i分别是对应的左右奇异向量,奇异值按从小到大或从大到小的顺序排列,则     

解大规模非对称线性方程组的Lanczos方法和精化Lanczos方法

引言科学工程计算的核心问题之一是数值求解大规模线性方程组,即给定n阶非奇异的非对1期贾仲孝等:解大规模非对称线性方程组的Lanczos方法和精化Lanczos方法称矩阵A和n维向量b,求一个。维向量x,使得Ax=b.(l)观察到该问题可以转化为     

Orthogonal projection method for eigenpair derivatives of large symmetric matrices

In this paper, orthogonal projection method (OPM) is proposed to calculate a few eigenpair derivatives of large real symmetric matrices, if the eigenvalues are simple. By projecting large matrix-valued functions to small subspaces dependent on parameter, linear systems of equations for eigenpair derivatives are greatly reduced from the original matrix size. Error bounds on the eigenpairs and their derivatives computed by OPM are established with trivial conditions, which is one of the main contributions of this paper. It is shown that if the deviations from small subspaces to invariant subspace corresponding to the desired eigenvalues and the derivatives of residuals of the computed eigenpairs tend to zero, then the computed eigenpairs and their derivatives converge to the desired eigenpairs and their derivatives. Next, a strategy to generate small subspaces in OPM is given based on the implicitly restarted Lanczos process (IRLP) for symmetric matrix-valued function. Convergence of the IRLP-based OPM is analysed in detail, which is the main emphasis of this method. Finally some numerical experiments are reported to show the efficiency of our method.     

基于交叉Gram矩阵的双侧H_2最优模型降阶方法

针对单输入单输出(SISO)线性时不变系统,提出了Grassmann流形上基于交叉Gram矩阵的双侧H2最优模型降阶方法。首先,将误差系统的H2范数通过交叉Gram矩阵表示,并且把它看成关于变换矩阵的代价函数。其次,引入Grassmann流形,将代价函数看作是定义在Grassmann流形上的非负实值函数。然后,在Grassmann流形上进行线性搜索,寻找使得代价函数尽可能小的一组变换矩阵。运用此方法对大规模SISO线性时不变系统进行降阶,可以得到精度较高的降阶系统。最后,数值算例验证了该算法的近似效果。     

On the evaluation of the eigenvalues of the finite differences laplacian over a hexagon

Fast numerical methods for the evaluation of the eigenvalues of the finite-differences laplacian over a regular hexagon are devised. At first we show how this eigenvalue problem can be splitted into 3 eigenvalue problems, of lower dimension (reduced by a factor of 1/6), for the discrete laplacian over a regular triangle with suitable boundary conditions. Then, expressing explicitely the eigenvalues and the eigenvectors of the discrete laplacian over a triangle in terms of the coefficients of the discrete Fourier transform, we show how to deal efficiently with each subproblem. In particular we show that each step of the shifted inverse power method, for the approximation of the eigenvalues, costs O(n²log n) arithmetic operations in a sequential model of computation, and O(log n) steps with n² processors in a parallel model of computation, where n is the number of the nodes on the edge of the hexagon. Similar estimates hold for the orthogonal iterations (subspaces iterations) method and for Lanczos method. This approach includes the deflation of the eigenvalues of the triangle from those of the hexagon. These results improve the methods given by Bauer and Reiss [1] allowing a higher precision in the approximation of the eigenvalues of the laplacian and reducing the computational cost, either in a sequential or in a parallel model of computation.     

Parallel eigenvalue calculation based on multiple shift–invert Lanczos and contour integral based spectral projection method

We discuss the possibility of using multiple shift–invert Lanczos and contour integral based spectral projection method to compute a relatively large number of eigenvalues of a large sparse and symmetric matrix on distributed memory parallel computers. The key to achieving high parallel efficiency in this type of computation is to divide the spectrum into several intervals in a way that leads to optimal use of computational resources. We discuss strategies for dividing the spectrum. Our strategies make use of an eigenvalue distribution profile that can be estimated through inertial counts and cubic spline fitting. Parallel sparse direct methods are used in both approaches. We use a simple cost model that describes the cost of computing k eigenvalues within a single interval in terms of the asymptotic cost of sparse matrix factorization and triangular substitutions. Several computational experiments are performed to demonstrate the effect of different spectrum division strategies on the overall performance of both multiple shift–invert Lanczos and the contour integral based method. We also show the parallel scalability of both approaches in the strong and weak scaling sense. In addition, we compare the performance of multiple shift–invert Lanczos and the contour integral based spectral projection method on a set of problems from density functional theory (DFT).     

A parallel Lanczos method for symmetric generalized eigenvalue problems

The Lanczos algorithm is a very effective method for finding extreme eigenvalues of symmetric matrices. In this paper, we present our parallel version of the Lanczos method for symmetric generalized eigenvalue problem, PLANSO. PLANSO is based on a sequential package called LANSO which implements the Lanczos algorithm with partial reorthogonalization. It is portable to all parallel machines that support MPI and it is easy to interface with most parallel computing packages. Through numerical experiments, we demonstrate that it achieves similar parallel efficiency as PARPACK, but uses considerably less time. Received: 21 January 1998 / Accepted: 10 June 1999     

An extended basis inexact shift–invert Lanczos for the efficient solution of large-scale generalized eigenproblems

This paper proposes a technique, based on the Inexact Shift–Invert Lanczos (ISIL) method with Inexact Jacobi Orthogonal Component Correction (IJOCC) refinement, and a preconditioned conjugate-gradient (PCG) linear solver with multilevel preconditioner, for finding several eigenvalues for generalized symmetric eigenproblems. Several eigenvalues are found by constructing (with the ISIL process) an extended projection basis. Presented results of numerical experiments confirm the technique can be effectively applied to challenging, large-scale problems characterized by very dense spectra, such as resonant cavities with spatial dimensions which are large with respect to wavelengths of the resonating electromagnetic fields. It is also shown that the proposed scheme based on inexact linear solves delivers superior performance, as compared to methods which rely on exact linear solves, indicating tremendous potential of the ‘inexact solve’ concept. Finally, the scheme which generates an extended projection basis is found to provide a cost-efficient alternative to classical deflation schemes when several eigenvalues are computed.