A novel Krylov subspace method for approximating Fréchet derivatives of large-scale matrix functions

Technical Summary: A Novel Krylov Subspace Method for Approximating Fréchet Derivatives of Large-Scale Matrix Functions

Problem Statement
The paper addresses the computational challenge of approximating the matrix-vector product $L_f(A, E)b$, where $L_f(A, E)$ is the Fréchet derivative of a large-scale matrix function $f(A)$ in the direction $E$, applied to a vector $b$. This task is critical in sensitivity analysis for network centrality measures and in gradient-based optimization problems involving matrix functions (e.g., fitting evolution equations).

Direct computation of the full $n^2 \times n^2$ matrix representation of the Fréchet derivative is infeasible for large $n$. While standard Krylov subspace methods exist for approximating $f(A)b$, their application to $L_f(A, E)b$ faces significant theoretical and practical hurdles. The standard approach utilizes the identity:
$$f\left(\begin{bmatrix} A & E \\ 0 & A \end{bmatrix}\right) \begin{bmatrix} 0 \\ b \end{bmatrix} = \begin{bmatrix} L_f(A, E)b \\ f(A)b \end{bmatrix}$$
This allows the use of standard algorithms (like Arnoldi) on the block triangular matrix $\mathcal{A} = \begin{bmatrix} A & E \\ 0 & A \end{bmatrix}$. However, the authors note that $\mathcal{A}$ typically possesses unfavorable spectral properties (e.g., Jordan blocks of size $\ge 2$ even if $A$ is symmetric), which impede convergence analysis and observed convergence. Alternative approaches, such as finite differences or complex step approximations combined with Arnoldi, suffer from numerical cancellation errors or the introduction of complex arithmetic for real data.

Methodology
The authors propose a novel modification of the Arnoldi algorithm designed to preserve the block triangular structure of the embedding matrix $\mathcal{A}$ during the subspace construction.

1. Structure-Preserving Basis Construction: Instead of applying Arnoldi directly to the full $2n \times 2n$ matrix $\mathcal{A}$, the method constructs a basis for the Krylov subspace $K_k(\mathcal{A}, [0, b]^T)$ by splitting the basis vectors into top and bottom $n$-dimensional components.
2. Modified Arnoldi Algorithm (Basic Version): The algorithm maintains two orthonormal bases, $U_k$ and $V_k$, corresponding to the top and bottom blocks, respectively. The compressed matrix $U_k^* \mathcal{A} U_k$ retains the block triangular form $\begin{bmatrix} U_k^* A U_k & U_k^* E V_k \\ 0 & V_k^* A V_k \end{bmatrix}$.
3. Separate Orthonormalization (Efficient Version): To avoid the computational inefficiency of re-orthogonalizing split vectors, the authors introduce a recursive procedure (Algorithm 2). This method constructs non-orthonormal Krylov bases of the form $\begin{bmatrix} U_i R_i \\ V_i \end{bmatrix}$, where $V_i$ and $U_i$ are orthonormal, and $R_i$ is strictly upper triangular. This approach performs Gram-Schmidt orthogonalization separately on the bottom components ($Av_{i-1}$) and the top components ($Aw_{i-1} + Ev_{i-1}$), updating the triangular factor $R_i$ to maintain the basis property without destroying orthogonality.
4. Convergence Analysis: Theoretical analysis (Theorem 2.2) demonstrates that the error in approximating $L_f(A, E)b$ is bounded by the best polynomial approximation of the scalar derivative $f'$ on the numerical range $W(A)$ of the original matrix $A$. This parallels existing results for $f(A)b$ but relies on $f'$ rather than $f$, and crucially, depends on the properties of $A$ rather than the ill-conditioned $\mathcal{A}$.

Key Contributions

• Algorithmic Innovation: The proposal of a structure-preserving modification to the Arnoldi algorithm that avoids the spectral pitfalls of the standard block embedding approach.
• Theoretical Justification: A convergence bound showing that the method's performance is governed by the approximation of $f'$ on the numerical range of $A$, rather than the numerical range of the enlarged matrix $\mathcal{A}$.
• Efficiency Improvements: The development of a "Separate Orthonormalization" variant that reduces computational overhead by eliminating redundant orthogonalizations and extra matrix-vector products, at the cost of modest additional memory ($O(nk)$).
• Handling Low-Rank and General Directions: The method is applicable to both low-rank directions (where specialized tensor methods exist) and general directions, bridging a gap in existing literature.

Numerical Results
The authors validate their method through two primary sets of experiments:

1. Network Centrality Sensitivity: The algorithm is applied to compute sensitivity measures (Total Network Communicability, Subgraph Centrality, Estrada Index) for various real-world networks (Air500, Autobahn, US Power Grid, etc.).
   • Convergence: In all tested cases, the Modified Arnoldi method converges faster than the standard Arnoldi method applied to $\mathcal{A}$, finite difference approaches, and complex step approaches. Its convergence rate is comparable to the Arnoldi method for approximating $f(A)b$.
   • Performance: The runtime of the proposed method is comparable to the standard Arnoldi method for $A$ and significantly faster than competing methods (e.g., De La Cruz Cabrera et al., Schweitzer).
   • Robustness: The method remains effective even when the direction matrix $E$ is not low-rank, a scenario where other specialized low-rank methods fail.
2. Parameter Fitting for Heat Equation: The method is used in a gradient descent optimization loop to fit a parameter $\sigma$ in a 2D parametric heat equation.
   • The algorithm successfully computes the matrix function action and its Fréchet derivative simultaneously, enabling efficient gradient evaluation without finite difference approximations.
   • The optimization converges to the true parameter value ($\sigma = 0.85$) with high accuracy.

Significance and Claims
The paper claims that the proposed modification effectively mitigates the theoretical difficulties associated with the unfavorable spectral properties of the block triangular embedding matrix. By preserving the block structure and linking convergence to the numerical range of the original matrix $A$, the method offers a robust and efficient alternative for sensitivity analysis and optimization involving matrix functions. The authors emphasize that while existing methods show convergence in practice, their theoretical bounds are often poor; the new method provides both improved practical convergence and a rigorous theoretical foundation based on the properties of $A$. The work is presented as a practical advancement for large-scale problems where full derivative computation is impossible, specifically targeting applications in network science and parameter estimation.