text:
\section{Introduction}
Krylov subspace methods (KSMs) are among the most powerful algorithms for computing approximations to $f(\vec{A})\vec{b}$ when $\vec{A}$ is a $n\times n$ real symmetric matrix and $f:\mathbb{R}\to\mathbb{R}$ is an arbitrary function.
Such methods construct an approximation to $f(\vec{A})\vec{b}$ that lies in the Krylov subspace
\begin{align*}
\mathcal{K}_k(\vec{A},\vec{b})
:= \operatorname{span}\{\vec{b}, \vec{A}\vec{b}, \ldots, \vec{A}^{k-1}\vec{b} \},
\end{align*}
and, as a result, only need to access $\vec{A}$ through matrix-vector products.
This means KSMs are well suited for large-scale computations where storing $\vec{A}$ in fast memory is infeasible.

In the special case that $f(x) = 1/(x-z)$ for some $z\in\mathbb{C}$, KSMs such as conjugate gradient (CG) \cite{hestenes_stiefel_52}, minimum residual (MINRES) \cite{paige_saunders_75}, and quasi-minimum residual (QMR) \cite{freund_92} are able to provide \emph{optimal} approximations to $f(\vec{A})\vec{b}$ from $\mathcal{K}_k(\vec{A},\vec{b})$ while storing just a few vectors of length $n$.
In some cases, these algorithms can be leveraged to approximate other matrix functions.
For instance, if $f$ is Stieltjes, then $f$ can be efficiently approximated by a rational function of the form $\sum_i \omega_i/(x-z_i)$ where $z_i < 0$.
Thus, if $\vec{A}$ is positive definite, then $f(\vec{A})\vec{b}$ can be approximated by using CG to compute each term.
Moreover, the shift invariance of Krylov subspace enables the same Krylov subspace to be used for all terms, so only $k$ matrix-vector products are needed in total to obtain approximations from $\mathcal{K}_k$
\cite{eshof_frommer_lippert_schilling_van_der_vorst_02,frommer_simoncini_08,frommer_simoncini_08,guttel_schweitzer_21}.
A similar algorithm involving MINRES can be used in the case of indefinite systems \cite{pleiss_jankowiak_eriksson_damle_gardner_20}.

For general functions $f$ the situation is murkier. General purpose KSMs, like the well-known Lanczos method for matrix function approximation (Lanczos-FA) \cite{saad_92}, cannot be guaranteed to return an optimal or new-optimal approximation to $f(\vec{A})\vec{b}$ from $\mathcal{K}_k$, except in the case of particular functions such as the inverse or exponential \cite{druskin_greenbaum_knizhnerman_98}.
Work on more general spectrum dependent bounds remains somewhat ad-hoc, even for the basic case of rational functions \cite{eshof_frommer_lippert_schilling_van_der_vorst_02,ilic_turner_simpson_09,frommer_guttel_schweitzer_14,frommer_guttel_schweitzer_14a,chen_greenbaum_musco_musco_21}.
Moreover, in terms of computational cost, to return an approximation to $f(\vec{A})\vec{b}$ from $\mathcal{K}_k$ methods for general functions like Lanczos-FA require either (i) storing $k$ vectors of length $n$, or
(ii) storing a constant number of vectors of length $n$, but increasing the number of matrix-vector products by a factor of two.

In this paper, we describe an algorithm that resolves these two main issues (lack of optimality and poor memory performance) for the imporant case of rational function. We call the method the Lanczos method for optimal rational matrix function approximation (Lanczos-OR).
When $f$ is a rational function,
Lanczos-OR outputs the \emph{optimal} (in a certain norm) approximation to $f(\vec{A})\vec b$ in the span of $\mathcal{K}_k(\vec{A},\vec{b})$ using slightly more than $k$ matrix-vector products.
We provide a practical implementation of Lanczos-OR that
only requires storing a number of vectors of length $n$ proportional to the degree of the denominator in the rational function.
Therefore, for a fixed rational function, the storage costs do not grow with the iteration $k$. The approach used to derive this implementation of Lanczos-OR can also be used for computing the Lanczos-FA approximations to rational matrix functions, avoiding storage costs growing with $k$ in that widely used method.

Beyond rational functions, we show that Lanczos-OR can be used to derive algorithms for other functions, including the matrix sign function, by discretizing an integral representation of the function.
The induced approximations are closely related to the corresponding approximations generated by the Lanczos method for matrix function approximation (Lanczos-FA).
Like the Lanczos-FA approximations, the Lanczos-OR induced approximations are not guaranteed to be optimal, but tend to perform well in practice.
Moreover, on problems where Lanczos-FA exhibits erratic behavior, the Lanczos-OR induced approximations tend to have nicer behavior.

Lanczos-OR is closely related to existing methods for linear systems.
In particular, if $f(x) = 1/(x-z)$, then, depending on the choice of $z$, the CG, MINRES, and QMR iterates can all be obtained as special cases of Lanczos-OR.
In fact, the Lanczos-OR based approach for the matrix sign function and rational matrix functions is closely related to using QMR to approximate solutions to certain pairs of linear systems.
More generally, the Lanczos-OR iterate is mathematically equivalent to an optimal Galerkin projection method as described in \cite[Section 4]{lopez_simoncini_06}, provided the denominator matrix is positive definite.
However, this method was mostly viewed as of theoretical interest since it could be used to help explain the behavior of Lanczos-FA.
On the other hand, we show that such approximations can be computed efficiently.
Our approach is also somewhat more general in that it works with \emph{any} rational function.

