\section{Introduction}
\label{sec:intro}

This file is documentation for the SIAM \LaTeX\ style, including how
to typeset the main document, the {\scshape Bib}\TeX\xspace\ file, and any supplementary
material. More information
about SIAM's editorial style can be found in the style manual, available
at \url{https://www.siam.org/journals/pdf/stylemanual.pdf}.
The major changes in the SIAM online-only class are summarized in \cref{sec:changes}.
The SIAM \LaTeX\@ files can be found at
\url{https://www.siam.org/journals/auth-info.php}. The files that
are distributed are given below.
\begin{itemize}
\item \texttt{siamonline220329.cls} (required): Main SIAM online-only \LaTeX\ class file.
\item \texttt{siamplain.bst} (required): Bibliographic style file for
  {\scshape Bib}\TeX\xspace.
\item \texttt{docsiamonline.tex}: Produces this documentation.
\item \texttt{references.bib}: {\scshape Bib}\TeX\xspace\ database for this
  documentation and examples.
\item \texttt{ex\_article.tex}: Template for article.
\item \texttt{ex\_supplement.tex}: Template for supplement.
\item \texttt{ex\_shared.tex}: Template for shared information for
  article and supplement.
\end{itemize}
To use these files, put \texttt{siamonline220329.cls} and
\texttt{siamplain.bst} in the directory with your
paper or, alternatively, into your \LaTeX\@ and {\scshape Bib}\TeX\xspace\@ paths,
respectively.
The outline of a SIAM \LaTeX\ article is shown in
\cref{ex:outline}. Templates are provided and discussed in more detail
in \cref{sec:template}.

\section{Class options}
\label{sec:class-options}

Class options can be included in the bracketed argument of the
command, separated by commas. The possible class options are:
\begin{itemize}
\item \code{review} --- Recommended for submitting your manuscript to
  a SIAM journal.  Adds line numbers as well as the statement ``This
  manuscript is for review purposes only'' to the bottom of each page.
\item \code{final} --- Turns off the black boxes that help authors
  identify lines that are too long. The final published version will
  have this option on.
\item \code{supplement} --- Specifies that the file is a supplement
  and not the main document, causing changes in the appearance of the
  title and numbering; see \cref{sec:supplement} for details.
\end{itemize}

\begin{example}[label={ex:outline},listing only,%
  listing options={style=siamlatex,{morekeywords=[1]{maketitle},
      morekeywords=[2]{siamonline220329}},}]%
  {Document outline}
\documentclass{siamonline220329}
\begin{document}
\maketitle
\end{document}
\end{example}

\section{Front matter}
\label{sec:front}
The title and author parts are formatted using the standard
\code{\title}, \code{\author}, and \code{\maketitle} commands as
described in Lamport \cite{La86}. The title and author should be
declared in the preamble. The title and author names are automatically
converted to uppercase in the document.
If there is more than one author, each additional author should be preceded by
the \code{\and} command.
The addresses and support acknowledgments are added via
\code{\thanks}. Each author's thanks should specify their address.
The support acknowledgment should be put in the title thanks,
unless specific support needs to be specified for individual authors,
in which case it should follow the author address.
The header for this file was produced by the code in \cref{ex:header}, including an
example of a shared footnote. Each thanks produces a footnote, so the
footnote of the second author is \#3.
The command
\code{\headers{title}{authors}} command, with
the title (possibly shortened to fit) and the authors' names, creates
the page headers, automatically converted to uppercase.

~

\examplefile[label={ex:header},listing only,%
listing options={style=siamlatex,%
deletetexcs={and,thanks,title,author},%
{moretexcs=[2]{and,thanks,title,author,maketitle,headers,email}}}
]{Title and authors in preamble}{tmp_\jobname_header.tex}

Following the author and title is the abstract, key words listing, and MSC codes,
designated using the \code{abstract}, \code{keywords}, and \code{MSCcodes}
environments.
Authors are responsible for providing MSC codes which can be found
on the AMS web site \cite{AMSMSC2010}. The abstract, keywords, and MSC codes
for this document are specified in \cref{ex:abstract}.

\examplefile[label={ex:abstract},%
before upper={\preamble{\bs newcommand\{\bs BibTeX\}\{\{\bs scshape Bib\}\bs TeX\bs xspace\}}},
listing only,%
listing options={style=siamlatex,%
{morekeywords=[2]{abstract,keywords,MSCcodes}}}
]{Abstract, keywords, and MSC codes}{tmp_\jobname_abstract.tex}

A more complete example, including a PDF supplement, that uses the
included files \texttt{ex\_article.tex}, \texttt{ex\_supplement.tex},
and \texttt{ex\_shared.tex} is discussed in
\cref{sec:template}. The example files can be used as a starting point
for producing a document.

\section{Cross references and hyperlinks}
\label{sec:cr+hyp}

SIAM now supports cross references and hyperlinks via the
\texttt{cleveref} and \texttt{hyperef} packages, which are loaded by
the class file.

\subsection{Cleveref}
\label{sec:cleveref}

SIAM strongly recommends using the commands provided by
the \texttt{cleveref} package for cross referencing.
The package is automatically loaded and already customized
to adhere to SIAM's style guidelines.
To create a cross reference, use the command \code{\cref} (inside
sentence) or \code{\Cref} (beginning of a sentence) in place of the
object name and \code{\ref}.
The \texttt{cleveref}
package enhances \LaTeX's cross-referencing features, allowing
the format of cross references to be determined automatically
according to the ``type" of cross reference (equation, section, etc.)
and the context in which the cross reference is used.
So, the package
\emph{automatically} inserts the object name as well as the
appropriate hyperlink; see \cref{ex:cref}.
It may require two \LaTeX\@ compilations for the references to show up
correctly.
Additional examples are shown in the sections below
for equations, tables, figures, sections, etc.

\begin{example}[label=ex:cref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref,ref}}}]{Advantage of using cleveref}
The normal way to get a cross reference with a hyperlink requires a
lot of typing: \hyperref[thm:mvt]{Theorem~\ref*{thm:mvt}}.
The \texttt{cleveref} package gets both the name and hyperlink
automatically using a single macro: \cref{thm:mvt}.
It also handles multiple references with the same macro, such as
\cref{thm:mvt,fig:pgfplots,fig:testfig}.
\end{example}

\subsection{Hyperef}
\label{sec:hyperef}

Hyperlinks are created with the \code{\href} and \code{\url} commands,
as shown in \cref{ex:href}.
SIAM has also defined the \code{\email} command, as shown in
\cref{ex:header}.

\begin{example}[label={ex:href},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{href,url}}}]{Creating hyperlinks}
The \href{https://www.siam.org}{SIAM homepage} has general information.
There are times when the author may want to specify
the location explicitly instead by using \url{https://www.siam.org}.
\end{example}

Note that homepage links via \code{\url} in the \code{\thanks} environment require
special formatting for the tilde (\string~) character. The formatting
is used in the template and shown in \cref{ex:shared}.

\section{Math and equations}
\label{sec:math}
Here we show some example equations, with numbering, and examples of
referencing the equations. SIAM now includes the package
\texttt{amsmath} by default, and we include some of its features as
well, although the reader should consult the package user manual for
further guidance \cite{amsmath,shortmath}.
 Several of the example are adapted
from Mittlebach and Goossen's guide to \LaTeX~\cite{MiGo04}.

\Cref{ex:textmath} is a straightforward example of inline mathematics
equations that does not use any special packages or features.

\begin{example}[label={ex:textmath},bicolor]{Inline math}
The following shows an example of math in text:
Let $S=[s_{ij}]$ ($1\leq i,j\leq n$) be a $(0,1,-1)$-matrix of order $n$.
\end{example}

In \cref{ex:bbm}, we show the recommended method for getting
blackboard fonts using the \texttt{amsfonts} package. This is not
loaded by default and must be included in the preamble.

\begin{example}[label={ex:bbm},bicolor,before upper={\preamble{\bs
      usepackage\{amsfonts\}}},%
listing options={style=siamlatex,%
{morekeywords=[2]{mathbb}}}]{Blackboard math}
Blackboard bold characters, such as $\mathbb{C}$ and $\mathbb{R}$,
should be created with the \texttt{amsfonts} package, although this
is not included by default.
\end{example}

\newpage
\Cref{ex:smallmatrix} shows the \code{smallmatrix} environment for an
inline matrix from the \texttt{amsmath} package, which is included by
default.

\begin{example}[label={ex:smallmatrix},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{smallmatrix}}}]{Inline matrix}
Matrices of no more than two rows appearing in text can be created
as shown in the next example:
$B = \bigl[ \begin{smallmatrix} B_{11} & B_{12} \\
  B_{21} & B_{22} \end{smallmatrix} \bigr]$.
\end{example}

Bigger matrices can be rendered with environments from
the \texttt{amsmath} package, such as \code{bmatrix} and
\code{pmatrix} used in
\cref{ex:matrices}.

\begin{example}[label={ex:matrices},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{bmatrix,pmatrix}}}]{Creating matrices}
Display matrices can be rendered using environments from \texttt{amsmath}:
\begin{equation}\label{eq:matrices}
  S=\begin{bmatrix}1&0\\0&0\end{bmatrix}
  \quad\text{and}\quad
  C=\begin{pmatrix}1&1&0\\1&1&0\\0&0&0\end{pmatrix}.
\end{equation}
\Cref{eq:matrices} shows some example matrices.
\end{example}

\Cref{ex:dmo} shows how to use the \code{\DeclareMathOperator} command
from the \texttt{amsopn} package to declare the \code{\Range} macro.
(This example also uses the \texttt{braket} package for the
\code{\set} macro, but this is not necessarily recommended by SIAM.)
\begin{example}[label={ex:dmo},%
before upper={\preamble{\bs usepackage\{braket,amsfonts,amsopn\}}\\
\noindent\preamble{\bs DeclareMathOperator\{\bs Range\}\{Range\}}},%
bicolor,%
listing options={style=siamlatex,%
{moretexcs=[2]{Range}}}
]{Declaring math operators}
An example of a math operator:
\begin{equation}\label{eq:range}
  \Range(A) = \set{ y \in \mathbb{R}^n | y = Ax }.
\end{equation}
\end{example}

\newpage
\Cref{ex:foo} shows how to use the \code{align} environment from
\texttt{amsmath} to easily align multiple equations.

\begin{example}[label={ex:foo},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{align}}}]{Aligned equations}
\Cref{eq:a,eq:b,eq:c} show three aligned equations.
\begin{align}
  f &= g, \label{eq:a} \\
  f' &= g', \quad\text{and} \label{eq:b} \\
  \mathcal{L}f &= \mathcal{L}g \label{eq:c}.
\end{align}
\end{example}

Another way to number a set of equations is the
\code{subequations} environment from \texttt{amsmath}, as shown in \cref{ex:aligned}.

\begin{example}[label={ex:aligned},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{subequations}}}]{Subequations}
We calculate the Fr\'{e}chet derivative of $F$ as follows:
\begin{subequations}
\begin{align}
  F'(U,V)(H,K)
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T} -
  P(H\Sigma V^{T} + U\Sigma K^{T})\rangle \label{eq:aa} \\
  &= \langle R(U,V),H\Sigma V^{T} + U\Sigma K^{T}\rangle
  \nonumber \\
  &= \langle R(U,V)V\Sigma^{T},H\rangle +
  \langle \Sigma^{T}U^{T}R(U,V),K^{T}\rangle. \label{eq:bb}
\end{align}
\end{subequations}
\Cref{eq:aa} is the first line, and \cref{eq:bb} is the last line.
\end{example}

For an equation split over multiple lines, \cref{ex:ml} shows the
usage of the \code{multline} environment provided by \texttt{amsmath}.
\begin{example}[label={ex:ml},bicolor,%
listing options={style=siamlatex,%
{morekeywords=[2]{multline}}}]{Equation split across lines}
We claim that the projection $g(U,V)$ is given by the pair of matrices:
\begin{multline} \label{eq:ml}
  g(U,V) = \biggl( \frac{R(U,V)V\Sigma^{T}U^{T}
    - U\Sigma V^{T}R(U,V)^{T}}{2}U,\\
  \frac{R(U,V)^{T}U\Sigma V^{T}-V \Sigma^{T}U^{T}R(U,V)}{2}V \biggr).
\end{multline}
\end{example}

\section{Theorem-like environments}
\label{sec:thm}
SIAM loads \texttt{ntheorem} package and uses it to define the following
theorem-like environments:
\code{theorem},
\code{lemma},
\code{corollary},
\code{definition}, and
\code{proposition}.
SIAM also defines a \code{proof} environment that automatically
inserts the symbol ``$\,\proofbox\,$'' at the end of any proof, even if it ends in an
equation environment. \emph{Note that the document may need to be
  compiled twice for the mark to appear.}
Some of the calculus examples were adapted from \cite{CalcI}.
\Cref{ex:theorem} shows usage of the \code{theorem} environment.
An optional argument can be used to name the theorem.
\Cref{ex:cor} illustrates a corollary, without a name, and the
proof environment.

\begin{example}[label=ex:theorem,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{theorem}}}]{Theorem}
\begin{theorem}[Mean Value Theorem]\label{thm:mvt}
  Suppose $f$ is a function that is continuous on the closed interval
  $[a,b]$.  and differentiable on the open interval $(a,b)$.
  Then there exists a number $c$ such that $a < c < b$ and
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a}.
  \end{displaymath}
  In other words, $f(b)-f(a) = f'(c)(b-a)$.
\end{theorem}
\end{example}

\begin{example}[label=ex:cor,bicolor,parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{corollary,proof}}}]%
{Corollary and proof}
\begin{corollary}
  Let $f(x)$ be continuous and differentiable everywhere. If $f(x)$
  has at least two roots, then $f'(x)$ must have at least one root.
\end{corollary}
\begin{proof}
  Let $a$ and $b$ be two distinct roots of $f$.
  By \cref{thm:mvt}, there exists a number $c$ such that
  \begin{displaymath}
    f'(c) = \frac{f(b)-f(a)}{b-a} = \frac{0-0}{b-a} = 0.
  \end{displaymath}
\end{proof}
\end{example}

\newpage
SIAM also defines commands to create your own theorem- and
remark-like environments:
\begin{itemize}
\item \code{newsiamthm} --- Small caps header, italized body.
\item \code{newsiamremark} --- Italics header, roman body.
\end{itemize}
Each command takes two arguments. The first is the environment name,
and the second is the name to show in the document. These commands should
be used instead of \code{\newtheorem}.
\Cref{ex:claim,ex:ref} shows how to use the commands above, including how to
specify the plural version for \texttt{cleveref} if it is unusual.

