Homotopy Method for the Large, Sparse, Real Nonsymmetric Eigenvalue Problem

A homotopy method to compute the eigenpairs, i.e., the eigenvectors and eigenvalues, of a given real matrix A1 is presented. From the eigenpairs of some real matrix A0, the eigenpairs of A(t) ≡ (1 − t)A0 + tA1 are followed at successive "times" from t = 0 to t = 1 using continuation. At t = 1, the eigenpairs of the desired matrix A1 are found. The following phenomena are present when following the eigenpairs of a general nonsymmetric matrix: • bifurcation, • ill conditioning due to nonorthogonal eigenvectors, • jumping of eigenpaths. These can present considerable computational difficulties. Since each eigenpair can be followed independently, this algorithm is ideal for concurrent computers. The homotopy method has the potential to compete with other algorithms for computing a few eigenvalues of large, sparse matrices. It may be a useful tool for determining the stability of a solution of a PDE. Some numerical results will be presented

Spectral methods in general relativistic astrophysics

We present spectral methods developed in our group to solve three-dimensional partial differential equations. The emphasis is put on equations arising from astrophysical problems in the framework of general relativity.Comment: 51 pages, elsart (Elsevier Preprint), 19 PostScript figures, submitted to Journal of Computational & Applied Mathematic

Short-recurrence Krylov subspace methods for the overlap Dirac operator at nonzero chemical potential

The overlap operator in lattice QCD requires the computation of the sign function of a matrix, which is non-Hermitian in the presence of a quark chemical potential. In previous work we introduced an Arnoldi-based Krylov subspace approximation, which uses long recurrences. Even after the deflation of critical eigenvalues, the low efficiency of the method restricts its application to small lattices. Here we propose new short-recurrence methods which strongly enhance the efficiency of the computational method. Using rational approximations to the sign function we introduce two variants, based on the restarted Arnoldi process and on the two-sided Lanczos method, respectively, which become very efficient when combined with multishift solvers. Alternatively, in the variant based on the two-sided Lanczos method the sign function can be evaluated directly. We present numerical results which compare the efficiencies of a restarted Arnoldi-based method and the direct two-sided Lanczos approximation for various lattice sizes. We also show that our new methods gain substantially when combined with deflation.Comment: 14 pages, 4 figures; as published in Comput. Phys. Commun., modified data in Figs. 2,3 and 4 for improved implementation of FOM algorithm, extended discussion of the algorithmic cos

Computational aspects of helicopter trim analysis and damping levels from Floquet theory

Helicopter trim settings of periodic initial state and control inputs are investigated for convergence of Newton iteration in computing the settings sequentially and in parallel. The trim analysis uses a shooting method and a weak version of two temporal finite element methods with displacement formulation and with mixed formulation of displacements and momenta. These three methods broadly represent two main approaches of trim analysis: adaptation of initial-value and finite element boundary-value codes to periodic boundary conditions, particularly for unstable and marginally stable systems. In each method, both the sequential and in-parallel schemes are used and the resulting nonlinear algebraic equations are solved by damped Newton iteration with an optimally selected damping parameter. The impact of damped Newton iteration, including earlier-observed divergence problems in trim analysis, is demonstrated by the maximum condition number of the Jacobian matrices of the iterative scheme and by virtual elimination of divergence. The advantages of the in-parallel scheme over the conventional sequential scheme are also demonstrated

Diagonalization- and Numerical Renormalization-Group-Based Methods for Interacting Quantum Systems

In these lecture notes, we present a pedagogical review of a number of related {\it numerically exact} approaches to quantum many-body problems. In particular, we focus on methods based on the exact diagonalization of the Hamiltonian matrix and on methods extending exact diagonalization using renormalization group ideas, i.e., Wilson's Numerical Renormalization Group (NRG) and White's Density Matrix Renormalization Group (DMRG). These methods are standard tools for the investigation of a variety of interacting quantum systems, especially low-dimensional quantum lattice models. We also survey extensions to the methods to calculate properties such as dynamical quantities and behavior at finite temperature, and discuss generalizations of the DMRG method to a wider variety of systems, such as classical models and quantum chemical problems. Finally, we briefly review some recent developments for obtaining a more general formulation of the DMRG in the context of matrix product states as well as recent progress in calculating the time evolution of quantum systems using the DMRG and the relationship of the foundations of the method with quantum information theory.Comment: 51 pages; lecture notes on numerically exact methods. Pedagogical review appearing in the proceedings of the "IX. Training Course in the Physics of Correlated Electron Systems and High-Tc Superconductors", Vietri sul Mare (Salerno, Italy, October 2004

A Singular Value Thresholding Algorithm for Matrix Completion

This paper introduces a novel algorithm to approximate the matrix with minimum nuclear norm among all matrices obeying a set of convex constraints. This problem may be understood as the convex relaxation of a rank minimization problem and arises in many important applications as in the task of recovering a large matrix from a small subset of its entries (the famous Netflix problem). Off-the-shelf algorithms such as interior point methods are not directly amenable to large problems of this kind with over a million unknown entries. This paper develops a simple first-order and easy-to-implement algorithm that is extremely efficient at addressing problems in which the optimal solution has low rank. The algorithm is iterative, produces a sequence of matrices {X^k,Y^k}, and at each step mainly performs a soft-thresholding operation on the singular values of the matrix Y^k. There are two remarkable features making this attractive for low-rank matrix completion problems. The first is that the soft-thresholding operation is applied to a sparse matrix; the second is that the rank of the iterates {X^k} is empirically nondecreasing. Both these facts allow the algorithm to make use of very minimal storage space and keep the computational cost of each iteration low. On the theoretical side, we provide a convergence analysis showing that the sequence of iterates converges. On the practical side, we provide numerical examples in which 1,000 × 1,000 matrices are recovered in less than a minute on a modest desktop computer. We also demonstrate that our approach is amenable to very large scale problems by recovering matrices of rank about 10 with nearly a billion unknowns from just about 0.4% of their sampled entries. Our methods are connected with the recent literature on linearized Bregman iterations for ℓ_1 minimization, and we develop a framework in which one can understand these algorithms in terms of well-known Lagrange multiplier algorithms

Fast computation of the Gauss hypergeometric function with all its parameters complex with application to the Poschl-Teller-Ginocchio potential wave functions

The fast computation of the Gauss hypergeometric function 2F1 with all its parameters complex is a difficult task. Although the 2F1 function verifies numerous analytical properties involving power series expansions whose implementation is apparently immediate, their use is thwarted by instabilities induced by cancellations between very large terms. Furthermore, small areas of the complex plane are inaccessible using only 2F1 power series formulas, thus rendering 2F1 evaluations impossible on a purely analytical basis. In order to solve these problems, a generalization of R.C. Forrey's transformation theory has been developed. The latter has been successful in treating the 2F1 function with real parameters. As in real case transformation theory, the large canceling terms occurring in 2F1 analytical formulas are rigorously dealt with, but by way of a new method, directly applicable to the complex plane. Taylor series expansions are employed to enter complex areas outside the domain of validity of power series analytical formulas. The proposed algorithm, however, becomes unstable in general when |a|,|b|,|c| are moderate or large. As a physical application, the calculation of the wave functions of the analytical Poschl-Teller-Ginocchio potential involving 2F1 evaluations is considered.Comment: 29 pages; accepted in Computer Physics Communication