\subsection{Notation}
For integers $a<a'$, the notation $a:a'$ indicates the set of integers $\{ a, a+1, \ldots, a'-1 \}$.
When used in a loop, $k=a:a'$ indicates that $k$ takes values $a, a+1, \ldots, a'-1$ sequentially.
We denote the complex conjugate of $x$ by $\overline{x}$.

Throughout, matrices are denoted by bold uppercase letters and vectors are denoted by bold lowercase letters.
We use \emph{zero-indexed numpy style slicing} to indicate entries of matrices and vectors.
Thus, the submatrix consisting of rows \( r \) through \( r'-1 \) and columns \( c \) through \( c'-1 \) is denoted by \( [\vec{B}]_{r:r',c:c'} \).
If any of these indices are equal to \( 0 \) or \( n \), they may be omitted and if \( r' = r+1  \) or \( c' = c+1 \), then we will simply write \( r \) or \( c \).
Thus, \( [\vec{B}]_{:,:2} \) denotes the first two columns of \( \vec{B} \), and \( [\vec{B}]_{3,:} \) denotes the fourth row (corresponding to index 3) of \( \vec{B} \).

Throughout, $\vec{A}$ will be a real symmetric matrix.
We denote the set of eigenvalues of $\vec{A}$ by $\Lambda(\vec{A})$ and define $\mathcal{I}(\vec{A}) := [\lambda_{\textup{min}}(\vec{A}), \lambda_{\textup{max}}(\vec{A})]$.
Without loss of generality, we assume that $\|\vec{b}\|_2 = 1$, where $\vec{b}$ is the vector to which the rational function is applied.

\section{The Lanczos algorithm and some basic Krylov subspace methods}
\label{sec:background}

KSMs for symmetric matrices  are often based on the Lanczos algorithm, which we describe in \cref{alg:lanczos}.
Given a symmetric matrix $\vec{A}$ and vector $\vec{b}$, the Lanczos algorithm run for $k$ iterations will construct an orthonormal basis \( \vec{Q} := [ \vec{q}_0 , \ldots , \vec{q}_{k-1} ] \) such that the first $j\leq k$ columns form a basis for the Krylov subspace
\begin{align*}
\mathcal{K}_j(\vec{A},\vec{b})
= \operatorname{span}( \vec{b}, \vec{A}\vec{b}, \ldots, \vec{A}^{j-1}\vec{b} )
= \{ p(\vec{A}) \vec{b} : \deg(p) < j \}.
\end{align*}
Moreover, the basis vectors satisfy a symmetric three term recurrence,
$$
\vec{A} \vec{Q} = \vec{Q} \vec{T} + \beta_{k-1} \vec{q}_{k} \vec{e}_k^\textup{\textsf{T}}.
$$
Above $\vec{e}_k$ is the $k^\text{th}$ standard basis vector and $\vec{T}$ is symmetric tridiagonal with  diagonals $(\alpha_0, \ldots, \alpha_{k-1})$ and off diagonals $(\beta_0, \ldots, \beta_{k-2})$  which are also computed by the algorithm. I.e.,
\begin{align*}
    {\vec{T}}
    := \begin{bmatrix}
        \alpha_0 & \beta_0 \\
        \beta_0 & \alpha_1 & \ddots \\
        & \ddots & \ddots & \beta_{k-2}\\
        & & \beta_{k-2} & \alpha_{k-1}
    \end{bmatrix}.
\end{align*}

Note that in exact arithmetic, the reorthogonalization step in line \ref{alg:reorth_line} of \cref{alg:lanczos}
is unnecessary as $\tilde{\vec{q}}_{j+1}$  is automatically orthogonal to $\{\vec{q}_i\}_{i=0}^{j-1}$.
However, in finite precision arithmetic this is no longer the case and the algorithm can behave very differently depending on whether or not reorthogonalization is used.
See \cref{sec:stability} for further discussion.

\begin{algorithm}
\caption{Lanczos}\label{alg:lanczos}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Procedure{Lanczos}{$\vec{A}, \vec{b}, k$}
\State \( \vec{q}_{-1} = \vec{0} \),
\( \beta_{-1} = 0 \),
\( \vec{q}_0  = \vec{b} / \| \vec{b} \|_2 \)
\For {\( j=0:k-1 \)}
    \State \( \tilde{\vec{q}}_{j+1} = \vec{A} \vec{q}_{j} - \beta_{j-1} \vec{q}_{j-1} \)
    \State \( \alpha_j = \langle \tilde{\vec{q}}_{j+1}, \vec{q}_j \rangle \)
    \State \( \tilde{\vec{q}}_{j+1} = \tilde{\vec{q}}_{j+1} - \alpha_j \vec{q}_j \)
        \State optionally, reorthogonalize \( \tilde{\vec{q}}_{j+1} \) against \( \{\vec{q}_i\}_{i=0}^{j-1} \) \label{alg:reorth_line}
    \State \( \beta_{j} = \| \tilde{\vec{q}}_{j+1} \|_2 \)
    \State \( \vec{q}_{j+1} = \tilde{\vec{q}}_{j+1} / \beta_{j} \)
\EndFor
\State \Return \( \{ (\vec{q}_j, \alpha_j, \beta_j ) \}_{j=0}^{k-1} \) and $\vec{q}_{k}$
\EndProcedure
\end{algorithmic}
\end{algorithm}

To simplify analysis, it will be useful to consider the recurrence that would obtained if the Lanczos algorithm were run to completion.
In exact arithmetic, for some $K\leq n$, $\beta_{K-1} = 0$ in which case the algorithm terminates.
Then the final basis \( \widehat{\vec{Q}} := [ \vec{q}_0 , \ldots , \vec{q}_{K-1} ]\) and symmetric tridiagonal $\widehat{\vec{T}}$ with  diagonals $(\alpha_0, \ldots, \alpha_{K-1})$ and off diagonals $(\beta_0, \ldots, \beta_{K-2})$
satisfy
 a three-term recurrence
\begin{align*}
    \vec{A} \widehat{\vec{Q}} = \widehat{\vec{Q}} \widehat{\vec{T}}.
\end{align*}

We emphasize that our algorithms \emph{do not} require Lanczos to be run to completion; the introduction of $\widehat{\vec{Q}}$ and $\widehat{\vec{T}}$ is for analysis purposes only.
We note that $\vec{Q} = [\widehat{\vec{Q}}]_{:,:k}$ and $\vec{T} =[\widehat{\vec{T}}]_{:k,:k}$.

Since the columns of $\widehat{\vec{Q}}$ are orthogonal, we have that
\begin{align*}
    \widehat{\vec{T}}
    = \widehat{\vec{Q}}^\textup{\textsf{T}} \vec{A} \widehat{\vec{Q}},
\end{align*}
from which we easily see that, after any number of iterations $k$, $\vec{T} = \vec{Q}^\textup{\textsf{T}} \vec{A} \vec{Q}$.
Note also that, for any shift $z\in\mathbb{C}$,
\begin{align*}
    (\vec{A} - z\vec{I}) \widehat{\vec{Q}} = \widehat{\vec{Q}} (\widehat{\vec{T}} - z\vec{I}).
\end{align*}
In other words, the Krylov subspaces generated by $(\vec{A},\vec{b})$ and $(\vec{A}-z\vec{I},\vec{b})$ coincide, and the associated tridiagonal matrices are easily related by a diagonal shift.
This shift invariance of Krylov subspaces is critical in several low-memory algorithms.

\subsection{The conjugate gradient and minimum residual algorithms}
\label{sec:linear_solvers}

CG, MINRES, and QMR are all algorithms for solving linear systems of equations $\vec{A}\vec{x} = \vec{b}$.
This is equivalent to approximating $\vec{A}^{-1} \vec{b}$; i.e., $f(\vec{A})\vec{b}$ where $f(x) = 1/x$.

It is well known that when $\vec{A}$ is positive definite, CG minimizes the $\vec{A}$-norm of the error over the Krylov subspace; i.e., the CG approximation $\mathsf{cg}_k$ is given by
\begin{align*}
    \mathsf{cg}_k := \argmin_{\vec{x}\in\mathcal{K}_k(\vec{A},\vec{b})}
    \| \vec{A}^{-1} \vec{b} - \vec{x} \|_{\vec{A}}.
\end{align*}
Recall that for a positive definite matrix $\vec{M}$, the norm $\|\vec{x}\|_{\vec{M}}$ denotes $\|\vec{x}\|_{\vec{M}} = \sqrt{\vec{x}^\textup{\textsf{T}}\vec{M}\vec{x}}$. Since $\vec{Q}$ is a basis for $\mathcal{K}_k(\vec{A},\vec{b})$ we can equivalently write
\begin{align*}
    \mathsf{cg}_k = \vec{Q} \argmin_{\vec{c}\in\mathbb{R}^k}
    \| \vec{A}^{1/2}(\vec{A}^{-1} \vec{b} - \vec{Q}\vec{c}) \|_2
    = \vec{Q} \argmin_{\vec{c}\in\mathbb{R}^k}
    \| \vec{A}^{-1/2} \vec{b} - \vec{A}^{1/2}\vec{Q}\vec{c} \|_2.
\end{align*}
The solution to this least squares problem is
\begin{align*}
    \mathsf{cg}_k
    = \vec{Q} ( \vec{Q}^\textup{\textsf{T}} \vec{A}^{1/2} \vec{A}^{1/2} \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} \vec{A}^{1/2} \vec{A}^{-1/2} \vec{b}
    = \vec{Q} \vec{T}^{-1} \vec{e}_0.
\end{align*}
Here we have used that $\vec{Q}^\textup{\textsf{T}}\vec{A}\vec{Q} = \vec{T}$ and, since we are assuming $\|\vec{b}\|_2=1$, $\vec{Q}^\textup{\textsf{T}} \vec{b} = \| \vec{b}\|_2\vec{e}_0 = \vec{e}_0$.

A common alternative to CG for indefinite systems is MINRES which minimizes the $\vec{A}^2$-norm of the error (the 2-norm of the residual) over the Krylov subspace; i.e., the MINRES approximation $\mathsf{mr}_k$ is given by
\begin{align*}
    \mathsf{mr}_k
    := \argmin_{\vec{x}\in\mathcal{K}_k(\vec{A},\vec{b})}
    \| \vec{A}^{-1} \vec{b} - \vec{x} \|_{\vec{A}^2}
    =
    \argmin_{\vec{x}\in\mathcal{K}_k(\vec{A},\vec{b})}
    \| \vec{b} - \vec{A}\vec{x} \|_{2}
\end{align*}
Now note that
\begin{align*}
    \vec{A} \vec{Q}
    = [\vec{A} \widehat{\vec{Q}}]_{:,:k}
    =[\widehat{\vec{Q}}\widehat{\vec{T}}]_{:,:k}
    =[\widehat{\vec{Q}}]_{:,:k+1} [\widehat{\vec{T}}]_{:k+1,:k}.
\end{align*}
Therefore, since $\widehat{\vec{T}}$ is symmetric and $\widehat{\vec{Q}}$ has orthonormal columns,
\begin{align*}
    \mathsf{mr}_k = \vec{Q} (\vec{Q}^\textup{\textsf{T}} \vec{A} \vec{A} \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} \vec{A} \vec{b}
    = \vec{Q} ([\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1,:k})^{-1} \vec{T} \vec{e}_0.
\end{align*}

More generally,  suppose $z$ is an arbitrary complex number.
Then a special case of the quasi-minimum residual method (QMR) \cite{freund_92} can be used to compute the optimal $(\vec{A}^2+|z|^2 \vec{I})$-norm
approximation to $(\vec{A}-z \vec{I})^{-1}\vec{b}$; i.e., the QMR approximation $\textsf{qmr}_k(z)$ is given by
\begin{align}
    \textsf{qmr}_k(z)
    &:= \argmin_{\vec{x}\in\mathcal{K}_k(\vec{A},\vec{b})}
    \| (\vec{A}-z \vec{I})^{-1} \vec{b} - \vec{x} \|_{(\vec{A}^2+|z|^2 \vec{I})}
    \nonumber\\&=
    \argmin_{\vec{x}\in\mathcal{K}_k(\vec{A},\vec{b})}
    \| (\vec{A}-\overline z \vec{I})^{1/2}(\vec{A}-z \vec{I})^{-1/2}\vec{b} - (\vec{A}-\overline z \vec{I})^{1/2}(\vec{A}-z \vec{I})^{1/2}\vec{x} \|_{2}.
\end{align}
Here we have used that $\vec{A}^2+|z|^2 \vec{I} = (\vec{A}-\overline z \vec{I})(\vec{A}-z \vec{I})$.
Next, using the shift invariance of Krylov subspace, we see that
\begin{align*}
    \textsf{qmr}_k(z)
    &:=  \vec{Q} (\vec{Q}^\textup{\textsf{T}} (\vec{A}-\overline z \vec{I})(\vec{A}-z \vec{I}) \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} (\vec{A}- z \vec{I})\vec{b}.
    \\&=  \vec{Q} ( \vec{Q}( \vec{A}^2+|z|^2 \vec{I})\vec{Q}^\textup{\textsf{T}})^{-1} \vec{Q}^\textup{\textsf{T}} (\vec{A}- \overline{z} \vec{I})\vec{b}
    \\&=  \vec{Q} ([\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1,:k}+|z|^2 \vec{I})^{-1} (\vec{T} - \overline{z} \vec{I})\vec{e}_0.
\end{align*}

Standard implementations of CG, MINRES, and QMR require storing just a constant number of vectors of length $n$.
This can be done by taking advantage of the tridiagonal structure of $\widehat{\vec{T}}$, and such implementations are closely related to solving $\vec{T}^{-1} \vec{e}_0$ by an $\vec{L}\vec{D}\vec{L}^\textup{\textsf{T}}$ factorization \cite{paige_saunders_75,freund_92,liesen_strakos_13,simonova_tichy_21}.
The low-memory implementations for Lanczos-OR and Lanczos-FA that we derive in \cref{sec:low_mem} are based on this idea.

\begin{remark}
Even when $\vec{A}$ is indefinite, we can use the above expression for $\mathsf{cg}_k$.
However, the resulting algorithm is not guaranteed to be optimal, and if $\vec{T}$ has an eigenvalue near to or at zero then $\vec{T}^{-1}$ will be poorly conditioned or even undefined.
\end{remark}

\subsection{Lanczos-FA}

When $f(x)$ is something other than $1/x$, then CG is generalized by the Lanczos method for matrix function approximation (Lanczos-FA), which produces an approximation to $f(\vec{A})\vec{b}$ using the information in $\vec{Q}$ and $\vec{T}$.
\begin{definition}
\label{def:lanczos_fa}
The Lanczos-FA iterate is defined as
\begin{align*}
    \textup{\textsf{lan-FA}}_k(f)
    :=  \vec{Q} f(\vec{T}) \vec{e}_0.
 \end{align*}
\end{definition}
As discussed, Lanczos-FA is not in general guaranteed to return an optimal approximation to $f(\vec{A})\vec{b}$ from the Krylov subspace. While some weaker convergence bounds are known\footnote{The error of Lanczos-FA after $k$ iterations can be related to the error of the best degree-$k$ uniform polynomial approximation to $f$ on the range of $\vec{A}$'s eigenvalues \cite{musco_musco_sidford_18,saad_92}.} and the method often performs well in practice, it can exhibit erratic behavior in certain cases, for instance when $f(x) = 1/\sqrt{x}$ or $\operatorname{sign}(x)$ and $\vec{A}$ is indefinite; see \cref{sec:comparison} for examples. Moreover, for general $f$, computing the Lanczos-FA iterate naively requires storing $\vec{Q}$, which amounts to storing $k$ length $n$ vectors. Unlike CG, the memory requirements thus grow with the number of iterations $k$.
\subsubsection{Two pass Lanczos-FA}\label{sec:2pass}
This storage cost can be avoided by incurring additional computational cost. In particular, to avoid storing $\vec{Q}$, we can use an implementation called two pass Lanczos-FA \cite{boricci_00a,frommer_simoncini_08}.
On the first pass, the tridiagonal matrix $\vec{T}$ is computed using the short-recurrence version of Lanczos; i.e., without storing all of $\vec{Q}$.
Once $\vec{T}$ has been computed, $f(\vec{T}) \vec{e}_0$ can be computed using $O(k^2)$ storage.
Lanczos is then run again and the product $\vec{Q} f(\vec{T}) \vec{e}_0$ is computed as the columns of $\vec{Q}$ become available.%
\footnote{Note that two-pass Lanczos may be well suited for high performance distributed memory computing where inner products can be significantly more expensive than structured matrix-vector products \cite{exascale}.
The second pass can avoid computing inner products as the values of the inner products are already stored in $\vec{T}$.}

\section{Optimal rational function approximation}

We now describe an optimal iterate for approximating $r(\vec{A})\vec{b}$ when $r(x) = M(x)/N(x)$ for polynomials $M$ and $N$.
We will describe a low-memory implementation in \cref{sec:low_mem} which can also be used to  compute the Lanczos-FA approximations to $r(\vec{A})\vec{b}$ without doubling the number of matrix-vector products.

\begin{definition}
\label{def:lanczos_or}
Let $r(x) = M(x) / N(x)$ be a rational function with $N(x)$ having leading coefficient one.
Write
\begin{align*}
    N(x) = \left( \prod_{i=1}^{q} (x-z_i) \right) \left(\prod_{i=1}^{q'} (x-{z'}_i)(x-\overline{z'}_i) \right)
\end{align*}
where $z_i$ are distinct and $z_i\neq \overline{z}_j$ for all $i,j$ with $j \neq i$
 Set $R(x) = \xi \prod_{i=1}^{q} (x-\overline{z}_i)^{\alpha_i}$ where $\alpha_i = 0$ if $z_i\in \mathbb{R} \setminus \mathcal{I}(\vec{A})$  and $\alpha_i = 1$ otherwise and $\xi$ is chosen so that $R(\lambda_{\textup{min}}(\vec A))N(\lambda_{\textup{min}}(\vec{A})) \geq 0$.
Define $\tilde{M}(x) = R(x) M(x)$ and $\tilde{N}(x) = R(x) N(x)$.
Then the Lanczos-OR iterate is defined as
\begin{align*}
    \textup{\textsf{lan-OR}}_k(M,N)
    :=  \vec{Q}  ([\tilde{N}(\widehat{\vec{T}})]_{:k,:k})^{-1} [\tilde{M}(\widehat{\vec{T}})]_{:k,:k} \vec{e}_0.
 \end{align*}
\end{definition}

\begin{remark}
Clearly $\tilde{M}(x) / \tilde{N}(x) = M(x)/N(x) =  r(x)$.
Indeed, we can simply view $R(x)$ as the minimum degree polynomial so that $\tilde{N}(x) = N(x)R(x)$ is real and non-negative over $\mathcal{I}(\vec{A})$.
In other words, $R(x)$ simply contains a factor $(x-z_i)$ for reach root $z_i$ of $N(x)$ which is contained in $\mathcal{I}(\vec{A})$ and has odd multiplicity.
\end{remark}

\iffalse
Indeed, suppose $M(x) = 1$, $N(x) = x$, and $\vec{A}$ is symmetric positive definite.
Then $R(x) = 1$ and
\begin{align*}
    \textup{\textsf{lan-OR}}(M,N)
    = \vec{Q} ([\widehat{\vec{T}}]_{:k,:k})^{-1} \vec{e}_0
    = \vec{Q} \vec{T}^{-1} \vec{e}_0
    = \mathsf{cg}_k.
\end{align*}
Alternately, suppose $\vec{A}$ is symmetric indefinite.
Then $R(x) = x$ and
\begin{align*}
    \textup{\textsf{lan-OR}}(M,N)
    &= \vec{Q} ([\widehat{\vec{T}}^{2}]_{:k,:k})^{-1}\vec{T} \vec{e}_0
    = \vec{Q} ([\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1:k})^{-1} \vec{T} \vec{e}_0
    = \mathsf{mr}_k.
\end{align*}
Finally, suppose $\vec{A}$ is symmetric, $M(x) = 1$, and $N(x) = (x-z)$ for some purely imaginary $z$ vb .
Then $R(x) = (x-\overline{z})$ and
\begin{align*}
    \textup{\textsf{lan-OR}}(M,N)
    &= \vec{Q} ([\widehat{\vec{T}}^{2}]_{:k,:k}+z^2\vec{I})^{-1} (\vec{T} - \overline{z} \vec{I}) \vec{e}_0
    \\&=  \vec{Q} ([\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1,:k}+|z|^2 \vec{I})^{-1} (\vec{T} - \overline{z} \vec{I})\vec{e}_0
    = \mathsf{qmr}_k(z).
\end{align*}
\fi

It is not hard to see that CG, MINRES, and the special case of QMR described in \cref{sec:linear_solvers} are all obtained as special cases of Lanczos-OR.
Specifically, if $r(x) = 1/x$ and $\vec{A}$ is positive definite, then $0\in \mathbb{R}\setminus \mathcal{I}(\vec{A})$ so $R(x) = 1$ and we obtain the CG iterate optimal in the $\vec{A}$-norm.
Alternately, if $r(x) = 1/x$ and $\vec{A}$ is indefinite, then $0\in \mathcal{I}(\vec{A})$ so $R(x) = x$ and we obtain the MINRES iterate optimal in the $\vec{A}^2$-norm.
Finally, if $r(x) = 1/(x-z)$ and $z\in \mathbb{C}\setminus \mathcal{I}(\vec{A})$ then $R(x) = (x+\overline{z})$ and we obtain the QMR iterate optimal in the $(\vec{A}^2 + |z|^2 \vec{I})$-norm.

In fact, Lanzcos-OR is \emph{always} optimal in a certain norm depending on $N$ and its relation to $\mathcal{I}(\vec{A})$.
Specifically, we have the following optimality guarantee, which we prove at the end of the section.
\begin{theorem}
\label{thm:lanczos_opt}
Let $r(x) = M(x)/N(x)$, $\tilde{M},\tilde{N}$ as in \cref{def:lanczos_or} and suppose $N(\lambda) \neq 0$ for all $\lambda \in \Lambda(\vec{A})$.
Define $\vec{H} = \tilde{N}(\vec{A})$.
Then $\vec{H}$ is positive definite and $\textup{\textsf{lan-OR}}_k(M,N)$ is the $\vec{H}$-norm optimal approximation to $r(\vec{A})\vec{b}$ from $\mathcal{K}_k(\vec{A},\vec{b})$.
\end{theorem}

The optimality of Lanczos-OR implies an a priori scalar polynomial error bound in terms of \emph{the eigenvalues of $\vec{A}$}.
This bound generalizes the well known minimax bounds for CG, MINRES, and QMR \cite{greenbaum_97}
\begin{corollary}
Let $r(x) = M(x)/N(x)$, $\tilde{M},\tilde{N}$ as in \cref{def:lanczos_or}, suppose $N(\lambda) \neq 0$ for all $\lambda \in \Lambda(\vec{A})$, and let $\vec{H}$ be as in \cref{thm:lanczos_opt}.
Then, with $\Lambda(\vec{A})$ denoting the set of eigenvalues of $\vec{A}$,
\begin{align*}
    \frac{\| \textup{\textsf{lan-OR}}_k(M,N) \|_{\vec{H}}}{\| \textup{\textsf{lan-OR}}_0(M,N) \|_{\vec{H}}}
    \leq \min_{\deg(p) < k} \max_{\lambda \in \Lambda(\vec{A})} | r(\lambda) - p(\lambda) |.
\end{align*}
\end{corollary}

\begin{proof}
By \cref{thm:lanczos_opt}, $\textup{\textsf{lan-OR}}_k(M,N)$ is the $\vec{H}$-norm optimal approximation to $r(\vec{A})\vec{b}$ from $\mathcal{K}_k(\vec{A},\vec{b})$.
Thus,
\begin{align*}
\| \textup{\textsf{lan-OR}}_k(M,N) \|_{\vec{H}}
= \min_{\vec{x}\in \mathcal{K}_k(\vec{A},\vec{b})} \| r(\vec{A})\vec{b} - \vec{x} \|_{\vec{H}}
= \min_{\deg(p) < k} \| r(\vec{A})\vec{b} - p(\vec{A}) \vec{b} \|_{\vec{H}}.
\end{align*}
Next, using the fact that $\vec{A}$ and $\vec{H}^2$ commute, we note that,
\begin{align*}
    \| r(\vec{A})\vec{b} - p(\vec{A}) \vec{b} \|_{\vec{H}}
    \leq \| r(\vec{A}) - p(\vec{A}) \|_2 \| \vec{b} \|_{\vec{H}}.
\end{align*}
Finally, using the definition of the spectral norm,
\begin{align*}
    \| r(\vec{A}) - p(\vec{A}) \|_2
    = \| (r-p)(\vec{A}) \|_2
    = \max_{\lambda \in \Lambda(\vec{A})} | r(\lambda) - p(\lambda) |.
\end{align*}
The result follows from the fact that $ \textup{\textsf{lan-OR}}_0(M,N)  = \vec{b}$
\end{proof}

While the optimally of the Lanczos-OR iterate is nice, it is not yet apparent that the iterate can be computed efficiently, since it involves the terms  $[\tilde{N}(\widehat{\vec{T}})]_{:k,:k}$ and $[\tilde{M}(\widehat{\vec{T}})]_{:k,:k}$, and computing $\widehat{\vec{T}}$ requires running Lanczos to completion.
The following lemma, which we prove at the end of the section, implies that these terms can be computed without much more information than is in $\vec{T}$.
This means the Lanczos-OR iterate can be computed with only a few more matrix vector products than required to compute the Lanczos-FA iterate.

\begin{lemma}
\label{thm:poly_That}
Suppose $p$ is a polynomial with $q := \deg(p) > 0$.
Then $[p(\widehat{\vec{T}})]_{:k,:k}$ can be computed using the coefficients generated by $k+\lfloor (q-1)/2 \rfloor$ iterations of Lanczos.
Morever, if $k':=k+\lfloor q/2 \rfloor$, then
\begin{align*}
    [p(\widehat{\vec{T}})]_{:k,:k} = [p([\widehat{\vec{T}}]_{:k',:k'})]_{:k,:k}.
\end{align*}
\end{lemma}

\begin{proof
It suffices to consider the case $p(x) = x^q$.
Let $\widehat{\vec{I}}_\ell$ be the $K\times \ell$ matrix with the $\ell\times \ell$ identity in the the top $\ell$ rows and zeros in the bottom $K-\ell$ rows.
Then, we have that
\begin{align*}
    \widehat{\vec{T}} \widehat{\vec{I}}_\ell
    = [\widehat{\vec{T}}]_{:,:\ell}
    = \widehat{\vec{I}}_{\ell+1} [\widehat{\vec{T}}]_{:\ell+1,:\ell}.
\end{align*}
Repeatedly applying this relation,
\begin{align*}
    \widehat{\vec{T}}^j \widehat{\vec{I}}_k
    =
    \widehat{\vec{I}}_{k+j}
    [\widehat{\vec{T}}]_{:k+j,k+j-1}
    \cdots
    [\widehat{\vec{T}}]_{:k+2,k+1}
    [\widehat{\vec{T}}]_{:k+1,k}
    = \widehat{\vec{I}}_{k+j} \vec{B}(k+j,k)
\end{align*}
where we have defined
\begin{align*}
    \vec{B}(k+j,k) := [\widehat{\vec{T}}]_{:k+j,k+j-1}
    \cdots
    [\widehat{\vec{T}}]_{:k+2,k+1}
    [\widehat{\vec{T}}]_{:k+1,k}.
\end{align*}
Therefore, since $(\widehat{\vec{I}}_{k+j})^\textup{\textsf{T}} \widehat{\vec{I}}_{k+j}$ is the $(k+j)\times (k+j)$ identity,
\begin{align*}
    [\widehat{\vec{T}}^{2j}]_{:k,:k}
    &= \vec{B}(k+j,k)^\textup{\textsf{T}} \vec{B}(k+j,k)
\intertext{and, since $(\widehat{\vec{I}}_{k+j-1})^\textup{\textsf{T}} \widehat{\vec{T}} \widehat{\vec{I}}_{k+j-1} = [\vec{T}]_{:k+j-1,:k+j-1}$,}
    [\widehat{\vec{T}}^{2j-1}]_{:k,:k}
    &= \vec{B}(k+j-1,k)^\textup{\textsf{T}}  [\widehat{\vec{T}}]_{:k+j-1,:k+j-1} \vec{B}(k+j-1,k).
\end{align*}
In both cases, we see that these expressions depend only on $[\widehat{\vec{T}}]_{:k+j,:k+j-1}$ which can be obtained using $k+j-1$ matrix-vector products.
The result follows by noting that $\lfloor (2j-1)/2 \rfloor = \lfloor(2j-1-1)/2 \rfloor = j-1$.

To complete the lemma, note that for $\ell+1\leq k'$,
\begin{align*}
    [\widehat{\vec{T}}]_{:k',:k'} \widetilde{\vec{I}}_\ell
    = [\widehat{\vec{T}}]_{:k',:\ell}
    = \widetilde{\vec{I}}_{\ell+1} [\widehat{\vec{T}}]_{:\ell+1,:\ell}
\end{align*}
where $\widetilde{\vec{I}}_\ell$ is defined like $\widehat{\vec{I}}_\ell$ but is of size $k' \times \ell$.
The same argument as above then implies that
\begin{align*}
    ([\widehat{\vec{T}}]_{:k',:k'})^j \widetilde{\vec{I}}_k
    =\widetilde{\vec{I}}_{k+j}\vec{B}(k+j,k)
\end{align*}
provided that \( k+j \leq k' \).
We therefore have that
\begin{align*}
    [\widehat{\vec{T}}^{q}]_{:k,:k} = [([\widehat{\vec{T}}]_{:k',:k'})^{q}]_{:k,:k}.
\end{align*}
\end{proof}

\iffalse
\begin{lemma}
Suppose $p$ is a polynomial with $q := \deg(p) > 0$.
Then $[p(\widehat{\vec{T}})]_{:k,:k}$ can be computed using the coefficients generated by $k+\lfloor (q-1)/2 \rfloor$ iterations of Lanczos.
\end{lemma}

\begin{proof
It suffices to consider the case $p(x) = x^q$.
Let $\widehat{\vec{I}}_\ell$ be the $K\times \ell$ matrix with the $\ell\times \ell$ identity in the the top $\ell$ rows and zeros in the bottom $K-\ell$ rows.
Then, we have that
\begin{align*}
    \widehat{\vec{T}} \widehat{\vec{I}}_\ell
    = [\widehat{\vec{T}}]_{:,:\ell}
    = \widehat{\vec{I}}_{\ell+1} [\widehat{\vec{T}}]_{:\ell+1,:\ell}.
\end{align*}
Repeatedly applying this relation,
\begin{align*}
    \widehat{\vec{T}}^j \widehat{\vec{I}}_k
    =
    \widehat{\vec{I}}_{k+j}
    [\widehat{\vec{T}}]_{:k+j,k+j-1}
    \cdots
    [\widehat{\vec{T}}]_{:k+2,k+1}
    [\widehat{\vec{T}}]_{:k+1,k}
    = \widehat{\vec{I}}_{k+j} \vec{B}(k+j,k)
\end{align*}
where we have defined
\begin{align*}
    \vec{B}(k+j,k) := [\widehat{\vec{T}}]_{:k+j,k+j-1}
    \cdots
    [\widehat{\vec{T}}]_{:k+2,k+1}
    [\widehat{\vec{T}}]_{:k+1,k}.
\end{align*}
Therefore, since $(\widehat{\vec{I}}_{k+j})^\textup{\textsf{T}} \widehat{\vec{I}}_{k+j}$ is the $(k+j)\times (k+j)$ identity,
\begin{align*}
    [\widehat{\vec{T}}^{2j}]_{:k,:k}
    &= \vec{B}(k+j,k)^\textup{\textsf{T}} \vec{B}(k+j,k)
\intertext{and, since $(\widehat{\vec{I}}_{k+j-1})^\textup{\textsf{T}} \widehat{\vec{T}} \widehat{\vec{I}}_{k+j-1} = [\vec{T}]_{:k+j-1}$,}
    [\widehat{\vec{T}}^{2j-1}]_{:k,:k}
    &= \vec{B}(k+j,k)^\textup{\textsf{T}}  [\widehat{\vec{T}}]_{:k+j-1,:k+j-1} \vec{B}(k+j,k).
\end{align*}
In both cases, we see that these expressions depend only on $[\widehat{\vec{T}}]_{:k+j,:k+j-1}$ which can be obtained using $k+j-1$ matrix-vector products.
The result follows by noting that $\lfloor (2j-1)/2 \rfloor = \lfloor(2j-1-1)/2 \rfloor = j-1$.
\end{proof}
\fi

\subsection{Proof of optimality}

The proof of \cref{thm:lanczos_opt} is straightforward and essentially follows our derivations of the methods in \cref{sec:linear_solvers}.
For clarity, we first state and prove a useful lemma

\iffalse
\begin{lemma}
    \label{thm:lanczos_exact_poly}
    Let $p$ be a polynomial with $\deg(p) < k$.
    Then
    \begin{align*}
    p(\vec{A})\vec{b} = \vec{Q} p( \vec{T} ) \vec{e}_0.
    \end{align*}
\end{lemma}

\begin{proof
Observe that
\begin{align*}
    p(\vec{A}) \vec{b}
    = p(\vec{A}) \widehat{\vec{Q}} \vec{e}_0
    = \widehat{\vec{Q}} p(\widehat{\vec{T}})
    \vec{e}_0.
\end{align*}
Since $\widehat{\vec{T}}$ is tridiagonal, for any $j > q$, $\vec{e}_j^\textup{\textsf{T}} \widehat{\vec{T}}^q \vec{e}_0 = 0$.
Thus, we see that entries $0$ through $k-1$ of $p( \hat{\vec{T}} ) \vec{e}_0$ are the same as
those of $p( \vec{T} ) \vec{e}_0$ and also that only columns $0$ through $k-1$ of $\widehat{\vec{Q}}$ contribute.
\end{proof}
\fi

\begin{lemma}
\label{thm:H_orth_projector}
    Suppose $\vec{H}$ is positive definite and $\vec{f}\in\mathbb{R}^n$ is an arbitrary vector.
    Then,
    \begin{align*}
        \vec{Q}( \vec{Q}^\textup{\textsf{T}} \vec{H} \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} \vec{H} \vec{f} = \argmin_{\vec{y} \in \mathcal{K}_k(\vec{A},\vec{b})} \| \vec{f} - \vec{y} \|_{\vec{H}}
    \end{align*}
\end{lemma}

\begin{proof
Since $\vec{y}\in\mathcal{K}_k(\vec{A},\vec{b})$, we have $\vec{y} = \vec{Q} \vec{c}$ for some vector $\vec{c}$.
Thus, we can consider the problem
\begin{align*}
    \argmin_{\vec{c}\in\mathbb{R}^k} \| \vec{f} - \vec{Q} \vec{c} \|_{\vec{H}}
    =
    \argmin_{\vec{c}\in\mathbb{R}^k} \| \vec{H}^{1/2} \vec{f} - \vec{H}^{1/2} \vec{Q} \vec{c} \|_2.
\end{align*}
But this is just a standard least squares problem which has solution
\begin{align*}
    \vec{c} = ( (\vec{H}^{1/2} \vec{Q})^\textup{\textsf{T}} (\vec{H}^{1/2} \vec{Q}) )^{-1} (\vec{H}^{1/2} \vec{Q})^\textup{\textsf{T}} (\vec{H}^{1/2} \vec{f}).
\end{align*}
The result follows by simplifying the expression for $\vec{Q}\vec{c}$.
\end{proof}
With these lemmas in place, we are prepared to prove the optimality of Lanczos-OR.

\begin{proof}[Proof of \cref{thm:lanczos_opt}]

We first prove that $\vec{H}$ is positive definite.
By \cref{def:lanczos_or}, we observe that $\tilde N(x)$ is defined as
\begin{align*}
    \tilde{N}(x) = \xi \left( \prod_{i=1}^{s} (x-z_i) \right) \left(\prod_{i=1}^{s'} (x-{z'}_i)(x-\overline{z'}_i) \right)
\end{align*}
for $z_i \in \mathbb{R} \setminus \mathcal{I}(\vec{A})$.
Clearly $\prod_{i=1}^{q'} (x-{z'}_i)(x-\overline{z'}_i) \geq 0$ for all $x \in \mathbb{R}$, and for $z_i \in \mathbb{R} \setminus \mathcal{I}(\vec{A})$, $(x-z_i)$ does not change signs over $\mathcal{I}(\vec{A})$.
The choice of $\xi$ ensures that $\tilde{N}(\lambda_{\textup{min}}(\vec{A}))>0$, and by assumption $N(\lambda) \neq 0$ for all $\lambda \in \Lambda(\vec{A})$ so that that $\tilde{N}(\lambda) \neq 0$ for all $\lambda\in\Lambda(\vec{A})$.
It follows that $\tilde{N}(\lambda) > 0$ for all $\lambda\in\Lambda(\vec{A})$; i.e. that $\vec{H} = \tilde{N}(\vec{A})$ is positive definite.

We now prove the Lanczos-OR iterate is optimal in $\vec{H}$-norm.
By \cref{thm:H_orth_projector} we have that the $\vec{H}$-norm optimal approximation of $r(\vec{A})\vec{b}$ is given by
\begin{align*}
    \vec{Q}(\vec{Q}^\textup{\textsf{T}} \vec{H} \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} \vec{H} r(\vec{A}) \vec{b}
\end{align*}
By definition, $\vec{H} = \tilde{N}(\vec{A})$ so $\vec{H}r(\vec{A}) = \tilde{M}(\vec{A})$.
\begin{align*}
    \vec{Q}^\textup{\textsf{T}} \vec{H} r(\vec{A}) \vec{b}
    = \vec{Q}^\textup{\textsf{T}} M(\vec{A}) \vec{b}
    = \vec{Q}^\textup{\textsf{T}} \widehat{\vec{Q}}\tilde{M}(\widehat{\vec{T}}) \widehat{\vec{Q}}^\textup{\textsf{T}} \vec{b}
    = [\tilde{M}(\widehat{\vec{T}})]_{:k,:k} \vec{e}_0.
\end{align*}
Then, since $\vec Q$ consists of the first $k$ columns $\widehat{\vec{Q}}$,
\begin{align*}
    \vec{Q}^\textup{\textsf{T}} \vec{A}^q \vec{Q} = [\widehat{\vec{Q}}^\textup{\textsf{T}} \vec{A}^q \widehat{\vec{Q}}]_{:k,:k}

    = [ \widehat{\vec{T}}^{q} ]_{:k,:k}
\end{align*}
so since $\vec{H}$ is a linear combination of powers of $\vec{A}$ we obtain
\begin{align*}
    \vec{Q}^\textup{\textsf{T}} \vec{H} \vec{Q}
    = \vec{Q}^\textup{\textsf{T}} \tilde{N}(\vec{A}) \vec{Q}
    =  [\tilde{N}(\widehat{\vec{T}}) ]_{:k,:k}.
\end{align*}
The result follows by combining the above expressions.
\end{proof}

\section{Algorithms for other matrix functions}

We now show that Lanczos-OR can be used to derive algorithms for other matrix functions, in particular, the matrix sign function.

\subsection{The matrix sign function}
\label{sec:sign_function}

A range of past works on Krylov subspace methods for matrix function approximation use integral representations of $f$ to provide error bounds/estimates and even to derive more advanced schemes such as restarted Lanczos-FA \cite{ilic_turner_simpson_09,frommer_guttel_schweitzer_14,frommer_guttel_schweitzer_14a,frommer_schweitzer_15,chen_greenbaum_musco_musco_21}.
In this section we show that Lanczos-OR yields a new algorithm for the matrix sign function.
The sign function has long been used in lattice quantum chromodynamics for studying matter interactions \cite{eshof_frommer_lippert_schilling_van_der_vorst_02}.
The sign function also induces a distribution function such that $\vec{b}^\textup{\textsf{T}} f(\vec{A})\vec{b}$ is equivalent to an integral of $f$ against this distribution function \cite{golub_meurant_09}.
As such, approximations to the sign function yield approximations to the spectrum of a matrix \cite{lin_saad_yang_16,braverman_krishnan_musco_21,chen_trogdon_ubaru_21}.

\iffalse
We begin by noting that, for any $x > 0$,
\begin{align*}
    \frac{1}{\sqrt{x}} = \frac{2}{\pi} \int_{0}^{\infty} \frac{1}{x + z^2} \d{z}.
\end{align*}
Thus, if $f(x) = 1/|x| = 1/\sqrt{x^2}$,
\begin{align*}
    f(\vec{A})\vec{b} = \frac{2}{\pi} \int_{0}^{\infty} (\vec{A}^2 + z^2\vec{I})^{-1} \vec{b} \: \d{z}.
\end{align*}
We can approximate $(\vec{A}^2 + z^2\vec{I})^{-1} \vec{b}$ by the Lanczos-OR iterate $\vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z^2\vec{I})^{-1} \vec{e}_0$, which is the optimal approximation over Krylov subspace in the $(\vec{A}^2 + z^2\vec{I})$-norm.
Noting that $[\widehat{\vec{T}}^2]_{:k,:k} = [\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1:k}$ gives the approximation
\begin{align*}
    \frac{2}{\pi} \int_{0}^{\infty}  \vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z^2\vec{I})^{-1} \vec{e}_0 \: \d{z}
    = \vec{Q} \left( [\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1:k} \right)^{-1/2} \vec{e}_0.
\end{align*}
Similarly,
\fi
We begin by noting that, for any $a>0$,
\begin{align*}
    \frac{1}{\sqrt{a}}
    = \frac{2}{\pi} \int_{0}^{\infty} \frac{1}{a + z^2} \d{z}.
\end{align*}
Thus, if $f(x) = \operatorname{sign}(x)=x/|x| = x/\sqrt{x^2}$, we have
\begin{align*}
    f(\vec{A})\vec{b} = \frac{2}{\pi} \int_{0}^{\infty} \vec{A}(\vec{A}^2 + z^2\vec{I})^{-1} \vec{b} \: \d{z}.
\end{align*}
The Lanczos-OR approximation to $\vec{A}(\vec{A}^2 + z^2\vec{I})^{-1} \vec{b}$ is $\vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z^2\vec{I})^{-1} \vec{T} \vec{e}_0$, which is optimal over Krylov subspace in the $(\vec{A}^2 + z^2\vec{I})$-norm.
This yields the approximation
\begin{align*}
    \frac{2}{\pi} \int_{0}^{\infty}  \vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z^2\vec{I})^{-1} \vec{T} \vec{e}_0 \: \d{z}
    &= \vec{Q} \left([\widehat{\vec{T}}^2]_{:k,:k} \right)^{-1/2} \vec{T}  \vec{e}_0.
\end{align*}
Thus, we can define the induced iterate as
\begin{align*}
\textup{\textsf{sign-OR}}_k := \vec{Q} \left([\widehat{\vec{T}}^2]_{:k,:k} \right)^{-1/2} \vec{T}  \vec{e}_0
=
\vec{Q} \left( [\widehat{\vec{T}}]_{:k,:k+1}[\widehat{\vec{T}}]_{:k+1,:k} \right)^{-1/2} \vec{T}  \vec{e}_0.
\end{align*}

\subsubsection{Relation to Lanczos-FA}

The Lanczos-OR and Lanczos-FA iterates for a given rational matrix function are clearly related.
In particular, $\tilde{N}(\vec{T})$ and $[\tilde{N}(\widehat{\vec{T}})]_{:k,:k}$ differ only in the bottom rightmost $(q-1)\times (q-1)$ principle submatrix, where $q = \deg(\tilde{N})$.
Using this fact, it can be shown that the Lanczos-OR and Lanczos-FA iterates ``tend to coalesce as convergence takes place'' \cite[Proposition 5.1]{lopez_simoncini_06}.
We now show that a similar phenomenon occurs with the induced Lanczos-OR approximation to the sign function and the Lanczos-FA approximation.

\begin{theorem}
\label{thm:sign_approx_compare}
The Lanczos-FA and induced Lanczos-OR approximations to the matrix sign function satisfy
\begin{align*}
    \| \textup{\textsf{lan-FA}}_k(\operatorname{sign}) - \textup{\textsf{sign-OR}}_k \|_2
    \leq
    (\beta_{k-1})^2

    \frac{\sigma_{\textup{max}}(\vec{T})}{2 \sigma_{\textup{min}}(\vec{T})^3},
\end{align*}
where $\sigma_{\textup{max}}(\vec{T})$ and $\sigma_{\textup{min}}(\vec{T})$ are the largest and smallest singular values of $\vec{T}$ respectively.
\end{theorem}

\begin{proof
Let ${N}(x) = x^2 + z^2$ and ${M}(x) = x$.
Note that ${N}(\vec{T}) = [{N}(\widehat{\vec{T}})]_{:k,:k} - \beta_{k-1}^2 \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}}$ so,
\begin{align*}
    {M}(\vec{T}) \vec{e}_0
    = ([{N}(\widehat{\vec{T}})]_{:k,:k} - \beta_{k-1}^2 \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}}){N}(\vec{T})^{-1} {M}(\vec{T}) \vec{e}_0.
\end{align*}
Thus, rearranging terms and multiplying by $([{N}(\vec{T})]_{:k,:k})^{-1}$ we find that
\begin{align*}
    &{N}(\vec{T})^{-1} {M}(\vec{T}) \vec{e}_0 - ([{N}(\widehat{\vec{T}})]_{:k,:k})^{-1}
    {N}(\vec{T})^{-1}{M}(\vec{T}) \vec{e}_0
    \\&\hspace{6em}= \beta_{k-1}^2 ([{N}(\widehat{\vec{T}})]_{:k,:k})^{-1} \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}} {N}(\vec{T})^{-1}{M}(\vec{T}) \vec{e}_0
\end{align*}
After left multiplying with $\vec{Q}$, this implies that
\begin{align*}
    \textup{\textsf{lan-FA}}_k(M,N) - \textup{\textsf{lan-OR}}_k(M,N)
    = \beta_{k-1}^2 \vec{Q} ([{N}(\widehat{\vec{T}})]_{:k,:k})^{-1} \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}} {M}(\vec{T}) \vec{e}_0.
\end{align*}
Now, suppose that $f(x) = \operatorname{sign}(x)$ and set $ \textsf{error}_k := \textup{\textsf{lan-FA}}_k(\operatorname{sign}) - \textup{\textsf{sign-OR}}_k$ as above.
Then, since the Lanczos-FA approximation can also be induced by an integral over $z\in[0,\infty)$, we have that,
\begin{align*}
    \textsf{error}_k
    = \beta_{k-1}^2
    \frac{2 }{\pi} \int_0^\infty
    \vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z^2 \vec{I})^{-1} \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}} (\vec{T}^2 + z^2 \vec{I})^{-1} \vec{T} \vec{e}_0
    \d{z}.
