On the achievement of high fidelity and scalability for large‐scale diagonalizations in grid‐based DFT simulations

Recent advance in high performance computing (HPC) resources has opened the possibility to expand the scope of density functional theory (DFT) simulations toward large and complex molecular systems. This work proposes a numerically robust method that enables scalable diagonalizations of large DFT Hamiltonian matrices, particularly with thousands of computing CPUs (cores) that are usual these days in terms of sizes of HPC resources. The well‐known Lanczos method is extensively refactorized to overcome its weakness for evaluation of multiple degenerate eigenpairs that is the substance of DFT simulations, where a multilevel parallelization is adopted for scalable simulations in as many cores as possible. With solid benchmark tests for realistic molecular systems, the fidelity of our method are validated against the locally optimal block preconditioned conjugated gradient (LOBPCG) method that is widely used to simulate electronic structures. Our method may waste computing resources for simulations of molecules whose degeneracy cannot be reasonably estimated. But, compared to LOBPCG method, it is fairly excellent in perspectives of both speed and scalability, and particularly has remarkably less (< 10%) sensitivity of performance to the random nature of initial basis vectors. As a promising candidate for solving electronic structures of highly degenerate systems, the proposed method can make a meaningful contribution to migrating DFT simulations toward extremely large computing environments that normally have more than several tens of thousands of computing cores.     

Numerical Experiences with New Truncated Newton Methods in Large Scale Unconstrained Optimization

Recently, in [12] a very general class oftruncated Newton methods has been proposed for solving large scale unconstrained optimization problems. In this work we present the results of an extensive numericalexperience obtained by different algorithms which belong to the preceding class. This numerical study, besides investigating which arethe best algorithmic choices of the proposed approach, clarifies some significant points which underlies every truncated Newton based algorithm.     

Krylov subspace methods, biorthogonal polynomials and Padé-type approximants

It is shown that Krylov subspace methods for solving systems of linear equations can be based on formal biorthogonal polynomials and on Padé-type and Padé approximants. New algorithms for their implementation are derived. This revised version was published online in June 2006 with corrections to the Cover Date.     

Model reduction using the Vorobyev moment problem

Given a nonsingular complex matrix and complex vectors v and w of length N, one may wish to estimate the quadratic form w ^* A ^− 1 v, where w ^* denotes the conjugate transpose of w. This problem appears in many applications, and Gene Golub was the key figure in its investigations for decades. He focused mainly on the case A Hermitian positive definite (HPD) and emphasized the relationship of the algebraically formulated problems with classical topics in analysis - moments, orthogonal polynomials and quadrature. The essence of his view can be found in his contribution Matrix Computations and the Theory of Moments, given at the International Congress of Mathematicians in Zürich in 1994. As in many other areas, Gene Golub has inspired a long list of coauthors for work on the problem, and our contribution can also be seen as a consequence of his lasting inspiration. In this paper we will consider a general mathematical concept of matching moments model reduction, which as well as its use in many other applications, is the basis for the development of various approaches for estimation of the quadratic form above. The idea of model reduction via matching moments is well known and widely used in approximation of dynamical systems, but it goes back to Stieltjes, with some preceding work done by Chebyshev and Heine. The algebraic moment matching problem can for A HPD be formulated as a variant of the Stieltjes moment problem, and can be solved using Gauss-Christoffel quadrature. Using the operator moment problem suggested by Vorobyev, we will generalize model reduction based on matching moments to the non-Hermitian case in a straightforward way. Unlike in the model reduction literature, the presented proofs follow directly from the construction of the Vorobyev moment problem. The work was supported by the GAAS grant IAA100300802 and by the Institutional Research Plan AV0Z10300504.     

Inexact Newton method via Lanczos decomposed technique for solving box-constrained nonlinear systems

This paper proposes an inexact Newton method via the Lanczos decomposed technique for solving the box-constrained nonlinear systems. An iterative direction is obtained by solving an affine scaling quadratic model with the Lanczos decomposed technique. By using the interior backtracking line search technique, an acceptable trial step length is found along this direction. The global convergence and the fast local convergence rate of the proposed algorithm are established under some reasonable conditions. Furthermore, the results of the numerical experiments show the effectiveness of the pro- posed algorithm.     

The geometry of PLS1 explained properly: 10 key notes on mathematical properties of and some alternative algorithmic approaches to PLS1 modelling

