We present a novel Krylov subspace method for approximating $L_f(A, E) \mathbf{b}$, the matrix-vector product of the Fréchet derivative $L_f(A, E)$ of a large-scale matrix function $f(A)$ in direction $E$, a
task that arises naturally in the sensitivity analysis of quantities involving matrix functions such as centrality measures for networks.
It also arises in the context of gradient-based methods for optimization problems that feature matrix functions, e.g., when fitting an evolution equation to an observed solution trajectory.
In principle, the well-known identity
\[
f\left( \begin{bmatrix}
A & E \\\\
0 & A
\end{bmatrix} \right) \begin{bmatrix}
0 \\\\
\mathbf{b}
\end{bmatrix} = \begin{bmatrix}
L_f(A, E) \mathbf{b} \\\\
f(A) \mathbf{b}
\end{bmatrix}
\]
allows one to directly apply any standard Krylov subspace method, such as the Arnoldi algorithm, to address this task.
However, this comes with the major disadvantage that the involved block-triangular matrix has unfavorable spectral properties, which impede the convergence analysis and, to a certain extent, also the observed convergence. To avoid these difficulties, we propose a novel modification of the Arnoldi algorithm that aims at better preserving the block-triangular structure. In turn, this allows one to bound the convergence of the modified method by the best polynomial approximation of the derivative $f^\prime$ on the numerical range of $A$. Several numerical experiments illustrate our findings.