\end{align*}
Note that $[\widehat{\vec{T}}^2]_{:k,:k} - \vec{T}^2 = \beta_{k-1}^2 \vec{e}_{k-1} \vec{e}_{k-1}^\textup{\textsf{T}}$ is positive semidefinite.
Therefore, using that $\sigma_{\textup{min}}([\widehat{\vec{T}}^2]_{:k,:k}) \geq \sigma_{\textup{min}}(\vec{T}^2) = \sigma_{\textup{min}}(\vec{T})^2$,
\begin{align*}
    \| \textsf{error}_k \|_2
    &
    = \beta_{k-1}^2
     \left\| \vec{Q} \left( \frac{2 }{\pi} \int_0^\infty
    ([\widehat{\vec{T}}^2]_{:k,:k} + z^2 \vec{I})^{-1} \vec{e}_{k-1}\vec{e}_{k-1}^\textup{\textsf{T}} (\vec{T}^2 + z^2 \vec{I})^{-1} \d{z} \right)  \vec{T} \vec{e}_0
     \right\|_2
    \\&
    \leq \beta_{k-1}^2
     \left( \frac{2 }{\pi} \int_0^\infty
    \| ([\widehat{\vec{T}}^2]_{:k,:k} + z^2 \vec{I})^{-1} \|_2 \| (\vec{T}^2 + z^2 \vec{I})^{-1} \|_2 \d{z} \right)  \| \vec{T} \|_2
    \\& \leq
    \beta_{k-1}^2
     \left( \frac{2 }{\pi} \int_0^\infty
    | (\sigma_{\textup{min}}(\vec{T})^2 + z^2 )^{-1} |~ | (\sigma_{\textup{min}}(\vec{T})^2 + z^2 )^{-1} | \d{z} \right)  \sigma_{\textup{max}}(\vec{T})
    \\&= \beta_{k-1}^2
    \frac{\sigma_{\textup{max}}(\vec{T})}{2 \sigma_{\textup{min}}(\vec{T})^3}.