\begin{example}[label=ex:claim,bicolor,%
before upper={\preamble{\bs newsiamthm\{claim\}\{Claim\}}\\
\noindent\preamble{\bs newsiamremark\{hypothesis\}\{Hypothesis\}}\\
\noindent\preamble{\bs crefname\{hypothesis\}\{Hypothesis\}\{Hypotheses\}}},%
parbox=false,%
listing options={style=siamlatex,%
{morekeywords=[2]{claim,proof,hypothesis}}}]{New theorem-like environment}
\begin{claim}\label{cl:constant}
  If $f'(x) = 0$ for all $x \in (a,b)$ then $f(x)$ is constant on $(a,b)$.
\end{claim}
\begin{hypothesis}\label{hyp1}
The function $f$ is continuously differentiable.
\end{hypothesis}
\begin{hypothesis}\label{hyp2}
The random variable is normally distributed.
\end{hypothesis}
\end{example}

\begin{example}[label=ex:ref,bicolor,listing options={style=siamlatex,%
{morekeywords=[2]{cref}}}]{References}
We can reference multiple types of objects with a single reference:
\cref{cl:constant,thm:mvt,hyp1}.
\end{example}

\section{Tables}
\label{sec:tab}
Table captions should go above the tables.
\Cref{ex:simpletable} shows the code to generate a
\cref{tab:simpletable}.
A more complicated example is shown in
\cref{ex:table}, which generates \cref{tab:KoMa14}. This
example uses subfloats via the \texttt{subfig} package, as well as
special column options from the \texttt{array} package.

\begin{tcbverbatimwrite}{tmp_\jobname_simpletable.tex}
\begin{table}[tbhp]
\footnotesize
  \caption{Example table}\label{tab:simpletable}
\begin{center}
  \begin{tabular}{|c|c|c|} \hline
   Species & \bf Mean & \bf Std.~Dev. \\ \hline
    1 & 3.4 & 1.2 \\
    2 & 5.4 & 0.6 \\ \hline
  \end{tabular}
\end{center}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:simpletable},%
listing only, listing options={style=siamlatex}]%
{Example table.}{tmp_\jobname_simpletable.tex}

\input{tmp_\jobname_simpletable.tex}

