Randomized Strong Recursive Skeletonization: Simultaneous compression and factorization of $\mathcal{H}$-matrices in the Black-Box Setting

Anna Yesypenko and Per-Gunnar Martinsson

The University of Texas at Austin

The hierarchical matrix ($\mathcal{H}^{2}$-matrix) formalism provides a way to reinterpret the Fast Multipole Method and related fast summation schemes in linear algebraic terms.
The idea is to tessellate a matrix into blocks in such as way that each block is either small or of numerically low rank; this enables the storage of the matrix and the application of it to a vector in linear or close to linear complexity.
A key motivation for the reformulation is to extend the range of dense matrices that can be represented.
Additionally, $\mathcal{H}^{2}$-matrices in principle also extend the range of operations that can be executed to include matrix inversion and factorization.
While such algorithms can be highly efficient for certain specialized formats (such as HBS/HSS matrices based on “weak admissibility”), inversion algorithms for general $\mathcal{H}^{2}$-matrices tend to be based on nested recursions and recompressions, making them challenging to implement efficiently.
An exception is the strong recursive skeletonization (SRS) algorithm by Minden, Ho, Damle, and Ying, which involves a simpler algorithmic flow.
However, SRS greatly increases the number of blocks of the matrix that need to be stored explicitly, leading to high memory requirements.
This manuscript presents the randomized strong recursive skeletonization (RSRS) algorithm, which is a reformulation of SRS that incorporates the randomized SVD (RSVD) to simultaneously compress and factorize an $\mathcal{H}^{2}$-matrix.
RSRS is a “black box” algorithm that interacts with the matrix to be compressed only via its action on vectors; this extends the range of the SRS algorithm (which relied on the “proxy source” compression technique) to include dense matrices that arise in sparse direct solvers.
Moreover, RSRS leads to dramatically simpler data structures, faster execution time, and easier parallelization.
RSRS is particularly effective when applied to geometries that have lower dimensionality than the ambient space (as when discretizing a boundary integral equation on a surface in $\mathbb{R}^{3}$), as it enables off-diagonal blocks to be compressed using their actual numerical ranks, rather than the rank imputed by a proxy surface.

§ INTRODUCTION

Dense matrices that arise in mathematical physics often have
internal structure that makes it possible to solve problems
involving elliptic operators in linear or close to linear time.
Early techniques such as the Fast Multipole Method exploited
mathematical properties of the kernel function directly to enable
the fast application of dense matrices to vectors. The
$\mathcal{H}^{2}$-matrix methodology followed, and reinterpreted the
FMM in linear algebraic terms, where a matrix is tessellated into
blocks in such as way that each block is either small or of
numerically low rank. This reinterpretation opened the door to
linear complexity algorithms for a wider range of algebraic
operations, including the matrix-matrix multiplication, and the
construction of invertible factorizations.

The manuscript describes the algorithm
Randomized Strong Recursive Skeletonization (RSRS)
for simultaneously compressing
and inverting an $\mathcal{H}^{2}$-matrix, given a means of applying the
matrix and its adjoint to vectors.
The precise problem formulation is this: Suppose that $\mtx{A}$ is
an $\mathcal{H}^{2}$-matrix, and that you are given a fast method for
applying $\mtx{A}$ and its adjoint $\mtx{A}^{*}$ to vectors.
We then seek to build two “test-matrices” $\mtx{\Omega}$ and $\mtx{\Psi}$ with the property that the $\mathcal{H}^{2}$ representation of $\mtx{A}^{-1}$ can be constructed from the set $\{ \mtx Y, \mtx Z, \mtx \Omega, \mtx \Psi\}$, where

\begin{equation*}
\mtx{Y} = \mtx{A}\mtx{\Omega}
\qquad\mbox{and}\qquad
\mtx{Z} = \mtx{A}^* \mtx{\Psi}.
\end{equation*}

RSRS is immediately applicable in a range of important environments.
First, it can be
used to derive a rank-structured representation of any integral operator
for which a fast matrix-vector multiplication algorithm,
such as the Fast Multipole Method <cit.>,
is available.
In this context, RSRS will directly output an invertible factorization
of the integral operator.
Second, it can greatly simplify algebraic operations involving products of
rank-structured or sparse matrices. As an illustration, the Dirichlet-to-Neumann (DtN)
operator for a bounded domain can often be constructed using boundary integral equation
techniques, with the overall action of the DtN operator calculated by applying an
iterative solver combined with a post-processing step. RSRS would let us directly
build and factorize the DtN operator.
Relatedly, the perhaps key application of rank-structured
matrix algebra is the acceleration of sparse direct solvers, as the
dense matrices that arise during LU factorization are often rank structured
In the course of such a solver, a typical operation would
be to form a Schur complement such as
$\mtx{S} = \mtx{B}(I,J)\mtx{B}(J,J)^{-1}\mtx{B}(J,I)$,
where $I$ and $J$ are two index sets.

If $\mtx{B}(J,J)$ is an $\mathcal{H}^{2}$-matrix, then $\mtx{B}(J,J)^{-1}$
can easily be applied to vectors via an LU factorization. If, additionally,
$\mtx{B}(I,J)$ and $\mtx{B}(J,I)$ are either sparse or rank structured,
then $\mtx{S}$ can easily be applied to a vector.
RSRS in this environment allows us to directly compute a
factorization of $\mtx{S}$, which is precisely what is needed in
order to move the overall LU factorization forwards.

The fact that RSRS simultaneously compresses and inverts the matrix is a key
feature of the algorithm.
Previous work on black-box randomized algorithms for compressing rank
structured matrices <cit.> has approached the tasks of compression and inversion as two separate computational stages, with one executed before the other.
This approach is very natural for rank structured formats based on “weak admissibility” (where all off-diagonal blocks of the matrix are treated as low rank) since in this case simple and exact inversion algorithms are available <cit.>.
In the “strong admissibility” case (where only interactions between well separated parts of the computational domain are compressed), LU factorization or inversion are much more challenging tasks that require repeated recompression of off-diagonal blocks as the algorithm proceeds <cit.>.
This complication has greatly limited the appeal of direct solvers based on rank structured matrices for problems in involving fully three dimension geometries where strong admissibility is essential.

In terms of prior work, RSRS draws on the randomized SVD (RSVD) for
low rank approximation of matrices <cit.>, and on earlier methods that apply the RSVD to the problem of reconstructing rank structured matrices from matrix-vector products
It in particular draws on ideas from the recent dissertation of James Levitt <cit.>.
Finally, it draws heavily from the prior work on strong recursive skeletonization, including the original work <cit.>, and the more recent <cit.>.

The manuscript is structured as follows:
Section <ref> surveys the key ideas underlying the RSVD and how they can be applied to simultenously compress all the off-diagonal blocks of a rank structured matrix. Section <ref> also introduces the core idea that is the main contribution of the present work.
Section <ref> introduces the interpolatory decomposition and shows how it can be used to solve simple linear systems involving rank deficiencies in their coefficient matrices.
Section <ref> introduces our interpretation of the Strong Recursive Skeletonization (SRS) algorithm.
Section <ref> describes Randomized Strong Recursive Skeletonization in full, and discusses refinements to the basic scheme.
Section <ref> presents numerical experiments that demonstrate how
RSRS performs in terms of speed, memory requirements, and precision.

§ RANDOMIZED COMPRESSION AND FACTORIZATION OF RANK STRUCTURED MATRICES

In this section, we highlight the key methodological contributions of the manuscript
for a simple model problem. For points $\{\pxx_{j}\}_{j=1}^{N}$ shown in Figure <ref>,
consider the $N\times N$ matrix $\mtx{A}$ with entries

\begin{equation} \label{eq:kernel_eval}
\mtx{A}(i,j) = k(\pxx_{i},\,\pxx_{j})
\end{equation}

given by a kernel function $k$, for instance $k(\pxx,\pyy) = \log|\pxx - \pyy|$. The singularity at $x=y$ can be
handled with an appropriate quadrature rule.
The matrix (<ref>) and its adjoint can be applied rapidly to vectors using the fast multipole method.
The objective is to recover an approximate representation of $\mtx A^{-1}$ using random sketches

\begin{equation} \label{eq:rand_sample_setting}
\underset{N \times p}{\mtx Y}\ =\ \underset{N \times N}{\mtx A}\ \underset{N \times p}{\mtx \Omega}, \qquad
\underset{N \times p}{\mtx Z}\ =\ \underset{N \times N}{\mtx A^*}\ \underset{N \times p} {\mtx \Psi}, \qquad \mtx \Omega, \mtx \Psi \sim \mathcal{N}(\mtx 0, \mtx I)
\end{equation}

using Gaussian test matrices $\mtx \Omega, \mtx \Psi$. We will discuss the needed number of samples later, but for now, consider that $p \ll N$.
The material in Section 2.1 follows <cit.>, while the techniques in Sections 2.2 and 2.3 are based on ideas in <cit.>.

[fill=none, thick] (, +*)̱ rectangle (+ *,̱ +*+̱ )̱;
plot [only marks, mark=*, mark size=0.5, domain=0:8, samples=700] ( ,rnd + 0.5 + 7.5);
at (-2.5*,̱ ++0.5*)̱ Points $\{\pxx_j\}_{j=1}^N$;

* 0.25
at (-1.5*,̱ ++0.5*)̱ Level 1;
[fill=none, thick] (, +*)̱ rectangle (+ *,̱ +*+̱ )̱;
in 1,...,

[thick] (+ 4**,̱ +*)̱ – (+ 4**,̱ +*+̱)̱;
[scale=0.8] at (+ 4**-̱ 2*,̱ +*+̱0.5*)̱ ${\xtmp}$;

* 0.5
at (-1.5*,̱ ++0.5*)̱ Level 2;
[fill=none, thick] (, +*)̱ rectangle (+ *,̱ +*+̱ )̱;
in 1,...,

[thick] (+ 2**,̱ +*)̱ – (+ 2**,̱ *+̱ +)̱;
[scale=0.8] at (+ 2**-̱ ,̱ +*+̱0.5*)̱ ${\xtmp}$;

at (-1.5*,̱ ++0.5*)̱ Level 3;
[fill=none, thick] (, +*)̱ rectangle (+ *,̱ +*+̱ )̱;
in 1,...,
[thick] (+ * ,̱ +*)̱ – (+ *,̱ +*+̱ )̱;
[scale=0.8] at (+ *-̱ 0.5*,̱ +*+̱0.5*)̱ ${\xtmp}$;