\end{align*}
\end{proof}

Since $|\beta_{k-1}|$ tends to decrease as the Lanczos method converges, this seemingly implies that the induced Lanczos-OR iterate and the Lanczos-FA iterate tend to converge in this limit.
However, recall that $\vec{T} = [\widehat{\vec{T}}]_{:k,:k}$ changes at each iteration $k$.
Thus, there is the difficulty that $\vec{T}$ may have an eigenvalue near zero, in which case the preceding bound could be useless.
However, it is known that $\vec{T}$ cannot have eigenvalues near zero in two successive iterations, assuming that the eigenvalues of $\vec{A}$ are not too close to zero.
Specifically,  \cite[Equation 3.10]{greenbaum_druskin_knizhnerman_99} asserts that
\begin{align*}
\max\{ \sigma_{\textup{min}}([\vec{T}]_{:k-1}),\sigma_{\textup{min}}([\vec{T}]_{:k,:k})\} > \frac{\sigma_{\textup{min}}(\vec{A})^2}{(2+\sqrt{3})\|\vec{A}\|_2}.
\end{align*}
Since $\beta_{k-1}$ has little to do with the minimum magnitude eigenvalue of $\vec{T}$ (recall that the Lanczos recurrence is shift invariant), we expect that the ``overall'' convergence of the induced Lanczos-OR iterate and the Lanczos-FA iterate will be similar as Lanczos converges.

This implies that existing spectrum dependent bounds for Lanczos-FA for the sign function \cite{chen_greenbaum_musco_musco_21} can be carried over to the Lanczos-OR based approximation.
More interestingly, it means that understanding the induced approximation provides a way of understanding Lanczos-FA.
Since Lanczos-FA often exhibits oscillatory behavior, bounds for the induced Lanczos-OR based approximation may be easier to obtain.

\subsection{Rational function approximation}
\label{sec:rational_func}

We can use a similar approach to derive (non-optimal) approximations to rational matrix functions $r(\vec{A})\vec{b}$.
In many settings, particularly when the rational function $r$ is used as a proxy for a function $f$, this approach is more natural than computing the Lanczos-OR approximation to the rational function directly.
The convergence of such methods is closely related to the quality of the (scalar) rational function approximation as well as the quality of the approximation to the rational matrix function.
Specifically, for any output $\mathsf{alg}(r)$ meant to approximate $r(\vec{A})\vec{b}$, we have the following bound:
\begin{align}
\| f(\vec{A}) \vec{b} - \mathsf{alg}(r) \|
&\leq
\nonumber
\| f(\vec{A}) \vec{b} - r(\vec{A}) \vec{b} \| +  \| r(\vec{A}) \vec{b} - \mathsf{alg}(r) \|
\\&\leq \nonumber
\| f(\vec{A}) - r(\vec{A}) \| \|\vec{b} \|_2 +  \| r(\vec{A}) \vec{b} - \mathsf{alg}(r) \|
\\&\leq \nonumber
\|\vec{b}\|_2  \max_{\lambda \in \Lambda(\vec{A})} | f(\lambda) - r(\lambda) | +  \| r(\vec{A}) \vec{b} - \mathsf{alg}(r) \|
\\& \leq \underbrace{\|\vec{b}\|_2  \max_{\lambda \in \mathcal{I}(\vec{A})} | f(\lambda) - r(\lambda) |}_{\text{approximation error}}
    + \underbrace{\| r(\vec{A}) \vec{b} - \mathsf{alg}(r) \|}_{\text{application error}}.
    \label{eqn:triangle_ineq}
\end{align}
In many cases, very good or even optimal scalar rational function approximations to a given function on a single interval are known or can be easily computed.
Thus, the approximation error term can typically be made small with a rational function of relatively low degree.

Of course, this bound is only meaningful if the approximation error term is small relative to the application error.
Indeed, we also have
\begin{align}
\| f(\vec{A}) \vec{b} - \mathsf{alg}(r) \|
&\geq
\big| \| f(\vec{A}) \vec{b} - r(\vec{A}) \vec{b} \| -  \| r(\vec{A}) \vec{b} - \mathsf{alg}(r) \| \big|.
\label{eqn:reverse_triangle_ineq}
\end{align}
This shows that the size of $\| f(\vec{A}) \vec{b} - \mathsf{alg}(r) \|$ is rouhgly the size of $\| f(\vec{A}) \vec{b} - r(\vec{A}) \vec{b} \|$ when $\mathsf{alg}(r)$ is a good approximation to $r(\vec{A}) \vec{b}$.

As we noted, rational function approximations commonly are obtained by discretizing an integral representation using a numerical quadrature approximation.
For instance, the matrix sign function may be approximated as
\begin{align*}
    r_q(\vec{A})\vec{b} = \sum_{i=1}^{q} \omega_i \vec{A} (\vec{A}^2 + z_i^2 \vec{I})^{-1} \vec{b}
\end{align*}
where $z_i$ and $\omega_i$ are appropriately chosen quadrature nodes and weights \cite{hale_higham_trefethen_08}.

We can of course write $r_q(x) = M_q(x) / N_q(x)$, so it's tempting to set $R_q(x) = 1$ and $\vec{H}_q = N_q(\vec{A})$ and then use Lanczos-OR to compute the $\vec{H}_q$-norm optimal approximation.
However, while $r_q(x)$ is convergent to $f(x)$ as $q\to\infty$, $N_q(x) := \prod_{i=1}^q (x^2
+ z_i^2 )$ is not convergent to any fixed function.
In fact $N_q(x)$ will increase in degree and $\vec{H}_q$ will be increasingly poorly conditioned.
This presents a numerical difficulty in computing the Lanczos-OR iterate in this limit.
More importantly, it is not clear that it is meaningful to approximate a function in this way.
Indeed, it seems reasonable to expect that, for fixed $k$, as $q\to\infty$,  our approximation should be convergent to something.
However, we cannot guarantee $\textup{\textsf{lan-OR}}_k(M_q,N_q)$ is convergent in this limit.