\begin{tcbverbatimwrite}{tmp_\jobname_table.tex}
\newcolumntype{R}{>{$}r<{$}}  %
\newcolumntype{V}[1]{>{[\;}*{#1}{R@{\;\;}}R<{\;]}} %
\begin{table}[tbhp]
\footnotesize
 \captionsetup{position=top}
  \caption{Example table adapted from Kolda and Mayo \rm{\cite{KoMa14}}.}\label{tab:KoMa14}
\begin{center}
  \subfloat[$\beta=1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} &
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
 718 & 11.3476 &  0.5544 &  0.3155 &  1.2018 &  0.0977 &   45 & 0.17 & 0.06 \\ \hline
 134 &  3.7394 &  0.2642 & -1.1056 &  0.2657 & -0.3160 &   31 & 0.12 & 0.05 \\ \hline
   4 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.21 & 0.10 \\ \hline
 \end{tabular}}

  \subfloat[$\beta=-1$]{
    \begin{tabular}{|r|R|V{3}|c|r@{\,$\pm$\,}l|} \hline
occ. & \multicolumn{1}{c|}{$\lambda$} & \multicolumn{4}{c|}{$\mathbf{x}$} &
fevals &  \multicolumn{2}{c|}{time (sec.)}\\ \hline
  72 & -1.1507 &  0.2291 &  0.6444 &  0.3540 & -0.8990 &   34 & 0.14 & 0.06 \\ \hline
 624 & -6.3985 &  0.1003 &  0.1840 &  0.5305 &  1.2438 &   48 & 0.19 & 0.08 \\ \hline
   2 & \multicolumn{6}{c|}{\emph{--- Failed to converge ---}} & 0.23 & 0.02 \\ \hline
 \end{tabular}}
\end{center}
\end{table}
\end{tcbverbatimwrite}

\examplefile[label={ex:table},%
before upper={\preamble[\scriptsize]{\bs usepackage\{array\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\bs usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example table with subtables.}{tmp_\jobname_table.tex}

\input{tmp_\jobname_table.tex}

\newpage
\section{Figures}
\label{sec:fig}

It is recommended that all figures be generated in high resolution.
In the past, SIAM has required encapsulated
postscript (EPS) format for final production. This is still an
acceptable format, but SIAM also now allows high-resolution
PDF, JPEG, and PNG figures.
If working with EPS images and using \texttt{pdflatex}, we
recommend the package \texttt{epstopdf} to automatically convert EPS
images to PDF for inclusion in PDF documents created by
\texttt{pdflatex}. Please note that \texttt{epstopdf} requires \texttt{texlive-font-utils} which may not be part of a standard \LaTeX\ installation.
\Cref{ex:fig} shows the code to generate \cref{fig:testfig}.  This
example uses the \texttt{graphicx} package for the
\code{\includegraphics} command.

\begin{tcbverbatimwrite}{tmp_\jobname_fig.tex}
\begin{figure}[tbhp]
  \centering
  \subfloat[$\epsilon_{\max}=5$]{\label{fig:a}\includegraphics{lexample_fig1}}
  \subfloat[$\epsilon_{\max}=0.5$]{\label{fig:b}\includegraphics{lexample_fig2}}
  \caption{Example figure using external image files.}
  \label{fig:testfig}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:fig},%
before upper={\preamble[\scriptsize]{\bs usepackage\{graphicx,epstopdf\}}\\[-0.4em]
\noindent\preamble[\scriptsize]{\bs usepackage[caption=false]\{subfig\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example figure with subfigures and external files}{tmp_\jobname_fig.tex}

\input{tmp_\jobname_fig.tex}

Another option for figures is a graphics-generator that is platform- and
format-independent.
PGF is a TeX macro package for generating  such graphics and works together with the most important TeX backend drivers, including pdftex and dvips.
The user-friedly syntax layer called TikZ.
Here we show an example using \texttt{PGFPLOTS}, useful for drawing high-quality
plots directly in \LaTeX.
\Cref{ex:data} and \cref{ex:pgfplots} shows the data and code, respectively, to generate \cref{fig:pgfplots}, adapted from \cite{pgfplots}.

\examplefile[label={ex:data},listing only,
listing options={style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example data file (data.dat)}{data.dat}

\begin{tcbverbatimwrite}{tmp_\jobname_tikz.tex}
\begin{figure}[tbhp]
  \centering
  \begin{tikzpicture}
    \begin{loglogaxis}[height=2.75in, grid=major,
      xlabel={Degrees of Freedom}, ylabel={$L_2$ Error},
      legend entries={$d=2$,$d=3$}]
      \addplot table [x=d2_dof,y=d2_l2_err] {data.dat};
      \addplot table [x=d3_dof,y=d3_l2_err] {data.dat};
    \end{loglogaxis}
  \end{tikzpicture}
  \caption{Example \texttt{PGFPLOTS} figure.}
  \label{fig:pgfplots}
\end{figure}
\end{tcbverbatimwrite}

\examplefile[label={ex:pgfplots},%
before upper={\preamble[\scriptsize]{\bs usepackage\{pgfplots\}}},%
listing only, listing options={%
  style=siamlatex}]%
{Example TikZ/PGF for platform-independent graphics.}{tmp_\jobname_tikz.tex}

\input{tmp_\jobname_tikz.tex}

\newpage
\section{Algorithms}
\label{sec:algs}

SIAM automatically includes the \texttt{algorithm} package in the
class definition. This provides the float environment.
Users have the choice of \texttt{algpseudocode},
\texttt{algorithmic}, and other packages for actually formatting the algorithm.
For example, \cref{alg:buildtree} is produced by the code in \cref{ex:alg}.
In order to reference lines within the algorithm, we need to tell the
\texttt{cleveref} package how to do the referencing, which is the
second line of \cref{ex:alg}. Then we can use the code
\code{\cref{line3}} to produce \cref{line3}.

\begin{tcbverbatimwrite}{tmp_\jobname_alg.tex}
\begin{algorithm}
\caption{Build tree}
\label{alg:buildtree}
\begin{algorithmic}[1]
\STATE{Define $P:=T:=\{ \{1\},\ldots,\{d\}$\}}
\WHILE{$\#P > 1$}
\STATE\label{line3}{Choose $C^\prime\in\mathcal{C}_p(P)$ with $C^\prime := \operatorname{argmin}_{C\in\mathcal{C}_p(P)} \varrho(C)$}
\STATE{Find an optimal partition tree $T_{C^\prime}$ }
\STATE{Update $P := (P{\setminus} C^\prime) \cup \{ \bigcup_{t\in C^\prime} t \}$}
\STATE{Update $T := T \cup \{ \bigcup_{t\in\tau} t : \tau\in T_{C^\prime}{\setminus} \mathcal{L}(T_{C^\prime})\}$}
\ENDWHILE
\RETURN $T$
\end{algorithmic}
\end{algorithm}
\end{tcbverbatimwrite}

\examplefile[float=htpb,label={ex:alg},%
before upper={\preamble[\scriptsize]{\bs usepackage\{algorithmic\}}\\[-0.4em]
\preamble[\scriptsize]{\bs Crefname\{ALC@unique\}\{Line\}\{Lines\}}},%
listing only, listing options={%
  style=siamlatex,basicstyle=\ttfamily\scriptsize}]%
{Example algorithm}{tmp_\jobname_alg.tex}

\input{tmp_\jobname_alg.tex}

\section{Sections}
\label{sec:sec}
Sections are denoted using standard \LaTeX\ section commands, i.e.,
\code{\section}, \code{\subsection}, etc.
If you wish to end the section title with something other that a period (the default), you have to add the command \code{\nopunct} at the end of the title.

Appendices are created with the normal sectioning commands, following
the command \code{
\section{Introduction}
\label{sec:introduction}

Kernel methods \cite{SS02} are a class of machine learning algorithms that are often used for prediction and clustering tasks.
Kernel methods require manipulating a positive semidefinite (psd) \emph{kernel matrix}  $\mat{A} \in \mathbb{R}^{N \times N}$.  The
entries $a_{ij} = K(\vec{x}^{(i)},\vec{x}^{(j)})$ of the kernel matrix
are determined by a positive definite \emph{kernel function} $K(\cdot,\cdot)$ that measures the similarity between pairs of data points~\cite[Def.~2.5]{SS02}.
The dimension $N$ of the kernel matrix is equal to the number of data points, which
can easily grow into the millions or higher.
The size of modern data sets poses a challenge for kernel methods: even generating all $N^2$ entries of $\mat{A}$ may be infeasible when $N \geq 10^5$.

To speed up kernel computations, a standard practice is to replace the kernel matrix $\mat{A}$ with a low-rank approximation.
Remarkably, in the common setting in which the $N \times N$ kernel matrix has rapidly decaying eigenvalues, it is possible to generate an accurate approximation using only $\mathcal{O}(N)$ entry evaluations and $\mathcal{O}(N)$ arithmetic operations.
Nevertheless, it has been formidable to design a practical and provably justified algorithm that achieves this goal.

The \emph{column Nystr\"om approximation}
\cite[\S19.2]{MT20} is the most widely used technique for constructing a low-rank psd approximation $\mat{\hat{A}}$ of a kernel matrix $\mat{A}$.
Given a set $\set{S} = \{s_1, \dots, s_k\} \subseteq \{1,\ldots,N\}$ of $k$ column indices, the associated Nystr{\"o}m approximation is the rank-$k$ psd matrix
\begin{equation} \label{eq:column_nystrom}
  \mat{\hat{A}} = \mat{A}(:,\set{S}) \mat{A}(\set{S},\set{S})^\dagger \mat{A}(\set{S},:) \in \mathbb{R}^{N \times N}.
\end{equation}
In this expression,  $\mat{A}(:,\set{S})$ is the submatrix with the selected columns, $\mat{A}(\set{S},:)$ is the submatrix with the selected rows,
and $\mat{A}(\set{S},\set{S})^\dagger$ is the Moore--Penrose pseudoinverse of the submatrix with the selected rows and columns.
It can be shown~\cite[Eq.~(1.2)]{lm00} that the matrix~\cref{eq:column_nystrom} is the best psd approximation of $\mat{A}$ such that $\operatorname{range}(\mat{\hat{A}}) \subseteq \operatorname{range}(\mat{A}(:, \set{S}))$ and the residual $\mat{A} - \mat{\hat{A}}$ is also psd.

Let us emphasize that the column Nystr{\"o}m approximation~\cref{eq:column_nystrom} only requires evaluating $kN$ entries of the matrix $\mat{A}$.  As a consequence, this approximation is particularly well-suited for matrices that must be accessed by evaluating entry-wise formulas, such as kernel matrices~\cite{SS02} and Hamiltonian matrices from quantum chemistry~\cite[Ch.~2]{szabo2012modern}.

From a mathematical point of view, algorithms for column Nystr{\"o}m decomposition~\cref{eq:column_nystrom} differ primarily in how they select the set $\set{S}$ of $k$ columns.
The challenge is to identify columns that capture as much action of the matrix as possible.
This paper studies a natural approach, which we call the \textit{randomly pivoted Cholesky} (\mbox{\textsc{RPCholesky}}\xspace) algorithm
(\cref{sec:partial-cholesky}).
Similar approaches have been used in the past to approximate a rectangular matrix by a subset of its rows \cite{DV06}
or to sample from a determinantal point process \cite{Pou20}, and a remark of \cite{MW17} mentions \mbox{\textsc{RPCholesky}}\xspace as a possibility for psd low-rank approximation.
However, \mbox{\textsc{RPCholesky}}\xspace has not been investigated as a serious strategy for psd matrix approximation before now.

Our work provides the first in-depth exploration of the empirical and theoretical behavior of \mbox{\textsc{RPCholesky}}\xspace.
Across a range of test problems, we find that \mbox{\textsc{RPCholesky}}\xspace matches or improves over alternative methods for selecting columns in the Nystr{\"o}m approximation.
To explain this surprising observation, we prove that \mbox{\textsc{RPCholesky}}\xspace achieves a rigorous error bound that is nearly optimal within the class of Nystr{\"o}m algorithms.
Altogether, these results make a strong case that \mbox{\textsc{RPCholesky}}\xspace is the method of choice for low-rank approximation of a psd matrix under a constraint on the number of entry evaluations.

\subsection{Randomly pivoted Cholesky}
\label{sec:partial-cholesky}

We can best understand \mbox{\textsc{RPCholesky}}\xspace as a member of a family of algorithms based on \textit{partial Cholesky decomposition}.
To that end, observe that the column Nystr\"om approximation~\cref{eq:column_nystrom} can be constructed incrementally as follows.
Beginning with $\mat{\hat{A}}^{(0)} \coloneqq \mat{0}$ and $\mat{A}^{(0)} \coloneqq \mat{A}$, the partial Cholesky decomposition
performs the following operations at each step $i = 1, 2, \ldots, k$:
\begin{enumerate}
\item Select a pivot index $s_i \in \{1, \dots, N\}$.
\item \label{item:cholesky_approximation}
Update the Nystr\"om approximation:
\begin{equation*}
\mat{\hat{A}}^{(i)} = \mat{\hat{A}}^{(i-1)} + \frac{\mat{A}^{(i-1)}(:,s_i)\mat{A}^{(i-1)}(s_i,:)}{a^{(i-1)}_{s_is_i}}.
\end{equation*}
\item \label{item:cholesky_residual}
Update the residual (Schur complement):
\begin{equation*}
\mat{A}^{(i)} = \mat{A}^{(i-1)} - \frac{\mat{A}^{(i-1)}(:,s_i)\mat{A}^{(i-1)}(s_i,:)}{a^{(i-1)}_{s_is_i}}.
\end{equation*}
\end{enumerate}
By induction, each approximation $\mat{\hat{A}}^{(i)}$ and each residual $\mat{A}^{(i)}$ remain psd.
After $k$ steps, the partial Cholesky decomposition $\mat{\hat{A}}^{(k)}$ coincides with the Nystr\"om approximation \cref{eq:column_nystrom} defined by the index set $\set{S} = \{ s_1,\ldots, s_k \}$.
With a careful implementation (\cref{alg:rpcholesky}), steps \ref{item:cholesky_approximation} and \ref{item:cholesky_residual} involve just $kN$ entry evaluations and $\mathcal{O}(k^2 N)$ arithmetic operations.

There are several standard methods for choosing the set $\set{S}$ of pivots.
A classical strategy is the \emph{greedy method} \cite{FS01}, also known as diagonal or complete pivoting \cite{Hig02}.   At each step $i$, the pivot $s_i$ is the index of the largest diagonal entry of the current residual matrix, with ties broken arbitrarily:
\begin{equation} \label{eq:greedy}
  s_i \in \argmax\big\{ a^{(i-1)}_{jj} : 1\le j \le N \big\}.
\end{equation}
Another common strategy is to select each pivot $s_i$ \emph{uniformly at random} from $\{1, \dots, N\}$, either with or without replacement \cite{KMT09,WS00}.
However, both the greedy method and the uniform sampling method can fail to find salient columns.
The greedy method cannot explore the small diagonal entries of the residual matrix, while the uniform method cannot take advantage of the largest diagonal entries (\cref{sec:greedy,sec:uniform}).

To circumvent this problem, we consider \mbox{\textsc{RPCholesky}}\xspace as an alternative strategy for selecting the pivots.
At each step $i$, this algorithm draws a random pivot $s_i$ according to the probability distribution
\begin{equation} \label{eq:rpc_sampling}
\prob \left\{ s_i = j \right\} = \frac{a^{(i-1)}_{jj}}{\tr \mat{A}^{(i-1)}}
\quad
\text{for each $j = 1, \dots, N$.}
\end{equation}
In other words, each pivot is sampled in proportion to the diagonal entries of the current residual matrix.
\mbox{\textsc{RPCholesky}}\xspace thus strikes a balance between exploring small diagonal entries and exploiting large ones.
\mbox{\textsc{RPCholesky}}\xspace was mentioned in \cite{MW17} as a possible strategy for psd low-rank approximation,
but we are not aware of any earnest investigation of its behavior.

\Cref{alg:rpcholesky_naive} describes a naive implementation of \mbox{\textsc{RPCholesky}}\xspace, while \cref{sec:speed_up} presents a more efficient implementation.

\begin{algorithm}[h!]
  \caption{\mbox{\textsc{RPCholesky}}\xspace: Naive implementation} \label{alg:rpcholesky_naive}
  \textbf{Input:} Psd matrix $\mat{A} \in \mathbb{R}^{N\times N}$, iteration count $k$ or tolerance $\eta$

  \textbf{Output:} Rank-$k$ psd approximation $\mat{\hat{A}}^{(k)}$
  \begin{algorithmic}
    \State Initialize $\mat{\hat{A}}^{(0)} \leftarrow \mat{0}_{N\times N}$ and $\mat{A}_0 \leftarrow \mat{A}$
    \For{$i$ = $1$ to $k$}
    \Comment{Alternately, run until $\tr \mat{A}^{(i)} \le \eta \tr \mat{A}$}
    \State Sample $s_i \sim \diag \mat{A}^{(i-1)} / \tr \mat{A}^{(i-1)}$ \Comment{Sample new pivot}
    \State $\mat{\hat{A}}^{(i)} \leftarrow \mat{\hat{A}}^{(i-1)} + \mat{A}^{(i-1)}(:,s_i)\mat{A}^{(i-1)}(s_i,:) / a^{(i-1)}_{s_is_i}$ \Comment{Update low-rank approximation}
    \State $\mat{A}^{(i)} \leftarrow \mat{A}^{(i-1)} - \mat{A}^{(i-1)}(:,s_i)\mat{A}^{(i-1)}(s_i,:) / a^{(i-1)}_{s_is_i}$ \Comment{Update Schur complement}
    \EndFor
  \end{algorithmic}
\end{algorithm}

\subsection{Numerical results}
\label{sec:compare}

Given that \mbox{\textsc{RPCholesky}}\xspace is a natural algorithm, we were perplexed by the lack of attention that it has received in the literature on randomized matrix computations.  Therefore, we undertook a suite of experiments to compare \mbox{\textsc{RPCholesky}}\xspace with other column Nystr{\"o}m approximations on a range of test problems.
Surprisingly, we discovered that \mbox{\textsc{RPCholesky}}\xspace reliably produces some of the best matrix approximations using the fewest number of entry evaluations.  This work is documented in \cref{sec:rpc}.

For now, let us offer a visually appealing example that is representative of our broader findings.
We aim to approximate the kernel matrices associated with two data sets:

\begin{enumerate}
    \item \textbf{Smile:} A Gaussian kernel matrix constructed from $10^4$ data points depicting a smile in $\mathbb{R}^2$.
    The eyes are constructed from just $10^2$ points, making them easy to miss for certain sampling methods.
    \item \textbf{Outliers:} A Gaussian kernel matrix constructed from a point cloud of $10^4$ points in $\mathbb{R}^{20}$.  Of these points, $50$ are outliers that are positioned very far away from the others.
\end{enumerate}

\begin{figure}[h!]
    \begin{subfigure}{.49\textwidth}
    \centering
    \includegraphics[width = \textwidth]{smile_trace}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.24\textwidth}
    \centering
    \includegraphics[width = \textwidth]{smile_uniform}
    Uniform
    \end{subfigure}
    \hfill
    \begin{subfigure}{.24\textwidth}
    \centering
    \includegraphics[width = \textwidth]{smile_rpc}
    \mbox{\textsc{RPCholesky}}\xspace
    \end{subfigure}
    \begin{subfigure}{.49\textwidth}
    \includegraphics[width = \textwidth]{outliers_trace}
    \end{subfigure}
    \hfill
    \begin{subfigure}{.24\textwidth}
    \centering
    \includegraphics[width = \textwidth]{outliers_greedy}
    Greedy
    \end{subfigure}
    \hfill
    \begin{subfigure}{.24\textwidth}
    \centering
    \includegraphics[width = \textwidth]{outliers_rpc}
    \mbox{\textsc{RPCholesky}}\xspace
    \end{subfigure}
    \caption{\textit{Left:} Mean relative trace-norm error $\tr (\mat{A} - \mat{\hat{A}}^{(k)}) \slash \tr \mat{A}$ and standard deviation error bars for different Nystr\"om methods.
    \textit{Right:} Selected pivots (colored stars) and data points (grey circles) for the \textbf{Smile} and \textbf{Outliers} data, with \textbf{Outliers} data projected onto a random plane.
    The greedy method chooses many pivots from the 50 outliers, while \mbox{\textsc{RPCholesky}}\xspace mostly avoids the outliers.
    \label{fig:comparison}}
\end{figure}

\Cref{fig:comparison} charts the accuracy of the column Nystr{\"o}m approximations obtained with different pivoting strategies.
The left-hand panels show the average relative trace-norm error
over $100$ independent trials.
The right-hand panels show the pivots selected by several Nystr\"om methods with $k = 40$.
For additional computational details, see~\cref{app:numerical_details}.

These tests bring to light the failure modes of uniform sampling and the greedy method.
Uniform sampling fails to select the pivots representing small, important features of the data,
such as the eyes in the \textbf{Smile} example.
The greedy method heavily emphasizes outliers,
leading to poor approximation accuracy in the \textbf{Outliers} example.
By contrast, \mbox{\textsc{RPCholesky}}\xspace avoids these failure modes and achieves high accuracy for both test problems.

For completeness, we also compare \mbox{\textsc{RPCholesky}}\xspace with two other column selection schemes for Nystr\"om approximation: \emph{determinantal point process} (DPP) sampling \cite{DM21a} and \emph{ridge leverage score} (RLS) sampling \cite{AM15,CLV17,MM17,RCCR18}.
These methods offer strong theoretical guarantees, yet existing implementations are complicated and require far more entry evaluations than \mbox{\textsc{RPCholesky}}\xspace to achieve the same approximation quality.
See \cref{sec:dpp,sec:rls} for full details about these approaches.

\Cref{fig:entries} displays the computational cost of different column Nystr\"om methods, as measured by the number of entry evaluations.
This is the total number of entry evaluations required to first form the sampling distribution and then generate the sampled columns.
Uniform sampling is the cheapest Nystr\"om method, requiring just $kN$ entry evaluations.
The greedy method and \mbox{\textsc{RPCholesky}}\xspace follow closely behind, requiring $(k+1)N$ entry evaluations.
RLS sampling requires roughly $3kN$ entry evaluations and DPP sampling (using the \texttt{vfx} sampler \cite{GPBV19}) requires between $80kN$ and $1000kN$ entry evaluations,
making these latter two methods comparatively very expensive.

\begin{figure}[h!]
    \centering
    \includegraphics[width = .49\textwidth]{entries_smile.png}
    \hfill
    \includegraphics[width = .49\textwidth]{entries_outliers.png}

    \caption{Entry evaluations
    for different Nystr\"om methods
    applied to the \textbf{Smile} (\textit{left}) and \textbf{Outliers} (\textit{right}) matrices.
    The \texttt{vfx} DPP sampler fails for $k\ge 60$ on the \textbf{Smile} matrix.}
    \label{fig:entries}
\end{figure}

In addition to being expensive, DPP sampling and RLS sampling exhibit other limitations.
For the $\textbf{Smile}$ kernel matrix with
$k \geq 60$,
DPP sampling using the \texttt{vfx} sampler
\cite{GPBV19} returns an error message due to the matrix $\mat{A}$ being close to low-rank and fails to produce any output.
It is necessary to use a slower
\texttt{GS} sampler \cite{GPBV19} based on a complete eigendecomposition of the matrix $\mat{A}$.
For the $\textbf{Smile}$ kernel matrix,
RLS sampling provides a low-quality approximation,
often failing to include data points from one or both eyes.
The poor RLS sampling results are consistent with the limited theoretical error bounds for this method
(see \cref{sec:rls}),
which necessitate using hundreds of columns to guarantee an accurate low-rank approximation for any psd matrix.

\subsection{Theoretical results}

After observing the excellent performance of the \mbox{\textsc{RPCholesky}}\xspace algorithm empirically,
we sought a more rigorous explanation.
This paper proves that
\mbox{\textsc{RPCholesky}}\xspace attains error bounds that are nearly optimal within the class of
column Nystr{\"o}m approximations.  \Cref{sec:theoretical-analysis} studies the theoretical performance
of \mbox{\textsc{RPCholesky}}\xspace, and \cref{sec:related-work} outlines results for competing algorithms.
We summarize the comparison here.

Let $\mat{A}$ be a psd matrix, and let $\mat{\hat{A}}$ be a column Nystr{\"o}m approximation~\cref{eq:column_nystrom} of $\mat{A}$.
We consider the approximation error as $\tr(\mat{A} - \mat{\hat{A}})$, which is variously known as the trace norm, nuclear norm, or Schatten 1-norm of the residual $\mat{A} - \mat{\hat{A}}$.
We compare against the error attained by a best rank-$r$ psd approximation $\lowrank{\mat{A}}_r$,
obtained from an $r$-truncated eigenvalue decomposition.

A randomized column Nystr{\"o}m approximation $\mat{\hat{A}}$ is called an \textit{$(r, \varepsilon)$-approximation} of $\mat{A}$ when
\begin{equation} \label{eq:1+eps}
\expect \tr( \mat{A} - \mat{\hat{A}} ) \leq (1 + \varepsilon) \cdot \tr( \mat{A} - \lowrank{\mat{A}}_r )
\end{equation}
for parameters $r \in \mathbb{N}$ and $\varepsilon > 0$.
The expectation averages over the random choices in the algorithm (e.g., the columns that are sampled).
The efficiency of a particular sampling scheme depends on how large we must take the number $k$ of columns.
Here is the key theoretical question:

\begin{question*}
How many columns $k$ are sufficient to guarantee that a randomized column Nystr{\"o}m method attains
an $(r, \varepsilon)$-approximation~\cref{eq:1+eps} for every $N \times N$ psd input matrix?
\end{question*}

To achieve an $(r, \varepsilon)$-approximation for the worst-case matrix,
a column Nystr{\"o}m approximation must use at least $k \geq r \slash \varepsilon$
columns.  See \cref{thm:near-optimal}\ref{item:dpp_c} for an example.
We will establish the remarkable fact that \mbox{\textsc{RPCholesky}}\xspace achieves an $(r,\varepsilon)$-approximation
for every psd matrix after accessing only a modest multiple of $r \slash \varepsilon$ columns.

\Cref{tab:comparison} presents the best available upper bounds on the number $k$ of columns for \mbox{\textsc{RPCholesky}}\xspace and other column Nystr{\"o}m approximation methods
to achieve~\cref{eq:1+eps}.
These bounds are expressed terms of $r$, $\varepsilon$, $N$, and the relative approximation error
\begin{equation*} \label{eq:relative_error}
    \eta \coloneqq \frac{\tr(\mat{A} - \lowrank{\mat{A}}_r)}{\tr \mat{A}}
    \in [0, 1].
\end{equation*}
\begin{table}[h!]
  \centering
  \renewcommand{\arraystretch}{1.1}
  \begin{tabular}{ccc}
    \toprule
    Method & Number $k$ of columns & Reference\\ \midrule
    Greedy method & $(1 - (1 + \varepsilon)\eta) N$ & \cref{thm:greedy} \\
    Uniform sampling* & $\frac{r-1}{\varepsilon \eta} + \frac{1}{\varepsilon}$ & \Cref{thm:uniform} \\
    DPP sampling & $\frac{r}{\varepsilon} + r - 1$ &

    \Cref{thm:near-optimal}
    \\
    RLS sampling** & $136 \frac{r}{\varepsilon} \log\mleft(\frac{3r}{\varepsilon \eta}\mright)$ & \Cref{cor:musco} \\
    \midrule
    \mbox{\textsc{RPCholesky}}\xspace & $\frac{r}{\varepsilon} + r \log_+ \mleft(\frac{1}{\varepsilon \eta}\mright)$
    & \cref{thm:main_bound}  \\
    \mbox{\textsc{RPCholesky}}\xspace & $\frac{r}{\varepsilon} + r + r \log_+ \mleft(\frac{2^r}{\varepsilon}\mright)$
    & \cref{thm:main_bound}  \\\bottomrule
  \end{tabular}
  \caption{Upper bounds on the number of columns for column Nystr{\"o}m approximation methods to produce an $(r,\varepsilon)$-approximation \cref{eq:1+eps}.
  *The result for uniform sampling assumes the diagonal entries of $\mat{A}$ are set to one.
  **The result for RLS sampling assumes that $\varepsilon < \frac{1}{4}$.}
  \label{tab:comparison}
\end{table}

DPP sampling has the strongest error bounds of all the methods in \cref{tab:comparison}.  To achieve an $(r, \varepsilon)$-approximation with DPP sampling,
it is sufficient to take the number $k$ of columns as
\begin{equation} \label{eq:cannot_improve}
    k \geq \frac{r}{\varepsilon} + r - 1.
\end{equation}
\mbox{\textsc{RPCholesky}}\xspace has the second strongest error bounds, requiring no more columns than
\begin{equation}
    \label{eq:can_improve}
    k \geq \frac{r}{\varepsilon} + r \log_+\Bigl(\frac{1}{\varepsilon \eta}\Bigr),
\end{equation}
where $\log_+(x) \coloneqq \max\{\log x, 0\}$ for $x > 0$.
The main difference between the DPP result \cref{eq:cannot_improve} and the \mbox{\textsc{RPCholesky}}\xspace result \cref{eq:can_improve}
is the multiplicative factor $\log(1/\eta)$ present in the latter.
However, because the relative error $\eta$ appears inside the logarithm, this factor has only a modest impact on the computational scaling.
Indeed, $\log (1/\eta) < 37$ when $\eta > 10^{-16}$, the machine precision in double-precision arithmetic.
See \cref{sec:theoretical-analysis} for additional \mbox{\textsc{RPCholesky}}\xspace error bounds,
which consider the limit as $\eta \rightarrow 0$.

The error bounds for \mbox{\textsc{RPCholesky}}\xspace are \emph{significantly stronger} than the bounds for the greedy method and uniform sampling.
For a worst-case matrix $\mat{A}$, the greedy method requires $\Theta(N)$ columns to approach the best rank-$r$ approximation error,
and the uniform sampling method requires $\Theta(r/\eta)$ columns (\cref{sec:greedy,sec:uniform}).
In contrast, \mbox{\textsc{RPCholesky}}\xspace requires a number of columns that is independent of the dimension $N$ and depends only logarithmically on the relative error $\eta$.
These results rigorously justify why \mbox{\textsc{RPCholesky}}\xspace does not exhibit the same failure modes as the greedy method and the uniform sampling method.

Our theoretical work confirms that \mbox{\textsc{RPCholesky}}\xspace achieves a given approximation error with roughly the same number of columns
as the best theoretical approach
for column Nystr{\"o}m approximation.
At the same time, \mbox{\textsc{RPCholesky}}\xspace uses fewer entry evaluations than other more complicated algorithms to attain this level of error.
Given the simplicity, efficiency, and reliability of \mbox{\textsc{RPCholesky}}\xspace, we believe that there is a strong
case for deploying this algorithm in applications.

\subsection{Plan for paper}

The rest of the paper is organized as follows.  \Cref{sec:rpc} describes a fast implementation of \mbox{\textsc{RPCholesky}}\xspace
and presents numerical demonstrations of \mbox{\textsc{RPCholesky}}\xspace in machine learning tasks.  \Cref{sec:theoretical-analysis} presents the theoretical analysis of \mbox{\textsc{RPCholesky}}\xspace,
while \cref{sec:related-work} discusses the performance of alternative methods.
\Cref{sec:conclusion} surveys future directions.

\subsection{Notation}

The symbol $\preceq$ denotes the \emph{psd order} on symmetric matrices: $\mat{H} \preceq \mat{A}$ if and only if $\mat{A} - \mat{H}$ is psd.
Submatrices are expressed using MATLAB notation, for example, $\mat{A}(:,i)$ represents the $i$th column of $\mat{A}$ and $\mat{A}(\set{S},:)$ denotes the submatrix of $\mat{A}$ with rows indexed by the set $\set{S}$.
\mbox{\textsc{RPCholesky}}\xspace and its analysis extend to \emph{complex-valued} psd matrices without significant modification, but for simplicity this work focuses on the real setting.
The transpose of $\mat{A}$ is denoted $\mat{A}^*$, which becomes the \textit{conjugate transpose} for complex matrices.

\section{Algorithms and experiments}
\label{sec:rpc}

This section gives efficient pseudocode for \mbox{\textsc{RPCholesky}}\xspace (\cref{sec:speed_up}) and presents a range of numerical examples (\cref{sec:krr,sec:spectral_clustering}).

\subsection{Speeding up \mbox{\textsc{RPCholesky}}\xspace}
\label{sec:speed_up}

To accelerate the naive implementation of \mbox{\textsc{RPCholesky}}\xspace in \cref{alg:rpcholesky_naive}, we can apply the following techniques:
\begin{itemize}
    \item The Nystr\"om approximation $\mat{\hat{A}}$ can be stored and updated in factored form $\mat{\hat{A}} = \mat{F}\mat{F}^*$ where $\mat{F} \in \mathbb{R}^{N\times k}$ is a thin matrix.
    \item Rather than forming and maintaining the residual matrix $\mat{A}^{(i)}$, it suffices to generate the single column $\mat{A}^{(i)}(:,s_i)$ at every iteration $i = 1, 2, \ldots, k$.
\end{itemize}
The resulting \mbox{\textsc{RPCholesky}}\xspace algorithm requires $(k+1)N$ entry evaluations and $2k^2N/3$ arithmetic operations; see \cref{alg:rpcholesky}.

\begin{algorithm}
  \caption{\mbox{\textsc{RPCholesky}}\xspace: Efficient implementation} \label{alg:rpcholesky}
  \textbf{Input:} Psd matrix $\mat{A} \in \mathbb{R}^{N\times N}$ and iteration count $k$ or tolerance $\eta$

  \textbf{Output:} Factor matrix $\mat{F} \in \mathbb{R}^{N\times k}$ defining Nystr\"om approximation $\mat{\hat{A}} = \mat{F}\mat{F}^*$ and selected pivots $\set{S}$
  \begin{algorithmic}
    \State Initialize $\mat{F} \leftarrow \mat{0}_{N\times k}$,
    $\set{S} \leftarrow \emptyset$, and $\vec{d} \leftarrow \diag \mat{A}$ \For{$i = 1$ to $k$} \Comment{Alternately, run until $\sum_{j=1}^N d_j \le \eta \tr \mat{A}$}
    \State Sample $s \sim \vec{d} / \sum_{j=1}^N d_j$
    \State $\set{S} \leftarrow \set{S} \cup \{s\}$
    \State $\vec{g} \leftarrow \mat{A}(:,s) - \mat{F}(:,1:i-1) \mat{F}(s,1:i-1)^*$
    \State $\mat{F}(:,i) \leftarrow \vec{g} / \sqrt{g_s}$
    \State $\vec{d} \leftarrow \vec{d} - |\mat{F}(:,i)|^2$ \Comment{$|\mat{F}(:,i)|^2$ is the vector with $j$th entry $|\mat{F}(j,i)|^2$}
    \EndFor
  \end{algorithmic}
\end{algorithm}

To provide an even greater speedup in some settings,
\mbox{\textsc{RPCholesky}}\xspace can be modified to sample multiple columns in a single iteration, as in \cref{alg:rpcholesky_blocking}.
In experiments, this \textit{blocked} version of \mbox{\textsc{RPCholesky}}\xspace produces similar quality approximations to standard \mbox{\textsc{RPCholesky}}\xspace, and it reduces the runtime when generating multiple columns of $\mat{A}$ is faster then generating columns one at a time.
For a discussion of these runtime considerations, see \cref{app:numerical_details}.

\begin{algorithm}[h!]
  \caption{\mbox{\textsc{RPCholesky}}\xspace: Blocked variant} \label{alg:rpcholesky_blocking}
  \textbf{Input:} Psd matrix $\mat{A} \in \mathbb{R}^{N\times N}$, iteration count $k$, and block size $B$

  \textbf{Output:} Factor matrix $\mat{F} \in \mathbb{R}^{N\times k}$ defining Nystr\"om approximation $\mat{\hat{A}} = \mat{F}\mat{F}^*$ and selected pivots $\set{S}$
  \begin{algorithmic}
    \State Initialize $\mat{F} \leftarrow \mat{0}_{N\times k}$,
    $\set{S} \leftarrow \emptyset$,
    $\vec{d} \leftarrow \diag \mat{A}$, and
    $i \leftarrow 0$
    \While{$i < k$}
    \State Sample $s_1,s_2,\ldots,s_{\min(B,k-i)} \stackrel{\mathrm{iid}}{\sim} \vec{d} / \sum_{j=1}^N d_j$
    \State $\set{U} \leftarrow \Call{Unique}{s_1,s_2,\ldots,s_{\min(B,k-i)}}$ \Comment{Identify unique elements}
    \State $\set{S} \leftarrow \set{S} \cup \set{U}$
    \State $\mat{G} \leftarrow \mat{A}(:,\set{U}) - \mat{F}(:,1:i) [\mat{F}(\set{U},1:i)]^*$
    \State $\mat{R} \leftarrow \Call{Chol}{\mat{G}(\set{U},:)}$ \Comment{$\mat{G}(\set{U},:) = \mat{R}^*\mat{R}$}
    \State $\mat{F}(:,i:i+|\set{U}|) \leftarrow \mat{G}\mat{R}^{-1}$
    \State $\vec{d} \leftarrow \vec{d} - \Call{SquaredRowNorms}{\mat{F}(:,i:i+|\set{U}|)}$
    \State $i \leftarrow i + |\set{U}|$
    \EndWhile
  \end{algorithmic}
\end{algorithm}

\subsection{Kernel ridge regression} \label{sec:krr}

One powerful application for \mbox{\textsc{RPCholesky}}\xspace is to accelerate \emph{kernel ridge regression} (KRR)~\cite[\S4.9.1]{SS02}.
KRR is a nonlinear extension of least-squares regression that approximates an unknown input--output map
using a positive definite kernel function $K : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$
and a set of input--output pairs $(\vec{x}^{(1)},y_1),\ldots,(\vec{x}^{(N)},y_N)$.
As output, KRR provides a prediction function of the form
\begin{equation}
    \label{eq:krr_ansatz}
    f(\cdot \, ; \, \vec{\beta}) \coloneqq \sum_{i=1}^N \beta_i K\mleft(\vec{x}^{(i)}, \cdot\mright),
\end{equation}
with the coefficient vector $\vec{\beta}$ chosen to minimize a regularized least-squares loss:
\begin{equation*} \label{eq:minimization}
    \min_{\vec{\beta} \in \mathbb{R}^N} \frac{1}{N}
    \sum_{j=1}^N \mleft\lvert
    f\mleft(\vec{x}^{(j)}\, ;\,\vec{\beta}\mright) - y_j
    \mright\rvert^2
    + \lambda\, \sum_{i,j = 1}^N \beta_i \beta_j K(\vec{x}^{(i)}, \vec{x}^{(j)}).
\end{equation*}
Explicitly,
the vector $\vec{\beta}$ is the solution to a linear system
\begin{equation} \label{eq:krr_linear_system}
    \vec{\beta} = (\mat{A} + \lambda N \, \mathbf{I})^{-1} \vec{y},
\end{equation}
where $\mat{A}$ is the $N\times N$ kernel matrix
and $\vec{y}$ is the vector of output values.

Directly computing the $\vec{\beta}$ coefficients by \cref{eq:krr_linear_system} would require solving a dense linear system at $\mathcal{O}(N^3)$ cost, but
\mbox{\textsc{RPCholesky}}\xspace provides a faster alternative.
First, we apply \mbox{\textsc{RPCholesky}}\xspace to the kernel matrix $\mat{A}$ to generate a random set of pivots $\set{S} = \{s_1, \ldots, s_k\}$.
Next, we identify a restricted prediction function
\begin{equation} \label{eq:krr_acc}
    \hat{f}\mleft(\cdot \, ; \, \vec{\hat{\beta}} \mright) = \sum_{i = 1}^k \hat{\beta}_i K\mleft(\vec{x}^{(s_i)}, \cdot\mright)
\end{equation}
that minimizes the loss
\begin{equation*} \label{eq:minimization2}
    \min_{\vec{\hat{\beta}}\in\mathbb{R}^k} \frac{1}{N}
    \sum_{j=1}^N \mleft\lvert  \hat{f}\mleft(\vec{x}^{(j)} \, ; \, \vec{\hat{\beta}} \mright) - y_j
    \mright\rvert^2
    + \lambda\, \sum_{i,j = 1}^k \hat{\beta}_i \hat{\beta}_j K(\vec{x}^{(s_i)}, \vec{x}^{(s_j)}).
\end{equation*}
The coefficient vector $\vec{\hat{\beta}} \in \mathbb{R}^k$ satisfies a smaller linear system
\begin{equation} \label{eq:krr_acc_2}
    \vec{\hat{\beta}} = \bigl( \mat{A}(\set{S},:)\mat{A}(:,\set{S})
    + \lambda N\, \mat{A}(\set{S},\set{S}) \bigr)^{-1}
    \mat{A}(\set{S}, :) \vec{y}.
\end{equation}
involving a $k \times k$ matrix, which makes it relatively inexpensive to solve.

\Cref{alg:krr} provides pseudocode for \mbox{\textsc{RPCholesky}}\xspace-accelerated KRR.
The most expensive steps are the \mbox{\textsc{RPCholesky}}\xspace step and evaluation of the matrix--matrix product $\mat{A}(\set{S},:) \mat{A}(:,\set{S})$; each requires $\mathcal{O}(k^2 N)$ operations.
After performing these steps, the coefficient vector $\vec{\hat{\beta}}$ can be computed in $\mathcal{O}(k^3)$ operations via a standard dense linear solver.
Evaluating the prediction function \cref{eq:krr_acc} from \mbox{\textsc{RPCholesky}}\xspace-accelerated KRR requires just $\mathcal{O}(k)$ operations,
which improves over the $\mathcal{O}(N)$ cost of evaluating \cref{eq:krr_ansatz}.

\begin{algorithm}[h!]
  \caption{\mbox{\textsc{RPCholesky}}\xspace-accelerated kernel ridge regression} \label{alg:krr}
  \textbf{Input:} Data points $\set{X} = \{\vec{x}^{(1)},\ldots,\vec{x}^{(N)}\} \subseteq \mathbb{R}^d$; output values $\vec{y}\in\mathbb{R}^N$; approximation rank $k$; regularization parameter $\lambda > 0$

  \textbf{Output:} Pivots $\set{S}$ and coefficients $\vec{\hat{\beta}}$ defining a prediction function $\hat{f}(\cdot)$ by \cref{eq:krr_acc}
  \begin{algorithmic}
    \State $\mat{A} \leftarrow \Call{KernelMatrix}{\set{X}}$
    \State $(\sim,\set{S}) \leftarrow \Call{RPCholesky}{\mat{A},k}$
    \State $\vec{\hat{\beta}} \leftarrow \Bigl(\mat{A}(\set{S},:)\mat{A}(:,\set{S}) + \lambda N \, \mat{A}(\set{S},\set{S})\Bigr)^{-1} \mat{A}(\set{S},:) \vec{y}$
  \end{algorithmic}
\end{algorithm}

To demonstrate the effectiveness of \mbox{\textsc{RPCholesky}}\xspace-accelerated KRR, we use \cref{alg:krr} to predict the energy of the highest occupied molecular orbital (HOMO) for small organic molecules from the QM9 data set \cite{RDRv14,RvBR12}.
The HOMO energy quantifies the electron-donating capacity of a molecule, and it is typically obtained by expensive first-principles calculations.
Estimating HOMO energies and related properties from data has been a major goal of chemical machine learning \cite{STR+19}.

To represent molecules as vectors for KRR, we apply the Coulomb matrix representation with $\ell_2$ row norm sorting
\cite[{\S}III.A]{STR+19}, a standard approach based on Coulomb repulsions between atomic nuclei.
We standardize the data and, following  \cite{STR+19}, evaluate the similarity between data points using the $\ell_1$ Laplace kernel:
\begin{equation} \label{eq:l1_laplace}
    K(\vec{x},\vec{y}) = \exp \mleft(  - \frac{1}{\sigma} \sum_{j=1}^d \mleft| x_j - y_j \mright| \mright).
\end{equation}
We divide the data into
100,000 data points for training and roughly 33,000 data points for testing and set the bandwidth $\sigma$ and ridge parameter $\lambda$ using cross-validation.
We measure out-of-sample prediction errors using
the symmetric mean absolute percentage error (SMAPE)
\begin{equation} \label{eq:smape}
    \operatorname{SMAPE} = \frac{1}{N_{\rm test}} \sum_{i=1}^{N_{\rm test}} \frac{\mleft| y^{\rm test}_i - \hat{f}\mleft(\vec{x}^{(i),\,{\rm test}} \, ; \, \vec{\hat{\beta}} \mright)  \mright|}{\mleft(\mleft| y^{\rm test}_i\mright| + \mleft| \hat{f}\mleft(\vec{x}^{(i),\,{\rm test}} \, ; \, \vec{\hat{\beta}} \mright) \mright| \mright)/2}.
\end{equation}

Forming and storing the full kernel matrix for the QM9 data set would require 40 GB and 4 trillion arithmetic operations (20,000$\times$ the operation count for the \textbf{Smile} matrix from \cref{sec:compare}), making direct KRR highly expensive.
Prior work \cite{STR+19} applied KRR to this data set using a random subsample of $N \leq \text{64,000}$ training points.
Here, \mbox{\textsc{RPCholesky}}\xspace allows us to use take advantage of more of the data ($N = \text{100,000}$ training points) at a reduced computational cost by using a small approximation rank $k\le 1000 \ll N$.

\begin{figure}[h!]
    \begin{subfigure}{0.49\textwidth}
    \centering
    \includegraphics[width = \textwidth]{chemistry_smape.png}
    \end{subfigure}
    \hfill
    \begin{subfigure}{0.49\textwidth}
    \centering
    \includegraphics[width = \textwidth]{chemistry_errors.png}
    \end{subfigure}
    \caption{
    \textit{Left:} Prediction error \cref{eq:smape} for different Nystr\"om algorithms.
    \textit{Right:} Relative trace-norm error
    for different Nystr\"om algorithms.
    }
    \label{fig:chemistry}
\end{figure}

\Cref{fig:chemistry} displays the results.
The prediction accuracies for \mbox{\textsc{RPCholesky}}\xspace, uniform sampling, and RLS sampling are very close, with uniform sampling being the slight favorite.
However, an examination shows that \textbf{uniform sampling has 18$\times$ higher prediction errors} than \mbox{\textsc{RPCholesky}}\xspace for the molecules in the data set with the largest number of atoms (29); see \cref{tab:big_molecules}.
This observation suggests that \mbox{\textsc{RPCholesky}}\xspace can be better than uniform sampling at representing less populated regions of data space, as seen earlier in the \textbf{Smile} example (\cref{sec:compare}).
The greedy method has notably worse performance than the three other methods, and we were unable to use DPP sampling because of the large values of $N = 10^5$ and $k = 10^3$.

\begin{table}[h!]
    \centering
    \caption{Prediction error \cref{eq:smape} for the largest molecules in the QM9 dataset}
    \label{tab:big_molecules}
    \begin{tabular}{rrcc} \toprule
        Compound \# & Composition & Uniform & \mbox{\textsc{RPCholesky}}\xspace \\ \midrule
        1,996 & \texttt{CC(C)CC1CO1} & 0.519 & 0.030 \\
        8,664 & \texttt{CC(CO)(C=O)C=O} & 0.582 & 0.016 \\
        13,812 & \texttt{CC(O)C1CC1CO} & 0.523 & 0.033 \\
        64,333 & \texttt{CC1(C)CCOC(=N)O1} & 0.525 & 0.027 \\
        81,711 & \texttt{OC1C2CC1OCCO2} & 0.576 & 0.016 \\
        109,816 & \texttt{CCC1C(O)C1CC\#C} & 0.514 & 0.023 \\
        118,229 & \texttt{COCC(C)OC(C)C} & 0.575 & 0.075 \\
        122,340 & \texttt{CC1C(CCCO)N1C} & 0.557 & 0.017 \\
        131,819 & \texttt{OCCCN1C=NC=N1} & 0.570 & 0.030 \\ \bottomrule
    \end{tabular}
\end{table}

\subsection{Spectral clustering} \label{sec:spectral_clustering}

We can also use \mbox{\textsc{RPCholesky}}\xspace to accelerate \emph{kernel spectral clustering} \cite{von07a}.
Kernel spectral clustering is a machine learning approach that uses a nonnegative, positive definite kernel function $K : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}_+$ to identify features in a data set $\vec{x}^{(1)},\ldots,\vec{x}^{(N)}$ and applies clustering based on those features.

In a standard formulation of kernel spectral clustering \cite{von07a}, the features are defined by the dominant left eigenvectors $\vec{v}^{(1)}, \ldots, \vec{v}^{(m)}$ of the matrix $\mat{D}^{-1} \mat{A}$.
In this expression, $\mat{A} \in \mathbb{R}^{N \times N}$ is the kernel matrix
and $\mat{D} \in \mathbb{R}^{N \times N}$ is the diagonal matrix that lists the row sums of $\mat{A}$.
The eigenvectors are normalized such that $\mleft(\vec{v}^{(i)}\mright)^*\mat{D}\vec{v}^{(j)} = \delta_{ij}$,
and each data point $\vec{x}^{(i)}$ is projected as
\begin{equation*} \label{eq:spectral_clustering}
    \vec{\tilde{x}}^{(i)} = \mleft(\vec{v}^{(1)}_i,\vec{v}^{(2)}_i,\ldots,\vec{v}^{(m)}_i \mright).
\end{equation*}
The projected data points $\vec{\tilde{x}}^{(1)}, \ldots, \vec{\tilde{x}}^{(N)}$ are clustered using any standard method,
such as k-means with $m$ centers.

Directly computing the eigendecomposition of $\mat{D}^{-1} \mat{A}$ would require $\mathcal{O}(N^3)$ operations.
However, \mbox{\textsc{RPCholesky}}\xspace provides a cheaper alternative (\cref{alg:spectral_clustering}) that involves just $\mathcal{O}(k^2 N)$ operations.
In the \mbox{\textsc{RPCholesky}}\xspace-accelerated approach, a low-rank approximation $\mat{\hat{A}}^{(k)}$ is used in place of the kernel matrix $\mat{A}$,
and the diagonal matrix $\mat{\hat{D}}$ containing the row sums of $\mat{\hat{A}}^{(k)}$ is used in place of $\mat{D}$.
The data is projected onto the dominant normalized left eigenvectors of $\mat{\hat{D}}^{-1} \mat{\hat{A}}^{(k)}$
and the projected data points are clustered using $k$-means.

\begin{algorithm}[h!]
  \caption{\mbox{\textsc{RPCholesky}}\xspace-accelerated spectral clustering} \label{alg:spectral_clustering}
  \textbf{Input:} Data points $\set{X} = \{\vec{x}^{(1)},\ldots,\vec{x}^{(N)}\} \subseteq \mathbb{R}^d$; number $m$ of clusters; approximation rank $k$

  \textbf{Output:} Partition of $\set{X}$ into $m$ clusters $\set{X} = \set{X}_1\cup \cdots \cup \set{X}_m$
  \begin{algorithmic}
    \State $\mat{A} \leftarrow \Call{KernelMatrix}{\set{X}}$
    \State $\mat{F} \leftarrow \Call{RPCholesky}{\mat{A},k}$
    \State $\mat{\hat{D}}\leftarrow \diag(\mat{F}(\mat{F}^* \mathbb{1}))$
    \State $\mat{G} \leftarrow \mat{\hat{D}}^{-1/2}\mat{F}$
    \State $(\mat{Q},\mat{R}) \leftarrow \Call{QR}{\mat{G}, \texttt{'econ'}}$
    \State $(\mat{U},\sim,\sim) \leftarrow \Call{SVD}{\mat{R}}$
    \State $\mat{V} \leftarrow \mat{\hat{D}}^{-1 \slash 2}\mat{Q}\mat{U}$
    \State $\set{\tilde{X}} \leftarrow \{ \mat{V}(i,1:m) : 1\le i \le |\set{X}| \}$
    \State $\set{X}_1,\ldots,\set{X}_m \leftarrow$ \Call{k-means}{$\set{\tilde{X}}$, $m$}
  \end{algorithmic}
\end{algorithm}

In \cref{fig:clustering}, we apply
\mbox{\textsc{RPCholesky}}\xspace-accelerated spectral clustering to separate the letters ``RPC'', written in papyrus font.
After reading $140$ of the $2187$ columns of the kernel matrix, \mbox{\textsc{RPCholesky}}\xspace reliably produces a perfect clustering as shown in \cref{fig:clustering} (top left).
In contrast, given the same approximation rank, RLS and uniform sampling frequently produce an incorrect clustering.
The fraction of perfect clusterings over $100$ independent trials is given at the bottom of \cref{fig:clustering}.
For this example, the greedy method outperforms all the randomized sampling methods, but we have already seen (\cref{sec:compare,sec:krr}) that the greedy method can be unreliable and has serious failure modes.

\begin{figure}[t]
    \centering
    \begin{minipage}{.48\textwidth}
        \centering
        \includegraphics[height=2.6cm]{clustered.png}

        \mbox{\textsc{RPCholesky}}\xspace
        \includegraphics[height=2.6cm]{rls_clustered_3.png}

        RLS
    \end{minipage}
    \begin{minipage}{.48\textwidth}
        \centering
        \includegraphics[height=2.6cm]{dpp_clustered.png}

        DPP

        \includegraphics[height=2.6cm]{uniform_clustered_1.png}

        Uniform
    \end{minipage}
    {\centering
    \includegraphics[width = .5\textwidth]{accuracies.png}
    }

    \caption{
    \textit{Top:} A correct clustering produced by \mbox{\textsc{RPCholesky}}\xspace along with selected incorrect clusterings produced by RLS, DPP, and uniform sampling with rank $k = 140$.
    \textit{Bottom:} Fraction of correct clusterings produced by Nystr\"om algorithms over $100$ independent trials.
    }
    \label{fig:clustering}
\end{figure}

\section{Theoretical analysis of \mbox{\textsc{RPCholesky}}\xspace}
\label{sec:theoretical-analysis}

Given the appealing computational profile of the \mbox{\textsc{RPCholesky}}\xspace algorithm,
we would like to understand when it is guaranteed to produce an accurate
low-rank approximation.
In this section, we establish the following new result:

\begin{theorem}[Randomly pivoted Cholesky] \label{thm:main_bound}
Fix $r \in \mathbb{N}$ and $\delta \in (0, 1)$, and let $\mat{A}$ be a psd matrix.
The column Nystr\"om approximation $\mat{\hat{A}}^{(k)}$ produced by \mbox{\textsc{RPCholesky}}\xspace (\cref{alg:rpcholesky})
attains the bound
\begin{equation}
    \expect \tr \big(\mat{A} - \mat{\hat{A}}^{(k)} \big)
    \leq \tr (\mat{A} - \lowrank{\mat{A}}_r)
    + \delta \tr \mat{A},
\end{equation}
provided that the number $k$ of columns satisfies
\begin{equation}
    k \geq \frac{r \eta}{\delta} + \min\mleft\{ r \log\mleft(\frac{1}{\delta}\mright), r + r \log_+ \mleft(\frac{2^r \eta}{\delta}\mright)\mright\}.
\end{equation}
  The relative error $\eta$ is defined by $\eta \coloneqq \tr(\mat{A} - \lowrank{\mat{A}}_r) / \tr(\mat{A})$.
  As usual, $\log_+(x) \coloneqq \max\{  \log x, 0 \}$ for $x > 0$.
\end{theorem}

The major takeaway from \cref{thm:main_bound} is that the number $k$ of columns depends \textit{logarithmically} on the tolerance $1/\delta$, as long as $\delta \geq \eta$.
The error bounds for \mbox{\textsc{RPCholesky}}\xspace become more powerful as $\eta \rightarrow 0$ and
\mbox{\textsc{RPCholesky}}\xspace achieves the theoretically optimal error
\begin{equation*}
\expect \tr \big(\mat{A} - \mat{\hat{A}}^{(k)} \big)
= \tr (\mat{A} - \lowrank{\mat{A}}_r)
\end{equation*}
with just $k = r$ columns when $\eta = 0$.

By setting $\delta = \varepsilon \eta$ for $\varepsilon > 0$,
\cref{thm:main_bound} guarantees that \mbox{\textsc{RPCholesky}}\xspace must produce
an $(r, \varepsilon)$-approximation as soon as the number $k$ of columns satisfies
\begin{equation*}
\label{eq:main_bound}
    k \geq \frac{r}{\varepsilon} + r \log_+\mleft( \frac{1}{\varepsilon \eta} \mright).
\end{equation*}
In practice, the logarithmic factor is typically a modest constant, so the bound
is comparable with the minimal cost of $k \geq \frac{r}{\varepsilon}$ columns.
Thus, \mbox{\textsc{RPCholesky}}\xspace improves on the greedy method, in which the number $k$
of columns may be proportional to the dimension $N$ in the worst case.
It also improves on the uniform sampling method, in which the number $k$ of columns can depend \textit{polynomially}
on the relative error $1 / \eta$.  See~\cref{sec:uniform,sec:greedy} for details.

Additionally, \cref{thm:main_bound} demonstrates that \mbox{\textsc{RPCholesky}}\xspace produces an $(r, \varepsilon)$-approximation
when the number $k$ of columns satisfies
\begin{equation}
\label{eq:alternative_bound}
    k \geq \frac{r}{\varepsilon} + r + r \log_+\mleft( \frac{2^r}{\varepsilon} \mright).
\end{equation}
This alternative error bound is significant because it completely eliminates the dependence on the relative error $\eta$.
When $r$ is small, \cref{eq:alternative_bound} is competitive with the cost of
$k \geq \frac{r}{\varepsilon} + r - 1$ columns for DPP sampling.  However, when $r$ is large,
\cref{eq:alternative_bound} deteriorates due to the quadratic dependence on $r$.

Although some error bounds for \mbox{\textsc{RPCholesky}}\xspace can be deduced from earlier work \cite[Thms.~4--5 \&  Prop.~1]{DV06},
\cref{thm:main_bound} gives bounds that are strictly sharper.

\subsection{Proof of Theorem~\ref{thm:main_bound}}

The proof of \cref{thm:main_bound} is based on the properties of the \emph{expected residual function}:
\begin{equation} \label{eq:Phi}
  \mat{\Phi}(\mat{A}) \coloneqq \expect \left[\mat{A}^{(1)}\: \middle| \: \mat{A}\right]
  = \mat{A} - \frac{\mat{A}^2}{\tr\mat{A}}.
\end{equation}
This function returns the expectation of the residual $\mat{A}^{(1)}$ after applying one step of \mbox{\textsc{RPCholesky}}\xspace
to $\mat{A}$.
The equality \eqref{eq:Phi} follows from a short computation using the sampling distribution~\cref{eq:rpc_sampling} and the definition of the residual (\cref{alg:rpcholesky_naive}).
Note that $\mat{\Phi}$ is defined on psd matrices of any dimension.

The first lemma describes some basic facts about the expected residual function $\mat{\Phi}$.
We postpone the proof to \cref{sec:concave_proof}.

\begin{lemma}[Expected residual] \label{lem:concave}
    The expected residual map $\mat{\Phi}$ defined in~\cref{eq:Phi} is positive, monotone, and concave with respect to the psd order. That is, for all psd $\mat{A},\mat{H}$ with the same dimensions,
    %
    \begin{gather}
    \mathbf{0} \preceq \mat{\Phi}(\mat{A}) \preceq \mat{\Phi}(\mat{A} + \mat{H}); \label{eq:monotonicity} \\
    \theta \mat{\Phi}(\mat{A}) + (1-\theta)\mat{\Phi}(\mat{H})
    \preceq \mat{\Phi}(\theta \mat{A} + (1-\theta) \mat{H})
    \quad\text{for all $\theta \in [0,1]$}. \label{eq:concavity}
    \end{gather}
\end{lemma}

The second lemma describes how the trace of the residual declines after multiple steps
of the \mbox{\textsc{RPCholesky}}\xspace procedure.  The proof appears in \cref{sec:contraction_proof}.

\begin{lemma}[Contraction rate] \label{lem:trace_bound}
Consider the $k$-fold composition of the expected residual~\cref{eq:Phi}:
\begin{equation*}
    \mat{\Phi}^{\circ k} \coloneqq \underbrace{\mat{\Phi} \circ \mat{\Phi} \circ \cdots \circ \mat{\Phi}}_{\textnormal{$k$ times}}
    \quad\text{for each $k \in \mathbb{N}$.}
\end{equation*}
Fix $r \in \mathbb{N}$.
For each psd matrix $\mat{H}$ and each $\Delta > 0$,
\begin{equation*}
    \tr \mat{\Phi}^{\circ k}(\mat{H}) \le \tr \mleft(\mat{H} - \lowrank{\mat{H}}_r\mright) + \Delta \tr \mat{H} \quad \text{when} \quad
    k \ge \frac{r \tr \mleft(\mat{H} - \lowrank{\mat{H}}_r\mright)}{\Delta \tr \mat{H}} + r \log_+ \mleft( \frac{1}{\Delta} \mright).
\end{equation*}
\end{lemma}

Last, we present a bound which shows that the error after $k$ steps of \mbox{\textsc{RPCholesky}}\xspace
is comparable with the error in the best rank-$k$ approximation.  The proof appears in \cref{sec:doubling_proof}.

\begin{lemma}[Error doubling] \label{lem:single_step}
  For each psd matrix $\mat{A}$, the residual matrix $\mat{A}^{(k)}$ after applying $k$ steps of \mbox{\textsc{RPCholesky}}\xspace satisfies
  %
  \begin{equation*} \label{eq:repeated_doubling}
    \expect \tr \mat{A}^{(k)} \le 2^k \tr(\mat{A} - \lowrank{\mat{A}}_k)
    \quad\text{for each $k \in \mathbb{N}$.}
  \end{equation*}
\end{lemma}

With these results at hand, we quickly establish the main error bound for \mbox{\textsc{RPCholesky}}\xspace.

\begin{proof}[Proof of \cref{thm:main_bound}]
  Let $\mat{A}$ be a psd matrix.  By the concavity~\cref{eq:concavity} of the expected residual map~\cref{eq:Phi} and the matrix version of Jensen's inequality~\cite[Thm.~4.16]{Car09}, the residual matrices satisfy
  %
  \begin{equation*}
    \expect \mat{A}^{(j)} = \expect\mleft[\expect \left[ \mat{A}^{(j)} \,\middle|\, \mat{A}^{(j-1)} \right]\mright] = \expect \mat{\Phi}(\mat{A}^{(j-1)}) \preceq \mat{\Phi}(\expect \mat{A}^{(j-1)})
    \quad\text{for each $j \in \{1, \dots, k\}$.}
  \end{equation*}
  %
  Next, by monotonicity~\cref{eq:monotonicity} of $\mat{\Phi}$ and the last display,
  %
  \begin{equation*}
  \label{eq:intermediate}
      \expect \mat{A}^{(k)}
      \preceq \mat{\Phi}(\expect \mat{A}^{(k-1)})
      \preceq \mat{\Phi} \circ \mat{\Phi}(\expect \mat{A}^{(k-2)})
      \preceq \cdots \preceq \mat{\Phi}^{\circ k}(\mat{A}).
  \end{equation*}

  Using the fact that the trace is linear and preserves the psd order,
  \begin{equation*}
      \expect \tr\big(\mat{A} - \mat{\hat{A}}^{(k)} \big) =
      \expect \tr \mat{A}^{(k)} \leq \tr \mat{\Phi}^{\circ k}(\mat{A}).
  \end{equation*}
  Next, we apply \cref{lem:trace_bound} with $\mat{H} = \mat{A}$ and $\Delta = \delta$ to see that
  \begin{equation*}
  \expect \tr\big(\mat{A} - \mat{\hat{A}}^{(k)}\big) \leq \tr \mleft(\mat{A} - \lowrank{\mat{A}}_r\mright)
  + \delta \tr \mat{A}
  \quad\text{when}\quad
  k \geq \frac{r \eta}{\delta} + r \log\left(\frac{1}{\delta}\right),
  \end{equation*}
  where we have recognized the relative error $\eta$.

  By a similar argument,
  \begin{equation*}
  \label{eq:k-r}
    \expect \tr\big(\mat{A} - \mat{\hat{A}}^{(k)}\big) = \expect \tr \mat{A}^{(k)} \leq \tr \mat{\Phi}^{\circ (k - r)}\big(\expect \mat{A}^{(r)}\big).
  \end{equation*}
  We can apply \cref{lem:trace_bound} with $\mat{H} = \expect \mat{A}^{(r)}$
  and $\Delta = \delta \cdot \tr \mat{A} \slash \expect \tr \mat{A}^{(r)}$
  to see that
  \begin{equation}
  \label{eq:combine1}
      \expect \tr\big(\mat{A} - \mat{\hat{A}}^{(k)}\big) \leq \tr\mleft( \expect \mat{A}^{(r)} - \lowrank{ \expect \mat{A}^{(r)} }_r \mright)
      + \delta \tr \mat{A}
   \end{equation}
   provided the number of columns satisfies
   \begin{equation}
   \label{eq:combine2}
      k-r \geq \frac{r \tr\mleft( \expect \mat{A}^{(r)} - \lowrank{ \expect \mat{A}^{(r)} }_r \mright)}
      {\delta \tr \mat{A}} + r \log\left(\frac{\expect \tr \mat{A}^{(r)}}{\delta \tr \mat{A}}\right).
  \end{equation}
  This bound can be simplified as follows.
  Observe that the (random) residual matrix $\mat{A}^{(r)}$ is a Schur complement of $\mat{A}$, so it satisfies $\expect \mat{A}^{(r)} \preceq \mat{A}$.
  By the Ky Fan variational principle~\cite[Thm.~8.17]{Zha11}, the best rank-$r$ approximation error in the trace norm is monotone with respect to the psd order, so
  \begin{equation*}
  \tr\mleft( \expect \mat{A}^{(r)} - \lowrank{ \expect \mat{A}^{(r)} }_r \mright)
    \leq \tr( \mat{A} - \lowrank{ \mat{A} }_r ).
  \end{equation*}
  Additionally, \cref{lem:single_step}) guarantees that $\expect \tr \mat{A}^{(r)} \leq 2^r \tr( \mat{A} - \lowrank{\mat{A}}_r)$.
  Using these facts,
  we can simplify \cref{eq:combine1}--\cref{eq:combine2} to show that
  \begin{equation*}
  \expect \tr\big(\mat{A} - \mat{\hat{A}}^{(k)}\big) \leq \tr \mleft(\mat{A} - \lowrank{\mat{A}}_r\mright)
  + \delta \tr \mat{A}
  \quad\text{when}\quad
  k \geq r + \frac{r \eta}{\delta} + r \log\left(\frac{2^r \eta }{\delta}\right),
  \end{equation*}
  This completes the proof of the second bound.
\end{proof}

\subsection{Proof of \cref{lem:concave}} \label{sec:concave_proof}
    Let $\mat{A}, \mat{H}$ be psd matrices, and recall the definition~\cref{eq:Phi} of the expected residual map $\mat{\Phi}$.
    First, to prove that $\mat{\Phi}$ is positive, note that
  \begin{equation}
  \mat{\Phi}(\mat{H}) = \left( \mathbf{I} - \frac{\mat{H}}{\tr \mat{H}} \right)\mat{H}  \succeq \mathbf{0}.
  \end{equation}
  Next, to establish concavity, we choose $\theta \in [0, 1]$ and make the calculation
  %
  \begin{equation*}
    \mat{\Phi}(\theta \mat{A} + (1 - \theta) \mat{H})
    - \theta \mat{\Phi}(\mat{A})
    - (1 - \theta) \mat{\Phi}(\mat{H})
    = \frac{\theta(1 - \theta)}{\theta \tr \mat{A} + (1 - \theta) \tr \mat{H}}
    \left( \sqrt{\frac{\tr \mat{H}}{\tr \mat{A}}} \mat{A} - \sqrt{\frac{\tr \mat{A}}{\tr \mat{H}}} \mat{H} \right)^2 \succeq \mathbf{0}.
  \end{equation*}
  %
  Last, to establish monotonicity, observe that $\mat{\Phi}$ is positive homogeneous; that is, $\mat{\Phi}(\tau \mat{A}) = \tau \mat{\Phi}(\mat{A})$ for $\tau \ge 0$.
  Invoking the concavity property~\cref{eq:concavity} with $\theta = 1/2$,
  %
  \begin{equation*}
    \mat{\Phi}(\mat{A} + \mat{H}) = 2 \, \mat{\Phi} \mleft( \frac{\mat{A}+\mat{H}}{2} \mright) \succeq \mat{\Phi}(\mat{A}) + \mat{\Phi}(\mat{H}) \succeq \mat{\Phi}(\mat{A}).
  \end{equation*}
  We have used the positivity of $\mat{\Phi}(\mat{H})$ in the last step. \hfill $\proofbox$

\subsection{Proof of \cref{lem:trace_bound}}
\label{sec:contraction_proof}

    First, we show that it suffices to consider the case of a diagonal matrix.  Let $\mat{H}$ be an $N \times N$ psd matrix
    with eigendecomposition $\mat{H} = \mat{V}\mat{\Lambda}\mat{V}^*$.  The definition~\cref{eq:Phi} of the expected residual map implies that
    %
    \begin{equation*}
      \mat{\Phi}(\mat{H}) = \mat{V} \mat{\Phi}(\mat{\Lambda}) \mat{V}^*.
    \end{equation*}
    %
    By iteration, the same relation holds with $\mat{\Phi}^{\circ k}$ in place of $\mat{\Phi}$.
    In particular, $\tr \mat{\Phi}^{\circ k}(\mat{H}) = \tr \mat{\Phi}^{\circ k}(\mat{\Lambda})$.
    Therefore, we may restrict our attention to the diagonal case where $\mat{H} = \mat{\Lambda}$.

    Second, we obtain an upper bound on $\tr \mat{\Phi}^{\circ k}(\mat{\Lambda})$.
    Because the map $\mat{\Lambda} \mapsto \tr \mat{\Phi}^{\circ k}(\mat{\Lambda})$
    is concave and invariant under permutations of its arguments,
    averaging together some of the eigenvalues $\lambda_{11}, \ldots, \lambda_{NN}$ can only increase the value of $\tr \mat{\Phi}^{\circ k}(\mat{\Lambda})$.
    Therefore, by introducing the function
    %
    \begin{equation*}
      f_k(a,b,r,q) \coloneqq \tr \mat{\Phi}^{\circ k}\Big( \diag\Big(\underbrace{\frac{a}{r},\ldots,\frac{a}{r}}_{\textnormal{$r$ times}},\underbrace{\frac{b}{q},\ldots,\frac{b}{q}}_{\textnormal{$q$ times}}\Big) \Big),
    \end{equation*}
    %
    we obtain the upper bound
    \begin{equation}
    \label{eq:combine_with}
      \tr \mat{\Phi}^{\circ k}(\mat{\Lambda})
      \leq f_k( \tr \lowrank{\mat{\Lambda}}_r,
      \tr \mleft( \mat{\Lambda} - \lowrank{\mat{\Lambda}}_r\mright), r, N - r).
    \end{equation}
    For further reference, we note that $f_k(a, b, r, q)$ is weakly increasing as a function of $q$.
    Indeed, for every tuple $(a, b, r, q)$,
    %
    \begin{align*}
        f_k(a, b, r, q)
        &= \tr \mat{\Phi}^{\circ k}\Big( \diag\Big(\underbrace{\frac{a}{r},\ldots,\frac{a}{r}}_{\textnormal{$r$ times}},\underbrace{\frac{b}{q},\ldots,\frac{b}{q}}_{\textnormal{$q$ times}}\Big)\Big) = \tr \mat{\Phi}^{\circ k}\Big( \diag\Big(\underbrace{\frac{a}{r},\ldots,\frac{a}{r}}_{\textnormal{$r$ times}},\underbrace{\frac{b}{q},\ldots,\frac{b}{q}}_{\textnormal{$q$ times}},0\Big)\Big) \\
          &\le \tr \mat{\Phi}^{\circ k}\Big( \diag\Big(\underbrace{\frac{a}{r},\ldots,\frac{a}{r}}_{\textnormal{$r$ times}},\underbrace{\frac{b}{q+1},\ldots,\frac{b}{q+1}}_{\textnormal{$q+1$ times}}\Big)\Big) = f_k(a,b,r,q+1).
    \end{align*}
    %
    We have exploited the fact that $\mat{\Phi}$ is defined for matrices of every dimension.

    Next, we derive a worst-case expression for the error $f_k(a, b, r, q)$.
    For each $k = 0, 1, 2, \ldots$,
    define the nonnegative quantities $a^{(k)}$ and $b^{(k)}$ via the relation
    %
    \begin{equation*}
      \mat{\Phi}^{\circ k} \Big( \diag \Big( \underbrace{\frac{a}{r},\ldots,\frac{a}{r}}_{\text{$r$ times}}, \underbrace{\frac{b}{q},\ldots,\frac{b}{q}}_{\text{$q$ times}} \Big) \Big)
      \eqqcolon \diag \Big( \underbrace{\frac{a^{(k)}}{r},\ldots,\frac{a^{(k)}}{r}}_{\text{$r$ times}}, \underbrace{\frac{b ^{(k)}}{q},\ldots,\frac{b^{(k)}}{q}}_{\text{$q$ times}} \Big).
    \end{equation*}
    %
    By the definition~\cref{eq:Phi} of the expected residual map $\mat{\Phi}$, the quantities $a^{(k)}$ and $b^{(k)}$ satisfy the recurrence relations
    %
    \begin{equation*}
      a^{(k)} - a^{(k-1)} = -\frac{ \mleft( a^{(k-1)} \mright)^2 }{r\mleft(a^{(k-1)} + b^{(k-1)}\mright)}
      \quad\text{and}\quad
      b^{(k)} - b^{(k-1)} = -\frac{ \mleft( b^{(k-1)} \mright)^2 }{q\mleft(a^{(k-1)} + b^{(k-1)}\mright)}
    \end{equation*}
    %
    with initial conditions $a^{(0)} = a$ and $b^{(0)} = b$.
    This construction guarantees $f_k(a, b, r, q) = a^{(k)} + b^{(k)}$.
    Additionally, the quantities $a^{(k)}$ and $b^{(k)}$ depend continuously on $q$, and they converge as $q \rightarrow \infty$
    to limiting values $\overline{a}^{(k)}$ and $\overline{b}^{(k)} \equiv b$, where the sequence $\overline{a}^{(k)}$ satisfies
    \begin{equation*}
      \overline{a}^{(k)} - \overline{a}^{(k-1)} = -\frac{ \mleft( \overline{a}^{(k-1)} \mright)^2 }{r\mleft(\overline{a}^{(k-1)} + b\mright)}
      \quad\text{with initial condition $\overline{a}^{(0)} = a$.}
    \end{equation*}
    It follows that
    \begin{equation*}
        f_k(a, b, r, q)
        \leq \overline{a}^{(k)} + b \leq a + b.
    \end{equation*}
    We have used the facts that $f_{k}(a,b,r,q)$ is increasing in $q$ and $\overline{a}^{(k)}$ is decreasing in $k$.

    All that remains is to determine how quickly $\overline{a}^{(k)}$ decreases as a function of $k$.
    To that end, we pass from discrete time to continuous time.
    At each instant $t = 0, 1, 2, \ldots$,
    the discrete-time process $\overline{a}^{(t)}$ is bounded from above
    by the continuous-time process $x(t)$ satisfying the ODE
    %
    \begin{equation*}
      \frac{{\mathrm{d}}}{{\mathrm{d}} t} x(t) = -\frac{x(t)^2}{r(x(t) + b)} \quad\text{with initial condition $x(0) = a$.}
    \end{equation*}
    The comparison between discrete- and continuous-time processes holds because $x \mapsto -x^2 \slash (rx + rb)$ is decreasing over the range $x \in (0, \infty)$.
    Next, assuming $\Delta (a + b) \leq a$, we can use separation of variables to solve for the time $t_{\star}$ at which $x(t_{\star}) = \Delta (a + b)$:
    \begin{equation*}
      t_{\star} =
      \int_{x=a}^{\Delta (a + b)}
      -\frac{r(x + b)}{x^2} \,{\mathrm{d}} x
      = rb \mleft( \frac{1}{\Delta (a + b)} - \frac{1}{a}\mright) + r \log \mleft( \frac{a}{\Delta (a + b)} \mright)
      \leq \frac{rb}{\Delta( a + b)}
      + r \log\left( \frac{1}{\Delta} \right).
    \end{equation*}
    %
    It follows (even for $\Delta (a + b) > a$) that
    \begin{equation*}
      \overline{a}^{(k)} \leq \Delta (a + b)
      \quad \textnormal{when} \quad
    k \geq
    \frac{rb}{\Delta (a + b)} + r \log_+\left( \frac{1}{\Delta} \right).
    \end{equation*}
    To conclude, substitute $a = \tr \lowrank{\mat{\Lambda}}_r$ and $b = \tr\mleft(\mat{\Lambda} - \lowrank{\mat{\Lambda}}_r\mright)$ and combine with \cref{eq:combine_with}.
    \hfill $\proofbox$
\subsection{Proof of \cref{lem:single_step}} \label{sec:doubling_proof}

  Let $\mat{P}$ denote the orthogonal projection onto the $k$ dominant eigenvectors of $\mat{A}$ and set $\mat{P}_\perp \coloneqq \mathbf{I} - \mat{P}$.
  Consider a hierarchical sampling procedure for the pivot index $s_1$:
  %
  \begin{itemize}
  \item With probability $\tr (\mat{A}\mat{P})/\tr \mat{A}$, the ``good'' event $\set{G}$ occurs and the pivot index is chosen as
    %
    \begin{equation*} \label{eq:samp1}
      \prob \mleft( s_1 = i \,\mid\, \set{G} \mright) = \frac{(\mat{A}\mat{P})_{ii}}{\tr (\mat{A}\mat{P})}.
    \end{equation*}
    %
  \item With probability $\tr(\mat{A}\mat{P}_\perp)/\tr \mat{A}$, the ``bad'' event $\set{B} = \set{G}^\set{c}$ occurs and the pivot is selected as
    %
    \begin{equation*} \label{eq:samp2}
      \prob \mleft( s_1 = i \,\mid\, \set{B} \mright) = \frac{(\mat{A}\mat{P}_\perp)_{ii}}{\tr (\mat{A}\mat{P}_\perp)}.
    \end{equation*}
  \end{itemize}
  %
  It suffices to analyze this two-step sampling procedure.  Indeed, a short calculation confirms that this procedure leads to the same sampling distribution for $s_1$ as in \mbox{\textsc{RPCholesky}}\xspace.

  On the bad event $\set{B}$, apply the crude bound
  %
  \begin{equation} \label{eq:combo1}
    \tr\big( \mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1} \big)
    = \sum_{j = k}^{N} \lambda_j(\mat{A}^{(1)})
    \le
    \lambda_k(\mat{A}) + \tr(\mat{A} - \lowrank{\mat{A}}_{k}).
  \end{equation}
  %
  This relation holds because $\mat{A}^{(1)} \preceq \mat{A}$, and the Weyl monotonocity principle \cite[Thm.~8.11]{Zha11} yields $\lambda_j(\mat{A}_1) \le \lambda_j(\mat{A})$ for $j=k,\ldots,N$.

  On the good event $\set{G}$, apply the Ky Fan variational principle \cite[Thm.~8.17]{Zha11} to write
  %
  \begin{equation*}
    \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)
    = \min \mleft\{ \sum_{j=k}^N \vec{u}_j^* \mat{A}^{(1)} \vec{u}_j^{\vphantom{*}} \,:\, \vec{u}_k,\ldots,\vec{u}_N \textrm{ orthonormal} \mright\}.
  \end{equation*}
  %
  Choose the vectors $\vec{u}_{k+1},\ldots,\vec{u}_N$ to be unit-norm eigenvectors $\vec{v}_{k+1}(\mat{A}),\ldots,\vec{v}_N(\mat{A})$ associated with the smallest eigenvalues.  Set $\vec{u}_k \coloneqq \mat{P} \mathbf{e}_{s_1} / \norm{\smash{\mat{P}\mathbf{e}_{s_1}}}$, which is orthonormal to the other vectors by the choice of $\mat{P}$.
  This gives the bound
  %
  \begin{equation} \label{eq:towards_good}
    \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)
    \leq \sum_{j=k+1}^N \vec{v}_j(\mat{A})^* \mat{A}^{(1)} \vec{v}_j(\mat{A})^{\vphantom{*}} + \frac{\mathbf{e}_{s_1}^* \mat{P}\mat{A}^{(1)}\mat{P}\mathbf{e}_{s_1}}{\mathbf{e}_{s_1}^*\mat{P}\mathbf{e}_{s_1}}
    \leq \tr(\mat{A} - \lowrank{\mat{A}}_k) + \frac{(\mat{P}\mat{A}^{(1)}\mat{P})_{s_1s_1}}{p_{s_1s_1}}.
  \end{equation}
  %
  Using the fact that $\mat{A}^{(1)} = \mat{A} - \mat{A} \mathbf{e}_{s_1}\mathbf{e}_{s_1}^* \mat{A} / a_{s_1s_1}$, it follows
  %
  \begin{equation*}
    (\mat{P}\mat{A}^{(1)}\mat{P})_{s_1s_1} = (\mat{A}\mat{P})_{s_1s_1} - \frac{((\mat{A}\mat{P})_{s_1s_1})^2}{a_{s_1s_1}}
    = \frac{(\mat{A}\mat{P})_{s_1s_1}(\mat{A}\mat{P}_\perp)_{s_1s_1}}{a_{s_1s_1}}.
  \end{equation*}
  %
  Therefore, taking the conditional expectation of both sides of \cref{eq:towards_good},
  %
  \begin{equation} \label{eq:combo2}
  \begin{split}
    \expect \left[ \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr) \, \middle|\, \set{G} \right]
    &\le \tr(\mat{A} - \lowrank{\mat{A}}_k)
    + \sum_{i=1}^N \frac{(\mat{A}\mat{P})_{ii}}{\tr (\mat{A}\mat{P})} \cdot \frac{(\mat{A}\mat{P})_{ii}(\mat{A}\mat{P}_\perp)_{ii}}{a_{ii} p_{ii}} \\
    &\le \mleft( 1 + \frac{\lambda_1(\mat{A})}{\tr(\mat{A}\mat{P})} \mright) \tr(\mat{A} - \lowrank{\mat{A}}_k).
    \end{split}
  \end{equation}
  %
  The last line uses the inequalities $(\mat{A}\mat{P})_{ii} \le a_{ii}$ and $(\mat{A}\mat{P})_{ii} \le \lambda_1(\mat{A}) p_{ii}$.

  Finally, combine the bound \cref{eq:combo1} for the bad event and the bound \cref{eq:combo2} for the good event to verify that
  %
  \begin{align*} \label{eq:single_step_1}
    \expect \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)
    &= \expect \left[ \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)  \,\middle|\, \set{G} \right] \cdot \prob(\set{G})
    + \expect \left[ \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)  \,\middle|\, \set{B} \right] \cdot \prob(\set{B}) \\
    &\leq \mleft( 1 + \frac{\lambda_1(\mat{A}) + \lambda_k(\mat{A})}{\tr \mat{A}} \mright)
    \tr(\mat{A} - \lowrank{\mat{A}}_{k}).
  \end{align*}
  Since $\lambda_1(\mat{A}) + \lambda_k(\mat{A}) \le \tr \mat{A}$, this bound can be weakened to give
  %
  \begin{equation*}
    \expect \tr\bigl(\mat{A}^{(1)} - \lowrank{\mat{A}^{(1)}}_{k-1}\bigr)
    \leq 2 \tr(\mat{A} - \lowrank{\mat{A}}_{k}).
  \end{equation*}
  %
  By iterating, we conclude that
  %
  \begin{equation*}
    \expect \tr \mat{A}^{(k)}
    \le 2 \expect \tr\bigl(\mat{A}^{(k-1)} - \bigl\llbracket \mat{A}^{(k-1)} \bigr \rrbracket_1 \bigr)
    \le 4 \expect \tr\bigl(\mat{A}^{(k-2)} - \bigl \llbracket \mat{A}^{(k-2)} \bigr \rrbracket_2 \bigr)
    \leq \cdots \leq
    2^k \tr(\mat{A} - \lowrank{\mat{A}}_k).
  \end{equation*}
  This estimate is the statement of the doubling bound. \hfill $\proofbox$

\section{Comparison with other methods}
\label{sec:related-work}

This section compares \mbox{\textsc{RPCholesky}}\xspace with other psd low-rank approximation methods.
We cover the greedy method (\cref{sec:greedy}), uniform sampling (\cref{sec:uniform}), DPP sampling (\cref{sec:dpp}), RLS sampling (\cref{sec:rls}), and
alternative algorithms that are not based on the column Nystr\"om approximation (\cref{sec:non-nystrom}).
For brevity, the proofs are deferred to \cref{app:more_proofs}.

\subsection{The greedy method}
\label{sec:greedy}

The \emph{greedy method} \cref{eq:greedy} is a column Nystr\"om approximation with a long history in numerical analysis \cite[\S10.3]{Hig02} under the name \textit{complete pivoting} or \textit{diagonal pivoting}.
Fine and Scheinberg \cite{FS01} later popularized the method's use for kernel computations.
Despite its popularity, the greedy method is known to fail dramatically when applied to some worst-case matrices $\mat{A}$, as discussed by Higham \cite[\S10.3.1]{Hig02}.
Further supporting this conclusion, the following theorem analyzes the large number of columns needed for the greedy method to provide a $(r, \varepsilon)$-approximation.
\begin{theorem}[Greedy method] \label{thm:greedy}
The greedy method has the following error properties:
\begin{enumerate}[label=(\alph*)]
    \item If $k \geq (1 - (1 + \varepsilon) \eta) N$,
    the Nystr\"om approximation $\mat{\hat{A}}^{(k)}$ produced by the greedy method is a $(r, \varepsilon)$-approximation \cref{eq:1+eps}. \label{item:greedy_a}
    \item For any $N \geq r \geq 1$, $\varepsilon > 0$, and $\eta \in [0, 1]$, there exists a psd matrix $\mat{A} \in \mathbb{R}^{N \times N}$ such that $\tr(\mat{A} - \lowrank{\mat{A}}_r) = \eta \tr \mat{A}$ and $\mat{\hat{A}}^{(k)}$ produced by the greedy method is not a $(r, \varepsilon)$-approximation for $k \leq
    \min\{\eta N, (1 - (1 + \varepsilon) \eta) N\} - 1$. \label{item:greedy_b}
\end{enumerate}
\end{theorem}
\begin{proof}
See \cref{app:more_proofs}.
\end{proof}

\cref{thm:greedy} shows that the greedy method requires accessing $\Omega(N)$ columns to make progress in approximating a worst-case matrix $\mat{A}$.
This potential for slow convergence is a major concern in applications, and \cref{sec:compare,sec:krr} provide examples of kernel matrices for which the greedy method produces worse-quality approximations than other Nystr\"om methods.
To address the deficiencies of the greedy method, Gu and Miranian \cite{GM04} introduced a \emph{strong rank-revealing} pivoting strategy which produces higher-quality approximations on worst-case examples.
However, this approach is sufficiently complicated that it has not seen widespread use.

\subsection{Uniform sampling}
\label{sec:uniform}

Another popular column Nystr\"om approximation scheme selects pivots \emph{uniformly at random}, either with or without replacement.
Uniform sampling without replacement leads to a larger number of unique columns selections, hence, a higher-quality Nystr\"om approximation \cite{KMT09}.
Uniform sampling can be effective when the dominant $r$ eigenvectors have mass that is spread out equally over all the coordinates (a property known as ``incoherence'') \cite{Git11}.
However, uniform sampling fails when the dominant eigenvectors are heavily concentrated on just a few coordinates.
This potential for bad behavior is reflected in the following worst-case error bounds:

\begin{theorem}[Uniform sampling]     \label{thm:uniform}
Suppose $\mat{A}$ is a psd matrix with ones on the diagonal,
and the pivots are sampled uniformly at random either with or without replacement.
\begin{enumerate}[label=(\alph*)]
    \item If $k \geq \frac{r-1}{\varepsilon \eta} + \frac{1}{\varepsilon}$, the Nystr\"om approximation $\mat{\hat{A}}^{(k)}$ produced by uniform sampling is a $(r, \varepsilon)$-approximation. \label{item:uniform_a}
    \item For any $r \geq 2$ and any $\varepsilon > 0$ and $\eta > 0$ satisfying $\sqrt{\frac{\varepsilon}{\varepsilon + 1}} + \eta < 1$, there exists a psd matrix $\mat{A}$ such that $\tr(\mat{A} - \lowrank{\mat{A}}_r) = \eta \tr \mat{A}$ and $\mat{\hat{A}}^{(k)}$ produced by uniform sampling is not a $(r, \varepsilon)$-approximation for $k < \frac{r-1}{\eta} (\sqrt{\varepsilon^{-1} + 1} - 1)^2$. \label{item:uniform_b}
\end{enumerate}
\end{theorem}
\begin{proof}
See \cref{app:more_proofs}.
\end{proof}

Although uniform sampling is typically applied to a kernel matrix $\mat{A}$ with ones on the diagonal,
a more general \emph{diagonal sampling} \cite{FKV04} method randomly selects the pivots with probabilities proportional to $\diag \mat{A}$.
Diagonal sampling attains exactly the same error bounds as uniform sampling,
as seen from the proof of \cref{thm:uniform}.
\mbox{\textsc{RPCholesky}}\xspace is equivalent to performing diagonal sampling iteratively on the residual matrix, and it leads to significantly stronger error bounds.

\subsection{Determinantal point process sampling}
\label{sec:dpp}

Determinantal point process (DPP) sampling \cite{DM21a} is a column Nystr\"om approximation method that randomly selects a pivot set $\set{S} \subseteq \{1,\ldots,N\}$ with cardinality $|\set{S}| = k$ according to the distribution
\begin{equation}
\label{eq:exact_dpp}
\prob \mleft\{ \{s_1, \ldots, s_k\} = \set{S} \mright\} = \frac{\det \mat{A}(\set{S}, \set{S})}
{\sum_{|\set{S}^{\prime}| = k}
\det \mat{A}(\set{S}^{\prime}, \set{S}^{\prime})}.
\end{equation}
DPP sampling
satisfies the following error bounds (see \cite[Thm.~1]{BW09}, \cite{GS12}):

\begin{theorem}[$k$-DPP sampling] \label{thm:near-optimal}
The Nystr\"om approximation produced by $k$-DPP sampling is nearly optimal in the following sense:
\begin{enumerate}[label=(\alph*)]
    \item If $k \geq r/\varepsilon + r - 1$,
    the approximation $\mat{\hat{A}}^{(k)}$ produced by a $k$-DPP is an $(r, \varepsilon)$-approximation. \label{item:dpp_a}
    \item For any $r \geq 1$ and $\varepsilon > 0$, there exists a psd matrix $\mat{A}$ such that $\mat{\hat{A}}^{(k)}$ produced by a $k$-DPP is not an $(r, \varepsilon)$-approximation for $k < r/\varepsilon + r-1$. \label{item:dpp_b}
    \item For any $r \geq 1$ and $\varepsilon > 0$, there exists a matrix $\mat{A}$ such that no rank-$k$ column Nystr\"om approximation is an $(r, \varepsilon)$-approximation for $k < r/\varepsilon$. \label{item:dpp_c}
\end{enumerate}
\end{theorem}
\begin{proof}
Parts \ref{item:dpp_a} and \ref{item:dpp_b} were essentially proved in \cite{BW09} and \cite{GS12}.
For part \ref{item:dpp_c}, see \cref{app:more_proofs}.
\end{proof}

The theory above demonstrates that DPP sampling has nearly optimal theoretical properties, but implementing DPP sampling for large $k$ values remains expensive relative to peer methods \cite{AGR16,DCV19}.
\mbox{\textsc{RPCholesky}}\xspace is much cheaper to perform and is identical to performing $1$-DPP sampling iteratively on the residual matrix.

\subsection{Ridge leverage score sampling}
\label{sec:rls}

Like DPP sampling, ridge leverage score (RLS) sampling \cite{AM15,CLV17,MM17,RCCR18} is a Nystr\"om approximation defined by a distribution on the pivot set $\set{S}$.
However, the sampling distribution is potentially more tractable than in DPP sampling since it only depends on a single matrix inverse rather than numerous determinants of $k\times k$ submatrices.
For any $\lambda > 0$, the vector of ridge leverage scores is defined as
\begin{equation*} \label{eq:rls}
    \vec{\ell}^{\lambda}
    = \diag (\mat{A} (\mat{A} + \lambda \mathbf{I})^{-1}),
\end{equation*}
and the vector of sampling probabilities is defined by \cite{MM17} to be
\begin{equation*}
    \vec{p} = \min\left\{1, 16 \,\vec{\ell}^{\lambda} \log\left(\frac
    {1}
    {\delta}
    \sum_{j=1}^N \vec{\ell}^{\lambda}_j \right)\right\}
\end{equation*}
where $\delta \in (0, 1/8)$ is a parameter called the ``failure probability''.
The coordinate set $\set{S} \subseteq \{1,\ldots,N\}$ is generated
by independently including each index $i$ with probability $\vec{p}_i$.
This form of RLS sampling has been analyzed by~\cite[Thm.~3]{MM17},
and their main result can be summarized as follows:
\begin{theorem}[Ridge leverage score sampling]
\label{thm:musco}
With probability at least $1 - \delta$,
the Nystr\"om approximation $\mat{\hat{A}}$ produced by RLS sampling satisfies
\begin{equation*}
\label{eq:musco1}
    \lVert \mat{A} - \mat{\hat{A}} \rVert \leq \lambda
\end{equation*}
and
\begin{equation*}
\label{eq:musco2}
    \rank \mat{\hat{A}} \leq 32 \mleft(r + \frac{\tr (\mat{A} - \lowrank{\mat{A}}_r)}{\lambda}\mright)\,
    \log\mleft( \frac{r}{\delta} + \frac{\tr (\mat{A} - \lowrank{\mat{A}}_r)}{\delta \lambda}\mright).
\end{equation*}
\end{theorem}

\Cref{thm:musco} is different from the other error bounds presented in this work since it guarantees control in the spectral norm, rather than the trace norm, and it introduces a random approximation rank and a failure probability $\delta$ not seen in the other bounds.
Due to these differences, the following calculation is needed to convert \cref{thm:musco} into a bound in the expected trace norm:
\begin{corollary}[Ridge leverage score $(r,\varepsilon)$-approximation] \label{cor:musco}
For any $\varepsilon < \frac{1}{4}$, let $\mat{\hat{A}}$ be the Nystr\"om approximation obtained from RLS sampling with
\begin{equation*}
    \lambda = \frac{\varepsilon}{2r} \tr (\mat{A} - \lowrank{\mat{A}}_r), \qquad \delta = \frac{\varepsilon \eta}{2},
\end{equation*}
and set
\begin{equation*}
    \mat{\hat{A}}^{(k)} =
    \begin{cases}
    \mat{\hat{A}}, & \rank \mat{\hat{A}} \leq k, \\
    \mat{0}, & \rank \mat{\hat{A}} > k.
    \end{cases}
\end{equation*}
Then,
$\mat{\hat{A}}^{(k)}$ is an $(r, \varepsilon)$-approximation if
$k \geq 136 \frac{r}{\varepsilon}
\log\Bigl(\frac{3r}{\varepsilon \eta}\Bigr)$.
\end{corollary}
\begin{proof}
See \cref{app:more_proofs}.
\end{proof}

\Cref{cor:musco} guarantees better worst-case performance for RLS sampling than for the greedy method or uniform sampling,
at least when $N$ and $\eta^{-1}$ are very large.
However, the available implementations of RLS sampling \cite{CLV17,MM17,RCCR18} are fairly complicated,
as they require selecting parameters $\lambda$ and $\delta$ and approximating the resulting RLS sampling distribution.
Because of these preprocessing steps, RLS sampling requires a higher number of entry evaluations than \mbox{\textsc{RPCholesky}}\xspace for a fixed approximation rank $k$.

\subsection{Algorithms that are not based on column Nystr\"om approximation} \label{sec:non-nystrom}

In addition to column Nystr\"om approximation, there are two other popular approaches for randomized low-rank approximation of psd matrices.
First, \emph{random embeddings} can produce a high-quality approximation (even higher quality than \mbox{\textsc{RPCholesky}}\xspace) if it is possible to compute matrix--vector products with the input matrix.
Yet, these methods may be too expensive to deploy in the entry-wise access model.
Random embedding methods include the single-pass Nystr\"om method \cite{TYUC17b}, randomized subspace iteration \cite{HMT11}, and randomized block Krylov iteration \cite{TW22}.

Second, \emph{random features} (RF) methods
can be used to approximate kernel matrices for many well-known kernel functions $K$ \cite{MT20,RR07}.
RF methods do not require any entry evaluations,
and the computations are highly parallelizable.
However, due to the dependence on Monte Carlo sampling, the accuracy is limited by the $1/\sqrt{k}$ convergence rate from the central limit theorem, regardless of the decay of the eigenvalues in the input matrix \cite{yang2012nystrom}.
In contrast, \mbox{\textsc{RPCholesky}}\xspace takes advantage of the eigenvalue decay, so it can attain much higher accuracy.

\section{Conclusion} \label{sec:conclusion}

This work has demonstrated the power of \mbox{\textsc{RPCholesky}}\xspace for the low-rank approximation of a positive semidefinite kernel matrix $\mat{A} \in \mathbb{R}^{N \times N}$.
\mbox{\textsc{RPCholesky}}\xspace can be programmed in just a few lines of code, and it accelerates many kernel algorithms such as kernel ridge regression and spectral clustering.
\mbox{\textsc{RPCholesky}}\xspace reduces the computational cost of these algorithms from $\mathcal{O}(N^3)$ operations to just
$\mathcal{O}(k^2 N)$ operations,
where the approximation rank $k$ can be much smaller than the matrix dimension $N$.

Numerical experiments suggest that \mbox{\textsc{RPCholesky}}\xspace dominates other column Nystr\"om approximation algorithms in terms of efficiency and memory footprint.
Given a fixed approximation rank $k$, \mbox{\textsc{RPCholesky}}\xspace requires a very small number of entry evaluations, just $(k + 1) N$.
\mbox{\textsc{RPCholesky}}\xspace typically produces approximations that match or improve on the greedy method, uniform sampling, RLS sampling, and DPP sampling.
Moreover, theoretical error bounds guarantee that \mbox{\textsc{RPCholesky}}\xspace converges nearly as fast as possible in the expected trace norm.

Taken as a whole, this work paves the way for greater use of \mbox{\textsc{RPCholesky}}\xspace.
The method can potentially be pushed even further, for example, through combinations with other accelerated methods for prediction, clustering, and other learning tasks.
The future is indeed bright for this simple yet surprisingly effective algorithm.

\section*{Acknowledgements}

We thank Mateo D{\'i}az, Zachary Frangella, Marc
Gilles, Eitan Levin, Eliza O'Reilly, and Jonathan Weare for helpful discussions and corrections.