Approximation Of Matrix Functions Research Articles

Low-Memory Krylov Subspace Methods for Optimal Rational Matrix Function Approximation

We describe a Lanczos-based algorithm for approximating the product of a rational matrix function with a vector. This algorithm, which we call the Lanczos method for optimal rational matrix function approximation (Lanczos-OR), returns the optimal approximation from a given Krylov subspace in a norm depending on the rational function's denominator, and can be computed using the information from a slightly larger Krylov subspace. We also provide a low-memory implementation which only requires storing a number of vectors proportional to the denominator degree of the rational function. Finally, we show that Lanczos-OR can be used to derive algorithms for computing other matrix functions, including the matrix sign function and quadrature based rational function approximations. In many cases, it improves on the approximation quality of prior approaches, including the standard Lanczos method, with little additional computational overhead.

Error Bounds for Lanczos-Based Matrix Function Approximation

We analyze the Lanczos method for matrix function approximation (Lanczos-FA), an iterative algorithm for computing $f(\mathbf{A}) \mathbf{b}$ when $\mathbf{A}$ is a Hermitian matrix and $\mathbf{b}$ is a given vector. Assuming that $f : \mathbb{C} \rightarrow \mathbb{C}$ is piecewise analytic, we give a framework, based on the Cauchy integral formula, which can be used to derive a priori and a posteriori error bounds for Lanczos-FA in terms of the error of Lanczos used to solve linear systems. Unlike many error bounds for Lanczos-FA, these bounds account for fine-grained properties of the spectrum of $\mathbf{A}$, such as clustered or isolated eigenvalues. Our results are derived assuming exact arithmetic, but we show that they are easily extended to finite precision computations using existing theory about the Lanczos algorithm in finite precision. We also provide generalized bounds for the Lanczos method used to approximate quadratic forms $\mathbf{b}^\textsf{H} f(\mathbf{A}) \mathbf{b}$, and demonstrate the effectiveness of our bounds with numerical experiments.

Some uses of the field of values in numerical analysis

In this expository paper we illustrate the role that the field of values (or numerical range) of a matrix plays in connection with certain problems of numerical analysis. These include the approximation of matrix functions and the convergence of preconditioned Krylov subspace methods for solving large systems of equations arising from the discretization of partial differential equations.

Inexact Arnoldi residual estimates and decay properties for functions of non-Hermitian matrices

This paper derives a priori residual-type bounds for the Arnoldi approximation of a matrix function together with a strategy for setting the iteration accuracies in the inexact Arnoldi approximation of matrix functions. Such results are based on the decay behavior of the entries of functions of banded matrices. Specifically, a priori decay bounds for the entries of functions of banded non-Hermitian matrices will be exploited, using Faber polynomial approximation. Numerical experiments illustrate the quality of the results.

Multiple orthogonal polynomials applied to matrix function evaluation

Multiple orthogonal polynomials generalize standard orthogonal polynomials by requiring orthogonality with respect to several inner products. This paper discusses an application to the approximation of matrix functions and presents quadrature rules that generalize the anti-Gauss rules proposed by Laurie.

Enhanced matrix function approximation

Matrix functions of the form $f(A)v$, where $A$ is a large symmetric matrix, $f$ is afunction, and $v\\ne 0$ is a vector, are commonly approximated by first applying a few,say $n$, steps of the symmetric Lanczos process to $A$ with the initial vector $v$ in order todetermine an orthogonal section of $A$. The latter is represented by a (small)$n\\times n$ tridiagonal matrix to which $f$ is applied. This approach uses the $n$ firstLanczos vectors provided by the Lanczos process. However, $n$ steps of the Lanczosprocess yield $n+1$ Lanczos vectors. This paper discusses how the $(n+1)$stLanczos vector can be used to improve the quality of the computed approximation of$f(A)v$. Also the approximation of expressions of the form $v^Tf(A)v$ is considered.

A rational Arnoldi process with applications