On the other hand, akin to our approach for approximating the integral described in the previous subsection, we can compute the \emph{term-wise optimal} approximations to each term in the sum representation of $r_q$ and output
\begin{align*}
    \sum_{i=1}^{q} \omega_i \vec{Q} ([\widehat{\vec{T}}^2]_{:k,:k} + z_i^2\vec{I})^{-1} \vec{T} \vec{e}_0.
\end{align*}
In this case, as $q\to\infty$, the approximation is convergent to the integral output.

Interestingly, this is exactly what would be obtained by using QMR to approximate each of the corresponding linear systems in the partial fractions decomposition.
\begin{lemma}
Suppose $z\in\mathbb{R}$ and let $M(x)= 1$ and $N(x) = x^2+z^2$.
Then
\begin{align*}
    \textup{\textsf{lan-OR}}_k(M,N)
    = \frac{1}{2z} \left( \mathsf{qmr}_k(-\bm{i} z) - \mathsf{qmr}_k(\bm{i} z) \right).
\end{align*}
\end{lemma}

\begin{proof}
This follows easily from the fact that both methods are optimal with respect to the same norm.
    Note that
    \begin{align*}
        \frac{1}{x^2+z^2}
        = \frac{1}{x-\bm{i} z} \frac{1}{x+\bm{i} z}
        = \frac{1}{2z} \left( \frac{1}{x-\bm{i} z} - \frac{1}{x+\bm{i} z} \right)
    \end{align*}
Let $\vec{P} = \vec{Q} (\vec{Q}^\textup{\textsf{T}} (\vec{A}^2 + z^2 \vec{I}) \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} (\vec{A}^2 + z^2 \vec{I})^{1/2}$.
Then
\begin{align*}
    \vec{P}(\vec{A}^2+z^2\vec{I})^{-1}\vec{A} \vec{b}
    &= \frac{1}{2z} \vec{P}\left( (\vec{A} - \bm{i} z \vec{I})^{-1} - (\vec{A} + \bm{i} z \vec{I})^{-1} \right) \vec{A} \vec{b}
    \\&= \frac{1}{2z} \left( \vec{P}(\vec{A} - \bm{i} z \vec{I})^{-1} \vec{A} \vec{b} - \frac{1}{2z}\vec{P} (\vec{A} + \bm{i} z \vec{I})^{-1} \vec{A} \vec{b} \right)
\end{align*}
But the left hand side is the Lanczos-OR iterate whereas the right hand side contains the QMR iterates.
\end{proof}

Whether it is better to use Lanczos-OR or QMR is somewhat unclear.
The Lanczos-OR based approach avoids the need for complex arithmetic, which simplifies implementation slightly.
However, QMR has been studied longer, and therefore is likely to have more practical low-memory implementations.

\iffalse

\subsection{Error analysis}

\Tyler{
What can we say about optimality? At each point $z$, the approximation is optimal in $\vec{A}^2$ norm up to a factor $\kappa(\vec{A}^2(\vec{A}^2+z^2\vec{I})^{-1})$.
This approaches $\vec{A}^2$ as $z\to\infty$.

However, as $z$ increases, the solution to the systems decreases in size and the systems become better conditioned and easier to solve.
Thus, in some sense, it's the behavior for $z$ near the origin which is the most important.

In fact,
\begin{align*}
    \frac{2}{\pi} \int_0^a \frac{x}{x^2+z^2}\d{z}
    = \frac{2}{\pi}\arctan \left( \frac{a}{x} \right).
\end{align*}
So maybe can discard upper part of integral or break integral into pieces.

the $\vec{A}^2$ optimal approximation is given by
\begin{align*}
     \frac{2}{\pi} \int_{0}^{\infty} \vec{Q} (\vec{Q}^\textup{\textsf{T}}\vec{A}^2 \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}} \vec{A}^2 \vec{A}(\vec{A}^2 + z^2\vec{I})^{-1} \vec{b} \: \d{z}.
\end{align*}
whereas our approximation is given by
\begin{align*}
     &\frac{2}{\pi} \int_{0}^{\infty} \vec{Q} (\vec{Q}^\textup{\textsf{T}}(\vec{A}^2+z^2\vec{I}) \vec{Q})^{-1} \vec{Q}^\textup{\textsf{T}}  (\vec{A}^2+z^2\vec{I}) \vec{A}(\vec{A}^2 + z^2\vec{I})^{-1} \vec{b} \: \d{z}
     \\=&
     \frac{2}{\pi} \int_{0}^{\infty} \vec{Q} (\vec{Q}^\textup{\textsf{T}}\vec{A}^2\vec{Q}+z^2\vec{I})^{-1} \vec{Q}^\textup{\textsf{T}}   \vec{A}\vec{b} \: \d{z}
\end{align*}

Maybe something like: when $z$ is small these are close, and when $z$ is big the contribution to the integral is small.

If think a whole paper just on the sign function and a fine grained comparison of different algorithms might be interesting.
I.e., trying to compare Lanzos-FA, this algorithm, the harmonic ritz values, etc.

I think there should be some MINRES-CG type comparions based on numerical experiments.
}

\fi

\section{Implementing Lanczos-OR using low memory}
\label{sec:low_mem}

We now describe a low-memory implementation of Lanczos-OR which is similar in spirit to CG, MINRES, and QMR.

For convenience, we will denote $\vec{M} := [\tilde{M}(\vec{T})]_{:k,:k}$ and $\vec{N} := [\tilde{N}(\widehat{\vec{T}})]_{:k,:k}$ so that the Lanczos-OR output is given by $\vec{Q} \vec{N}^{-1} \vec{M} \vec{e}_0$.
Then, at a high level, our approach is to:
\begin{itemize}
    \item Take one iteration of Lanczos to generate one more column of $\widehat{\vec{Q}}$ and $\widehat{\vec{T}}$
    \item Compute one more column of $\vec{M}$ and $\vec{N}$
    \item Compute one more factor of $\vec{L}^{-1} = \vec{L}_{k-1}\cdots \vec{L}_1 \vec{L}_0$ and one more entry of $\vec{D}$ where $\vec{L}$ and $\vec{D}$ are defined by the LDL factorization $\vec{N} =\vec{L}\vec{D}\vec{L}^\textup{\textsf{T}}$
    \item Compute one more term of the sum:
    \begin{align*}
        \vec{Q}\vec{N}^{-1}\vec{M}\vec{e}_0
        = \vec{Q} \vec{L}^{-1} \vec{D}^{-1} \vec{L}^{-\textup{\textsf{T}}} \vec{M} \vec{e}_0
        = \sum_{i=0}^{k-1} \frac{[\vec{L}^{-\textup{\textsf{T}}}\vec{M}\vec{e}_0]_i}{[\vec{D}]_{i,i}} [\vec{Q}\vec{L}^{-1}]_{:,i}
    \end{align*}
\end{itemize}

There are two critical observations which must be made in order to see that this gives a memory-efficient implementation.
The first is that, since $\widehat{\vec{T}}$ is tridiagonal, $\vec{M}$, $\vec{N}$, and therefore $\vec{L}$ are all of half-bandwidth $q := \max(\deg(\tilde{M}),\deg(\tilde{N}))$.
This means that it is possible to compute the entries of $\vec{D}$ and the factors of $\vec{L}^{-1} = \vec{L}_{k-1}\cdots \vec{L}_1 \vec{L}_0$ one by one as we get the entries of $\widehat{\vec{T}}$.
The second is that because $\vec{L}$ is of bandwidth $q$, we can compute $[\vec{Q}\vec{L}^{-1}]_{:,i}$ without saving all of $\vec{Q}$.
More specifically, $[\vec{L}^{-1} \vec{M} \vec{e}_0]_i$ and $[\vec{Q} \vec{L}^{-1}]_{:,i}$ can respectively be computed from $\vec{L}_{j-1}\cdots \vec{L}_1 \vec{L}_0 \vec{M} \vec{e}_0$ and $\vec{Q} \vec{L}_0^\textup{\textsf{T}} \vec{L}_1^\textup{\textsf{T}} \cdots \vec{L}_{k-1}^\textup{\textsf{T}}$ and can therefore be maintained iteratively as the factors of $\vec{L}^{-1}$ are computed.
Moreover, because of the banded structure of the factors $\vec{L}_i$, we need only maintain a sliding window of the columns of $\vec{Q}\vec{L}^{-1}$ which will allow us to access the relevant columns when we need them and discard them afterwards.

For clarity of exposition, we only describe how to compute $\vec{M}$ and $\vec{N}$ in the case that $\tilde{M}$ and $\tilde{N}$ are degree at most two.
The rest of the subroutines are fully described for any degree.
We describe the syntax used in our pseudocode descriptions of streaming algorithms in detail in \cref{sec:streaming_notation}.
The syntax follows Python and other object oriented languages closely.

\subsection{Computing LDL factorization}

For the time being, we will assume that we can sequentially access the rows of $\vec{M}$ and $\vec{N}$.
Our first step is to compute an LDL factorization of $\vec{N}$.
A LDL factorization can be computed via a symmetrized version of Gaussian elimination and is guaranteed to exist if \( \vec{N} \) is positive definite \cite{higham_02}.
Gaussian elimination can be viewed as transforming the starting matrix \( \vec{N}_0 = \vec{N} \) to a diagonal matrix \( \vec{N}_{k-1} = \vec{D} \) via a sequence of row and column operations
\begin{align*}
    \vec{N}_{i+1} = \vec{L}_i \vec{N}_i \vec{L}_i^\textup{\textsf{T}}
\end{align*}
where
\begin{align*}
    \vec{L}_i :=
    \vec{I}_k + \vec{l}_i\vec{e}_i^\textup{\textsf{T}}
    ,&&
    \vec{l}_i :=
        \bigg[ \underbrace{\vphantom{\bigg|}0 , \cdots , 0}_{i+1\text{ zeros}} , -\frac{[\vec{N}_{i}]_{i+1,i}}{[\vec{N}_i]_{i,i}}  , \cdots , -\frac{[\vec{N}_i]_{k-1,i}}{[\vec{N}_i]_{i,i}}  \bigg]^\textup{\textsf{T}}.
\end{align*}
Note that the entries of \( \vec{L}_i \) are chosen to introduce zeros to the \( i \)-th row and column of \( \vec{N}_i \) such that \( [\vec{N}_{i+1}]_{:i+1,:i+1} \) is diagonal.
Therefore, if the algorithm terminates successfully, we will have obtained a factorization
\begin{align*}
    \vec{D} = (\vec{L}_{k-1} \cdots \vec{L}_1 \vec{L}_0) \vec{N} (\vec{L}_0^\textup{\textsf{T}} \vec{L}_1^\textup{\textsf{T}} \cdots \vec{L}_{n-1}^\textup{\textsf{T}})
\end{align*}
where \( \vec{D} \) is diagonal and each \( \vec{L}_i \) is unit lower triangular.
To obtain the factorization \( \vec{N} = \vec{L} \vec{D} \vec{L}^{\textup{\textsf{T}}} \), simply define  \( \vec{L} := (\vec{L}_{k-1} \cdots \vec{L}_1 \vec{L}_0)^{-1} \) and note that
\begin{align*}
    \vec{L} = \vec{I}_k - \sum_{i=0}^{k-1} \vec{l}_i \vec{e}_i^\textup{\textsf{T}}.
\end{align*}
We remark that that $\vec{l}_{k-1}$ is the zeros vector and is only included in sums for ease of indexing later on.
For further details on LDL factorizations, we refer readers to \cite{higham_02}.

To implement a LDL factorization, observe that the procedure above defines a recurrence
\begin{align*}
    [\vec{D}]_{j,j} &= [\vec{N}]_{j,j} - \sum_{\ell=0}^{j-1} [\vec{L}_{j,\ell}]^2 [\vec{D}]_{\ell,\ell}
    \\
    [\vec{L}]_{i,j} &= \frac{1}{[\vec{D}]_{j,j}} \left( [\vec{N}]_{i,j} - \sum_{\ell=0}^{j-1}[\vec{L}]_{j,\ell} [\vec{L}]_{i,\ell} [\vec{D}]_{\ell,\ell} \right)
    ,&& i > j.
\end{align*}
We therefore have \cref{alg:LDL}.
\begin{algorithm}[H]
\caption{LDL factorization}\label{alg:LDL}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
    \Procedure{LDL}{$\vec{N}$}
    \For{\(j=0:k\)}
    \State \( [\vec{D}]_{j,j} = [\vec{N}]_{j,j} - \sum_{\ell=0}^{j-1} [\vec{L}_{j,\ell}]^2 [\vec{D}]_{\ell,\ell} \)
    \State \( [\vec{L}]_{j,j} = 1 \)
    \For{\(i=j+1:k\)}
    \State \( [\vec{L}]_{i,j} = (1/[\vec{D}]_{j,j})( [\vec{N}]_{i,j} - \sum_{\ell=0}^{j-1}[\vec{L}]_{j,\ell} [\vec{L}]_{i,\ell} [\vec{D}]_{\ell,\ell} ) \)
    \EndFor
    \EndFor
    \State \Return \( \vec{L}, \vec{D} \)
\EndProcedure
\end{algorithmic}
\end{algorithm}

\subsubsection{Streaming version}

The fact that $\vec{L}$ has the same half bandwidth as $\vec{N}$ allows a more efficient implementation of \cref{alg:LDL} where terms which are known to be zero are not computed and only the important diagonals of $\vec{L}$ are stored.
Moreover, \cref{alg:LDL} already only accesses $\vec{N}$ one column at a time so it is easily converted to a streaming algorithm.
Making these changes gives the implementation \cref{alg:streaming_LDL} which is fed a stream of the columns of $\vec{N}$ in order, as shown in \cref{fig:streaming_LDL}.
Here the diagonal of $\vec{D}$ is stored as $\texttt{d}$ and the $(j+1)$-st diagonal of $\vec{L}$ is stored as $[\texttt{L}]_{j,:}$.
Thus, $\vec{L}_{i,j} = [\texttt{L}]_{i-j-1,j}$ as long as $i-j \in \{0,1,\ldots, q\}$.

