Source-linked AI summary
A Krylov subspace algorithm for evaluating the phi-functions appearing in exponential integrators
Jitse Niesen, Will M. Wright
TL;DR
Exponential integrators repeatedly require costly actions of matrix exponentials and phi-functions for large sparse systems. The paper develops phipm, which combines Arnoldi or Lanczos Krylov projection with time-stepping and adaptivity in both step size and subspace dimension. Comparisons show that adapting both dimensions improves efficiency, often substantially, although stability behavior and boundary-condition issues remain in the reported setting.
Problem
Evaluating linear combinations of matrix exponentials and phi-functions repeatedly and accurately is a major computational burden in exponential integrators.
Method
phipm uses Arnoldi or Lanczos Krylov projection, augmented-matrix exponentials, time-stepping, and adaptive step sizes and Krylov dimensions.
Results
Adapting both the Krylov-subspace dimension and time-step length significantly increases overall efficiency, often producing very large computational gains.
Takeaways & Limitations
The phipm routine provides a reusable efficient kernel for computing linear combinations of phi-functions in exponential integrators.
Takeaways & Limitations
The reported Heston experiment shows zig-zagging step sizes and variable error estimates that may indicate stability problems, while its boundary and initial conditions are non-matching.
Abstract
from arXiv · showhide
We develop an algorithm for computing the solution of a large system of linear ordinary differential equations (ODEs) with polynomial inhomogeneity. This is equivalent to computing the action of a certain matrix function on the vector representing the initial condition. The matrix function is a linear combination of the matrix exponential and other functions related to the exponential (the so-called phi-functions). Such computations are the major computational burden in the implementation of exponential integrators, which can solve general ODEs. Our approach is to compute the action of the matrix function by constructing a Krylov subspace using Arnoldi or Lanczos iteration and projecting the function on this subspace. This is combined with time-stepping to prevent the Krylov subspace from growing too large. The algorithm is fully adaptive: it varies both the size of the time steps and the dimension of the Krylov subspace to reach the required accuracy. We implement this algorithm in the Matlab function phipm and we give instructions on how to obtain and use this function. Various numerical experiments show that the phipm function is often significantly more efficient than the state-of-the-art.
1. Introduction
Exponential integrators target ODEs split into stiff linear and non-stiff nonlinear parts, but their computational cost is dominated by evaluating matrix exponentials and related phi-functions. The paper develops an adaptive Krylov-subspace solver combining time-stepping with variable subspace dimensions.
- Exponential integrators split ODEs into stiff linear and non-stiff nonlinear parts and use matrix exponentials and phi-functions within the integrator.
- Evaluating phi-functions is the dominant computational cost in exponential integrators.
- Large sparse matrices arising from semidiscretized PDEs make Krylov methods effective for evaluating matrix-function actions without forming the matrix function.
- Earlier software and methods applied Krylov approximations and time-stepping to the matrix exponential, first phi-functions, and polynomial inhomogeneities.
- The proposed solver combines time-stepping with adaptive Krylov dimensions, excludes small dense matrices, and is intended as a reusable kernel for exponential integrators.
- Rational Krylov, Leja-point, Laplace-inversion, and contour-integration approaches are identified as alternatives outside this paper’s scope.
2. The ϕ-functions
Exponential-integrator stages reduce to linear combinations of phi-functions acting on vectors, so the paper focuses on evaluating these expressions efficiently and accurately for large sparse matrices. A linear-ODE result with polynomial inhomogeneity motivates their role and the later time-stepping procedure.
- The scalar phi-function definition is extended to matrices using available matrix-function definitions.
- Every exponential-integrator stage can be expressed as a linear combination of phi-functions acting on vectors.
- The matrix A is often a large sparse Jacobian or an approximation to one, while the polynomial degree p is typically below five.
- The paper addresses repeated, accurate evaluation of these matrix-function expressions, which can be reused across different exponential integrators.
- A lemma gives the exact solution of linear differential equations with polynomial inhomogeneity and provides background for the time-stepping evaluation procedure.
3. The algorithm
The algorithm projects phi-function actions for large matrices onto Krylov subspaces, computes the reduced functions through an augmented matrix exponential, and adapts both error control and computational dimensions. Time-stepping limits subspace growth, while experiments motivate adaptive choices but also reveal unresolved stability and precision questions.
- 3.1. The basic method.: Arnoldi constructs an orthonormal Krylov basis because successive vectors Ajv are nearly dependent; for symmetric A, Lanczos produces a tridiagonal reduced matrix.
- 3.1. The basic method.: Krylov projection approximates ϕp(A)v by βVmϕp(Hm)e1, replacing an n-by-n matrix-function action with a smaller m-by-m computation.Here β = ∥v∥, and Hm is the projected Hessenberg matrix.
- 3.1. The basic method.: The augmented matrix ˆHm allows ϕp(Hm)e1 to be obtained from the top m entries of the last column of exp(ˆHm), computed with Padé approximation, scaling, and squaring.
- 3.1. The basic method.: Computing only the needed exponential column or exploiting Hessenberg structure could reduce cost, but these savings are left for future algorithmic work.
- 3.2. Error estimation.: The Krylov error estimate is computed from a series term already available from Arnoldi and is also used as a corrector for adaptivity.
- 3.3. Time-stepping.: Time-stepping is used when large spectral radius would otherwise require a large Krylov dimension, while the implementation adapts step size and subspace dimension.
- 3.4. Adaptivity.: For finite step sizes, the error follows a heuristic order q smaller than the asymptotic order m, and experiments confirm proportionality to τ^(m+1) as τ approaches zero.
4. Numerical experiments
The experiments evaluate phipm on a discretized Heston PDE system and compare it with Crank–Nicolson, two ADI schemes, and ode15s. They also examine phipm's adaptive Krylov dimension, time stepping, and error estimation.
- Heston equation in financial mathematics: The comparison includes Crank–Nicolson, first- and second-order ADI schemes, ode15s with its Jacobian supplied, and phipm.Crank–Nicolson and the ADI methods use decreasing step sizes, while phipm and ode15s use geometrically decreasing tolerances.
- Heston equation in financial mathematics: The Heston test discretizes a European call-option pricing problem into a sparse system with 5,100 degrees of freedom and 44,800 nonzero elements.The grid has 100 points in the asset-price direction and 51 in the variance direction.
- Heston equation in financial mathematics: Crank–Nicolson outperforms the ADI methods, while phipm is most efficient at an accuracy of around 10^-6.The reported comparison uses maximum error against CPU time; similar behavior is reported for four additional parameter sets, with the crossover near 10^-6.
- Heston equation in financial mathematics: During integration at tolerance 10^-4, phipm varies both the Krylov subspace dimension and step size, with rejected steps explicitly identified.These variations are presented as evidence that the algorithm's adaptivity is effective.
- Heston equation in financial mathematics: At tolerance 10^-4, the estimated error is always larger than the actual error for the accepted steps in this experiment.The authors nevertheless note zig-zagging step sizes and substantial step-to-step variation in the estimated error as a possible stability concern.
- Heston equation in financial mathematics: The runtime comparison favors methods implemented in optimized low-level routines over phipm's native Matlab computations as a relevant implementation consideration.The authors state that ode15s, Crank–Nicolson, and the ADI methods use highly optimized low-level linear algebra, whereas phipm uses native Matlab.
Experiment 1.
Experiment 1 compares phiv, phip, and phipm with expv on four sparse matrices, measuring average speedup and relative error. The tests cover symmetric and unsymmetric matrices, including large discretized PDE systems.
- The experiment evaluates matrix-exponential and phi-function routines on four sparse matrices using specified vectors, times, and tolerances.The matrices include Harwell–Boeing and GHS indef examples, with dimensions ranging from 900 to 392,257.
- The third example tests a forward computation followed by a reverse computation, e^-tA e^tA b0, against the known solution b0.This provides the only example with a known exact solution and uses Tol = 10^-14.
- The final example evaluates e^tA b0 + tϕ1(tA)b1 on a large symmetric sparse Helmholtz-equation matrix.For this test, only phiv, phip, and phipm are compared.
- Average speedups and relative errors are measured against expv for the first three comparisons and against phiv for the final comparison.Timings are averaged over 100 runs and summarized in Table 1.
Experiment 2.
Experiment 2 isolates the efficiency gained by allowing the Krylov subspace dimension to vary. It compares phipm with phip on the same sparse-matrix test set and summarizes speedups and errors in Table 2.
- The comparison evaluates ϕ0(tA)b0 + tϕ1(tA)b1 + ··· + t4ϕ4(tA)b4 with phip and phipm.All operand vectors are set to the all-ones vector, using the same matrix, time, and tolerance settings as Experiment 1 except for the stated restriction.
- The implementations are identical except that phipm can vary the Krylov subspace dimension m.This design directly gauges the effect of adaptive subspace size.
- The comparison uses four large sparse matrices from the University of Florida sparse matrix collection, with the gr example restricted to the forward computation.
- Table 2 reports phipm’s average speedup relative to phip together with relative errors.
Discussion of the results.
The comparisons report that phipm is more efficient in every tested case, and that jointly adapting Krylov dimension and time-step length improves overall efficiency.
- phipm is more efficient than the comparison codes in all tested cases, sometimes by a considerable margin.
- Adapting both the Krylov subspace dimension and time-step length significantly increases overall efficiency.Because phipm may be called repeatedly throughout exponential-integrator simulations, the paper reports that this can produce large overall computational gains.
5. Conclusion and future work
The conclusion presents phipm as an efficient adaptive routine for linear combinations of phi-functions, while identifying stability, implementation, and method-comparison work for the future.
- Conclusion: phipm computes the action of linear combinations of ϕ-functions on operand vectors by combining time stepping with adaptive Krylov-subspace sizing.
- Conclusion: The routine is positioned as an extension of codes available in expokit and Mathematica.
- Future work: The paper notes that the ϕ-functions are building blocks of exponential integrators and that a lower-level implementation would be useful.Implementation of exponential integrators using phipm is also identified as ongoing work.
- Future work: Planned code improvements include investigating stability concerns associated with error estimates and computing the reduced-matrix ϕ-function through augmentation.The authors also suggest exploiting symmetry, Hessenberg structure, and methods such as rational Chebyshev approximants.
- Future work: Future work includes GPU-based parallelization, Krylov methods for option-pricing problems, and comparisons with alternative methods.The alternatives include restricted-denominator rational Krylov, real Leja, quadrature, and contour-integration approaches.