The insight from, and conclusions of this paper motivate efficient and numerically robust ‘new’ variants of algorithms for solving the single response partial least squares regression (PLS1) problem. Prototype MATLAB code for these variants are included in the Appendix. The analysis of and conclusions regarding PLS1 modelling are based on a rich and nontrivial application of numerous key concepts from elementary linear algebra. The investigation starts with a simple analysis of the nonlinear iterative partial least squares (NIPALS) PLS1 algorithm variant computing orthonormal scores and weights. A rigorous interpretation of the squared P ‐loadings as the variable‐wise explained sum of squares is presented. We show that the orthonormal row‐subspace basis of W ‐weights can be found from a recurrence equation. Consequently, the NIPALS deflation steps of the centered predictor matrix can be replaced by a corresponding sequence of Gram–Schmidt steps that compute the orthonormal column‐subspace basis of T ‐scores from the associated non‐orthogonal scores. The transitions between the non‐orthogonal and orthonormal scores and weights (illustrated by an easy‐to‐grasp commutative diagram), respectively, are both given by QR factorizations of the non‐orthogonal matrices. The properties of singular value decomposition combined with the mappings between the alternative representations of the PLS1 ‘truncated’ X data (including P ^t W ) are taken to justify an invariance principle to distinguish between the PLS1 truncation alternatives. The fundamental orthogonal truncation of PLS1 is illustrated by a Lanczos bidiagonalization type of algorithm where the predictor matrix deflation is required to be different from the standard NIPALS deflation. A mathematical argument concluding the PLS1 inconsistency debate (published in 2009 in this journal) is also presented. Copyright © 2014 John Wiley & Sons, Ltd.     

Fast generalized cross validation using Krylov subspace methods

The task of fitting smoothing spline surfaces to meteorological data such as temperature or rainfall observations is computationally intensive. The generalized cross validation (GCV) smoothing algorithm, if implemented using direct matrix techniques, is O(n ³) computationally, and memory requirements are O(n ²). Thus, for data sets larger than a few hundred observations, the algorithm is prohibitively slow. The core of the algorithm consists of solving series of shifted linear systems, and iterative techniques have been used to lower the computational complexity and facilitate implementation on a variety of supercomputer architectures. For large data sets though, the execution time is still quite high. In this paper we describe a Lanczos based approach that avoids explicitly solving the linear systems and dramatically reduces the amount of time required to fit surfaces to sets of data.      

Implicitly restarted projection algorithm for solving optimization problems

An algorithm for solving the problem of minimizing a quadratic function subject to ellipsoidal constraints is introduced. This algorithm is based on the impHcitly restarted Lanczos method to construct a basis for the Krylov subspace in conjunction with a model trust region strategy to choose the step. The trial step is computed on the small dimensional subspace that lies inside the trust region.

One of the main advantages of this algorithm is the way that the Krylov subspace is terminated. We introduce a terminationcondition that allows the gradient to be decreased on that subspace.

A convergence theory for this algorithm is presented. It is shown that this algorithm is globally convergent and it shouldcope quite well with large scale minimization problems. This theory is sufficiently general that it holds for any algorithm that projects the problem on a lower dimensional subspace.     

Error estimates for large-scale ill-posed problems

The computation of an approximate solution of linear discrete ill-posed problems with contaminated data is delicate due to the possibility of severe error propagation. Tikhonov regularization seeks to reduce the sensitivity of the computed solution to errors in the data by replacing the given ill-posed problem by a nearby problem, whose solution is less sensitive to perturbation. This regularization method requires that a suitable value of the regularization parameter be chosen. Recently, Brezinski et al. (Numer Algorithms 49, 2008) described new approaches to estimate the error in approximate solutions of linear systems of equations and applied these estimates to determine a suitable value of the regularization parameter in Tikhonov regularization when the approximate solution is computed with the aid of the singular value decomposition. This paper discusses applications of these and related error estimates to the solution of large-scale ill-posed problems when approximate solutions are computed by Tikhonov regularization based on partial Lanczos bidiagonalization of the matrix. The connection between partial Lanczos bidiagonalization and Gauss quadrature is utilized to determine inexpensive bounds for a family of error estimates. In memory of Gene H. Golub. This work was supported by MIUR under the PRIN grant no. 2006017542-003 and by the University of Cagliari.