\begin{algorithm}[H]
\caption{Streaming LDL}\label{alg:streaming_LDL}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Class{streaming-LDL}{$q,k$}
\State{stream: $[\vec{N}]_{0,0:q+1}, [\vec{N}]_{1,1:q+2}, \ldots, [\vec{N}]_{k-1,k-1:q+k-1}$}
\State $\texttt{L} = \textsc{zeros}(q,k)$
\State $\texttt{d} = \textsc{zeros}(k)$
\State $\texttt{j}\gets0$
\Procedure{read-stream}{$\vec{n}$}
\State $[\texttt{d}]_{\texttt{j}} \gets [\vec{n}]_0 - \sum_{\ell=\max(0,\texttt{j}-q)}^{\texttt{j}-1} [\texttt{L}]_{\texttt{j}-\ell-1,\ell}^2 [\texttt{d}]_{\ell}$
\For{$i=\texttt{j+1} : \min(\texttt{j}-q+1,n)$}
\State $[\texttt{L}]_{i-\texttt{j}-1,\texttt{j}} \gets (1/[\texttt{d}]_{\texttt{j}})([\vec{n}]_{i-\texttt{j}} - \sum_{\ell=\max(0,i-q)}^{i-1} [\texttt{L}]_{i-\ell-1,\ell} [\texttt{L}]_{\texttt{j}-\ell-1,\ell} [\texttt{d}]_{\ell} )$
\EndFor
\State $\texttt{j}\gets \texttt{j}+1$
\EndProcedure
\EndClass
\end{algorithmic}
\end{algorithm}

\begin{figure}
    \centering
    \begin{subfigure}{.48\textwidth}\centering
    \includegraphics[height=4.5cm]{imgs/streaming_LDL.pdf}
    \caption{Pattern for $\vec{N}$ in \cref{alg:streaming_LDL}.}
    \label{fig:streaming_LDL}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.48\textwidth}\centering
    \includegraphics[height=4.5cm]{imgs/streaming_tridiag_square.pdf}
    \caption{Pattern for $\vec{T}$ in \cref{alg:streaming_tridiag_square}.}
    \label{fig:streaming_tridiag_square}
    \end{subfigure}
    \\
    \begin{subfigure}{\textwidth}\centering
    \includegraphics[height=4.5cm]{imgs/streaming_binv_Q.pdf}
    \hfill
    \includegraphics[height=4.5cm]{imgs/streaming_binv_L.pdf}
    \hfill
    \includegraphics[height=4.5cm]{imgs/streaming_binv_d.pdf}
    \hfill
    \includegraphics[height=4.5cm]{imgs/streaming_binv_M.pdf}
    \caption{Pattern for $\vec{Q}$, $\vec{L}$, $\vec{d}$, $\vec{M}$ in \cref{alg:streaming_banded_prod}.}
    \label{fig:streaming_banded_prod}
    \end{subfigure}
    \caption{Access patterns for inputs to streaming functions used in low-memory implementations of Lanczos-OR and Lanczos-FA.
    Indices indicate what information should be streamed into the algorithm at the given iteration.}
\end{figure}

\subsection{Inverting the LDL factorization}

Once we have computed a a factorization \( \vec{N} = \vec{L}\vec{D}\vec{L}^\textup{\textsf{T}} \), we can easily evaluate $\vec{Q}\vec{L}^{-1}\vec{D}^{-1} \vec{L}^{-\textup{\textsf{T}}}\vec{M}\vec{e}_1$ using the fact that \( \vec{L}^{-1} = \vec{L}_{k-1} \cdots \vec{L}_1 \vec{L}_0 \).
Moreover, because the \( \vec{L}_j \) can be computed one at a time, there is hope that we can derive a memory efficient implementation.

Towards this end, define \( \vec{y}_j := \vec{L}_{j-1} \cdots \vec{L}_1 \vec{L}_0 \vec{M}\vec{e}_0 \) and \( \vec{X}_j := \vec{Q} \vec{L}_0^\textup{\textsf{T}} \vec{L}_1^\textup{\textsf{T}} \cdots \vec{L}_{j-1}^\textup{\textsf{T}}  \).
Then, setting \( \vec{y}_{0} = \vec{M} \vec{e}_1 \) we have that
\begin{align*}
    \vec{y}_{j+1}
    = \vec{L}_{j} \vec{y}_{j}
    = (\vec{I} + \vec{l}_j \vec{e}_j^\textup{\textsf{T}} ) \vec{y}_{j}
    = \vec{y}_{j} +  (\vec{e}_{j}^\textup{\textsf{T}} \vec{y}_{j}) \vec{l}_j.
\end{align*}
Similarly, setting \( \vec{X}_{0} = \vec{Q} \) we have that
\begin{align*}
    \vec{X}_{j+1}
    = \vec{X}_{j} \vec{L}_j^\textup{\textsf{T}}
    = \vec{X}_{j} ( \vec{I} + \vec{e}_j \vec{l}_j^\textup{\textsf{T}})
    = \vec{X}_{j} + \vec{X}_{j} \vec{e}_j \vec{l}_j^\textup{\textsf{T}}.
\end{align*}
Then \( \vec{Q}\vec{L}^{-1}\vec{D}^{-1} \vec{L}^{-\textup{\textsf{T}}}\vec{M}\vec{e}_1 = \vec{X}_{k} \vec{D}^{-1} \vec{y}_{k} \) can be computed accessing \( \vec{L} \), and therefore \( \vec{N} \), column by column.

\subsubsection{Streaming version}
\label{sec:streaming_banded_solve}

Recall that $[\vec{l}_i]_{:,\ell}$ is zero if $\ell \leq i$ or $\ell > i +q$.
Since $[\vec{l}_i]_{:i}$ is zero, we have
\begin{align*}
    [\vec{y}_j]_j
    &= [\vec{y}_{j} +  (\vec{e}_{j}^\textup{\textsf{T}} \vec{y}_{j}) \vec{l}_j]_j
    = [\vec{y}_{j+1}]_j
    =\cdots
    = [\vec{y}_{k}]_j
\intertext{and}
    [\vec{X}_{j}]_{:,j}
    &= [\vec{X}_{j} + \vec{X}_{j} \vec{e}_j \vec{l}_j^\textup{\textsf{T}}]_{:,j}
    = [\vec{X}_{j+1}]_{:,j}
    = \cdots = [\vec{X}_{k}]_{:,j}.
\end{align*}
We therefore have that
\begin{align*}
    \vec{X}_{k} \vec{D}^{-1} \vec{y}_{k}
    = \sum_{j=0}^{k-1} \frac{[\vec{y}_{k}]_j}{[\vec{D}]_{j,j}} [\vec{X}_{k}]_{:,j}
    = \sum_{j=0}^{k-1} \frac{[\vec{y}_{j}]_j}{[\vec{D}]_{j,j}} [\vec{X}_{j}]_{:,j}.
\end{align*}
Similarly, since $[\vec{l}_i]_{i+q+1:}$ is zero,
\begin{align*}
    [\vec{y}_j]_{j+q:}
    &= [\vec{y}_{j-1} + (\vec{e}_{j-1}^\textup{\textsf{T}}\vec{y}_{j-1}) \vec{l}_{j-1}]_{j+q:}
    = [\vec{y}_{j-1}]_{j+q:}
    = \cdots
    = [\vec{y}_0]_{j+q:}
\intertext{and}
    [\vec{X}_{j}]_{:,j+q:}
    &= [\vec{X}_{j-1}+\vec{X}_{j-1}\vec{e}_{j-1}\vec{l}_{j-1}^\textup{\textsf{T}}]_{:,j+q:}
    = [\vec{X}_{j-1}]_{:,j+q:}
    =\cdots
    = [\vec{X}_{0}]_{:,j+q:}.
\end{align*}
By definition, $\vec{y}_0 = \vec{M}\vec{e}_0$ and $\vec{X}_0 = \vec{Q}$.
Thus, we see that it is not necessary to know the later columns of $\vec{X}_j$ immediately.

We can define a streaming algorithm by maintaining only the relevant portions of the $\vec{X}_i$ and $\vec{y}_i$.
Towards this end, define the length \( q+1 \) vector \( \bar{\vec{y}}_j := [\vec{y}_j]_{j:j+q+1} \) the \( n\times (q+1) \) matrix \( \bar{\vec{X}}_j = [\vec{X}_j]_{:,j:j+q+1} \).
Using the above observations, we see that these quantities can be maintained by the recurrences
\begin{align*}
    \bar{\vec{y}}_j &=
    \begin{bmatrix}
        [\bar{\vec{y}}_{j-1}]_{1:}\\0
    \end{bmatrix}
    +  [\bar{\vec{y}}_{j-1}]_0 [\vec{l}_j]_{j+1:j+q+1}
    \intertext{and}
    [\bar{\vec{X}}_j]_{:,:q}
    &= [\bar{\vec{X}}_{j-1}]_{:,1:}
        + ( [\bar{\vec{X}}_{j-1}]_{:,1} ) ([\vec{l}_j]_{j+1:j+q+1})^\textup{\textsf{T}}.
    \\ [\bar{\vec{X}}_j]_{:,q} &= [\vec{Q}]_{j+q}
\end{align*}
Note then that,
\begin{align*}
    \vec{X}_{k-1} \vec{D}^{-1} \vec{y}_{k-1}
    = \sum_{j=0}^{k-1}  \frac{[\bar{\vec{y}}_j]_{1}}{[\vec{D}]_{j+1,j+1}} [\bar{\vec{X}}_j]_{:,1}.
\end{align*}

This results in \cref{alg:streaming_banded_prod} whose streaming pattern is outlined in \cref{fig:streaming_banded_prod}.

\begin{algorithm}[H]
\caption{Streaming banded product}\label{alg:streaming_banded_prod}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Class{streaming-banded-prod}{$n,k,q$}
\State{stream:$ $}
\State $\texttt{X\_} \gets \textsc{zeros}(n,q+1)$
\State $\texttt{y\_} \gets \textsc{zeros}(q+1)$
\State $\texttt{out} \gets \textsc{zeros}(n)$
\State $\texttt{j}\gets-1$
\Procedure{read-stream}{$\vec{v},\vec{l},d,\vec{y}_0$}
\If{$\texttt{j}=-1$}
\State $[\texttt{X\_}]_{:,:q} = \vec{v}$
\Else
\If{$\texttt{i}=-1$}
\State $\texttt{y\_} \gets \vec{y}_0$
\EndIf
\State $\texttt{out} \gets \texttt{out} + ([\texttt{y\_}]_{0}/d)[\texttt{X\_}]_{:,0}$
\State $[\texttt{y\_}]_{:q} \gets [\texttt{y\_}]_{1:} - [\texttt{y\_}]_{0} \vec{l}$
\State $[\texttt{y\_}]_{-1} \gets 0$
\State $[\texttt{X\_}]_{:,-1} \gets \vec{v}$
\State $[\texttt{X\_}]_{:,:q} \gets [\texttt{X\_}]_{:,1:} + [\texttt{X\_}]_{:,0} \vec{l}^\textup{\textsf{T}}$
\EndIf
\State $\texttt{j}\gets \texttt{j}+1$
\EndProcedure
\EndClass
\end{algorithmic}
\end{algorithm}

\begin{algorithm}[H]
\caption{Streaming banded inverse}\label{alg:streaming_banded_inv}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Class{streaming-banded-inv}{$n,k,q$}
\State{stream:$ $}
\State $\texttt{LDL} \gets \textsc{streaming-LDL}(k,q)$
\State $\texttt{Q0} \gets \textsc{zeros}(n,q)$
\State $\texttt{j}\gets0$
\Procedure{read-stream}{$\vec{q},\vec{n},\vec{y}_0$}
\If{$\texttt{j}<q$}
\State $[\texttt{Q0}]_{:,j} \gets \vec{v}$
\If{$\texttt{j}=q-1$}
    \State $\texttt{b-prod} \gets \textsc{streaming-banded-prod}(n,k,q)$
    \State $\texttt{b-prod.}\textsc{read-stream}(\texttt{V0},\texttt{none},\texttt{none},\texttt{none})$
\EndIf
\Else
\State $\texttt{LDL.}\textsc{read-stream}(\vec{n})$
\State $\texttt{b-inv.}\textsc{read-stream}(\vec{q},-[\texttt{LDL.}\texttt{L}]_{:,j-q},[\texttt{LDL.}\texttt{d}]_{j-q},\vec{y}_0)$
\EndIf
\State $\texttt{j}\gets \texttt{j}+1$
\EndProcedure
\EndClass
\end{algorithmic}
\end{algorithm}

\subsection{Computing polynomials in $\vec{T}$}

The last major remaining piece is to construct $\vec{M} = \tilde{M}(\vec{T})$ and $\vec{N} = [\tilde{N}(\widehat{\vec{T}})]_{:k,:k}$.
Recall that we have assumed $\tilde{M}$ and $\tilde{N}$ are of degree at most two for convenience.
In iteration $\ell$ of Lanczos, we obtain $\alpha_{\ell}$ and $\beta_{\ell}$.
Observe that $\vec{T}^2$ is symmetric and that, defining $\beta_{-1} = \beta_{k} = 0$, the lower triangle is given by
\begin{align*}
    [\vec{T}^2]_{i,j}
    =
    \begin{cases}
    \beta_{j-1}^2 + \alpha_j^2 + \beta_j^2 & j=i \\
    (\alpha_j + \alpha_{j+1}) \beta_i & j=i-1 \\
    \beta_j \beta_{j+1} & j = i-2 \\
    0 & \text{o.w.}
    \end{cases}
\end{align*}

We can use this to implement the streaming algorithm, \cref{alg:streaming_tridiag_square}, for computing the entries of $\vec{T}^2$.
Rather than being fed the entire tridiagonal matrix $\vec{T}$, \cref{alg:streaming_tridiag_square} is fed a stream of the columns of $\vec{T}$ in order, as shown in \cref{fig:streaming_tridiag_square}.
The algorithm respectively stores the $j$-th diagonals of $\vec{T}$ and $\vec{T}^2$  as $[\texttt{T}]_{j,:}$ and $[\texttt{Tp2}]_{j,:}$.

\begin{algorithm}[H]
\caption{Streaming tridiagonal square}\label{alg:streaming_tridiag_square}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Class{streaming-tridiagonal-square}{$k$}
\State{stream: $(\alpha_0, \beta_0), \ldots, (\alpha_{k-1},\beta_{k-1})$}
\State $\texttt{T} \gets \textsc{zeros}(2,k)$
\State $\texttt{Tp2} \gets \textsc{zeros}(3,k)$
\State $\texttt{j}\gets0$
\Procedure{read-stream}{$\alpha,\beta$}
\State $[\texttt{T}]_{0,\texttt{j}} = \alpha$
\State $[\texttt{T}]_{1,\texttt{j}} = \beta$
\If{$\texttt{i}=0$}
\State $[\texttt{Tp2}]_{0,\texttt{j}} \gets [\texttt{T}]_{0,\texttt{j}}^2+[\texttt{T}]_{1,\texttt{j}}^2$
\Else
\State $[\texttt{Tp2}]_{0,\texttt{j}} \gets [\texttt{T}]_{0,\texttt{j}}^2+[\texttt{T}]_{1,\texttt{j}}^2+[\texttt{T}]_{1,\texttt{j}-1}^2$
\State $[\texttt{Tp2}]_{1,\texttt{j}} \gets ([\texttt{T}]_{0,\texttt{j}}+[\texttt{T}]_{0,\texttt{j}-1})[\texttt{T}]_{1,\texttt{j}-1}$
\State $[\texttt{Tp2}]_{2,\texttt{j}} \gets [\texttt{T}]_{1,\texttt{j}}[\texttt{T}]_{1,\texttt{j}-1}$
\EndIf
\State $\texttt{j}\gets \texttt{j}+1$
\EndProcedure
\EndClass
\end{algorithmic}
\end{algorithm}