The hierarchical tessellation of boxes for a quasi one dimensional model problem.
Boxes on the same level which share an edge or corner are called neighbors. Two boxes which are not neighbors are called well-separated.
A box b has a far-field consisting of all boxes which are well-separated from b. As an example, box 9 has neighbors {8,9,10} and far-field {11,…, 15}.

The methodology uses that $\mtx A$ has certain rank deficiencies to recover a sparse factorization of $\mtx A^{-1}$ efficiently.
Consider the tesselation of points into boxes $8, \dots, {15}$ on Level 3 of Figure <ref>, where each box has $m:=N/8$ points.
The matrix $\mtx A$ has the following property for interactions between far-field blocks

\begin{equation}\label{e:low_rank_far}
\underset{m \times m} {\mA_{i,j}}\
\approx\
\underset{m \times k} {\mtx Q_{i\phantom{j}}}\ \underset{k \times k} {\mtA_{ij}}\ \underset{k \times m}{\mtx U_j^*}\end{equation}

where index set $I_i$ and $I_j$ correspond to the indices of points of well-separated colleague boxes $i$ and $j$, respectively.
The low rank property of far-field interactions in (<ref>) is depicted in Figure <ref>.
The blocks corresponding to interactions between near neighbors are stored densely.

In the next sections, we discuss methods for recovering the bases $\mtx Q_i$ and $\mtx U_j$, which we call block nullification,
in the black-box setting of (<ref>).
We also discuss techniques for recovering block-sparse interactions (e.g. interactions between near neighbors),
which we call block extraction.
First, to provide context, we review standard methods for recovering globally low-rank factors from matrix-vector products.

§.§ Review of randomized sketching for a low rank matrix

Suppose we would like to compute a low-rank approximation to the matrix $\mtx A \in \mathbb{R}^{m \times n}$
where the rank $k$ is approximate and known apriori. Randomized low rank compression
provides a powerful set of techniques for accomplishing this task via the action of $\mtx A$ and its adjoint on a small number of vectors.

A low rank approximation to $\mtx A$ can be computed in two stages. First, we would like to compute $\mtx Q \in \mathbb{R}^{m \times k}$ with orthonormal columns for which
$\mtx A \approx \mtx Q \mtx Q^* \mtx A$. Once $\mtx Q$ is known, then we can form the matrix $\mtx B = \mtx A^* \mtx Q$ and the low-rank factorization of $\mtx A$ takes the form

\begin{equation}\label{eq:qb}
\underset{m \times n}{\mtx A} = \underset{m \times k} {\mtx Q}\ \underset{k \times n}{\mtx B \vphantom{Q}}.
\end{equation}

To find an approximate basis for the range of $\mtx A$ (e.g. find $\mtx Q$ such that $\mtx A \approx \mtx Q \mtx Q^* \mtx A$), we generate a randomized sketch of $\mtx A$

\begin{equation}
\underset{m \times (k+l)}{\mtx Y}\ =\ \underset{m \times n \vphantom{()}}{\mtx A}\ \underset{n \times (k+l)}{\mtx \Omega}, \qquad \mtx \Omega \sim \mathcal{N}(\mtx 0,\mtx I),
\end{equation}

where $l$ is a small oversampling parameter (e.g. $l=5$).
The sketch $\mtx Y$ approximately spans the column space of $\mtx A$, and desired orthonormal basis $\mtx Q$
can be computed by running a procedure to orthonormalize $\mtx Y$

\begin{equation}
\underset{m \times (k+l)} {\mtx Q} = \texttt{orth}\big(\mtx Y\big), \qquad \text{where}\ \mtx Y = \mtx Q \mtx R
\end{equation}

When the rank $k$ is not known apriori, there are adaptive algorithms available which build a low rank representation
of $\mtx A$ from successive sketches. Randomized sketching methods can be used to construct a wide range of decompositions,
including the interpolative decomposition which we introduce in Section <ref>.
These algorithms are especially useful in the black-box setting of (<ref>) because the matrix $\mtx A$ is only accessed through its action on vectors.
In the next section, we discuss a modification of
randomized low rank matrices which can be used to recover the orthogonal bases for the row and column blocks of (<ref>).

§.§ Block Nullification

In the black-box setting of (<ref>), we would like to recover the bases
$\mtx Q_i,\mtx U_i \in \mathbb{R}^{m \times k}$ for all blocks $i= 8, \dots, 15$ as defined in (<ref>).
Consider that we do this using some modification of randomized sketching methods.
Let the test and sketch matrices be tesselated according to the tree decomposition in Figure <ref> on Level 3,
so that

\begin{equation}\label{eq:dense_samples}
\mtx Y = \begin{pmatrix}\mtx Y_{8}\\ \mtx Y_{9}\\\vdots\\ \mtx Y_{15}\end{pmatrix},\qquad \mtx \Omega = \begin{pmatrix}\mtx \Omega_{8}\\ \mtx \Omega_{9}\\\vdots\\ \mtx \Omega_{15}\end{pmatrix}, \qquad \mtx Y_i \in \mathbb{R}^{m \times p},\ \mtx \Omega_i \in \mathbb{R}^{m \times p}
\end{equation}

and similar tesselations for $\mtx Z, \mtx \Psi$.

The structure of $\mtx A$ has dense interactions between near neighbor blocks, c.f. Figure <ref>,
which complicate the use of randomized sketching techniques.
Ideally, the structure of the
test matrices $\mtx \Omega$ would correspond to the sparsity pattern of the low rank blocks that we need to sample.
Consider a `structured' test matrix $\mtx \Omega' \in \mathbb{R}^{N \times k}$ for sampling the far-field of block $i = 9$, where
$\mtx \Omega'_{8}, \mtx \Omega'_{9}, \mtx \Omega'_{10} = \mtx 0$ and the other blocks are Gaussian random matrices.
Then sketching to produce $\mtx Y' = \mtx A \mtx \Omega'$, extracting the $I_9$ block, and post-processing would yield a basis for $\mtx Q_9$:

\begin{equation}\label{eq:orth_Omegaprime}
\underset{m \times k}{\mtx U_9}\ =\ \texttt{orth}(\underset{m \times k}{\mtx Y'_9}), \qquad \text{where}\
\underset{N \times k}{\mtx Y'}\ =\ \underset{N \times N}{\mtx A}\ \underset{N \times k}{\mtx \Omega'}
\end{equation}

To recover the basis for $i=12$, another structured test matrix $\mtx \Omega'' \in \mathbb{R}^{N \times k}$ can be designed
with zero blocks $\mtx \Omega''_{11}, \mtx \Omega''_{12}, \mtx \Omega''_{13} = \mtx 0$ and procedure in (<ref>) can be repeated
for the sketch $\mtx Y'' = \mtx A \mtx \Omega''$. Using these structured test matrices, tailored for sampling the far field of every box, would require $\frac k m N$ total samples
and $k^2 N$ post-processing cost.

Block nullification accomplishes the same aim with much fewer samples and slightly increased post-processing costs. Let us again consider sampling
the far-field of box $i=9$, and consider that we construct a matrix $\mtx N'$ from the fully dense test matrix $\mtx \Omega$ defined in (<ref>)

\begin{equation}
\underset{p \times (p-3m)} {\mtx N'} = \texttt{null} \begin{pmatrix} \underset{m \times p}{\mtx \Omega_{8}}\\ {\mtx \Omega_{9}}\\ {\mtx \Omega_{10}}\end{pmatrix}
\end{equation}

where the operation $\texttt{null}$ gives an orthogonal basis for the nullspace of matrix.
Suppose that $p=(3m+k)$, then $\mtx N'$ is $k$-dimensional with high probability. Multiplying $\mtx N'$ on the right of $\mtx \Omega$ gives

\begin{equation}
\underset{N \times p}{\mtx \Omega}\ \underset{ p \times k}{\mtx N'} =
\begin{pmatrix} \underset{m \times p} {\mtx \Omega_8}\ \underset{p \times k} {\mtx N'}\\
\mtx 0 \\ \mtx 0 \\ \mtx 0 \\ \mtx \Omega_{11} \mtx N'\\ \vdots \\ \mtx \Omega_{15} \mtx N'\end{pmatrix}\
\begin{pmatrix} \underset{m \times k} {\mtx \Omega'_8}\\
\mtx 0 \\ \mtx 0 \\ \mtx 0 \\ \mtx \Omega'_{11}\\ \vdots \\ \mtx \Omega'_{15}\end{pmatrix}:= \mtx \Omega' \in \mathbb{R}^{N \times k},
\end{equation}

where $\mtx \Omega'$ is the `structured' Gaussian random matrix discussed earlier.
The blocks $\mtx \Omega'_9, \mtx \Omega'_{10}, \mtx \Omega'_{11}$ are zero matrices by construction.
Likewise, the structured random matrix $\mtx \Omega''$ for sampling the far field for $i=12$
can be constructed in a similar fashion.

\begin{equation} \label{eq:null_12}
\underset{N \times p} {\mtx \Omega}\ \underset{p \times k}{ \mtx N''} := \underset{N \times k}{\mtx \Omega''}, \text{where}\ \underset{p \times (p-3m)} {\mtx N''} = \texttt{null} \begin{pmatrix} \mtx \Omega_{11}\\ \mtx \Omega_{12}\\ \mtx \Omega_{13}\end{pmatrix}
\end{equation}

To recover $\mtx Q_{9}$ as described in (<ref>), we need the $I_9$ block of the sketch $\mtx Y'$, which we can construct by post-processing
$\mtx \Omega$