SummaryThe rational Arnoldi process is a popular method for the computation of a few eigenvalues of a large non‐Hermitian matrix and for the approximation of matrix functions. The method is particularly attractive when the rational functions that determine the process have only few distinct poles , because then few factorizations of matrices of the form A − zjI have to be computed. We discuss recursion relations for orthogonal bases of rational Krylov subspaces determined by rational functions with few distinct poles. These recursion formulas yield a new implementation of the rational Arnoldi process. Applications of the rational Arnoldi process to the approximation of matrix functions as well as to the computation of eigenvalues and pseudospectra of A are described. The new implementation is compared to several available implementations.

Convergence rates for inverse-free rational approximation of matrix functions

This article deduces geometric convergence rates for approximating matrix functions via inverse-free rational Krylov methods. In applications one frequently encounters matrix functions such as the matrix exponential or matrix logarithm; often the matrix under consideration is too large to compute the matrix function directly, and Krylov subspace methods are used to determine a reduced problem. If many evaluations of a matrix function of the form f(A)v with a large matrix A are required, then it may be advantageous to determine a reduced problem using rational Krylov subspaces. These methods may give more accurate approximations of f(A)v with subspaces of smaller dimension than standard Krylov subspace methods. Unfortunately, the system solves required to construct an orthogonal basis for a rational Krylov subspace may create numerical difficulties and/or require excessive computing time. This paper investigates a novel approach to determine an orthogonal basis of an approximation of a rational Krylov subspace of (small) dimension from a standard orthogonal Krylov subspace basis of larger dimension. The approximation error will depend on properties of the matrix A and on the dimension of the original standard Krylov subspace. We show that our inverse-free method for approximating the rational Krylov subspace converges geometrically (for increasing dimension of the standard Krylov subspace) to a rational Krylov subspace. The convergence rate may be used to predict the dimension of the standard Krylov subspace necessary to obtain a certain accuracy in the approximation. Computed examples illustrate the theory developed.

Solution of nonlinear time-dependent PDEs through componentwise approximation of matrix functions

Exponential propagation iterative (EPI) methods provide an efficient approach to the solution of large stiff systems of ODEs, compared to standard integrators. However, the bulk of the computational effort in these methods is due to products of matrix functions and vectors, which can become very costly at high resolution due to an increase in the number of Krylov projection steps needed to maintain accuracy. In this paper, it is proposed to modify EPI methods by using Krylov subspace spectral (KSS) methods, instead of standard Krylov projection methods, to compute products of matrix functions and vectors. Numerical experiments demonstrate that this modification causes the number of Krylov projection steps to become bounded independently of the grid size, thus dramatically improving efficiency and scalability. As a result, for each test problem featured, as the total number of grid points increases, the growth in computation time is just below linear, while other methods achieved this only on selected test problems or not at all.

New block quadrature rules for the approximation of matrix functions

Golub and Meurant have shown how to use the symmetric block Lanczos algorithm to compute block Gauss quadrature rules for the approximation of certain matrix functions. We describe new block quadrature rules that can be computed by the symmetric or nonsymmetric block Lanczos algorithms and yield higher accuracy than standard block Gauss rules after the same number of steps of the symmetric or nonsymmetric block Lanczos algorithms. The new rules are block generalizations of the generalized averaged Gauss rules introduced by Spalević. Applications to network analysis are presented.

Rational Krylov approximation of matrix functions: Numerical methods and optimal pole selection

AbstractMatrix functions are a central topic of linear algebra, and problems of their numerical approximation appear increasingly often in scientific computing. We review various rational Krylov methods for the computation of large‐scale matrix functions. Emphasis is put on the rational Arnoldi method and variants thereof, namely, the extended Krylov subspace method and the shift‐and‐invert Arnoldi method, but we also discuss the nonorthogonal generalized Leja point (or PAIN) method. The issue of optimal pole selection for rational Krylov methods applied for approximating the resolvent and exponential function, and functions of Markov type, is treated in some detail. (© 2013 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)

Efficient mixed rational and polynomial approximation of matrix functions