Since we maintain the columns of $\vec{T}^2$ with \cref{alg:streaming_tridiag_square}, we can easily compute $\vec{M}$ and $\vec{N}$ using \cref{alg:get_poly}.

\begin{algorithm}[H]
\caption{Get polynomial of tridiagonal matrix}\label{alg:get_poly}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Procedure{get-poly}{$P,\texttt{STp2},k,j$}
\State $a,b,c = P(0), P'(0), P''(0)$
\State $\vec{p} \gets \textsc{zeros}(3)$
\State $[\vec{p}]_{:3} \gets a [\texttt{STp2.}\texttt{Tp2}]_{:,j}$
\State $[\vec{p}]_{:2} \gets b [\texttt{STp2.}\texttt{T}]_{:,j}$
\State $[\vec{p}]_{:1} \gets c $
\EndProcedure
\end{algorithmic}
\end{algorithm}

\subsection{Putting it all together}

With this algorithm in place, putting everything together is straightforward, and the full implementation is shown in \cref{alg:rational_lanczos}.

\begin{algorithm}[H]
\caption{Streaming banded rational inverse}\label{alg:rational_lanczos}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Class{banded-rational}{$n,k,\tilde{M},\tilde{N}$}
\State $\texttt{b-inv} \gets \textsc{banded-inv}(n,k,2)$
\State $\texttt{STp2} \gets \textsc{streaming-tridiagonal-square}(k)$
\State $\texttt{j}\gets0$
\Procedure{read-stream}{$\vec{q},\alpha,\beta$}
\If{$\texttt{j}<k$}
    \State $\texttt{STp2.}\textsc{read-stream}(\alpha,\beta)$
    \State $\texttt{b-inv.}\textsc{read-stream}($
    \State \hspace{1em} $\vec{q},$
    \State \hspace{1em} $\textsc{getpoly}(\tilde{N},\texttt{STp2},k,\texttt{j}-1)~\textbf{if}~\texttt{j}\geq 2~\textbf{else}~\texttt{none},$
    \State \hspace{1em} $\textsc{getpoly}(\tilde{M},\texttt{STp2},k,\texttt{j}-1)~\textbf{if}~\texttt{j}= 2~\textbf{else}~\texttt{none},$
    \State$)$
\EndIf
\State $\texttt{LDL.}\textsc{read-stream}(\vec{n})$
\State $\texttt{j}\gets \texttt{j}+1$
\EndProcedure
\Procedure{finish-up}{$ $}
\For{$i=k:k+2$}
    \State $\texttt{b-inv.}\textsc{read-stream}( \texttt{none}, \textsc{getpoly}(\tilde{N},\texttt{STp2},k,\texttt{j}-1), \texttt{none} )$
\EndFor
\EndProcedure
\Procedure{get-output}{$ $}
\State \Return \texttt{b-inv.}\texttt{b-prod.}\texttt{out}
\EndProcedure
\EndClass
\end{algorithmic}
\end{algorithm}

This can be incorporated into any Lanczos implementation and used to compute the Lanczoz-OR iterates.
For concreteness, we show this with the implementation of Lanczos from \cref{alg:lanczos}.
We call the resulting implementation Lanczos-OR-lm.

\begin{algorithm}[H]
\caption{Lanczos-OR-lm}
 \label{alg:main}
\fontsize{10}{10}\selectfont
\begin{algorithmic}[1]
\Procedure{Lanczos-OR-lm}{$\vec{A},\vec{b},k,M,N$}
\State \( \vec{q}_{-1} = \vec{0} \),
\( \beta_{-1} = 0 \),
\( \vec{q}_0  = \vec{b} / \| \vec{b} \| \)
\State Set $\tilde{M}$ and $\tilde{N}$ as in \cref{def:lanczos_or}
\State $\texttt{lam-lm} \gets \textsc{banded-rational}(n,k,\tilde{M},\tilde{N})$
\For {\( j=0:k \)}
    \State \( \tilde{\vec{q}}_{j+1} = \vec{A} \vec{q}_{j} - \beta_{j-1} \vec{q}_{j-1} \)
    \State \( \alpha_j = \langle \tilde{\vec{q}}_{j+1}, \vec{q}_j \rangle \)
    \State \( \tilde{\vec{q}}_{j+1} = \tilde{\vec{q}}_{j+1} - \alpha_j \vec{q}_j \)
        \State optionally, reorthogonalize \( \tilde{\vec{q}}_{j+1} \) against \( \{\vec{q}_i\}_{i=0}^{j-1} \)
    \State \( \beta_{j} = \| \tilde{\vec{q}}_{j+1} \| \)
    \State \( \vec{q}_{j+1} = \tilde{\vec{q}}_{j+1} / \beta_{j} \)

    \State $\texttt{lam-lm.}\textsc{read-stream}(\vec{q}_{j},\alpha_{j},\beta_{j})$

\EndFor
\State $\texttt{lam-lm.}\textsc{finish-up}()$
\EndProcedure
\end{algorithmic}
\end{algorithm}

We can easily obtain an implementation of Lanczos-FA, which we call Lanczos-FA-lm, by replacing $\beta_{k-1}$ with $0$ in the final iteration of the loop.

\subsection{Some comments on implementation}

Our main goal is to describe how to implement Lanczos-FA and Lanczos-OR in a way that requires $k$ matrix-vector products and $O(n)$ storage, when $M$ and $N$ are each at most degree two. As mentioned, the approach can be extended to any constant degree.
There are a range of improvements to our implementation which may be useful in practice.

First, the amount of storage used can be reduced somewhat.
Indeed, the implementation described above saves $\vec{T}$, $\vec{T}^2$, $\vec{L}$, and $\vec{d}$, but only accesses a sliding window of these quantities.
We have chosen to save them for convenience since they require only $O(k)$ storage.
However, storing only the relevant information from these quantities would result in an implementation with storage costs independent of the number of iterations $k$.
In this vein, a practical implementation would likely determine $k$ adaptively by monitoring the residual or other measures of the error.

Improvements to the number of vectors of length $n$ may be possible as well, although we expect these would limited to constant factors as well. For example, storage could possibly be reduced somewhat by incorporating the Lanczos iteration more explicitly with the inversion of the LDL facorization, much like the classical Hestenes and Stiefel implementation of CG \cite{hestenes_stiefel_52}.

\subsection{Lanczos-FA-lm and Lanczos-OR-lm in finite precision arithmetic}
\label{sec:stability}

As with other short-recurrence based Krylov subspace methods, the behavior of Lanczos-FA-lm and Lanzos-OR-lm in finite precision arithmetic may be different than in exact arithmetic.
Fortunately, quite a bit is known about the standard implementation of Lanczos \cite{paige_71,paige_76,paige_80,greenbaum_89,musco_musco_sidford_18}, and we have stated Lanczos-FA-lm and Lanczos-OR-lm in terms of this implementation.
Knowledge about the standard implementation of Lanczos carries over to Lanczos-FA.
For instance, assuming $\vec{Q} f(\vec{T}) \vec{e}_0$ is computed accurately from the output of the standard Lanczos algorithm, then many error bounds for Lanczos-FA are still applicable \cite{musco_musco_sidford_18,chen_greenbaum_musco_musco_21}.
It is more or less clear that Lanczos-FA-lm and Lanczos-OR-lm will accurately compute the expressions $\vec{Q} N(\vec{T})^{-1} M(\vec{T}) \vec{e}_0$ and $\vec{Q}  ([\tilde{N}(\widehat{\vec{T}})]_{:k,:k})^{-1} \tilde{M}(\vec{T}) \vec{e}_0$ provided that $M(\vec{T})$ $N(\vec{T})$, $[\tilde{N}(\widehat{\vec{T}})]_{:k,:k}$ are reasonably well conditioned.
Indeed, in practice solving linear systems by symmetric Gaussian elimination is accurate; see for instance \cite[Chapter 10]{higham_02}.
Thus, such bounds and techniques can be applied to Lanczos-FA-lm and Lanczos-OR-lm.

\iffalse
However, our implementation is equivalent, even in finite precision arithmetic, to \begin{itemize}
    \item running Lanczos to obtain $\vec{Q}$ and $\vec{T}$
    \item computing $\vec{M}$ and $\vec{N}$ via certain standard approachs
    \item computing an LDL factorizaton of $\vec{N}$ via a certain standard approach
    \item computing $\vec{Q} \vec{L}^{-1} \vec{D}^{-1} \vec{L}^{-\textup{\textsf{T}}}\vec{M}\vec{e}_0$ in an essentially standard way.
\end{itemize}
\fi

\section{Numerical experiments and comparison to related algorithms}
\label{sec:comparison}

We now provide several examples which illustrate various aspects of the convergence properties of Lanczos-OR and Lanczos-OR based algorithms, and show when these new methods can outperform more standard techniques like the classic Lanczos-FA.

\subsection{The matrix sign function}

As we noted in \cref{sec:sign_function}, Lanczos-OR can be used to obtain an approximation to the matrix sign function.
A related approach, which interpolates the sign function at the so called ``harmonic Ritz values'', is described in \cite[Section 4.3]{eshof_frommer_lippert_schilling_van_der_vorst_02}.
The harmonic Ritz values are characterized by the generalized eigenvalue problem
\begin{align*}
     [\widehat{\vec{T}}^2]_{:k,:k} \vec{y} = \theta \vec{T} \vec{y}
\end{align*}
and are closely related to MINRES in the sense that MINRES produces a polynomial interpolating $1/x$ at the harmonic Ritz values \cite{paige_parlett_vandervorst_95}.

\begin{example}

We construct a matrix with $400$ eigenvalues, 100 of which are the negatives of the values of a model problem \cite{strakos_91,strakos_greenbaum_92} with parameters $\kappa = 10^2$, $\rho=0.9$, and $n=100$ and 300 of which are the values of a model problem with parameters $\kappa = 10^3$, $\rho=0.8$, $n=300$.
Here, the model problem eigenvalues are given by
\begin{align*}
    \lambda_1 = 1
    ,&& \lambda_n = \kappa
    ,&& \lambda_i = \lambda_1 + \left( \frac{i-1}{n-1} \right) \cdot (\kappa -1) \cdot \rho^{n-i}, && i=2,\ldots,n-1.
\end{align*}

We compute the Lanczos-OR induced approximation, the Lanczos-FA approximation, the harmonic Ritz value based approximation from \cite{eshof_frommer_lippert_schilling_van_der_vorst_02}, and the optimal $\vec{A}^2$-norm approximation to the matrix sign function.
The results are shown in \cref{fig:sign_opt}.
In all cases, we use the Lanczos algorithm with full reorthogonalization.
Because eigenvalues of $\vec{T}$ may be near to zero, Lanczos-FA exhibits oscillatory behavior
On the other hand, the Lanczos-OR based approach and the harmonic Ritz value based approach have much smoother convergence.
Note that the Lanczos-OR induced approximation is not optimal, although it seems to perform close to optimally after a few iterations.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{imgs/sign.pdf}
    \caption{Comparison of $\vec{A}^2$-norm errors for approximating $\operatorname{sign}(\vec{A})\vec{b}$ (normalized by $\|\vec{b}\|_{\vec{A}^2}$).
    \emph{Legend}: Lanczos-OR induced approximation ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_solid.pdf}}}),
    interpolation at harmonic Ritz values ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dot.pdf}}}),
    Lanczos-FA ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dash_thin.pdf}}}),
    $\vec{A}^2$-norm optimal ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dashdot_thin.pdf}}}).
    \emph{Left}: $\vec{A}^2$-norm of error.
    \emph{Right}: optimality ratio.
    \emph{Remark}: Lanczos-OR exhibits smoother convergence than Lanczos-FA and is nearly optimal.
    }
    \label{fig:sign_opt}
\end{figure}

\end{example}

\begin{example}

In this example, we show the spectrum approximations induced by the algorithms from the previous example.
We now set $\vec{A}$ to be a diagonal matrix with $1000$ eigenvalues set to the quantiles of a Chi-squared distribution with parameters $\alpha=1$ and $\beta=10$.
We set $k=10$ and consider approximations to the function $c\mapsto \vec{b}^\textup{\textsf{T}} \ensuremath{\mathbbm{1}} [\vec{A} \leq c] \vec{b}$ for a range of values $c$.
Here $\ensuremath{\mathbbm{1}}[x\leq c] = (1-\operatorname{sign}(x-c))/2$ is one if $x\leq c$ and zero otherwise.
We pick $\vec{b}$ as a unit vector with equal projection onto each eigencomponent so that $\vec{b}^\textup{\textsf{T}} \ensuremath{\mathbbm{1}} [\vec{A} \leq c] \vec{b}$ gives the fraction of eigenvalues of $\vec{A}$ below $c$.
In the $n\to\infty$ limit, this function will converge pointwise to the cumulative distribution of a Chi-squared random distribution with parameters $\alpha=1$ and $\beta=10$.
The results are shown in \cref{fig:spec_approx}.

Note that the Lanczos-FA based approach is piecewise constant with jumps at each eigenvalue of $\vec{T}$.
On the other hand, the harmonic Ritz value and Lanczos-OR based approaches produce continuous approximations to the spectrum.
In this particular example, the spectrum of $\vec{A}$ is near to a smooth limiting density, so the harmonic Ritz value and Lanczos-OR based approaches seem to produce better approximations.

In general it is not possible to pick $\vec{b}$ with equal projection onto each eigencomponent since the eigenvectors of $\vec{A}$ are unknown.
However, by choosing $\vec{b}$ from a suitable distribution, it can be guaranteed that $\vec{b}$ has roughly equal projection onto each eigencomponent.
In this case, the Lanczos based approach above is referred to as stochastic Lanczos quadrature \cite{chen_trogdon_ubaru_21},

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{imgs/spec_approx.pdf}
    \caption{Comparison of spectrum approximations.
        \emph{Legend}: Lanczos-OR induced approximation ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_solid.pdf}}}),
        Lanczos-FA ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dash_thin.pdf}}}),
        harmonic Ritz values based approximation ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dashdot_thin.pdf}}}),
        limiting density
        ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dot_thin.pdf}}}).
        .
        \emph{Remark}: The Lanczos-OR and Harmonic Ritz value based approaches provide smooth approximations which match the smooth limiting density much better than the piecewise constant approximation computed with Lanczos-FA.
        }
    \label{fig:spec_approx}
\end{figure}
\end{example}

\subsection{Rational matrix functions}

We now illustrate the effectiveness of the Lanczos-OR based approach to approximating rational matrix functions described in \cref{sec:rational_func}.
Then we compare an existing low-memory approach, called multishift CG, to the analogous approaches based on Lanczos-OR-lm and Lanczos-FA-lm.