\begin{equation}
\underset{N \times k}{\mtx Y'}\ =\ \underset{N \times N}{\mtx A} \underset{N \times k}{\mtx \Omega'}\ =\
\underset{N \times N} {\mtx A} \left( { \underset{N \times p} {\mtx \Omega}\ \underset{ p \times k}{\mtx N'}} \right)\ =\
\underset{N \times p}{\mtx Y}\ \underset{p \times k}{\mtx N'}.
\end{equation}

Note that because we only need the $I_9$  block, the full sketch $\mtx Y'$ does not need to be formed explicitly. Instead, $\mtx Y'_9$ can be constructed as

\begin{equation}
\underset{m \times k}{\mtx Y'_9}\ =\
\underset{m \times p}{\mtx Y_9}\ \underset{p \times k}{\mtx N'},\qquad \text{then}\ \mtx Q_9 = \texttt{orth} \left( {\mtx Y'_9} \right).
\end{equation}

Likewise, the basis $\mtx Q_{12}$ can be recovered by computing $\mtx N''$ as in (<ref>) and orthonormalizing the sketch
$\mtx Y''_{12} = \mtx Y_{12} \mtx N''$.
Block nullification allows us to recover all bases $\mtx Q_8, \dots, \mtx Q_{15}$ using only $p=(3m+k)$ samples of $\mtx A$
and $\bigO\left({(m^2 + m^2k) N } \right)$ post-processing cost.
The big-O notation hides constants related to the geometry which are discussed in a later section.
The bases $\mtx U_8, \dots, \mtx U_{15}$ can be recovered using a similar procedure for post-processing $\mtx Z$ and $\mtx \Psi$.

[fill=gray, thick] (,)̱ rectangle (+*,̱*+̱)̱;
in 1,...,
[thick,white,fill=darkgray] (+ *-̱ * +̱ ,̱ + * +̱ )̱ rectangle (+ *-̱ * ,̱ + * )̱;
in 2,...,
[thick,white,fill=darkgray] (+ *-̱ * +̱ 2*,̱+ * +̱ )̱ rectangle (+*-̱ * +̱ ,̱ +* )̱;
in 2,...,
[thick,white,fill=darkgray] (+*-̱ * +̱ ,̱ +* )̱ rectangle (+*-̱ * ,̱ +* -̱ )̱;

in 1,...,
[thick, white] (,+*)̱ – (+2*,+*)̱;
[fill=none, thick] (,+)̱ rectangle (+*,̱+*+̱)̱;

[fill=gray, thick] (+0,+)̱ rectangle (+*,̱+*+̱)̱;

in 1,...,
[thick, white,fill=darkgray] (+ *-̱ * +̱ ,̱ +* +̱ )̱ rectangle (+*-̱ * ,̱ +* )̱;
in 2,...,
[thick, white,fill=darkgray] (+*-̱ * +̱ 2*,̱ +* +̱ )̱ rectangle (+*-̱ * +̱ ,̱ +* )̱;
in 2,...,
[thick, white, fill=darkgray] (+*-̱ * +̱ ,̱ +* )̱ rectangle (+*-̱ * ,̱ +* -̱ )̱;

in 1,...,
[thick, white] (+ *,̱+)̱ – (+ * ,̱+2*+)̱;

[fill=none, thick] (+0,+)̱ rectangle (+*,̱+*+̱)̱;

Setting for block nullification.

[thick,fill=none] (,)̱ rectangle (+*,̱*+̱)̱;
in 1,...,
[thick,white,fill=darkgray] (+ *-̱ * +̱ ,̱ + * +̱ )̱ rectangle (+ *-̱ * ,̱ + * )̱;
in 2,...,
[thick,white,fill=darkgray] (+ *-̱ * +̱ 2*,̱+ * +̱ )̱ rectangle (+*-̱ * +̱ ,̱ +* )̱;
in 2,...,
[thick,white,fill=darkgray] (+*-̱ * +̱ ,̱ +* )̱ rectangle (+*-̱ * ,̱ +* -̱ )̱;

in 1,...,
[thick, white] (,+*)̱ – (+2*,+*)̱;
[fill=none, thick] (,+)̱ rectangle (+*,̱+*+̱)̱;

Setting for block extraction.

[fill=darkgray] (0,)̱ rectangle (,̱2*)̱;
at (4*,̱1.5*)̱ full rank;
[fill=gray] (0,3*)̱ rectangle (,̱ 4*)̱;
at (4*,̱3.5*)̱ low rank;

The figure depicts a rank-structured matrix $\mtx A$ in two settings.
Figure <ref> shows a dense matrix $\mtx A$ that uses low rank bases for each row and column block.
Blocks which are well-separated, according to the geometry in Figure <ref>,
are compressed using uniform bases, c.f. (<ref>) for a formal definition.
Figure <ref> shows a special case of a rank-structured matrix, where the well-separated
blocks have rank exactly 0.

§.§ Block Extraction

In addition to low-rank factors, the small dense blocks corresponding to interactions between adjacent neighbors need to be recovered as well.
For this task, we describe a randomized sampling technique which is useful for block-sparse matrices.
The setting for using this technique appropriately in the context of rank-structured matrices will be described in a later section.
For now, we consider a special case of a rank structured matrix $\mtx A$ when the off-diagonal rank is exactly 0.
Then, $\mtx A$ is block-sparse with the structure shown in Figure <ref>.
The technique known as block extraction recovers the dense blocks $\mtx A_{8,8},\mtx A_{8,9}, \mtx A_{9,8}$, etc. by post-processing
$\{\mtx Y,\mtx Z, \mtx \Omega, \mtx \Psi\}$ from samples (<ref>).

Because $\mtx A$ is sparse, all the blocks of $\mtx A$ can be recovered using only $3m$ samples by constructing $\mtx \Omega$ with special sparse structure
(e.g. carefully placed identity matrices).
For general geometries, this approach requires solving a combinatorial graph problem to find the sparsity pattern of $\mtx \Omega$
to recover $\mtx A$ from black-box samples.
Block extraction uses a different approach, which is very similar in spirit to block nullification.
Let us again design a `structured' test matrix $\mtx \Omega'$
to recover the blocks $\mtx A_{9,8}, \mtx A_{9,9}, \mtx A_{9,10}$ of the following form

\begin{equation}\label{eq:extract_sparse}
\mtx \Omega'\ =\
\begin{pmatrix}
\underset{m \times m}{\mtx I}\\
& \underset{m \times m}{\mtx I}\\
&& \underset{m \times m}{\mtx I}\\
& \underset{m \times 3m}{\mtx \Omega'_{11}}\\
& \underset{m \times 3m}{\mtx \Omega'_{15}}
\end{pmatrix}
\end{equation}

Note that the entries on blocks $i=11,\dots,15$ can be arbitrary because $\mtx A$ is block-sparse.
Then applying $\mtx \Omega'$ to form the sketch $\mtx Y'$ and extracting the $I_9$ block gives

\begin{equation}
\begin{pmatrix} \mtx A_{98,}& \mtx A_{9,9}& \mtx A_{9,10}\end{pmatrix}\ =\ \underset{m \times 3m} {\mtx Y'_{9}}, \qquad \text{where}\
\underset{N \times 3m}{\mtx Y'}\ =\ \underset{N \times N}{\mtx A}\ \underset{N \times 3m}{\mtx \Omega'}.
\end{equation}

Consider that we draw $p$ samples
using dense Gaussian random test matrices, as in (<ref>) and that we compute a matrix $\mtx M$ from the fully
dense test matrix $\mtx \Omega$ as

\begin{equation}\label{eq:pinv_omega}
\underset{p \times 3m}{\mtx M} = \begin{pmatrix} \underset{m \times p}{\mtx \Omega_8}\\ \mtx \Omega_9\\ \mtx \Omega_{10} \end{pmatrix}^{\dagger},
\end{equation}

where $\dagger$ denotes the right pseudo-inverse. With high probability, the Gaussian random matrix in (<ref>) has linearly independent rows.
Then applying $\mtx M$ to the right of Gaussian test matrix $\mtx \Omega$ gives

\begin{equation}
\underset{N \times p}{\mtx \Omega}\ \underset{p \times 3m}{ \mtx M}\ =\ \underset{N \times 3m} {\mtx \Omega'},
\end{equation}

where $\mtx \Omega'$ has the sparsity pattern of (<ref>). Note that $\mtx \Omega_{11}' = \mtx \Omega_{11} \mtx M$
is not a Gaussian random matrix any longer (e.g. $\mtx \Omega'_{11} \not \sim \mathcal N(\mtx 0, \mtx I)$).
This, however, does not matter for the purposes of extracting sparse blocks because $\mtx A_{9,11} = \mtx 0$.
To construct the sample $\mtx Y'_{9}$, we simply compute $\mtx M$ from $\mtx \Omega$ using (<ref>), then compute

\begin{equation}
\underset{m \times 3m}{\mtx Y'_9}\ =\ \underset{m \times p}{\mtx Y_9}\ \underset{p \times 3m}{\mtx M}, \qquad \text{where}\
\underset{N \times 3m}{\mtx Y'}\ =\ \underset{N \times N}{\mtx A}\ \underset{N \times 3m}{\mtx \Omega'}\ =\
\underset{N \times N}{\mtx A} \left({\underset{N \times p}{\mtx \Omega}\ \underset{p \times 3m}{\mtx M} } \right)
\end{equation}

Note that all of $\mtx Y'$ does not need to be computed, only the $I_9$ block $\mtx Y'_9$ is needed.
Block extraction  allows us to recover all sparse blocks of $\mtx A$ using only $p=(3m+l)$ samples, where $l$ is a small oversampling parameter,
and $\bigO\left({m^2 N} \right)$ post-processing cost.
Similar techniques are often useful for extracting sparse blocks by sampling the adjoint operator and post-processing $\mtx Z$ and $\mtx \Psi$.

§.§ Factorizing $\mtx A$ using Randomized Sampling Techniques

The objective of this work is to compute an invertible factorization of $\mtx A$.
The details of this machinery are discussed in subsequent sections.
For now, assume the factorization takes the following form

\begin{equation}\label{eq:randsample_fact}
\mtx A = \mtx V_{8} \dots \mtx V_{15}\ \widetilde{\mtx A}\ \mtx W_{15} \dots \mtx W_{8},
\end{equation}

where $\widetilde{\mtx A}$ is block-diagonal. The matrices $\mtx V_i, \mtx W_i$ for $i=8,\dots 15$
diagonalize the row interactions and the column interactions of the block $i$.
In the subsequent sections, we describe the operators $\mtx V_i$ and $\mtx W_i$, but for this discussion,
it is important to note that they are sparse and easy to invert.

Consider the following simpler decomposition of $\mtx A$, which only includes the diagonalization operators
for block $i=8$

\begin{equation}\label{eq:rand_recover15}
{\mtx A}\ =\ \mtx V_{8}\ \hat{\mtx A}\ \mtx W_{8},
\end{equation}

where $\hat{\mtx A}$ is defined in terms of the factorization of (<ref>).
In later sections, we will describe how $\mtx V_{8}$ and $\mtx W_{8}$ can be recovered using only
a small number of random sketches (<ref>)
of the operator $\mtx A$ and its adjoint by post-processing $\{{\mtx Y}, {\mtx \Omega}, {\mtx Z}, {\mtx \Psi}\}$.

Suppose that $\mtx V_{8}$ and $\mtx W_{8}$ have been recovered and that we would like to use similar techniques
to recover the operators $\mtx V_{9}$ and $\mtx W_{9}$. Then, it would be natural to compute random sketches
of $\hat{\mtx A}$

\begin{equation*}
\hat{\mtx Y} = \hat{\mtx A} \hat{\mtx \Omega},\quad \hat{\mtx Z} =\hat{\mtx A}^* \hat{\mtx \Psi}
\end{equation*}

and post-process $\{\hat{\mtx Y}, \hat{\mtx \Omega}, \hat{\mtx Z}, \hat{\mtx \Psi}\}$.
A key observation of this manuscript, is that it is not necessary to draw sketches of $\hat{\mtx A}$
because the matrix can be written in terms of $\mtx A$ as

\begin{equation*}
\hat{\mtx A} = \mtx V_8^{-1} \mtx A \mtx W_8^{-1}.
\end{equation*}

Instead, test and sketch matrices can be updated $\{{\mtx Y}, {\mtx \Omega}, {\mtx Z}, {\mtx \Psi}\} \Rightarrow\{\hat{\mtx Y}, \hat{\mtx \Omega}, \hat{\mtx Z}, \hat{\mtx \Psi}\}$
using the following formulas

\begin{equation}
\begin{split}
\hat{\mtx Y} = \mtx V_{8}\ \mtx Y,\qquad \hat{\mtx \Omega} = \mtx W_{8}^{-1} \mtx \Omega\\
\hat{\mtx Z} = \mtx W_{8}^*\ \mtx Z,\qquad \hat{\mtx \Psi} = \mtx V_{8}^{-*} \mtx \Psi\\
\end{split}
\end{equation}

This observation allows us to sample the operator and its adjoint as in (<ref>), then recover the factorization of $\mtx A$
by reusing the random sketches drawn initially.

§ THE INTERPOLATIVE DECOMPOSTION AND RECURSIVE SKELETONIZATION

In this section, we review preliminaries which are useful in understanding the strong recursive skeletonization algorithm.
First, we discuss Gaussian elimination and block elimination matrices in Section <ref>.
Then, we discuss interpolative decompositions (ID) for low rank compression in Section <ref>, which use a subset
of the rows or columns of the original matrix as a basis.
Section <ref> describes how the interpolative decomposition can be used to represent off-diagonal
blocks in a rank-structured $\mtx A$ to efficiently compute a direct solver $\mtx A^{-1}$.

§.§ Gaussian elimination and block elimination matrices

Consider a matrix of the form

\begin{equation*}
\mtx{A}
\left(\begin{array}{rrr}
\mtx{A}_{11} & \mtx{A}_{12} & \mtx{0}      \\
\mtx{A}_{21} & \mtx{A}_{22} & \mtx{A}_{23} \\
\mtx{0}      & \mtx{A}_{32} & \mtx{A}_{33}
\end{array}\right).
\end{equation*}

If $\mtx{A}_{11}$ is non-singular, we can “decouple” it from the other blocks via
one step of block Gaussian elimination. We express this mathematically through
a factorization

\begin{equation}
\label{eq:basicdecomp}
\mtx{A}
\left(\begin{array}{ccc}
\mtx{A}_{11} & \mtx{A}_{12} & \mtx{0}      \\
\mtx{A}_{21} & \mtx{A}_{22} & \mtx{A}_{23} \\
\mtx{0}      & \mtx{A}_{32} & \mtx{A}_{33}
\end{array}\right)
\mtx{L}
\left(\begin{array}{ccc}
\mtx{A}_{11} & \mtx{0}      & \mtx{0}      \\
\mtx{0}      & \mtx{S}_{22} & \mtx{A}_{23} \\
\mtx{0}      & \mtx{A}_{32} & \mtx{A}_{33}
\end{array}\right)
\mtx{U}
\end{equation}

where $\mtx{L}$ and $\mtx{U}$ are “block elimination matrices” of the form

\begin{equation}
\label{eq:defLU}
\mtx{L} =
\left(\begin{array}{ccc}
\mtx{I} & \mtx{0} & \mtx{0}  \\
\mtx{A}_{21}\mtx{A}_{11}^{-1} & \mtx{I} & \mtx{0}  \\
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right)
\qquad\mbox{and}\qquad
\mtx{U} =
\left(\begin{array}{ccc}
\mtx{I} & \mtx{A}_{11}^{-1}\mtx{A}_{12} & \mtx{0}  \\
\mtx{0} & \mtx{I} & \mtx{0}  \\
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right),
\end{equation}

and where $\mtx{S}_{22}$ is the
“Schur complement”

\begin{equation*}
\mtx{S}_{22} = \mtx{A}_{22} - \mtx{A}_{21}\mtx{A}_{11}^{-1}\mtx{A}_{12}.
\end{equation*}

For future reference, let us introduce a function  that builds the elimination matrices.
To be precise,

\begin{equation*}
\texttt{elim}(\mtx{B},I,J,n)
\end{equation*}

is the $n\times n$ identity matrix, except that the matrix $\mtx{B}$ has been inserted in the block identified by
the index vector $I$ and $J$. In other words, in (<ref>),

\begin{equation*}
\mtx{L} = \texttt{elim}(\mtx{A}_{21}\mtx{A}_{11}^{-1},I_{2},I_{1},n)
\qquad\mbox{and}\qquad
\mtx{U} = \texttt{elim}(\mtx{A}_{11}^{-1}\mtx{A}_{12},I_{1},I_{2},n).
\end{equation*}

Note that block-elimination matrices of the form $\texttt{elim}$ are simple to invert by toggling the sign of the off-diagonal block, e.g.
\begin{equation*}
\mtx L^{-1} = \texttt{elim}(-\mtx{A}_{21}\mtx{A}_{11}^{-1},I_{2},I_{1},n)
\end{equation*}

When $\mtx{A}$ is symmetric and positive definite (spd),
it is preferable to compute its Cholesky factorization.
We would then first factorize the $\mtx{A}_{11}$ block as

\begin{equation*}
\mtx{A}_{11} = \mtx{C}_{11}\mtx{C}_{11}^{*}
\end{equation*}

where $\mtx{C}_{11}$ is lower triangular. Then instead of (<ref>), we would form the triangular

\begin{equation*}
\left(\begin{array}{rrr}
\mtx{A}_{11} & \mtx{A}_{12} & \mtx{0}      \\
\mtx{A}_{21} & \mtx{A}_{22} & \mtx{A}_{23} \\
\mtx{0}      & \mtx{A}_{32} & \mtx{A}_{33}
\end{array}\right)
\left(\begin{array}{rrr}
\mtx{C}_{11} & \mtx{0} & \mtx{0}  \\
\mtx{A}_{21}\mtx{C}_{11}^{-*} & \mtx{I} & \mtx{0}  \\
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right)
\left(\begin{array}{rrr}
\mtx{I} & \mtx{0}      & \mtx{0}      \\
\mtx{0} & \mtx{S}_{22} & \mtx{A}_{23} \\
\mtx{0} & \mtx{A}_{32} & \mtx{A}_{33}
\end{array}\right)
\left(\begin{array}{rrr}
\mtx{C}_{11}^{*} & \mtx{C}_{11}^{-1}\mtx{A}_{12} & \mtx{0}  \\
\mtx{0} & \mtx{I} & \mtx{0}  \\
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right).
\end{equation*}

When $\mtx{A}$ is not spd, one could form analogous triangular factorizations using either the (partially pivoted) LU decomposition of $\mtx{A}_{11}$, or potentially an LDL factorization.
In this manuscript, we will focus on block diagonal factorizations, but all techniques described can easily be adapted to form triangular factorizations instead.

§.§ The interpolatory decomposition

Let $\mtx{B}$ be an $m\times n$ matrix of rank $k$.
The interpolatory decomposition (ID) of $\mtx{B}$ is a low rank factorization where a subset
of the columns/rows are used to span its column/row space.
For instance, in the column ID, we pick a set of $k$ linearly independent columns identified through
the “skeleton” index vector $J_{\rm s}$, and collect the remaining indices in the “residual” index vector $J_{\rm r}$,
and set

\begin{equation*}
J = [J_{\rm r},\,J_{\rm s}],
\qquad
\mtx{B}_{\rm r} = \mtx{B}(:,J_{\rm r}),
\qquad
\mtx{B}_{\rm s} = \mtx{B}(:,J_{\rm s}),
\end{equation*}

Since $\mtx{B}$ has rank $k$, there exists a matrix $\mtx{T}_{\rm sr}$ of size $k \times (n-k)$ such that

\begin{equation*}
\mtx{B}_{\rm r} =  \mtx{B}_{\rm s} \mtx{T}_{\rm sr}.
\end{equation*}

This allows us to factor the matrix $\mtx{B}$ as

\begin{equation}\label{eq:id_factorization}
\mtx{B}(:,J) =
\bigl(\mtx{B}_{\rm r},\,\mtx{B}_{\rm s}\bigr) =
\bigl(\mtx{B}_{\rm s}\mtx{T}_{\rm sr},\,\mtx{B}_{\rm s}\bigr) =
\bigl(\mtx{0},\,\mtx{B}_{\rm s}\bigr)
\left(\begin{array}{rr}
\mtx{I} & \mtx{0} \\
\mtx{T}_{\rm sr} & \mtx{I}
\end{array}\right).
\end{equation}

Interpolatory decompositions can also be computed of matrices that are only of approximate low rank.
The errors induced can in theory be significantly larger than those resulting from the optimal low rank decomposition obtained by truncating a singular value decomposition. However, in practice the error tends to be modest as long as the singular values of the input matrix decay at a decent rate.
For numerical stability, we would like the matrix
\mtx{I} & \mtx{0} \\
\mtx{T}_{\rm sr} & \mtx{I}
\end{array}\right)$
to be as well-conditioned as possible, but in practice tends to mean that we want the entries of $\mtx{T}_{\rm sr}$ to be small. It has been demonstrated that it is always possible to pick the set $J_{\rm s}$ so that every entry of $\mtx{T}_{\rm sr}$ has modulus bounded by one, and practical algorithms that ensure that the entries are of modest size are known.
See <cit.> for a detailed discussion of different algorithms for computing the ID in practice.

For completeness, let us also describe the row ID, where a set of linearly independent rows of $\mtx{B}$ would be identified by an index vector $I_{\rm s}$. Collecting the remaining row indices in the index vector $I_{\rm r}$, we then get the factorization

\mtx{B}([I_{\rm r},I_{\rm s}],:)
\left(\begin{array}{rr}
\mtx{I} & \mtx{T}_{\rm rs} \\
\mtx{0} & \mtx{I}
\end{array}\right)
\left(\begin{array}{r}
\mtx{0} \\ \mtx{B}_{\rm s}
\end{array}\right),

where $\mtx{B}_{\rm rs}$ is the $(m-k)\times k$ matrix such that
$\mtx{B}(I_{\rm r},:) = \mtx{T}_{\rm rs}\mtx{B}(I_{\rm s},:)$.

§.§ Classical skeletonization (weak admissibility)

We in this section describe how the interpolatory decomposition that we introduced in Section <ref> can be used to exploit rank deficiencies in off-diagonal blocks to map a given linear system to another “equivalent one” one with a smaller coefficient matrix. The technique described has previously been published as a fast direct solver for boundary integral equations, and is now widely known as “recursive skeletonization”
<cit.>. It is based on “weak admissibility”, and forms the basis for the more complex algorithm based on strong admissibility that is the main focus of this manuscript.

As a model problem, let us consider a matrix $\mtx{A}$ that has been partitioned into $2\times 2$ blocks as

\begin{equation}
\label{eq:Abf}
\mtx{A} =
\left(\begin{array}{r|r}
\mtx{A}_{\rm bb} & \mtx{A}_{\rm ab} \\ \hline
\mtx{A}_{\rm fb} & \mtx{A}_{\rm ff}
\end{array}\right),
\end{equation}

where the off-diagonal blocks $\mtx{A}_{\rm fb}$ and $\mtx{A}_{\rm bf}$ each have low rank. Let $k$ denote the rank of these blocks, and let $m$ and $n$ represent the dimensions of the blocks $\mtx{A}_{\rm bb}$ and $\mtx{A}_{\rm ff}$, respectively. We will construct block elimination matrices that map $\mtx{A}$ to a block diagonal matrix where one block has size $(m-k) \times (m-k)$, and other has size $(k+n)\times (k+n)$. Let $I_{\rm b}$ and $I_{\rm f}$ denote index vectors that identify the blocks, so that

(1:n) = I_{\rm b} \cup I_{\rm f}.

(For now, our discussion is one purely of linear algebra, but we use notation that will serve us well when applying these ideas to solve potential problems in Section <ref>.)
The first step of the factorization process is to form the row ID of $\mtx{A}_{\rm af}$ and the column ID of $\mtx{A}_{\rm fa}$ to obtain factorizations

\begin{equation}
\label{eq:af_rowskel}
\mtx{A}_{\rm bf}([I_{\rm r},I_{\rm s}],:) =
\left(\begin{array}{r} \mtx{A}_{\rm rf} \\ \mtx{A}_{\rm sf}\end{array}\right) =
\left(\begin{array}{r} \mtx{T}_{\rm rs}\mtx{A}_{\rm sf} \\ \mtx{A}_{\rm sf}\end{array}\right)
\end{equation}

\begin{equation}
\label{eq:fa_colskel}
\mtx{A}_{\rm fb}(:,[J_{\rm r},J_{\rm s}]) =
\bigl[ \mtx{A}_{\rm fr},\quad \mtx{A}_{\rm fs}\bigr] =
\bigl[ \mtx{A}_{\rm fs}\mtx{T}_{\rm sr},\quad \mtx{A}_{\rm fs}\bigr].
\end{equation}

Let us reorder the indices within the current box to have the residual indices go first:

I = [I_{\rm r},\,I_{\rm s},\,I_{\rm f}],
\qquad\mbox{and}\qquad
J = [J_{\rm r},\,J_{\rm s},\,I_{\rm f}].

This results in the tessellation

\begin{equation}
\label{eq:victess}
\mtx{A}(I,J)
\left(\begin{array}{cc|c}
\mtx{A}_{\rm rr} & \mtx{A}_{\rm rs} & \mtx{A}_{\rm rf} \\
\mtx{A}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sf} \\ \hline
\mtx{A}_{\rm fr} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm ff}
\end{array}\right)
\left(\begin{array}{cc|c}
\mtx{A}_{\rm rr} & \mtx{A}_{\rm rs} & \mtx{T}_{\rm rs}\mtx{A}_{\rm sf} \\
\mtx{A}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sf} \\ \hline
\mtx{A}_{\rm fs}\mtx{T}_{\rm sr} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm ff}
\end{array}\right).
\end{equation}

Next, we factorize the matrix to introduce zero blocks in the “fr” and “rf” locations:

\begin{equation}
\label{eq:factor1}
\mtx{A}(I,J)
\left(\begin{array}{cc|c}
\mtx{I} & \mtx{T}_{\rm rs} & \mtx{0} \\
\mtx{0} & \mtx{I} & \mtx{0} \\ \hline
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right)
\left(\begin{array}{cc|c}
\mtx{X}_{\rm rr} & \mtx{X}_{\rm rs} & \mtx{0} \\
\mtx{X}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sf} \\ \hline
\mtx{0} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm ff}
\end{array}\right)
\left(\begin{array}{cc|c}
\mtx{I} & \mtx{0} & \mtx{0} \\
\mtx{T}_{\rm sr} & \mtx{I} & \mtx{0} \\ \hline
\mtx{0} & \mtx{0} & \mtx{I}
\end{array}\right).
\end{equation}

Recognizing that the outside factors in (<ref>) are both block elimination matrices, we define

\mtx{E} := \texttt{elim}(\mtx{T}_{\rm rs},I_{\rm r},I_{\rm s},n)
\qquad\mbox{and}\qquad
\mtx{F} := \texttt{elim}(\mtx{T}_{\rm sr},J_{\rm s},J_{\rm r},n).

so that we can rewrite (<ref>) as

\begin{equation}
\label{eq:factor2}
\mtx{E}
\left(\begin{array}{rr|r}
\mtx{X}_{\rm rr} & \mtx{X}_{\rm rs} & \mtx{0} \\
\mtx{X}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sf} \\ \hline
\mtx{0} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm ff}
\end{array}\right)
\mtx{F}.
\end{equation}

We next perform one step of block Gaussian elimination to decouple $\mtx{X}_{\rm rr}$, so that

\mtx{A}(I,J)
\mtx{E}
\mtx{L}
\left(\begin{array}{rr|r}
\mtx{X}_{\rm rr} & \mtx{0} & \mtx{0} \\
\mtx{0} & \mtx{X}_{\rm ss} & \mtx{A}_{\rm sf} \\ \hline
\mtx{0} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm ff}
\end{array}\right)
\mtx{U}
\mtx{F},

where $\mtx{L}$ and $\mtx{U}$ are the block elimination matrices

\mtx{L} = \texttt{elim}(\mtx{X}_{\rm sr}\mtx{X}_{\rm rr}^{-1},I_{\rm s},I_{\rm r},n)
\qquad\mbox{and}\qquad
\mtx{U} = \texttt{elim}(\mtx{X}_{\rm rr}^{-1}\mtx{X}_{\rm rs},J_{\rm r},J_{\rm s},n).

§ STRONG RECURSIVE SKELETONIZATION

In this section, we review the strong recursive skeletonization (SRS) algorithm of <cit.>
that is used to compute an invertible factorization of a dense matrix involving a kernel matrix whose kernel represents physical interactions between a set of points in two or three dimensions. The decomposition is related to classical skeletonization as introduced in Section <ref>, however, the
skeleletonization is performed with “strong admissibility” condition, where only well-separated boxes are compressed as low-rank.

First, we discuss a hierarchical tree data structure which gives useful terminology for discussing the algorithm.
In Section <ref>, we discuss the strong recursive skeletonization procedure for a single box of the tree,
and in Section <ref>, we describe how the algorithm can be applied recursively in a traversal of the tree.
Section <ref> discusses how the factorization can be computed efficiently
in circumstances where matrix entries of $\mtx A$ can be accessed easily, setting the stage for computing the factorization in the black-box setting.

In the introduction of the ideas, we focus on the simple quasi-one dimensional domain of Figure <ref>,
though the ideas are directly applicable to more compilicated domains, as shown in Figure <ref>.

§.§ Hierarchical tree structure

SRS relies on a hierarchical partitioning of given points $\{\pxx_j\}_{j=1}^N$
into an quad-tree or oct-tree. For the purposes of this discussion, we restrict our attention to
uniform (e.g. fully populated) trees.
Formally, we define a tree $\tree$, in which
each box $b$ is associated with a subset of the given points.
Initally, all points belong to the root node.
Letting $d$ denote the dimension, we split the root node into $2^d$ children
nodes, and then split each child again until the size of each node is below
some given threshold $m$.
We refer to a node with children as a parent node,
and a node with no children as a leaf node.
The depth of a node is defined as its distance from the root node,
and level $\ell$ of the tree is defined as the set of nodes with depth $\ell$,
so that level 0 consists of only the root node,
level 1 consists of the $2^d$ children of the root node,
and so on.
The levels of the tree represent successively finer partitions of the points.
The depth of the tree is defined as the maximum node depth,
denoted by $L \approx \log_2 \left( N / m \right)$.
See Figure <ref> for a hierarchical decomposition of points into a uniform quadtree.

§.§ Strong skeletonization for a single box

Let us consider again the simple model problem on a quasi one dimensional domain in Figure <ref> that we introduced in Section <ref>,
where $\mtx A$ is given by the evaluation of (<ref>) on a set of points $\{\pxx_i\}_{i=1}^N$.
Consider the tesselating the given points of Figure <ref>(a) into boxes $8, \dots, 15$ as shown in Figure <ref>(b).
The objective, in this section, is to describe the process of diagonalizing the interactions between a target box b
and the rest of the points. Once this machinery is in place, the procedure can be applied successively to all boxes on a level.

A rectangular domain contains the points $\{\pxx_{j}\}_{j=1}^{N}$,
marked with blue.

A tessellation of points into boxes $8, \dots 15$.
The colors indicate the index vectors defined in Section <ref>, with
magenta for $I_{\rm a}$,
green for $I_{\rm n}$, and
blue for $I_{\rm f}$.

Consider the target box 8, for concreteness. First, we would like to compress the interactions between box b and its far field. To this end, we split the set of points into three sets: The box itself which contains all points in the box to be compressed (box 8 in this case), the “near field set” which contains points in boxes directly adjacent to the active box (box 9 in this case), and the “far field set” which contains everything else (boxes 10 through 15 in this case). These sets are identified through three index vectors:

$I_{\rm b}$        The “target box”. (While compressing box 8, $I_{\rm b} = I_{8}$.)

$I_{\rm n}$        The “near field indices”. (While compressing box 8, $I_{\rm n} = I_{9}$.)

$I_{\rm f}$        The “far field indices”. (While compressing box 8, $I_{\rm f} = [I_{10},\,I_{11},\,\dots,\,I_{15}]$.)

Figure <ref>(b) illustrates these definitions.
Setting $I = J = [I_{\rm b},\,I_{\rm n},\,I_{\rm f}]$, we now obtain a partition of
$\mtx{A}$ into the blocks

\begin{equation}
\label{eq:Abnf}
\mtx{A}(I,J)
\left(\begin{array}{r|rr}
\mtx{A}_{\rm bb} & \mtx{A}_{\rm bn}  & \mtx{A}_{\rm bf} \\ \hline
\mtx{A}_{\rm nb} & \mtx{A}_{\rm nn}  & \mtx{A}_{\rm nf} \\
\mtx{A}_{\rm fb} & \mtx{A}_{\rm fn}  & \mtx{A}_{\rm ff}
\end{array}\right).
\end{equation}

(Comparing (<ref>) to the tessellation (<ref>) that is used in classical skeletonization, we see that the near field has been added. The purpose is to reduce the numerical ranks that arise.)
We now invoke the well-known fact from the literature on fast summation schemes that all interactions between the active box and its far-field has low numerical rank. (To be precise, the rank is bounded by $O(\log(1/\varepsilon))$ where $\varepsilon$ is the requested precision. This bound is independent of how many points are actually in the box.) Linear algebraically, this means that $\mtx{A}_{\rm bf}$ and $\mtx{A}_{\rm fb}$ are both of numerically low rank. We form the row and the column IDs, respectively,  cf. (<ref>) and (<ref>),

\begin{equation}\label{eq:Abf}
\mtx{A}_{\rm bf}([I_{\rm r},I_{\rm s}],:) =
\left(\begin{array}{r} \mtx{A}_{\rm rf} \\ \mtx{A}_{\rm sf}\end{array}\right) =
\left(\begin{array}{r} \mtx{T}_{\rm rs}\mtx{A}_{\rm sf} \\ \mtx{A}_{\rm sf}\end{array}\right)
\end{equation}

\begin{equation}\label{eq:Afb}
\mtx{A}_{\rm fb}(:,[J_{\rm r},J_{\rm s}]) =
\bigl[ \mtx{A}_{\rm fr},\quad \mtx{A}_{\rm fs}\bigr] =
\bigl[ \mtx{A}_{\rm fs}\mtx{T}_{\rm sr},\quad \mtx{A}_{\rm fs}\bigr].
\end{equation}

These factorizations partition the nodes in $I_{\rm b}$ into a
set of “skeleton nodes” that represent all interactions between the
box and the outside world, and a set of “residual nodes” that we will
decouple from the rest of the system.
To keep track of this, we update $I$ and $J$ to put the residual nodes first:

I = [I_{\rm r},\,I_{\rm s},\,I_{\rm n},\,I_{\rm f}],
\qquad\mbox{and}\qquad
J = [J_{\rm r},\,J_{\rm s},\,I_{\rm n},\,I_{\rm f}].

The corresponding tessellation of $\mtx{A}$ becomes, cf. (<ref>),

\begin{equation}
\label{eq:A_fourtess}
\mtx{A}(I,J)
\left(\begin{array}{rr|rr}
\mtx{A}_{\rm rr} & \mtx{A}_{\rm rs} & \mtx{A}_{\rm rn}  & \mtx{A}_{\rm rf} \\
\mtx{A}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sn}  & \mtx{A}_{\rm sf} \\ \hline
\mtx{A}_{\rm nr} & \mtx{A}_{\rm ns} & \mtx{A}_{\rm nn}  & \mtx{A}_{\rm nf} \\
\mtx{A}_{\rm fr} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm fn}  & \mtx{A}_{\rm ff}
\end{array}\right)
\left(\begin{array}{rr|rr}
\mtx{A}_{\rm rr} & \mtx{A}_{\rm rs} & \mtx{A}_{\rm rn}  & \mtx{T}_{\rm rs}\mtx{A}_{\rm sf} \\
\mtx{A}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sn}  & \mtx{A}_{\rm sf} \\ \hline
\mtx{A}_{\rm nr} & \mtx{A}_{\rm ns} & \mtx{A}_{\rm nn}  & \mtx{A}_{\rm nf} \\
\mtx{A}_{\rm fs}\mtx{T}_{\rm sr} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm fn}  & \mtx{A}_{\rm ff}
\end{array}\right)
\end{equation}

Eliminating the “fr” and the “rf” blocks from (<ref>), we next get

\begin{equation}\label{eq:mod_fr}
\mtx{A}(I,J)
\mtx{E}
\left(\begin{array}{rr|rr}
\mtx{X}_{\rm rr} & \mtx{X}_{\rm rs} & \mtx{X}_{\rm rn}  & \mtx{0} \\
\mtx{X}_{\rm sr} & \mtx{A}_{\rm ss} & \mtx{A}_{\rm sn}  & \mtx{A}_{\rm sf} \\ \hline
\mtx{X}_{\rm nr} & \mtx{A}_{\rm ns} & \mtx{A}_{\rm nn}  & \mtx{A}_{\rm nf} \\
\mtx{0} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm fn}  & \mtx{A}_{\rm ff}
\end{array}\right)
\mtx{F}
\end{equation}

\mtx{E} = \texttt{elim}(\mtx{T}_{\rm rs},I_{\rm r},I_{\rm s},n)
\qquad\mbox{and}\qquad
\mtx{F} = \texttt{elim}(\mtx{T}_{\rm sr},J_{\rm s},J_{\rm r},n).

Then block Gaussian elimination yields the factorization

\begin{equation}
\label{eq:rr_eliminated}
\mtx{A}(I,J)
\mtx{E}
\mtx{L}{
\left(\begin{array}{rr|rr}
\mtx{X}_{\rm rr} & \mtx{0} & \mtx{0}  & \mtx{0} \\
\mtx{0} & \mtx{X}_{\rm ss} & \mtx{X}_{\rm sn}  & \mtx{A}_{\rm sf} \\ \hline
\mtx{0} & \mtx{X}_{\rm ns} & \mtx{X}_{\rm nn}  & \mtx{A}_{\rm nf} \\
\mtx{0} & \mtx{A}_{\rm fs} & \mtx{A}_{\rm fn}  & \mtx{A}_{\rm ff}
\end{array}\right)}
\mtx{U}
\mtx{F}
\end{equation}

where $\mtx{L}$ and $\mtx{U}$ are the block elimination matrices

\begin{equation}\label{eq:LU_elim_srs}
\mtx{L} = \texttt{elim}(\mtx{X}_{\rm sr}\mtx{X}_{\rm rr}^{-1},I_{\rm s},I_{\rm r},n)
\qquad\mbox{and}\qquad
\mtx{U} = \texttt{elim}(\mtx{X}_{\rm rr}^{-1}\mtx{X}_{\rm rs},J_{\rm r},J_{\rm s},n).
\end{equation}

The equation (<ref>) successfully decouples the residual nodes of target box b from the rest of the problem.
Alternatively, the system can be written as

\begin{equation}\label{eq:Zmat}
\widetilde{\mtx A} = \mtx V^{-1} \mtx A(I,J) \mtx W^{-1}, \qquad \mtx V:= \mtx E \mtx L, \qquad \mtx W:=\mtx U \mtx F
\end{equation}

where the matrices $\mtx V$ and $\mtx W$
are defined as the left and right diagonalization operators, respectively, and the matrix $\widetilde{\mtx A}$
is written explicitly in the equation (<ref>).
The redundant nodes of the target box are decoupled in equation (<ref>), and the algorithm proceeds by repeating a similar
procedure to the other boxes of the tree, but now for the modified system $\widetilde{\mtx A}$.

§.§ Recursive Algorithm

In the previous section, we described the procedure for decoupling the redundant nodes of a target box from the rest of the system.
Now, we describe how this procedure can be applied successively to all boxes on a level.
Consider again the simple model of a quasi-one dimensional geometry in Figure <ref>.
For target box 8, the left and right diagonalization operators $\mtx V_8$ and $\mtx W_8$ have the form described in Section <ref>.
Applying these operators to the appropriate permutation of $\mtx A$ yields the matrix $\widetilde{\mtx A}$ defined in (<ref>).
For the sparsity pattern of $\widetilde{\mtx A}$, see the leftmost panel of Figure <ref>, where the modified entries of $\widetilde{\mtx A}$
are highlighted in red.
The procedure can be repeated for target box $9$ by computing left and right diagonalization operators $\mtx V_9$ and $\mtx W_9$.
Note that instead of computing skeletons for the original matrix $\mtx A$ as in equation (<ref>),
we instead need to compute skeletons which sparsify the modified matrix $\widetilde{\mtx A}$.
Applying these techniques to the boxes $8, \dots, 15$ on Level 3 successively, we can compute a factorization, which (for an appropriate permutation)
has the following form

\begin{equation}\label{eq:single_level_decomp}
\begin{pmatrix}
\widetilde{\mtx A}_{\rm r_8}\\
& \widetilde{\mtx A}_{\rm r_9}\\
&& \ddots\\
&&& \widetilde{\mtx A}_{\rm r_{15}}\\
&&&&\widetilde{\mtx A}_{\rm[s_8, \dots, s_{15}]}
\end{pmatrix}
= \mtx V^{-1}_{15} \cdots \mtx{V}^{-1}_{8}\
\mtx{A}\
\mtx W^{-1}_{15} \cdots \mtx{W}^{-1}_{8},
\end{equation}

where $\widetilde{\mtx A}$ is a matrix of modified interactions.
The redundant nodes $I_{\rm r}$ for every box $8,\dots 15$ are decoupled from the rest of the system.
The remaining nodes $I_{\rm s}$ of skeleton nodes for every box interact densely in the dense block $\widetilde{\mtx A}_{\rm[s_8, \dots, s_{15}]}$.
In the process of computing the single-level factorization (<ref>), the modified matrix of decoupled interactions
need to be stored and updated. See Figure <ref> for
the sparsity pattern of decoupled interactions at various stages of computing the single-level factorization (<ref>).
Typically, strong recursive skeletonization is applied in circumstances where entries
of $\mtx A$ can be accessed in $\bigO(1)$ time, and the entire matrix of decoupled interactions does not need to be stored explicitly. Then, only the
modified entries shown in red of Figure <ref> need to be stored at every stage of the computation.

Since the residual degrees of freedom have been decoupled from the problem in (<ref>), the task that remains
is to solve the surviving linear system on the “skeleton nodes”. To illustrate the structure of this problem, we show the coefficient matrix after the residual nodes have been dropped in Figure <ref>(a).
For intermediate size problems, it is often possible to directly compute the LU decomposition of the linear system that connects the skeleton nodes, since the ranks tend to be small when strong admissibility it used. However, for large scale problems, there will be too many skeleton nodes surviving. In this case, we can continue the skeletonization process in a multilevel fashion. In order to reintroduce rank deficiencies in the off-diagonal blocks, we merge the boxes by twos to form larger boxes, as shown in <ref>(b). The idea is now to simply repeat the skeletonization process outlined in Section <ref>. This results in the further sparsified coefficient matrix shown in Figure <ref>(c).

After diagonalizing box 8
After diagonalizing box 11
After diagonalizing box 15

Illustration of the combined compression and
factorization process described in Section
Upper row:
The computational domain after steps 1, 4, and 8.
The nodes are colored so that
$I_{\rm s}$ is red,
$I_{\rm r}$ is cyan,
$I_{\rm n}$ is green, and
$I_{\rm f}$ is blue.
(The gray nodes mark points that were identified as “residual” nodes
in the previous steps. These points have dropped out of the computation
and play no active role.)
Lower row:
The sparsity pattern of the matrix $\widetilde{\mtx{A}}$ after the residual nodes in
Box 8/11/15 have been decoupled. White blocks are zeros,
red blocks are entries that got modified, and black blocks are entries that
have not been modified.

(a)        (b)        (c)

Illustration of the merging of boxes described in Section <ref>. (a) On top, the points remaining active after compression of the first level has completed. Below, the corresponding coefficient matrix, with entries that have been modified shown in red. (b) Boxes on the finest level have been merged by pairs to create larger boxes and reintroduce rank deficiencies in the off-diagonal blocks. (c) The geometry and the coefficient matrix after compression at the coarser level has completed.

We will shortly formalize the description of the recursive skeletonization
process, but let us first briefly show what changes occur when the domain is “truly” two dimensional, as in Figure <ref>(a). It is now natural to organize the domain into a quadtree of boxes, rather than a binary tree, as shown in Figures <ref>(b) and <ref>(c). Other than that, the scheme proceeds just as it does for the quasi one dimensional example we discussed earlier. After all boxes at the finest level have been compressed, the remaining skeleton points are shown on top in Figure <ref>(a) in red, while the residual points are shown in gray. The matrix of interactions between the skeleton points is shown below. Observe that there are now many more updated (red) boxes than there were for the quasi one dimensional domain (as many as 25 in some rows). In order to reintroduce rank deficiencies in the off-diagonal blocks, we now merge boxes by sets of four to yield level 2, as illustrated in Figure <ref>(b). Once compression has completed at level 2, very few points remain as skeleton nodes, as we see on top in Figure <ref>(c), with the very sparse coefficient matrix shown below.

Computational domain
Ordering of boxes on level 3
Ordering of boxes on level 2

The square computational domain described in Section <ref>. (a) The original domain, with the points $\{\pxx_{i}\}_{i=1}^{N}$ shown in blue. (b) The boxes at the finest level. (c) The merged boxes at the next coarser level.

Active points after level 3:

Merge boxes by fours:

Active points after level 2:

(a)        (b)        (c)

Matrix after level 3:

Matrix retessellated:

Matrix after level 2:

(b)        (d)        (f)

The analog of Figure <ref> when the computational domain is the box shown in Figure <ref>. Observe that many more blocks get updated for a true two-dimensional domain — up to 25 boxes per row need to be explicitly stored.

To formalize the description of SRS, consider that the boxes are diagonalized in an upward traversal through the tree in order
$1, \dots, n$. Then the decomposition takes the form

\begin{equation}\label{eq:multi_level_decomp}
\widetilde{\mtx A}_{\rm diag} := \begin{pmatrix}
\widetilde{\mtx A}_{\rm r_1}\\
& \widetilde{\mtx A}_{\rm r_2}\\
&& \ddots\\
&&& \widetilde{\mtx A}_{\rm r_{n}}\\
&&&&\widetilde{\mtx A}_{\rm[S]}
\end{pmatrix}
= \mtx V^{-1}_{1} \cdots \mtx{V}^{-1}_{n}\
\mtx{A}\
\mtx W^{-1}_{n} \cdots \mtx{W}^{-1}_{1},
\end{equation}

where $[\rm S]$ is the set of all remaining skeleton nodes at the time of terminating the algorithm (typically on level 2 of the tree).
Because $\widetilde{\mtx A}_{\rm diag}$ is block-diagonal, it is easy to invert. The factorization (<ref>)
can also be used for a sparse factoriation of $\mtx A$.
In the next section, we discuss techniques which allow for the efficient calculation of elimination matrices $\mtx E$ and $\mtx F$ of Section <ref>
in the setting of boundary integral equations.

§.§ How to compress the far-field interactions

The aim is to compress the two off-diagonal blocks $\mtx{A}_{\rm b \, f}$ and $\mtx{A}_{\rm f,\ b}$ using the column ID.
Notice that the computational cost would be $\bigO(N)$ if the full matrix is formed, which turns out to be unnecessary,
particularly when $\mtx A$ arises from the discretization of a boundary integral equation.

The proxy method introduces points $\left\{\pxx^{(\rm proxy)}_i \right\}_{i=1}^{n_p}$ which replicate the effect of far-field points $I_{\rm f}$. Instead of using the matrix $\mtx A_{\rm b,f}$, a smaller matrix with $\bigO(1)$ rows and columns is formed and compressed

\begin{equation}  \label{e:proxy}
\mtx{A}_{\rm proxy,\rm b} \qquad \text{where}
\left(\mtx A_{\rm proxy, \rm b} \right)_{i,j} =
k {\left( \pxx^{(proxy)}_i,  \pxx[I_{\rm b_j}] \right)}.
\end{equation}

With the appropriate choice of proxy surface, the matrix $\mtx A_{\rm proxy, b}$ spans the row space of $\mtx A_{f,b}$, and likewise, $\mtx A_{b,\rm proxy}$ spans the column space of $\mtx A_{b,f}$ for any point distribution in the far-field. Therefore, the indices $I_{b}=I_{r} \cup I_{s}$ and interpolative matrix $\mtx T$ computed by forming and factorizing the matrix (<ref>) will satisfy (<ref>).
Because the matrix $\widetilde{\mtx A}$ of diagonalized interactions does modify entries of $\mtx A$, the proxy surface needs to be placed appropriately
to compute the SRS factorization; this is discussed in further detail in <cit.>.

For a box b, skeleton indices $I_{\rm s} \subseteq I_{\rm b}$ (shown in red) can be computed efficiently using
a proxy surface $\pxx_{\rm proxy}$, which are fictitious points replicating the effect of the
far-field points $I_{\rm f}$.

§ RANDOMIZED STRONG RECURSIVE SKELETONIZATION

In this section, we describe computing the SRS factorization described in Section <ref> in the black-box setting
where the matrix $\mtx A$ is only accessed through its action on vectors.
The algorithm relies on the randomized sampling machinery first introduced in Section <ref>.
Recall that our setting is the following

\begin{equation} \label{eq:rand_sample_setting_again}
\underset{N \times p}{\mtx Y}\ =\ \underset{N \times N}{\mtx A}\ \underset{N \times p}{\mtx \Omega}, \qquad
\underset{N \times p}{\mtx Z}\ =\ \underset{N \times N}{\mtx A^*}\ \underset{N \times p} {\mtx \Psi}, \qquad \mtx \Omega, \mtx \Psi \sim \mathcal{N}(\mtx 0, \mtx I)
\end{equation}

where the matrix $\mtx A$ and its adjoint are sketched with Gaussian random matrices.

As in our description of SRS, we first describe the process of diagonalizing the interactions of a target box b with respect to the rest
of the problem. Recall the notation used in Section <ref> for a target box $I_{\rm b}$, its near-field indices $I_{\rm n}$
and its far-field indices $I_{\rm f}$.

In order to compute the skeletonization operators $\mtx E$ and $\mtx F$, we need to compute the column ID of the sketches

\begin{equation}\label{eq:sketch_far}
\underset{m \times k}{\mtx Y'_{\rm b}} = \mtx A_{\rm bf}\ \mtx \Omega_{\rm f},\qquad \underset{m \times k}{\mtx Z'_{\rm b}} = \mtx A_{\rm fb}^*\ \mtx \Psi_{\rm f},
\end{equation}

respectively. The samples $\mtx Y$ and $\mtx Z$ in (<ref>), however, include the contributions of the near field blocks $I_{\rm n}$.
The contribution of $\mtx \Omega_{\rm n}$ and $\mtx \Psi_{\rm n}$, however, can be “nullified”
using block nullification as described in Section <ref> and sketches (<ref>) can be computed as

\begin{equation}\label{eq:sketch_blocknull}
\underset{m \times k}{\mtx Y'_{\rm b}}\ =\ \underset{m \times p}{\mtx Y_b}\ \underset{p \times k} {\texttt{null} \left(\mtx \Omega_{\rm n} \right)},\qquad
\underset{m \times k}{\mtx Z'_{\rm b}}\ =\ \underset{m \times p}{\mtx Z_{\rm b}}\ \underset{p \times k}{\texttt{null}\left(\mtx \Psi_{\rm n} \right)},
\end{equation}

The sparsified matrix $\hat{\mtx A} = \mtx E^{-1} \mtx A \mtx F^{-1}$ has decoupled interactions between $I_{\rm r}$ and $I_{\rm f}$, and
the interactions between $I_{\rm r}$ and the rest of the system can be entirely decoupled using block-elimination matrices $\mtx L, \mtx U$
as defined in <ref>. Ideally, we would like to extract the modified interactions between $\mtx X_{\rm r, \rm n}$ and $\mtx X_{\rm n,r}$
as well as $\mtx X_{\rm r, s}$ and $\mtx X_{\rm s,r}$ in equation (<ref>).
The samples $\{\mtx Y, \mtx \Omega, \mtx Z, \mtx \Psi\} \Rightarrow \{\hat{\mtx Y}, \hat{\mtx \Omega}, \hat{\mtx Z}, \hat{\mtx \Psi}\}$
can be updated to instead be samples  of $\hat{\mtx A}$.

\begin{equation}
\begin{split}
\hat{\mtx Y} = \mtx E\ \mtx Y,\qquad \hat{\mtx \Omega} = \mtx F^{-1} \mtx \Omega\\
\hat{\mtx Z} = \mtx F^* \mtx Z,\qquad \hat{\mtx \Psi} = \mtx E^{-*} \mtx \Psi\\
\end{split}
\end{equation}

Then the modified interactions $\mtx X_{\rm r,n}$ and $\mtx X_{\rm n,r}$ can be extracted using block extraction, as discussed in Section <ref>
using the formulas

\begin{equation}\label{eq:sketch_blocknull}
{\mtx X_{\rm r,n }}\ =\ {\hat{\mtx Y}_{\rm r}}\ \hat{\mtx \Omega}_{\rm n}^{\dagger},\qquad
{\mtx X_{\rm n, r }^*}\ =\ {\hat{\mtx Z}_{\rm r}}\ \hat{\mtx \Psi}_{\rm n}^{\dagger}.
\end{equation}

After computing the block elimination matrices $\mtx L$ and $\mtx U$, the random sketches can again to updated to be instead
of the matrix $\widetilde{\mtx A} = \mtx L^{-1} \hat{\mtx A} \mtx U^{-1}$ by updating the sketches $\{\hat{\mtx Y}, \hat{\mtx \Omega}, \hat{\mtx Z}, \hat{\mtx \Psi}\}$
to produce sketches $\{\widetilde{\mtx Y}, \widetilde{\mtx \Omega}, \widetilde{\mtx Z}, \widetilde{\mtx \Psi}\}$.

\begin{equation}
\begin{split}
\widetilde{\mtx Y} = \mtx L\ \hat{\mtx Y},\qquad \widetilde{\mtx \Omega} = \mtx U^{-1} \hat{\mtx \Omega}\\
\widetilde{\mtx Z} = \mtx U^* \hat{\mtx Z},\qquad \widetilde{\mtx \Psi} = \mtx L^{-*} \hat{\mtx \Psi}\\
\end{split}
\end{equation}

After completing this process for a single box, the procedure can be repeated for a sequence of boxes.
Over the coarse of the algorithm, the samples do not need to be redrawn, instead they are updated using the computed factorization at every stage.
Because the algorithm is multi-level, compression errors do propogate from one level to the next. To handle this issue,
the absolute tolerance for the compression stage is successively relaxed at every level of the computation.
The matrix $\widetilde{\mtx A}$ in equation (<ref>) is block-diagonal to some accuracy which may deviate slighly from the
desired compression tolerance. To estimate the extent to which our compression is successful for a particular level,
we use randomized sampling techniques (e.g. block nullification) to estimate the extent to which $\widetilde{\mtx A}_{\rm r,f}$ and $\widetilde{\mtx A}_{\rm f,r}$
deviate from desired compression tolerance, where in this case $I_{\rm f} = (1:N) \setminus I_r$, for every box.

§ NUMERICAL EXPERIMENTS

In this section, we illustrate the performance of the proposed algorithm.
For the example in Section <ref>, we consider a sparse Schur complement, which arises
in the context of sparse direct solvers. Sparse Schur complements are a composition of sparse
matrix operations; they can be applied fast to vectors, however, matrix entries are challenging to access efficiently.

$N$       number of points

$m$       leaf size of tree

$T_{\rm factorize}$     time for Algorithm <ref>

$M$     memory needed for $\mtx V_i$ and $\mtx W_i$ for all boxes
atol     absolute compression tolerance at the leaf level

relerr     relative maximum error (on a subset of points)

errsolve     relative maximum error (on a subset of points)

Notation for reported numerical results.

The RSRS algorithm is operates on a uniform quad-tree $\mathcal T$ which is partitioned to have at most $m$ points in
the leaf boxes. The user provides a parameter atol which dictates the absolute tolerance of compression
at the leaf level; the compression rank $k$ is chosen for each box adaptively.
The success of the
algorithm is measured by accessing the relative error of the computed factorization, as well as
the error in the computed inverse. For the computed invertible factorization $\mtx K\approx\mtx A$, we report

\begin{equation}\label{eq:error}
\rm{relerr} = \frac{\|\mtx A - \mtx K \|_2}{\| \mtx A \|_2}, \qquad \rm{errsolve} = \| \mtx I - \mtx K^{-1} \mtx A\|_2.
\end{equation}

Typical choices of user parameter atol may be $10^{-5}$ or even as large as $10^{-2}$; the resulting errors (<ref>)
in the factorization are often much smaller, as we will demonstrate in the numerical results.
Table <ref> summarizes the notation used to report the numerical results.

§.§ 3D Sparse Direct Solvers

Consider a boundary value problem of the form

\begin{equation}
\label{eq:bvp}
\left\{\begin{aligned}
\mathcal -\Delta u(x) =&\ f(x),\qquad&x \in \Omega,\\
u(x) =&\ g(x),\qquad&x \in \Gamma,
\end{aligned}\right.
\end{equation}

where $\mathcal A$ is a second order elliptic differential operator,
and $\Omega$ is a rectangular domain in three dimensions with boundary $\Gamma$.
Upon discretizing (<ref>) with 2nd order finite differences, one obtains a linear system

\begin{equation*}
\label{eq:Au=f}
\mtx{A}\mtx{u} = \mtx{f},
\end{equation*}

involving a sparse coefficient matrix.
Sparse direct solvers produce an invertible factorization of $\mtx A$ by leveraging sparsity, as well as $\mathcal H$-matrix structure when appropriate
Often, it is useful to decompose the domain for the purposes of parallelizing the computation.
The typical choice of domain decomposition is a hierarchical quad-tree or oct-tree, but a decomposition into thin slab subdomains is also
possible and advantageous for ease of parallelization.
This domain decomposition is employed in SlabLU, which is a simplified two-level sparse direct solver <cit.>.
The corresponding domain decomposition is shown in Figure <ref>.

Domain decomposition used in SlabLU. The even-numbered nodes correspond to the nodes interior to each subdomain. The odd-numbered nodes correspond to interfaces between slabs. The slab partitioning is chosen so that interactions between slab interiors are zero. The slabs have width of $b$ points.

To be precise, SlabLU uses a decomposition of the domain into elongated “slab” subdomains of width $b$, where the odd-numbered nodes correspond to slab interface nodes and the even-numbered nodes correspond to slab interior nodes.
With the slab decomposition, the linear system (<ref>) has the block form

\begin{equation}
\label{eq:dell1}
\begin{bmatrix}
\mtx { A}_{11} & \mtx {A}_{12} & \mtx 0 & \mtx 0 & \mtx 0 & \dots \\
\mtx { A}_{21} & \mtx {A}_{22} & \mtx {A}_{23} & \mtx 0 &  \mtx 0 & \dots \\
\mtx 0 & \mtx {A}_{32} & \mtx {A}_{33} & \mtx {A}_{34} &   \mtx 0 & \dots\\
\mtx{0} & \mtx 0 & \mtx {A}_{43} & \mtx {A}_{44} & \mtx {A}_{45} & \dots\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots
\end{bmatrix}
\begin{bmatrix} \mtx {u}_1\\\mtx {u}_2\\ \mtx {u}_3 \\ \mtx {u}_4 \\ \vdots \end{bmatrix} =
\begin{bmatrix} \mtx {f}_1\\ \mtx {f}_2\\ \mtx {f}_3\\ \mtx {f}_4 \\ \vdots        \end{bmatrix}.
\end{equation}

The system (<ref>) can be factorized by first eliminating the nodes internal to each slab are eliminated in parallel by computing sparse direct factorizations $\mtx A_{22}^{-1}, \mtx A_{44}^{-1}, \dots$ in parallel.
The remaining system has sparse Schur complements of the form

\begin{equation}
\label{eq:T11}
\mtx{T}_{11} = \mtx A_{11} - \mtx{A}_{ 12}\ \mtx{A}_{ 22}^{-1}\ \mtx{A}_{ 21}
\end{equation}

Owing to sparsity, $\mtx T_{11}$ and its transpose can be rapidly applied to random vectors

\begin{equation}
\underset{n^2 \times s}{\mtx Y} = \underset{n^2 \times n^2}{\mtx T_{11}}\ \underset{n^2 \times s}{\mtx \Omega}, \qquad \underset{n^2 \times s}{ \mtx Z} = \underset{n^2 \times n^2}{\mtx T_{11}^*}\ \underset{n^2 \times s}{\mtx \Psi}
\end{equation}

Figure <ref> reports experiments compressing and factorizing $\mtx T_{11}$ using the RSRS algorithm which only
accesses the matrix through its action on vectors.
The PDE is Poisson, and the slab thickness is fixed as $b=10$ for various problem sizes $N = n^2$.

Time to compute $\mtx T_{11}^{-1}$ scales linearly.

Memory to store $\mtx T_{11}^{-1}$ scales linearly.

Number of samples needed is independent of $N$.

The computed factorization is accurate and does not deteriorate with increasing $N$.

The figures above show the effectiveness of using RSRS on a sparse Schur complement $\mtx T_{11}$ which arises
while computing sparse direct solvers for discretized PDEs. The number of samples needed to construct a factorization $\mtx T_{11}^{-1}$
is constant as $N$ increases and instead depends only on the desired user tolerance as well as the leaf size $m$.
For atol=5e-2 and 5e-4, the leaf size $m$ is set as $m=200$. For atol=5e-6, the leaf size $m$ is set as $m=350$.

§ CONCLUSIONS

The manuscript introduces the algorithm
Randomized Strong Recursive Skeletonization (RSRS)
for simultaneously compressing
and inverting an $\mathcal{H}^{2}$-matrix, given a means of applying the
matrix and its adjoint to vectors.
RSRS is immediately applicable in a range of important environments, and
the numerical results demonstrate the effectiveness of the approach on Schur complements which arise
while factorizing a sparse direct solver.
It extends the range of the SRS algorithm using tailored approaches for randomized linear algebra applied to rank structured

§.§ Acknowledgments

The work reported was supported by the Office of Naval Research (N00014-18-1-2354),
by the National Science Foundation (DMS-1952735, DMS-2012606, and DMS-2313434),
and by the Department of Energy ASCR (DE-SC0022251).