This paper presents an efficient method for computing approximations for general matrix functions based on mixed rational and polynomial approximations. A method to obtain this kind of approximation from rational approximations is given, reaching the highest efficiency when transforming nondiagonal rational approximations with a higher numerator degree than the denominator degree. Then, the proposed mixed rational and polynomial approximation can be successfully applied for matrix functions which have any type of rational approximation, such as Padé, Chebyshev, etc., with maximum efficiency for higher numerator degrees than the denominator degrees. The efficiency of the mixed rational and polynomial approximation is compared with the best existing evaluating schemes for general polynomial and rational approximations, providing greater theoretical accuracy with the same cost in terms of matrix multiplications. It is well known that diagonal rational approximants are generally more accurate than the corresponding nondiagonal rational approximants which have the same computational cost. Using the proposed mixed approximation we show that the above statement is no longer true, and nondiagonal rational approximants are in fact generally more accurate than the corresponding diagonal rational approximants with the same cost.

Superlinear convergence of the rational Arnoldi method for the approximation of matrix functions

A superlinear convergence bound for rational Arnoldi approximations to functions of matrices is derived. This bound generalizes the well-known superlinear convergence bound for the conjugate gradient method to more general functions with finite singularities and to rational Krylov spaces. A constrained equilibrium problem from potential theory is used to characterize a max-min quotient of a nodal rational function underlying the rational Arnoldi approximation, where an additional external field is required for taking into account the poles of the rational Krylov space. The resulting convergence bound is illustrated at several numerical examples, in particular, the convergence of the extended Krylov method for the matrix square root.

A restarted Lanczos approximation to functions of a symmetric matrix

In this paper the error term in the 'first' method of Eiermann & Ernst (2006, SIAM J. Numer. Anal., 44, 2481-2504) for restarting the Lanczos approximation of the matrix-vector product f (A)b, where A ∈ ℝ n×n is symmetric, is re-derived and expressed as an explicit partial fraction expansion of the divided differences. The partial fraction representation makes the new variant slightly more stable (albeit still unstable) than the former method because it requires fewer finite-difference evaluations. We then present an error bound for the restarted Lanczos approximation of f (A)b for symmetric positive-definite A when f is in a particular class of completely monotone functions and illustrate for some important matrix function applications the usefulness of these bounds for terminating the restart process once the desired accuracy in the matrix function approximation has been achieved. Finally, in an attempt to overcome the inherent instability of our restart procedure, we propose a simple heuristic that identifies when to halt the iterations.

Construction of Data-Sparse $\mathcal{H}^2$-Matrices by Hierarchical Compression

Hierarchical matrices ($\mathcal{H}$-matrices) provide an elegant approach to handling large densely populated matrices: The matrix is split into a hierarchy of blocks, and each block is approximated by a low-rank matrix in factorized form. It has been demonstrated that this representation can be used to treat integral and partial differential equations, solve matrix equations from the field of control theory, and evaluate matrix functions efficiently. $\mathcal{H}^2$-matrices use a refined representation that employs a multilevel structure in order to reduce the storage requirements of hierarchical matrices. It has been shown that $\mathcal{H}^2$-matrices can significantly reduce storage requirements for large problems, in particular when combined with modern error control schemes. Until now, all algorithms for constructing an efficient approximation of a general matrix by an $\mathcal{H}^2$-matrix required a representation of the entire original matrix to be kept in storage; therefore the storage requirements of $\mathcal{H}^2$-matrix algorithms could be far larger than those of the final approximation. This paper presents a new approach that allows us to construct an $\mathcal{H}^2$-matrix without storing the entire original matrix. The central idea is to approximate submatrices and combine them by an efficient new algorithm to form approximations of larger matrices until the entire matrix has been treated. Using this new approach, many $\mathcal{H}$-matrix algorithms can be “refitted” easily to compute results in the more efficient representation. Possible applications include efficient matrix arithmetics for the construction of preconditioners, the approximation of matrix functions or solutions of matrix equations, or efficient compression schemes based on the popular cross approximation algorithms.