Throughout, this section, we will assume that $r$ is a rational function of the form
\begin{align}
\label{eqn:rat_form}
r(x) = \sum_{i=1}^{m} \frac{A_ix^2+B_ix+C_i}{a_ix^2+b_ix+c_i}.
\end{align}
This is relatively general since any real valued rational function $r:\mathbb{R}\to\mathbb{R}$ with numerator degree smaller than denominator degree and only simple poles can be written in this form (in fact, this would be true even if $A_i=0$).
A range of rational functions of this form appear naturally; for instance by a quadrature approximation to a Cauchy integral formula representation of $f$ \cite{hale_higham_trefethen_08}.
Similar rational functions are seen in \cite{eshof_frommer_lippert_schilling_van_der_vorst_02,frommer_simoncini_09}.

It is clear that $r(\vec{A})\vec{b}$ has the form
\begin{align*}
    r(\vec{A})\vec{b} = \sum_{i=1}^{m} (A_i\vec{A}^2 + B_i\vec{A} + C_i\vec{I}) \vec{x}_i
\end{align*}
where $\vec{x}_i$ is obtained by solving the linear system of equations $(a_i\vec{A}^2 + b_i \vec{A} + c_i \vec{I}) \vec{x}_i = \vec{b}$.

\begin{example}

In this example, we use the same spectrum as in the first example.
However, rather than approximating the sign function directly, we instead use Lanczos-OR to approximate each term of a proxy rational function of the form \cref{eqn:rat_form}.
In particular, we consider the best uniform approximation\footnote{Note that the eigenvalues of $\vec{A}$ live in $[-10^2,-1]\cup[1,10^3]$, so we could have used an asymmetric approximation to the sign function.
This would reduce the degree of the rational function required to obtain an approximation of given accuracy, but the qualitative behavior of Lanczos-OR-lm would not change substantially.} of degree $(39,40)$ to the sign function on $[-10^3,1]\cup[1,10^3]$.
Such an approximation is due to Zolotarev \cite{zolotarev_77}, an can be derived from the more well known Zolotarev approximation to the inverse square root function on $[1,10^6]$.
Our implementation follows the partial fractions implementation in the Rational Krylov Toolbox \cite{berljafa_elsworth_guttel_20} and involves computing the sum of $20$ terms of degree $(1,2)$.
The results are shown in \cref{fig:proxy_rat}.

\begin{figure}
    \label{fig:proxy_rat}
    \centering
    \includegraphics[width=\textwidth]{imgs/Lanczos-OR-lm_delay.pdf}
    \caption{$\vec{A}^2$-norm error in Lanczos-OR-lm based rational approximation to matrix sign function (normalized by $\|\vec{b}\|_{\vec{A}^2}$).
        \emph{Legend}:
        Lanczos-OR-lm based approximation of matrix sign function without/with reorthogonalization
        ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_solid.pdf}}}/{\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_solid_thin.pdf}}}),
        Lanczos-OR-lm based approximation of proxy rational matrix function without/with reorthogonalization
        ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dash.pdf}}}/{\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dash_thin.pdf}}}),
        Infinity norm error of proxy rational function approximation ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dot_thin.pdf}}}).
        \emph{Remark}:
        The convergence of Lanczos-OR to the sign function matches the convergence to the proxy rational function until the error is of the order of the error in the proxy rational function.
        }
\end{figure}
\end{example}

At least while while the application error for the Lanczos-OR approximation to the proxy rational matrix function is large relative to the approximation error, then as seen in \cref{eqn:triangle_ineq}, the error in approximating the matrix sign function is similar to the error in approximating the proxy rational matrix function.
However, as seen in \cref{eqn:reverse_triangle_ineq}, the final accuracy of approximating the matrix sign function is limited by the quality of the scalar approximation.

We also note that it really only makes sense to use Lanczos-OR-lm with a short-recurrence version of Lanczos, in which case the effects of a perturbed Lanczos recurrence are prevalent.
In particular, the algorithm encounters a delay of convergence as compared to what would happen with reorthogonalization.
This is because the example problem's spectrum has many outlying eigenvalues, so the Lanczos algorithm quickly loses orthogonality and begins to find ``ghost eigenvalues'' \cite{meurant_strakos_06,liesen_strakos_13}.

\subsubsection{Comparison of Lanczos-OR-lm, Lanczos-FA-lm, and multishift CG/MINRES}

In certain cases, the shift invariance of Krylov subspace can be used to simultaneously compute all of the $\vec{x}_i$ using the same number of matrix-vector products as would be required to approximate a single $\vec{x}_i$.
Specifically, suppose $(a_i\vec{A}^2 + b_i \vec{A} + c_i \vec{I})^{-1} \vec{b}$ can be written as $\vec{B}+z_i \vec{I}$ for all $i=1,\ldots, m$.
Then $\mathcal{K}_k(\vec{B}-z_i\vec{I},\vec{b})) = \mathcal{K}_k(\vec{B},\vec{b})$, so by constructing a single Krylov subspace $\mathcal{K}_k(\vec{B},\vec{b})$, one can compute all of the $\vec{x}_i$ and therefore $r(\vec{A})\vec{b}$.
The resulting algorithms are typically called multishift-CG or multishift-MINRES \cite{eshof_frommer_lippert_schilling_van_der_vorst_02,frommer_simoncini_08,frommer_simoncini_08,guttel_schweitzer_21,pleiss_jankowiak_eriksson_damle_gardner_20}.
However, such an approach only works when $a_i = 0$ or $a_i = a$ and $b_i = b$, and in the latter case, matrix-vector products with $\vec{B}$ require two matrix-vector products with $\vec{A}$ and convergence depends only on the properties of $\vec{A}^2$ rather than $\vec{A}$.

Instead, we might apply Lanczos-FA-lm to compute individual terms of $r(\vec{A})\vec{b}$.
However, if $(a_i\vec{A}^2 + b_i \vec{A} + c_i \vec{I})$ is indefinite, then Lanczos-FA may exhibit oscillatory behavior due to eigenvalues of $(a_i\vec{T}^2 + b_i \vec{T} + c_i \vec{I})$ near zero.
This may result in a breakdown of Lanczos-FA-lm similar to the breakdown which may be encountered by standard implementations of CG on indefinite linear systems.
Lanczos-OR-lm avoids such issues.

\iffalse
For the moment, suppose that $a_i = 0$ for all $i$; i.e., that the rational function has only real poles.
This is the case when $f$ is a Stieltjes function and $r$ is derived from a quadrature approximation.
Then, defining $h_z(x) := 1/(x-z)$, we see that,
\begin{align*}
    \textup{\textsf{lan-FA}}_k(h_0 ; \vec{A} - a_i \vec{I}, \vec{b} )
    = \textup{\textsf{lan-FA}}_k(h_{a_i} ; \vec{A}, \vec{b} ).
\end{align*}
Critically, the same Krylov subspace, $\mathcal{K}_k(\vec{A},\vec{b})$, can be used for the CG approximations to all of the $\vec{x}_i$.
This can be implemented into an algorithm which uses just $k$ matrix-vector products with $\vec{A}$ and a few vectors per term in the rational function.
We refer to the resulting algorithm as multishift CG and remark that it is equivalent to Lanczos-FA applied to the corresponding rational matrix function in exact arithmetic.

Now suppose $a_i = a \neq 0$ and $b_i = b$ are constant for all $i$; i.e., all of the complex poles have the same real part.
This is the case in the Zoloterov approximations to the absolute value and sign functions \cite{berljafa_elsworth_guttel_20}.
Then $a\vec{A}^2+b \vec{A} + c_i\vec{I}$ can be factored as $a(\vec{A}+\frac{b}{2a} \vec{I}) + (c_i-\frac{b^2}{4a})\vec{I}$, so multishift CG can be used with $(\vec{A}+\frac{b}{2a}\vec{I})^2$ to obtain all of the $\vec{x}_i$ efficiently.
Note that the resulting algorithm is no longer equivalent to Lanczos-FA as after $k$ matrix-vector products with $\vec{A}$; indeed, the $\vec{x}_i$ live in the space $\mathcal{K}_{k/2}((\vec{A}+\frac{b}{2a}\vec{I})^2,\vec{b})$ since matrix-vector products with $(\vec{A}+\frac{b}{2a}\vec{I})^2$ require two matrix-vector products with $\vec{A}$.

If we require two matrix-products per iteration, one may instead try to use two-pass Lanczos-FA.
This has the advantage of being usable with any function $f$.
However, if $\vec{A}+\frac{b}{2a}\vec{I}$ is indefinite and $c_i$ is small, then Lanczos-FA may exhibit oscillatory behavior due to the appearance of Ritz values near $-b/2a$.
However, Lanczos-OR does not encounter such oscillations due to it's the optimally condition.

\fi

\iffalse
The fact that Lanczos-FA-lm is equivalent to Lanczos-FA is important.
In the low-memory setting, is is typically infeasible to store $\vec{A}$, and therefore $(\vec{A}+\frac{b}{2}\vec{I})^2$.
As such, products with $(\vec{A}+\frac{b}{2}\vec{I})^2$ must be computed by successive products with $\vec{A}+\frac{b}{2}\vec{I}$.
Thus, if only $k$ matrix-vector products with $\vec{A}$ are allowed, the number of iterations of CG is $k/2$ and the $\vec{y}_i$ are contained in $\mathcal{K}_{k/2}((\vec{A}+\frac{b}{2}\vec{I})^2,\vec{b})$.
It's not hard to see that $\mathcal{K}_{k/2}((\vec{A}+\frac{b}{2}\vec{I})^2,\vec{b})$ is a proper subspace of $\mathcal{K}_{k}(\vec{A}+\frac{b}{2}\vec{I},\vec{b}) = \mathcal{K}_k[\vec{A},\vec{b}]$.
The implications of this fact are discussed further in \cref{sec:msCG_comparison}.
\fi

The following example highlights some of the tradeoffs between the algorithms:

\begin{example}
\label{ex:lanczos_msCG_comparison}
We construct several test problems by placing eigenvalues uniformly throughout the specified intervals.
In all cases, $\vec{b}$ has uniform weight onto each eigencomponent.
The outputs are computed using standard Lanczos, but we note that the spectrum and number of iterations are such that the behavior is quite similar to if full reorthgonalization were used.
In particular, orthogonality is not lost since no Ritz value converges.
The results of our experiments are shown in \cref{fig:msCG_squared}.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{imgs/lanczos_msCG_squared.pdf}
    \caption{Comparison of $(\vec{A}^2+c\vec{I})$-norm errors for CG and Lanczos-FA for computing $(\vec{A}^{2}+c\vec{I})^{-1}\vec{b}$ (normalized by $\|\vec{b}\|_{\vec{A}^2+c\vec{I}}$).
    Here CG works with $\vec{A}^2+c\vec{I}$ and requires two matrix-vector products per iteration whereas Lanczos-FA works with $\vec{A}$ and requires just one.
    \emph{Legend}: Lanczos-OR ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_solid.pdf}}}),
    Lanczos-FA ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dash_thin.pdf}}}),
    CG on squared system ({\protect\raisebox{0mm}{\protect\includegraphics[scale=.7]{imgs/legend/blue_dashdot_thin.pdf}}}).
    \emph{Top left}: eigenvalues on $[1,10]$, $c=0.05$.
    \emph{Top right}: eigenvalues on $[-1.5,-1]\cup[1,10]$, $c=0.05$.
    \emph{Bottom left}: eigenvalues on $[-10,-1]\cup[1,10]$, $c=0.05$.
    \emph{Bottom right}: eigenvalues on $[-10,-1]\cup[1,10]$, $c=0$.
    \emph{Remark}: Lanczos-OR converges without oscillations while automatically matching the rate of convergence of the better of Lanczos-FA and CG on the squared system.
    }
    \label{fig:msCG_squared}
\end{figure}
\end{example}
In the first three examples, we consider approximations to $r(x) = 1/(x^2+0.05)$ with eigenvalues spaced with increments of $0.005$ throughout $[1,10]$, $[-1.5,-1]\cup[1,10]$, and $[-10,-1]\cup[1,10]$ respectively.
For all these examples, the condition number of $\vec{A}^2+0.05\vec{I}$ is roughly 100 and the eigenvalues of $\vec{A}^2+0.05 \vec{I}$ fill out the interval $[1,100.05]$.
As such we observe that multishift CG converges at a rate (in terms of matrix products with $\vec{A}$) of roughly $\exp(-k/\sqrt{\kappa(\vec{A}^2)}) = \exp(-k/\sqrt{100})$ on all of the examples.

In the first example, $\vec{A}$ is positive definite.
Here Lanczos-FA and Lanczos-OR converge similarly to CG on $\vec{A}$ at a rate of roughly $\exp(-2k/\sqrt{10})$, where $k$ is the number of matrix-vector products with $\vec{A}$.

In the next example $\vec{A}$ is indefinite.
The convergence of CG is unchanged, because CG acts on $\vec{A}^2 + c\vec{I}$, it is unable to ``see'' the asymmetry in the eigenvalues of $\vec{A}$.
While the convergence of Lanczos-FA and Lanczos-OR is slowed considerably, both methods converges more quickly than CG due to the asymmetry in the intervals to the left and the right of the origin.
The convergence of these methods is at a rate of roughly $\exp(-k/\sqrt{15})$, although the exact rate is more complicated to compute \cite{fischer_96,schiefermayr_11}.
We also note the the emergence of oscillations in the error curve of Lanczos-FA.

In the third example, the asymmetry in the eigenvalue distribution about the origin is removed, and Lanczos-FA and Lanczos-OR converge at a rate very similar to that of multishift CG.
Note that Lanczos-FA displays larger oscillations, since the symmetry of the eigenvalue distribution of $\vec{A}$ ensures that $\vec{T}$ has an eigenvalue at zero whenever $k$ is odd.
However, the size of the oscillations is regularized by the fact that $c>0$.

In the final example, we use the same eigenvalue distribution as the third example, but now apply the function $r(x) = x^{-2}$.
Here CG and Lanczos-OR behave essentially the same, but the behavior of Lanczos-FA becomes far more oscillatory.
Indeed, the lack of the regularizing term $c\vec{I}$ means that $r(\vec{T}) = \vec{T}^{-2}$ is not even defined when $\vec{T}$ has an eigenvalue at zero.
Lanczos-FA-lm will break down in such settings, as the LDL factorization of $\vec{T}^2$ is not well defined.
Even in less extreme situations, the LDL factorization may become inaccurate.

\section{Conclusion}

We have introduced an optimal algorithm for computing approximations to rational matrix functions over Krylov subspace.
This is accompanied by an low-memory implementation.
The particular case of denominator degree two is likely the most important as it can be used in conjunction with the rational function approximations for a range of functions, including the matrix sign function.

There are a range of interesting directions for future work.
A few of the most interesting are summarized here:
\begin{itemize}
    \item Can we provide a sharper comparison between the Lanczos-OR and Lanczos-FA approximations to the matrix sign function?
    \item Is the induced algorithm for the matrix sign function nearly optimal/can we derive simple spectrum dependent bounds?
    \item For what other functions can we use Lanczos-OR to induce a new algorithm?
     \item How does Lanczos-OR generalize ``harmonic Ritz values'', and can this perspective provide any insight into Lanczos-FA?
\end{itemize}

redpajama_set_name: RedPajamaArXiv