Analytic approximation of matrix functions in [formula omitted

We consider the problem of approximation of matrix functions of class L p on the unit circle by matrix functions analytic in the unit disk in the norm of L p , 2 ≤ p < ∞ . For an m × n matrix function Φ in L p , we consider the Hankel operator H Φ : H q ( C n ) → H − 2 ( C m ) , 1 / p + 1 / q = 1 / 2 . It turns out that the space of m × n matrix functions in L p splits into two subclasses: the set of respectable matrix functions and the set of weird matrix functions. If Φ is respectable, then its distance to the set of analytic matrix functions is equal to the norm of H Φ . For weird matrix functions, to obtain the distance formula, we consider Hankel operators defined on spaces of matrix functions. We also describe the set of p -badly approximable matrix functions in terms of special factorizations and give a parametrization formula for all best analytic approximants in the norm of L p . Finally, we introduce the notion of p -superoptimal approximation and prove the uniqueness of a p -superoptimal approximant for rational matrix functions.

Analytic approximation of matrix functions and dual extremal functions

We study the question of the existence of a dual extremal function for a bounded matrix function on the unit circle in connection with the problem of approximation by analytic matrix functions. We characterize the class of matrix functions, for which a dual extremal function exists in terms of the existence of a maximizing vector of the corresponding Hankel operator and in terms of certain special factorizations that involve thematic matrix functions.

A Numerical Solution Using an Adaptively Preconditioned Lanczos Method for a Class of Linear Systems Related with the Fractional Poisson Equation

This study considers the solution of a class of linear systems related with the fractional Poisson equation (FPE) (−∇2)α/2φ=g(x,y) with nonhomogeneous boundary conditions on a bounded domain. A numerical approximation to FPE is derived using a matrix representation of the Laplacian to generate a linear system of equations with its matrix A raised to the fractional power α/2. The solution of the linear system then requires the action of the matrix function f(A)=A−α/2 on a vector b. For large, sparse, and symmetric positive definite matrices, the Lanczos approximation generates f(A)b≈β0Vmf(Tm)e1. This method works well when both the analytic grade of A with respect to b and the residual for the linear system are sufficiently small. Memory constraints often require restarting the Lanczos decomposition; however this is not straightforward in the context of matrix function approximation. In this paper, we use the idea of thick-restart and adaptive preconditioning for solving linear systems to improve convergence of the Lanczos approximation. We give an error bound for the new method and illustrate its role in solving FPE. Numerical results are provided to gauge the performance of the proposed method relative to exact analytic solutions.

Approximating functions of a large sparse positive definite matrix using a spectral splitting method

The computation of functions of large sparse matrices f(A) is an important topic in numerical linear algebra and finds application in many fields of applied mathematics and statistics. In previous research we considered ? matrices with compact spectrum ?( A ) ? [a,b] and proposed low degree matrix polynomial approximations p( A ) such that e = ?f( A ) ? p( A ) ? was small on the spectral interval, where the extreme eigenvalues a and b were calculated using Krylov subspace approximation. For the class of matrices examined, the thick restarted Lanczos scheme enabled rapid convergence to the extreme eigenvalues and these Ritz values were used to construct cubic near-minimax Chebyshev least squares approximations of the desired matrix functions. There is a good balance between accuracy and efficiency for this approximation method. The aim of the present study is to extend the previously developed matrix function approximation technique to enable ? matrices with a wider spectrum to be treated using a novel splitting of ?( A ). In this case, the decomposition of f( A ) as a sum of a 'singular' part and a 'regular' part is investigated. To perform the split a projector onto the singular part is here constructed using Krylov subspace approximation. Numerical results for a representative large sparse positive definite matrix appear promising.

QRT: A QR-Based Tridiagonalization Algorithm for Nonsymmetric Matrices

The stable similarity reduction of a nonsymmetric square matrix to tridiagonal form has been a long-standing problem in numerical linear algebra. The biorthogonal Lanczos process is in principle a candidate method for this task, but in practice it is confined to sparse matrices and is restarted periodically because roundoff errors affect its three-term recurrence scheme and degrade the biorthogonality after a few steps. This adds to its vulnerability to serious breakdowns or near-breakdowns, the handling of which involves recovery strategies such as the look-ahead technique, which needs a careful implementation to produce a block-tridiagonal form with unpredictable block sizes. Other candidate methods, geared generally towards full matrices, rely on elementary similarity transformations that are prone to numerical instabilities. Such concomitant difficulties have hampered finding a satisfactory solution to the problem for either sparse or full matrices. This study focuses primarily on full matrices. After outlining earlier tridiagonalization algorithms from within a general framework, we present a new elimination technique combining orthogonal similarity transformations that are stable. We also discuss heuristics to circumvent breakdowns. Applications of this study include eigenvalue calculation and the approximation of matrix functions.