# Block-encoding dense and full-rank kernels using hierarchical matrices:
applications in quantum numerical linear algebra

Quynh T. Nguyen  Department of Electrical Engineering and Computer Science,
Massachusetts Institute of Technology, USA  Department of Physics,
Massachusetts Institute of Technology, USA  Bobak T. Kiani  Department of
Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, USA  Research Laboratory of Electronics, Massachusetts Institute
of Technology, USA  Seth Lloyd  Research Laboratory of Electronics,
Massachusetts Institute of Technology, USA  Department of Mechanical
Engineering, Massachusetts Institute of Technology, USA  Turing Inc.,
Cambridge, MA, USA

###### Abstract

Many quantum algorithms for numerical linear algebra assume black-box access
to a block-encoding of the matrix of interest, which is a strong assumption
when the matrix is not sparse. Kernel matrices, which arise from discretizing
a kernel function $k(x,x^{\prime})$, have a variety of applications in
mathematics and engineering. They are generally dense and full-rank.
Classically, the celebrated fast multipole method performs matrix
multiplication on kernel matrices of dimension $N$ in time almost linear in
$N$ by using the linear algebraic framework of hierarchical matrices. In light
of this success, we propose a block-encoding scheme of the hierarchical matrix
structure on a quantum computer. When applied to many physical kernel
matrices, our method can improve the runtime of solving quantum linear systems
of dimension $N$ to $O(\kappa\operatorname{polylog}(\frac{N}{\varepsilon}))$,
where $\kappa$ and $\varepsilon$ are the condition number and error bound of
the matrix operation. This runtime is near-optimal and, in terms of $N$,
exponentially improves over prior quantum linear systems algorithms in the
case of dense and full-rank kernel matrices. We discuss possible applications
of our methodology in solving integral equations and accelerating computations
in N-body problems.

##  1 Introduction

One of the most promising applications of quantum information processing is in
solving linear algebraic problem in high dimensions efficiently. Since the
Harrow-Hassidim-Lloyd (HHL) algorithm was first developed for solving linear
systems of equations [1], a number of improvements have been proposed.
Significantly, the breakthrough algorithms in [2] and [3] solved sparse linear
systems of equations in time polylogarithmic in the system size $N$ and the
target error bound $\varepsilon$. When these algorithms are applied to dense
matrices, however, their complexity scales linearly in $N$ (see Table 1).
Using a quantum random access memory (QRAM) data structure [4], the work of
[5] improved the runtime to $\widetilde{O}(\sqrt{N})$. Classically, quantum-
inspired algorithms [6, 7] use an analog of QRAM called sample-query access
and apply sketching techniques [8] to solve linear systems with runtimes
independent of $N$ but scaling poorly with the rank $k$ of the matrix as
$\widetilde{O}(k^{6})$ (see Appendix A).

In fact, achieving a runtime logarithmic in $N$ for dense and full-rank
matrices is generally challenging unless there are specific symmetries or
structures inherent in the matrix. For example, quantum algorithms have been
developed to achieve polylogarithmic complexity in $N$ for Toepliz systems
[9], Hankel matrices [10], and linear group convolutions [11] which all
feature the group quantum Fourier transforms (QFT) [12] to implement
operations as sparse matrices in the Fourier basis. This approach of using
group QFTs, however, only applies for certain groups and requires the matrix
entries to be _exactly_ uniformly sampled from a generating function, thus
inapplicable to more general settings (e.g., matrix entries corresponding to
atoms on a lattice vibrating around their equilibrium positions).

Kernel matrices are a class of matrices whose entries are calculated as a
function of a kernel $k(x,x^{\prime})$. They have applications in mathematical
physics [13], engineering [14] and machine learning [15], where they are often
obtained by sampling or discretizing a smooth kernel $k(x,x^{\prime})$ that
reflects the physical model at hand. In general, kernel matrices are dense and
full-rank. Hierarchical matrices ($\mathcal{H}$-matrices) [16, 17] are a
classical framework for efficiently performing matrix operations on commonly
found kernel matrices by hierarchically splitting the “space" dimension of the
matrix entries. In terms of the matrix dimension $N$, classically performing
matrix-vector multiplication with $\mathcal{H}$-matrices takes time
$O(N\operatorname{polylog}N)$ which is a significant improvement over the
generic $O(N^{2})$ time. This $\mathcal{H}$-matrix structure lies at the heart
of the celebrated fast multipole method [13]. This improvement intuitively
arises because kernel matrices often decay or change slowly along entries that
are farther from the diagonals of the matrix. $\mathcal{H}$-matrices
progressively approximate the matrix in larger low-rank blocks for farther
off-diagonal entries to efficiently implement kernel matrices.

We present an efficient quantum algorithm to implement $\mathcal{H}$-matrices
in the block-encoding framework, which embeds matrices inside a unitary
operator that can be queried multiple times. This construction, which is at
the core of many quantum algorithms, allows one to conveniently perform matrix
computations on a quantum computer [18, 19, 20, 21]. As many recently proposed
algorithms use block-encodings as a black box, one of the major challenges for
quantum numerical linear algebra in this framework is to efficiently block-
encode matrices so that the asymptotic runtimes are efficient even when the
block-encoding complexity is included. Our $\mathcal{H}$-matrix block-encoding
procedure, combined with existing quantum linear systems algorithms, _e.g.,_
adiabatic solver in [22], provides near-optimal runtime in solving quantum
linear systems for many typical kernels that are neither sparse nor low-rank,
e.g., $k(x,x^{\prime})=\|x-x^{\prime}\|^{-1}$ (see Table 1 for comparisons to
prior work). Apart from the $\mathcal{H}$-matrix block-encoding, we also
provide in the Appendix a list of common kernel matrices that can be
efficiently block encoded using existing techniques.

We proceed as follows. First, we overview classical hierarchical matrices in
Section 2 and the quantum block-encoding framework in Section 3. We then
present an algorithm for efficiently block-encoding hierarchical matrices and
apply this to kernel matrices in Section 4. In Section 4.4, we generalize our
methodology to construct block encodings for a more general class of matrices
whose structure is derived from hierarchical matrices. Finally, we discuss
applications of our algorithm in calculating $N$-body gravitational potentials
and solving integral equations in Section 5.

| Complexity |  |
---|---|---|---
Result |  Forward $\widetilde{O}(\cdot)$ |  Inversion $\widetilde{O}(\cdot)$ | Input model | Algorithm
Classical [16] | $N$ | $N$ | Classical |  Classical $\mathcal{H}$-matrix
GST [6]/SM [7] | - | $\frac{N^{6}\kappa^{8}}{\varepsilon^{2}}$ | Sample query | Classical sketching
CKS [2]/GSLW [3] | $N\kappa$ | $N\kappa$ |  Oracle accessa |  Lemma 5.2b
KP [4]/WZP [5] | $\sqrt{N}\kappa$ | $\sqrt{N}\kappa$ |  QRAMc | Lemma 5.2
This paper | $\operatorname{polylog}(N)\kappa$ | $\operatorname{polylog}(N)\kappa$ | Oracle access |  Lemma 4.3, Lemma 5.2

a Naive block-encoding using oracle access via Lemma 48 or 49 in [3] (also see
Lemma 3.2)
b Lemma 5.2 uses the adiabatic solver in [22] which has optimal runtime
dependence on $\kappa$
c This input model is strictly more powerful than the oracle access model

Table 1: Runtimes of previous methods versus our method for matrix
multiplication (forward) and solving linear systems (inversion) of dense and
full-rank kernel matrices, e.g., $k(x,x^{\prime})=\|x-x^{\prime}\|^{-1}$. For
purposes of comparison, we assume the operator norm of the matrix is bounded.
The $\widetilde{O}$ notation hides terms that grow as
$\operatorname{polylog}(\frac{N}{\varepsilon})$. Here, $\kappa$ is the
condition number of the matrix operation and $\varepsilon$ is the error bound
on the output. For quantum algorithms, we assume the input vector is given as
a quantum state and the runtime is equal to the circuit depth times the number
of queries needed to obtain the output as a quantum state. We note that
classical $\mathcal{H}$-matrix literature and this work focus on kernel
matrices, whereas the quantum algorithms we compare to are “general-purpose”.
For these general-purpose algorithms, runtimes are calculated assuming they
are applied directly on kernel matrices. More details about prior works can be
found in Appendix A.

##  2 Background in $\mathcal{H}$-matrices and hierarchical splitting

To observe why hierarchical matrices are useful, consider a physical model
which requires summation or integration over spacetime. The $N$-body
gravitational force problem [23], for example, requires evaluating $N$
potentials of the form

$\Phi\left(\mathbf{x}_{j}\right)=\sum_{i=1\atop i\neq
j}^{N}\frac{m_{i}}{\|\mathbf{x}_{j}-\mathbf{x}_{i}\|^{2}},\hskip
14.22636pt1\leq j\leq N$ (1)

where $m_{i}$ and $\mathbf{x}_{i}$ denote the mass and location of particle
$i$. Summing or integrating over all particles is equivalent to performing
matrix-vector multiplication, generally requiring $O(N^{2})$ operations unless
the underlying matrix is sparse or has special properties. As we will show,
when particles are placed in a favorable arrangement (e.g., approximately
uniformly), this matrix-vector multiplication can be approximately performed
in time $O(N\operatorname{polylog}N)$.

Hierarchical matrices ($\mathcal{H}$-matrices) [16, 17, 24] are a powerful
algebraic framework for accelerating matrix computations involving
displacement kernels
$k(\mathbf{x},\mathbf{x}^{\prime})=k(\|\mathbf{x}-\mathbf{x}^{\prime}\|)$
which satisfy certain properties outlined below. The key idea is to (i) split
the kernel matrix into hierarchically off-diagonal blocks and (ii) approximate
each block by a low-rank matrix. This low rank approximation is justified
because off-diagonal blocks represent distant interactions that are small or
slowly changing with distance. We now detail the two ingredients of
$\mathcal{H}$-matrices, namely the hierarchical splitting and the off-diagonal
low-rank approximation.

Figure 1: (A) Structure of a hierarchical splitting for a 1D array of 16
particles (gray circles), which are close to uniformly distributed. The
accumulative potential at the destination location (star, which can also be
one of the sources) is decomposed into contributions from clusters of source
particles whose sizes grow with the distance to the destination. (B) Cluster
size decreases exponentially with the level of the hierarchical splitting. (C)
The associated kernel matrix of this 16-particle system takes the form of an
$\mathcal{H}$-matrix, assuming that the destination locations are the same as
the source locations. This matrix is decomposed into hierarchically structured
blocks using the hierarchical splitting in (A), where each block is
numerically low-rank.

###  2.1 Hierarchical splitting into admissible blocks

Hierarchical splitting partitions a kernel matrix into so-called “admissible”
blocks where entries representing more distant interactions are grouped
together into larger blocks. Consider a system of $N$ particles, whose
pairwise interaction is described by the kernel $k(x,x^{\prime})$. The kernel
matrix of this system is defined as

$\mathbf{K}=(k(x_{i},x_{j}))_{i,j=0}^{N-1},$ (2)

where $x_{i}$ is the location of particle $i$. For ease of analysis, we assume
the particles are uniformly distributed and choose $x_{i}=i/N,\text{ for any
}i\in\\{0,1,\ldots,N-1\\}$ and $N=2^{L}$. For example, this is the case when
the kernel matrix is obtained by discretizing an integral operator over the
domain $[0,1]$ (see [16]). In general, a kernel matrix need not be square and
the source locations $x_{i}$ need not be 1-dimensional nor equally spaced and
we will show how to generalize the current setup to these cases in later
sections. We refer the interested readers to [25, 26, 17, 27, 28] for further
details.

A set of $m$ particles located at $\\{x_{j_{1}},\ldots,x_{j_{m}}\\}$ is
denoted as a _cluster_ $\sigma$, storing the list of individual indices
$\sigma=\\{j_{1},\ldots,j_{m}\\}$.

######  Definition 2.1 (Admissible blocks).

Let $\mathbf{K}=(k(x_{i},x_{j}))_{i,j=0}^{N-1}$ and
$\mathcal{N}=\\{0,1,\ldots,N-1\\}$. Let $\mathcal{P}(\mathcal{N})$ be a
partition of $\mathcal{N}$, i.e.,
$\mathcal{N}=\bigcup_{\sigma\in\mathcal{P}(\mathcal{N})}\sigma$. For a cluster
$\sigma\in\mathcal{P}(\mathcal{N})$, define the center of $\sigma$ as
$c_{\sigma}=\operatorname{mean}\\{x_{i}|i\in\sigma\\}$ and the radius of
$\sigma$ as $r_{\sigma}=\max_{i\in\sigma}|x_{i}-c_{\sigma}|.$ Furthermore, for
$\sigma,\rho\in\mathcal{P}(\mathcal{N})$, define the distance between the two
clusters $\sigma$ and $\rho$ as
$\operatorname{dist}(\sigma,\rho)=\min_{i\in\sigma,j\in\rho}|x_{i}-x_{j}|.$
Then, the matrix block
$\mathbf{K}^{\sigma,\rho}=(k(x_{i},x_{j}))_{i\in\sigma,j\in\rho}$ is called
admissible if
$\eta\cdot\max\\{r_{\sigma},r_{\rho}\\}\leq\operatorname{dist}(\sigma,\rho)$,
where $\eta=2$.

Informally, a block is admissible if the distance between its two
corresponding clusters is larger than their radii (other values of $\eta$ can
give rise to different forms of the hierarchical splitting, see [25, 16]). We
now split the matrix $\mathbf{K}$ into a hierarchy of admissible blocks
consisting of $L=\log_{2}N$ levels of decreasing radii. At level $\ell$,
define the following partition of $\mathcal{N}=\\{0,1,\ldots,N-1\\}$:

$\mathcal{P}^{(\ell)}=\\{\sigma^{(\ell)}_{I}|I\in\\{0,\ldots,2^{\ell}-1\\}\\},$
(3)

where

$\sigma^{(\ell)}_{I}=\\{i\in\mathcal{N}|x_{i}\in[I/2^{\ell},(I+1)/2^{\ell})\\}.$
(4)

For example, $\mathcal{P}^{(0)}=\\{\mathcal{N}\\}$,
$\mathcal{P}^{(L)}=\\{\\{i\\}|i\in\mathcal{N}\\}$, and $\mathcal{P}^{(\ell)}$
for $\ell\in\\{2,3,4\\}$ are shown in Figure 1B. Note that admissible blocks
only exist for $\ell\geq 2$. We now construct a hierarchical splitting of our
kernel matrix $\mathbf{K}$ into admissible blocks starting from the coarsest
level $\ell=2$. First, we split $\mathbf{K}$ into admissible and inadmissible
blocks in $\mathcal{P}^{(2)}$

$\mathbf{K}=\mathbf{K}^{(2)}+\mathbf{K}^{(2)}_{inadmissible}.$

Next, we split $\mathbf{K}^{(2)}_{inadmissible}$ into admissible and
inadmissible blocks in $\mathcal{P}^{(3)}$ (excluding the admissible blocks of
$\mathcal{P}^{(3)}$ that are already covered in $\mathbf{K}^{(2)}$)

$\mathbf{K}^{(2)}_{inadmissible}=\mathbf{K}^{(3)}+\mathbf{K}^{(3)}_{inadmissible}.$

Recursively applying the decomposition
$\mathbf{K}^{(\ell-1)}_{inadmissible}=\mathbf{K}^{(\ell)}+\mathbf{K}^{(\ell)}_{inadmissible}$
up to level $L$, we obtain a hierarchically structured representation of
$\mathbf{K}$

$\mathbf{K}=\sum_{\ell=2}^{L}\mathbf{K}^{(\ell)}+\mathbf{K}_{ad},$ (5)

where $\mathbf{K}_{ad}=\mathbf{K}^{(L)}_{inadmissible}$ is a tridiagonal
matrix which represents the “adjacent" interaction of neighboring particles
(see Figure 1C). Observe that the matrices $\mathbf{K}^{(\ell)}$ are block-
sparse; that is, each block-row or block-column of $\mathbf{K}^{(\ell)}$
contains at most three non-zero admissible blocks. Furthermore, these blocks
represent “far-field” interactions which can be approximated by low-rank
matrices if the kernel $k(x,x^{\prime})$ satisfies some additional properties
discussed in the next section.

In Section D.3, we describe a generalization of the hierarchical splitting to
higher dimensions, where the two key properties still hold, i.e., there are
only $\log N$ levels and each level $\mathbf{K}^{(\ell)}$ is block-sparse.

###  2.2 Low-rank approximation of admissible blocks

For far-field interactions in admissible blocks, the kernel $k(x,x^{\prime})$
can be well approximated by a truncated Taylor series if it satisfies the
following _asymptotic smoothness_ condition:

$\displaystyle|\partial_{x^{\prime}}^{q}k(x,x^{\prime})|$ $\displaystyle\leq
Cq!|x-x^{\prime}|^{-q}\hskip 5.69046pt(\forall q\in\mathbb{N}),$ (6)

where $C$ is a constant. For instance, the log kernel
$k(x,x^{\prime})=\log(|x-x^{\prime}|)$ directly satisfies the above. This
condition enables approximating any admissible block
$\mathbf{K}^{\sigma,\rho}$ by a rank-$p$ matrix

$\mathbf{K}^{\sigma,\rho}\approx\mathbf{\widetilde{K}}^{\sigma,\rho}=\bm{\Psi}^{\sigma,\rho}\mathbf{D}\left(\bm{\Phi}^{\rho}\right)^{\dagger},$
(7)

where $\mathbf{D}=\operatorname{Diag}(1/q!)_{q\in P}\in\mathbb{R}^{p\times
p}$, $P=\\{0,1,\ldots,p-1\\}$ and

$\displaystyle\bm{\Psi}^{\sigma,\rho}=(\partial_{x^{\prime}}^{q}k(x_{i},c_{\rho}))_{i\in\sigma,q\in
P}\in\mathbb{R}^{|\sigma|\times p},$ (8)
$\displaystyle\bm{\Phi}^{\rho}=((x_{j}^{\prime}-c_{\rho})^{q})_{j\in\rho,q\in
P}\in\mathbb{R}^{|\rho|\times p}.$ (9)

Here, $p$ regulates the precision of the approximation (typically
$p=\operatorname{polylog}(N/\varepsilon)$ for an error bound $\varepsilon$). A
detailed example and error analysis of this far-field approximation can be
found in Appendix D. In addition, we note that conditions other than that of
Equation 6 can also be used to justify the approximation (see [29, 30, 31]) of
admissible blocks by low-rank kernels
$k(x,x^{\prime})=\sum_{q=0}^{p-1}\psi_{q}(x)\varphi_{q}(x^{\prime})$.

Combining this low-rank approximation and the hierarchical splitting forms a
complete $\mathcal{H}$-matrix. In Appendix D, we show that matrix-vector
multiplication with an $\mathcal{H}$-matrix can be performed in time $O(Np\log
N)$. In addition, other $\mathcal{H}$-matrix operations such as matrix-matrix
multiplication and matrix inversion run in time $O(Np^{2}\log^{2}N)$, we refer
the reader to [25, 16, 24] for further details. The rank $p$ can be chosen at
each level to further optimize computational complexity versus approximation
error [26].

##  3 Block-encoding arithmetic

Throughout this paper, we employ the block-encoding framework to perform
matrix computations on a quantum computer [18]. Here, an $s$-qubit (non-
unitary) operator $A\in\mathbb{C}^{2^{s}\times 2^{s}}$ (let $N=2^{s}$) is
embedded in a unitary operator $U\in\mathbb{C}^{2^{(s+a)}\times 2^{(s+a)}}$
using $a$ ancilla qubits such that the top left block of $U$ is proportional
to $A$:

$U=\begin{pmatrix}A/\alpha&\cdot\\\ \cdot&\cdot\end{pmatrix},$ (10)

where $\alpha$ is the _normalization factor_ of the block-encoding. Formally,
the $(s+a)$-qubit unitary $U$ is an $(\alpha,a,\varepsilon)$-block-encoding of
$A$ if

$\left\|A-\alpha(\bra{0^{a}}\otimes I_{s})U(\ket{0^{a}}\otimes
I_{s})\right\|\leq\varepsilon,$ (11)

where $\|\cdot\|$ denotes the operator norm of a matrix and $I_{s}$ is the
identity operator on $s$ qubits. In other words, applying the unitary $U$ to a
quantum state $\ket{\psi}$ and post-selecting on the measurement outcome
$\ket{0^{a}}$ on the ancilla qubits is equivalent to applying the operation
$A$ on $\ket{\psi}$:

$U\ket{0^{a}}\ket{\psi}=\frac{1}{\alpha}\ket{0^{a}}A\ket{\psi}+\ket{G},$ (12)

where $\ket{G}$ is a garbage state that is orthogonal to the subspace
$\ket{0^{a}}$, i.e., $(\bra{0^{a}}\otimes I_{s})\ket{G}=0$. The probability of
successfully post-selecting $\ket{0^{a}}$ is thus equal to
$\|A\ket{\psi}\|^{2}\alpha^{-2}\leq(\|A\|/\alpha)^{2}$. Observe that
$\alpha\geq\|A\|$ is required to embed $A$ inside a unitary matrix. Further,
since the probability of success of the post-selection step depends on the
_ratio_ $\|A\|/\alpha$, a block-encoding $U$ is _optimal_ if
$\alpha=\Theta(\|A\|)$. More generally, if the ratio $\alpha/\|A\|$ grows
logarithmically with the size of the matrix $N$, the probability of success
also only changes logarithmically with $N$. More generally, as we will see
some examples in Section 5, the complexity of quantum algorithms [3] in the
block-encoding framework depends on the ratio $\alpha/\|A\|$. This important
observation motivates the following definition of a “good” block-encoding.

######  Definition 3.1 (Good block-encoding).

Given a matrix $A\in\mathbb{C}^{N\times N}$, an
$(\alpha,a,\varepsilon)$-block-encoding $U$ of $A$ is good if and only if the
ratio $\alpha/\|A\|=O(\operatorname{polylog}(N))$, and additionally,
$a=O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ and the circuit depth of
$U$ is $O(\operatorname{polylog}(\frac{N}{\varepsilon}))$. The block-encoding
is called _optimal_ if $\alpha=\Theta(\|A\|)$.

Throughout this paper, we assume that matrices are normalized to have the
largest entry at most $1$, and we analyze the optimality of block-encodings by
comparing the normalization factor $\alpha$ to the operator norm $\|A\|$.

Before proceeding to block-encode a hierarchical matrix, we provide various
methods for block-encoding the individual matrix blocks in the hierarchical
splitting. Depending on the locations of the matrix blocks and how the values
of the kernel entries can be accessed by a quantum computer, there are various
different ways one can block-encode a given matrix. Previous work has proposed
efficient means to block-encode sparse matrices, purified density matrices,
etc., as well as transformations of these matrices in the block-encoding
framework [3, 32, 18]. We list some previously established results in this
framework in Appendix B, among which we particularly use Lemma B.5 to block-
encode sparse matrices with oracle access to entries, Lemma B.3 to linearly
combine block-encoded matrices, and Lemma B.4 to multiply block-encoded
matrices. Later, we will apply and combine these individual block-encoding
Lemmas to form a block-encoding of the complete hierarchical matrix.

First, we provide a version of Lemma B.5 applied to dense matrices, which we
refer to as the “naive” block-encoding approach since it is unaware of the
inherent structure of the block-encoded matrix.

######  Lemma 3.2 (Naive block-encoding of dense matrices with oracle access).

Let $A\in\mathbb{C}^{N\times N}$ (where $N=2^{s}$) and let
$\hat{a}\geq\max_{i,j}|a_{ij}|$. Suppose the following oracle is provided

$\mathcal{O}_{A}:\ket{i}\ket{j}\ket{0^{b}}\rightarrow\ket{i}\ket{j}\ket{\tilde{a}_{ij}},$

where $0\leq i,j<N$ and $\tilde{a}_{ij}$ is the (exact) $b$-qubit description
of $a_{ij}/\hat{a}$. Then one can implement a
$(N\hat{a},s+1,\varepsilon)$-block-encoding of $A$ with two uses of
$\mathcal{O}_{A}$, $O(\operatorname{polylog}(\frac{\hat{a}N}{\varepsilon}))$
one- and two-qubit gates and
$O(b,\operatorname{polylog}(\frac{\hat{a}N}{\varepsilon}))$ extra qubits
(which are discarded before the post-selection step).

Assuming $\hat{a}=1$, the normalization factor of this block-encoding is $N$,
which is suboptimal if $\|A\|\ll N$, e.g., this situation arises when the
matrix $A$ has many small entries. This suboptimal block-encoding could lead
to an exponentially small probability of success in the post-selection step.
Therefore, a structure-aware block-encoding procedure is needed to achieve an
optimal runtime in terms of $N$.

To directly implement the $\mathcal{H}$-matrices outlined in Section 2, we
need procedures to block-encode block-sparse and low-rank matrices which we
provide in the Lemmas below (all proofs are deferred to Appendix B). First we
define state preparation pairs, which are frequently used in this study to
prepare the coefficients in a linear combination of matrices.

######  Definition 3.3 (State preparation pair [3]).

Let $y\in\mathbb{C}^{m}$ and $\|y\|_{1}\leq\beta$, the pair of unitaries
$(P_{L},P_{R})$ is called a $(\beta,n,\varepsilon)$-state-preparation-pair of
$y$ if $P_{L}\ket{0^{n}}=\sum_{j=0}^{2^{n}-1}c_{j}\ket{j}$ and
$P_{R}\ket{0^{n}}=\sum_{j=1}^{2^{n}-1}d_{j}\ket{j}$ such that
$\sum_{j=0}^{m-1}\left|\beta c_{j}^{*}d_{j}-y_{j}\right|\leq\varepsilon_{1}$
and $c_{j}^{*}d_{j}=0$ for any $j\in\\{m,\ldots,2^{n}-1\\}$.

Next, we show how to implement a block-encoding of low-rank matrices.

######  Lemma 3.4 (Block-encoding of low-rank operators with state preparation
unitaries, inspired by Lemma 1 of [33]).

Let
$A=\sum_{i=0}^{p-1}\sigma_{i}\mathbf{u}_{i}\mathbf{v}_{i}^{\dagger}\in\mathbb{C}^{2^{s}\times
2^{s}}$ where $\|\mathbf{u}_{i}\|=\|\mathbf{v}_{i}\|=1$. Let $r=\lceil\log
p\rceil$ and $\sum_{i=0}^{p-1}|\sigma_{i}|\leq\beta$. Suppose the following
$(r+s)$-qubit unitaries are provided:

$\displaystyle G_{L}:$
$\displaystyle\ket{i}\ket{0^{s}}\rightarrow\ket{i}\ket{\mathbf{u}_{i}},$
$\displaystyle G_{R}:$
$\displaystyle\ket{i}\ket{0^{s}}\rightarrow\ket{i}\ket{\mathbf{v}_{i}},$

where $0\leq i<p$ and $\ket{\mathbf{u}_{i}}$ $(\ket{\mathbf{v}_{i}})$ is the
quantum state whose amplitudes are the entries of $\mathbf{u}_{i}$
$(\mathbf{v}_{i})$. Let $(P_{L},P_{R})$ be a $(\beta,r,\varepsilon)$-state-
preparation-pair for the vector $[\sigma_{0},\ldots,\sigma_{p-1}]$. Then, one
can construct a $(\beta,r+s,\varepsilon)$-block-encoding of $A$ with one use
of each of $G_{L},G_{R}^{\dagger},P_{L}^{\dagger},P_{R}$, and a
$\operatorname{SWAP}$ gate on two $s$-qubit registers.

The above lemma provides one method to block-encode the approximately low-rank
matrices described in Section 2, but as we will see, is not the only option
one has at their disposal. When the rank of a matrix is small and bounded, the
time needed to prepare unitaries $P_{L},P_{R}$ is typically negligible; e.g.,
they can be efficiently implemented with an oracle providing access to
$\sigma_{i}$ (and $\beta$ is close to the sum of the singular values). In
addition, the unitaries $G_{R},G_{L}$ can be constructed using existing state
preparation methods based on the structure of the entries in the singular
vectors (e.g., [34] utilized efficiently integrable distributions, [4] used a
QRAM data structure). In Section B.3, we also present an approach to
efficiently construct $G_{R},G_{L}$ from oracles providing access to the
entries of $\mathbf{u}_{i},\mathbf{v}_{i}$ for cases where the entries change
slowly (Lemma B.8).

Next, the following lemma block-encodes a block-sparse matrix, whose blocks
are also block-encoded.

######  Lemma 3.5 (Block-encoding of block-sparse matrices).

Let $A=\sum_{i,j=0}^{2^{t}-1}\ket{i}\bra{j}\otimes A^{ij}$ be a $d_{r}$-row-
block-sparse and $d_{c}$-column-block-sparse matrix, where each $A^{ij}$ is an
$s$-qubit operator. Let $U^{ij}$ be an $(\alpha_{ij},a,\varepsilon)$-block-
encoding of $A^{ij}$. Suppose that we have the access to the following
$2(t+1)$-qubit oracles

$\displaystyle\mathcal{O}_{r}:$
$\displaystyle\ket{i}\ket{k}\rightarrow\ket{i}\ket{r_{ik}},$
$\displaystyle\mathcal{O}_{c}:$
$\displaystyle\ket{l}\ket{j}\rightarrow\ket{c_{lj}}\ket{j},$

where $0\leq i<2^{t},0\leq k<d_{r}$, $0\leq j<2^{t},0\leq l<d_{c}$, $r_{ik}$
is the index for the $k$-th non-zero block of the $i$-th block-row of A, or if
there are less than $k$ non-zero blocks, then $r_{ik}=k+2^{t}$, and similarly
$c_{lj}$ is the index for the $l$-th non-zero block of the $j$-th block-column
of $A$, or if there are less than $l$ non-zero blocks, then $c_{lj}=l+2^{t}$.
Additionally, suppose the following oracle is provided:

$\displaystyle\mathcal{O}_{\alpha}:\ket{i}\ket{j}\ket{z}\rightarrow\ket{i}\ket{j}\ket{z\oplus\tilde{\alpha}_{ij}}$

where $0\leq i,j<2^{t}$ and $\tilde{\alpha}_{ij}$ is the $b$-qubit description
of $\alpha_{ij}$ (if $i$ or $j$ is out of range then $\tilde{\alpha}_{ij}=0$),
and let the $(2t+2+s+a)$-qubit unitary $W=\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j})\otimes U^{ij}+\left(I-\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j})\right)\otimes I_{s+a}$. Let
$\hat{\alpha}=\max_{i,j}\alpha_{ij}$. Then one can implement an
$(\hat{\alpha}\sqrt{d_{r}d_{c}},t+a+3,2\sqrt{d_{r}d_{c}}\varepsilon)$-block-
encoding of $A$, with one use of each of $\mathcal{O}_{r},\mathcal{O}_{c},$
and $W$, two uses of $\mathcal{O}_{\alpha}$,
$O(t+\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ additional
one- and two-qubit gates, and
$O(b,\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ added ancilla
qubits.

######  Remark 3.6.

Lemma 3.5 is efficient when the matrix is sparse over blocks (where each block
may be dense), which is not in general the same as the matrix being sparse. It
generalizes Lemma 48 of [3] (also see Lemma B.5) for block-encoding matrices
with oracle access to entries. In addition, if $\alpha_{ij}$ are the same for
all blocks, we do not need the oracle $\mathcal{O}_{\alpha}$, hence some
ancilla qubits and gates can be saved (see proof in Section B.4).

##  4 Block-encoding of kernel matrices via hierarchical splitting

In this section, we apply the block-encoding techniques delineated in Section
3 to block-encode $\mathcal{H}$-matrices. Our goal is to obtain good block-
encodings (Definition 3.1) of these matrices where the ratio between the
normalization factor and the operator norm is at most
$O(\operatorname{polylog}(N))$. We provide an optimal block-encoding procedure
for kernel matrices arising from a variety of polynomially decaying and
exponentially decaying kernels. We show that using only the hierarchical
splitting _without subsequent low-rank approximation_ gives an optimal block-
encoding for these kernel matrices, enabling matrix transformations (e.g.,
multiplication, inversion, polynomial transformations) with optimal query and
gate complexities. Then, we provide a block-encoding procedure for general
$\mathcal{H}$-matrices whose admissible blocks are approximated by low-rank
matrices and discuss the applicability of our procedure on other variants of
hierarchical splitting.

###  4.1 Polynomially decaying kernels

We use the hierarchical splitting described in Section 2 to implement an
optimal block-encoding procedure for polynomially decaying kernels, which take
the following form

$k(x,x^{\prime})=\begin{cases}C&\text{ if }x=x^{\prime}\\\
\|x-x^{\prime}\|^{-p}&\text{ otherwise}\end{cases},$ (13)

where $p>0$ and $|C|\leq 1$ represents self-interaction terms (e.g., $C=0$ in
Coulomb point source interactions). For concreteness, we consider the problem
of computing pairwise interactions in a system of $N$ particles. Particularly,
consider a 1D system of $N$ particles $\\{(x_{i},m_{i})\\}_{i=0}^{N-1}$, where
$x_{i}$ and $m_{i}$ are the location and mass of the particle $i$,
respectively. The potential at $x_{i}$ is the sum over the pairwise
interactions

$\Phi(x_{i})=\sum_{j=0}^{N-1}m_{j}k(x_{i},x_{j}).$ (14)

For simplicity and ease of analysis, we consider a system of equally distanced
particles in which $x_{j}=j$ and $m_{j}=1$ for any $j\in[N]$, where $N=2^{L}$.
Note that we choose the domain to be $[1,N]$ instead of $[0,1]$ as considered
in Section 2, so that the maximum entry in the kernel matrix is $1$. Had the
particles been placed in the region $[0,1]$, we could simply downscale the
matrix by dividing it by $N^{p}$ so that the entries are bounded by $1$.
Recall from Section 3 that doing so does not change the block-encoding and
algorithmic runtime as these depend on the ratio between the operator norm and
the normalization factor. This ratio is unchanged under global rescaling of
the matrix because the normalization factor will rescale accordingly in the
block-encoding procedure.

To access the kernel function on a quantum computer, we assume the following
oracle is given

$\mathcal{O}_{k}:\ket{i}\ket{j}\ket{0}^{\otimes
b}\rightarrow\ket{i}\ket{j}\ket{\tilde{k}(x_{i},x_{j})},$ (15)

where $\tilde{k}(x_{i},x_{j})$ is the $b$-qubit description of
$k(x_{i},x_{j})$. In fact, if the function $k(x_{i},x_{j})$ can be computed
via an efficient classical circuit, one can construct a comparably efficient
quantum circuit for $\mathcal{O}_{k}$ [35]. In particular, this is the case
when the distribution of particles is known and efficiently computable (e.g.
$x_{j}$ are equally spaced points in the discretization of integral
operators).

We optimally block-encode the kernel matrix
$\mathbf{K}=(k(x_{i},x_{j}))_{i,j=0}^{N-1}$ up to error $\varepsilon$ with a
normalization factor of $\Theta(\|\mathbf{K}\|)$ by using only two queries to
$\mathcal{O}_{k}$ and $O(\operatorname{polylog}(N,\frac{1}{\varepsilon}))$
additional resources (i.e., one- and two-qubit gates and ancilla qubits). At a
high level, our method proceeds as follows. First, we decompose $\mathbf{K}$
into a hierarchical structure of admissible blocks as described in Section 2

$\mathbf{K}=\sum_{\ell=2}^{L}\mathbf{K}^{(\ell)}+\mathbf{K}_{ad}.$ (16)

Note that this hierarchical decomposition can also be performed in higher
dimensional systems (see Section D.3). Next, for each admissible block (which
is a dense matrix) in the level $\mathbf{K}^{(\ell)}$, we block-encode it by
using the naive procedure for dense matrices (Lemma 3.2). Then, we use the
procedure for block-sparse matrices (Lemma 3.5) to form a block-encoding of
$\mathbf{K}^{(\ell)}$. Finally, we sum over hierarchical levels to obtain a
block-encoding of $\mathbf{K}$.

Specifically, at level $\ell$ of the hierarchy, each admissible block has
dimension $2^{L-\ell}$. In addition, the minimum pairwise distance between
particles in an admissible block in $\mathbf{K}^{(\ell)}$ is
$d_{\min}^{(\ell)}=2^{L-\ell}$, so the maximum entry of an admissible block
$\mathbf{K}^{\sigma,\rho}$ in $\mathbf{K}^{(\ell)}$ is bounded by
$k^{(\ell)}_{\max}=d_{\min}^{-p}=2^{-(L-\ell)p}$. It follows that an
admissible block at level $\ell$ can be
$(\alpha_{\ell},L-\ell+1,\frac{\varepsilon}{3})$-block-encoded via the naive
procedure for dense matrices (Lemma 3.2), where
$\alpha_{\ell}=2^{(L-\ell)(1-p)}$, using two queries to $\mathcal{O}_{k}$ and
$O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ additional one- and two-
qubit gates. Note that in using Lemma 3.2 to block-encode a block
$\mathbf{K}^{\sigma,\rho}$ at level $\ell$, one actually needs to query
$k(x_{i},x_{j})/k^{(\ell)}_{\max}$ rather than $k(x_{i},x_{j})$ itself. This
is readily achieved from the oracle $\mathcal{O}_{k}$ by conditioning on a
register that indexes the level $\ell$ and the block indices $\sigma,\rho$
(see Section C.1 for details). From here, each $\mathbf{K}^{(\ell)}$, which is
a block-sparse matrix of block-sparsity 3, can be
$(3\alpha_{\ell},L+3,\varepsilon)$-block-encoded via Lemma 3.5 (particularly
Remark 3.6). Observe that while $\mathbf{K}^{(\ell)}$ is a $2^{L}\times 2^{L}$
matrix, the normalization factor of this block-encoding is only
$3\alpha_{\ell}$, which is better than the normalization factor of $N$ that
one would get when naively applying Lemma 3.2 to block-encode
$\mathbf{K}^{(\ell)}$. In addition, we still only need two queries to
$\mathcal{O}_{k}$ and $O(\operatorname{polylog}(\frac{N}{\varepsilon}))$
additional gates because all blocks $\mathbf{K}^{\sigma,\rho}$ are constructed
in parallel in Lemma 3.5 (see detailed analysis in Section C.1).

The adjacent interaction part $\mathbf{K}_{ad}$, which is a tridiagonal
matrix, can be $(3,L+3,\varepsilon)$-block-encoded with the procedure for
sparse matrices (Lemma B.5), using two queries to $\mathcal{O}_{k}$ and
$O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ one- and two-qubit gates.
Finally, we use the procedure for linearly combining block-encoded matrices
(Lemma B.3) to obtain a block-encoding of
$\mathbf{K}=\sum_{\ell=2}^{L}\mathbf{K}^{(\ell)}+\mathbf{K}_{ad}$ as
$U=\sum_{\ell=2}^{L}3\alpha_{\ell}U^{(\ell)}+3U_{ad}$. This step requires a
$(\log L)$-qubit state preparation pair (Definition 3.3) which prepares the
coefficients in the linear combination. We neglect the error in implementing
this state preparation pair since it can be performed with $\log L=\log\log N$
qubit operations (see Section C.1). This entire procedure results in a
$(\alpha,L+\log L+3,\alpha\varepsilon)$-block-encoding of the _exact_ kernel
matrix $\mathbf{K}$, where for $L=\log N$,

$\displaystyle\alpha$ $\displaystyle=3+3\sum_{\ell=2}^{L}2^{(L-\ell)(1-p)}$
(17)
$\displaystyle=\begin{cases}3(1+\frac{2^{1-p}-N^{1-p}}{2^{1-p}-4^{1-p}})\leq\Theta(1)&\text{
if }p>1\\\ 3\log N&\text{ if }p=1\\\
3(1+\frac{N^{1-p}-2^{1-p}}{4^{1-p}-2^{1-p}})&\text{ if }p<1\\\ \end{cases}.$

As before, the linear combination step does not significantly increase the
circuit complexity since all levels can be constructed in parallel with only
two queries to the oracle $\mathcal{O}_{k}$. See Section C.1 for more detailed
analysis of the entire block-encoding procedure.

######  Remark 4.1 (Optimality).

The block-encoding of the polynomially decaying kernels in Equation 13 is
optimal in the normalization factor $\alpha$. Indeed, a lowerbound on the
operator norm of $\mathbf{K}$ can be obtained by evaluating the norm of the
vector $\mathbf{K}\ket{\mathbf{1}}$, where $\ket{\mathbf{1}}$ is the
normalized all-ones state. Observe that,
$\|\mathbf{K}\|\geq\|\mathbf{K}\ket{\mathbf{1}}\|\geq\int_{1}^{N}x^{-p}dx$.
Thus, for $p\neq 1$, $\|\mathbf{K}\|\geq\frac{N^{1-p}-1}{1-p}$ and for $p=1$,
$\|\mathbf{K}\|\geq\ln N$. Clearly, the existence of our block-encoding
procedure also implies that these (up to a constant factor) are an upper bound
for $\|\mathbf{K}\|$.

######  Remark 4.2 (Generalized polynomially decaying kernel).

Since the above analysis only relies on the maximum entry in admissible
blocks, it still holds for the case of non-uniform (bounded) masses $m_{i}$ or
more general kernels of the form
$k(x,x^{\prime})=|x-x^{\prime}|^{-p}G(x,x^{\prime})$, where $G(x,x^{\prime})$
is a bounded function such as
$G(x,x^{\prime})=\frac{1}{\sqrt{1+\epsilon|x-x^{\prime}|^{2}}}$,
$G(x,x^{\prime})=e^{-i|x-x^{\prime}|}$, etc. In these cases, we simply scale
the kernel matrix such that the maximum entry is $1$ (if needed) and apply the
above analysis. This generalized class includes a variety of kernels, ranging
from power-law interactions to common Green’s functions [36].

The above block-encoding procedure similarly applies for higher-dimensional
settings where coordinates $x_{j}$ are in a $d$-dimensional space with $d$
constant. In Section D.3, we show that there are still only $\log N$
hierarchical levels in this higher dimensional setting and that each level
$\mathbf{K}^{(\ell)}$ is still block-sparse. The two key features of the
hierarchical splitting that allow for optimal block-encodings are maintained,
and the normalization factor is still $\alpha=\Theta(\|\mathbf{K}\|)$. The
above optimal block-encoding construction and remarks are summarized in the
following lemma.

######  Lemma 4.3 (Optimal block-encoding of polynomially decaying kernel
matrix using $\mathcal{H}$-matrix).

Let $k(x,x^{\prime})=\|x-x^{\prime}\|^{-p}G(x,x^{\prime})$, where $p>0$ and
$G(x,x^{\prime})$ is a bounded function, i.e. $|G(x,x^{\prime})|<C$. Assume
there exists an efficient quantum circuit $\mathcal{O}_{k}$ that performs
$\mathcal{O}_{k}:\ket{i}\ket{j}\ket{z}\rightarrow\ket{i}\ket{j}\ket{z\oplus\tilde{k}(x_{i},x_{j})},$
where $\tilde{k}(x_{i},x_{j})$ is a qubit string description of
$k(x_{i},x_{j})$. Then the matrix $\mathbf{K}=(k(x_{i},x_{j}))_{i,j=0}^{N-1}$
can be $(\alpha,\operatorname{polylog}N,\varepsilon)$-block-encoded with
$\alpha=\Theta(\|\mathbf{K}\|)$ using the hierarchical matrix splitting. That
is, the block-encoding is optimal according to Definition 3.1. This block-
encoding can be constructed with two queries to $\mathcal{O}_{k}$ and
$O(\operatorname{polylog}\frac{N}{\varepsilon})$ additional one- and two-qubit
gates.

Our block-encoding procedure is optimal for common kernel matrices, which are
neither sparse nor low-rank, thus expanding the scope of application of
quantum computers in the block-encoding framework. In Section 5, for example,
we discuss applications of this block-encoding in $N$-body force computations
and solving quantum linear systems.

######  Remark 4.4.

When sources are not close to uniformly distributed, one must first apply
classical adaptive methods [27, 28] which discretize the domain to form a
finer mesh, such that the number of nonzero-mass sources in each mesh site is
$O(1)$. The dimension of the kernel matrix in this method depends on the
inverse of the smallest pairwise distances between the actual sources.
Therefore, the runtime of the block-encoding remains polylogarithmic in the
number of particles $N$ so long as the minimum pairwise distance is not
exponentially small (see Section C.3).

###  4.2 Exponentially decaying kernels

The procedure used in the previous section readily generalizes to
exponentially decaying kernels. Consider, for example, the following kernel

$k(x,x^{\prime})=e^{-\|x-x^{\prime}\|^{q}}P(\|x-x^{\prime}\|),$ (18)

where $q>0$ and $P(\|x-x^{\prime}\|)$ is a function growing polynomially or
slower. Specifically, $|P(y)|\leq O(y^{k})$. Typical kernels of this form
include the Laplace and Gaussian kernels, which are commonly used in machine
learning [15].

As before, we first hierarchically decompose the matrix into admissible
blocks. For each admissible block $\mathbf{K}^{\sigma,\rho}$ at level $\ell$,
the minimum pairwise distance between particles is $d_{\min}=2^{L-\ell}$ and
therefore the maximum entry in the block can be bounded as

$\displaystyle\max_{(i,j)\in(\sigma,\rho)}|k(x_{i},x_{j})|$
$\displaystyle\leq\max_{(i,j)\in(\sigma,\rho)}\frac{|x_{i}-x_{j}|^{k}}{\sum_{t=0}^{\infty}|x_{i}-x_{j}|^{qt}/t!}$
(19)
$\displaystyle\leq\max_{(i,j)\in(\sigma,\rho)}\frac{t!}{|x_{i}-x_{j}|^{qt-k}}$
$\displaystyle\leq t!2^{-(L-\ell)(qt-k)},$

where we choose an integer $t$ that satisfies $qt-k>1$.

Thus, an admissible block $\mathbf{K}^{\sigma,\rho}$ at level $\ell$ can be
$(\alpha_{\ell},L-\ell+3,\varepsilon/3)$-block-encoded via the naive procedure
for dense matrices (Lemma 3.2), where
$\alpha_{\ell}=t!2^{(L-\ell)(1-(qt-k))}$. From here, an efficient block-
encoding of $\mathbf{K}$ can be obtained via the same procedure described in
the previous section, where $p$ is now replaced by $qt-k>1$. It can be
verified that the normalization factor of this block-encoding is $\Theta(1)$,
which is optimal since the maximum entry of $\mathbf{K}$ is of magnitude
$\Theta(1)$.

We note, however, that one can make use of a sparsification approach, in which
the $\mathbf{K}$ is approximated as a sparse band matrix, and apply the block-
encoding procedure for a sparse matrix (Lemma B.5). This approach yields a
block-encoding with normalization factor logarithmic in the sparsification
error (see Section 5.3).

###  4.3 Block-encoding of general $\mathcal{H}$-matrices

In the setting most consistent with that of classical $\mathcal{H}$-matrix
literature, one may be provided oracle access to the low-rank vectors in
$\bm{\Psi}^{\sigma,\rho}$ and $\bm{\Phi}^{\rho}$ which approximate an
admissible block
$\mathbf{\widetilde{K}}^{\sigma,\rho}=\bm{\Psi}^{\sigma,\rho}\mathbf{D}\left(\bm{\Phi}^{\rho}\right)^{\dagger}$
(as seen in Equation 7). Such a situation may also arise when one wants to
apply an arbitrary $\mathcal{H}$-matrix which has no connection to any kernel
function, e.g., see [37]. In this case, we can use the block-encoding
procedure for low-rank matrices (Lemma 3.4) to block-encode
$\mathbf{\widetilde{K}}^{\sigma,\rho}$.

After block-encoding the admissible blocks in each level of the hierarchical
splitting, we follow the same procedure in previous sections to obtain the
entire kernel matrix. Namely we apply the procedure for block-sparse matrices
(Lemma 3.5) to construct a block-encoding of $\mathbf{K}^{(\ell)}$ for each
level $\ell$ of the hierarchy; then we sum over all levels using the procedure
for linear combination of block-encoded matrices (Lemma B.3) to obtain a
block-encoding of $\mathbf{K}$. The normalization factor of the block-encoding
of $\mathbf{K}$ is bounded as

$\alpha=O\left(\sum_{\ell=2}^{\log N}\alpha_{\ell}\right),$ (20)

where $\alpha_{\ell}$ is normalization factor of the admissible blocks
$\mathbf{K}^{\sigma,\rho}$ at level $\ell$ (for ease of analysis we assume
admissible blocks at the same level have the same block-encoding normalization
factor). As before, assuming the maximum entry of $\mathbf{K}$ is 1, the naive
block-encoding (Lemma 3.2) yields a block-encoding with a normalization factor
of $N$, the dimension of the matrix. Whereas, $\mathcal{H}$-matrix block-
encoding can yield an exponentially better normalization factor over the naive
approach when $\alpha_{\ell}$ are $O(\operatorname{polylog}(N))$. This is
_not_ the case for some kernels studied in classical $\mathcal{H}$-matrix
literature [25, 13] such as the log kernel
$k(x,x^{\prime})=\log(|x-x^{\prime}|)$ and the multiquadric kernel
$k(x,x^{\prime})=\sqrt{c+|x-x^{\prime}|^{2}}$. In fact, we show in Section C.6
that these particular kernels and most polynomially _growing_ kernels can
already be optimally block-encoded using the naive approach of Lemma 3.2. The
naive approach works intuitively because these kernels are numerically low-
rank and their top singular vectors are nearly uniform.

###  4.4 Variants of hierarchical splitting

The form of the hierarchical splitting can be slightly modified for other
types of decaying kernels. For example, consider kernels of the form
$k(x,x^{\prime})=\frac{1}{|(x-x^{\prime})+c|^{p}}$ on the domain
$\Omega=[0,N-1)$ and $-N<c<N$ is an integer. In this case, we can use a
“shifted” hierarchical splitting to obtain an optimal block-encoding of the
kernel matrix $\mathbf{K}=(|i-j+c|^{-p})_{i,j=0}^{N-1}$, in which the
hierarchy is shifted in the horizontal direction (along the rows). Similarly,
for kernels of the form $k(x,x^{\prime})=\frac{1}{(|x-x^{\prime}|-c)^{p}}$
($0\leq c<N$), we can apply a shifted hierarchical splitting along the skew-
diagonal direction. We illustrate this approach in Section C.2.

In Section C.4, we also describe a more general block-encoding that exploits
the structure of entry magnitudes in a given matrix even if the entries of the
admissible blocks are not neighboring (in fact, Definition 2.1 does not
require the entries of an admissible block to be contiguous). At a high level,
we use a pointer oracle which, conditioned on the level $\ell$, points to the
entries whose magnitudes are in the range $[2^{-(\ell+1)},2^{-\ell})$. When
this oracle can be efficiently implemented, this block-encoding procedure
might find applications in other classes of dense matrices.

######  Proposition 4.5 (Block-encoding using generalized hierarchical
splitting).

Let $A\in\mathbb{R}^{N\times N}$ be a Hermitian matrix with non-negative
entries which are provided via an oracle
$\mathcal{O}_{A}:\ket{i}\ket{j}\ket{z}\rightarrow\ket{i}\ket{j}\ket{z\oplus\tilde{a}_{ij}}$,
where $\tilde{a}_{ij}$ is an exact bit string description of $a_{ij}$. Suppose
$N=2^{L}$ and $2^{-L}\leq a_{ij}\leq 1$. For any $0\leq\ell<L=\log N$ and a
column $j$, define its level-$\ell$ row-index subset to be
$I_{\ell}(j)=\\{i|2^{-(\ell+1)}<a_{ij}\leq 2^{-\ell}\\}$ and let
$n_{\ell}(j)=|I_{\ell}(j)|$. Suppose there exist $n_{\ell}$ and $\gamma$ such
that $\gamma n_{\ell}\leq n_{\ell}(j)\leq n_{\ell}$ for any $j$. Furthermore,
suppose there exists an _invertible_ function $f_{j}(\ell,k)$ which returns
the row index of the $k$-th element (according to some fixed ordering, e.g.
from top to bottom in the column) in $I_{\ell}(j)$. Then, one can construct a
block-encoding of $A$ with a normalization factor of at most $\frac{2\log
N}{\gamma}\|A\|$ using two queries to $\mathcal{O}_{A}$,
$O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ elementary gates, and
$O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ ancilla qubits.

Intuitively, $\gamma$ measures how balanced the magnitudes of the matrix
entries are across rows and columns. If
$\frac{1}{\gamma}=\operatorname{polylog}(N)$, this yields a good block-
encoding according to Definition 3.1. For instance, $\frac{1}{\gamma}$ is a
constant for kernel matrices. In the proof of 4.5, we also present an
efficient state preparation procedure using methods of $\mathcal{H}$-matrix
splitting in Section C.4, which could be of independent interest.

##  5 Applications

In this section, we discuss applications of our quantum hierarchical block-
encoding method in two commonly studied settings: (i) gravitational potential
calculation [23] and (ii) solving integral equations [13]. The first problem
requires performing matrix-vector multiplication on a quantum computer, while
the second problem is equivalent to solving a quantum linear system. We also
compare our procedure to the sparsification approach and address potential
issues with large condition numbers.

We first state the following lemmas, which allow one to apply a block-encoded
matrix or its inverse to a quantum state. Other polynomial transformations of
block-encoded matrices can be performed by the quantum singular value
transform (Theorem 17 of [3]).

######  Lemma 5.1 (Applying block-encoded matrices).

If $U$ is an $(\alpha,a,\varepsilon)$-block-encoding of the matrix
$A\in\mathbb{C}^{2^{s}\times 2^{s}}$ and $U_{\mathbf{x}}$ is a unitary that
prepares the $s$-qubit quantum state $\ket{\mathbf{x}}$, then it takes
$O\left(\frac{\alpha}{\|A\|}\kappa\right)$ queries to $U$ and $U_{\mathbf{x}}$
to obtain a quantum state $\ket{\mathbf{y}}$ such that
$\|\ket{\mathbf{y}}-\ket{A\mathbf{x}}\|\leq\varepsilon$, where $\kappa$ is the
condition number of $A$.

###### Proof.

Applying the operator $U(I_{a}\otimes U_{\mathbf{x}})$ on the state
$\ket{0^{a}}\otimes\ket{0^{s}}$ we obtain

$\ket{0^{a}}\frac{A}{\alpha}\ket{\mathbf{x}}+\ket{\text{garbage}},$ (21)

where $\ket{\text{garbage}}$ is orthogonal to the subspace $\ket{0^{a}}$. By
post-selecting the subspace $\ket{0^{a}}$, we obtain the desired quantum state
$\ket{A\mathbf{x}}$. This step succeeds with a probability of
$\left\|\frac{A\ket{\mathbf{x}}}{\alpha}\right\|^{2}=\Omega((\frac{\|A\|}{\alpha\kappa})^{2})$.
Using amplitude amplification [38], this can be improved to
$\Omega(\frac{\|A\|}{\alpha\kappa})$, giving us the specified query
complexity. ∎

######  Lemma 5.2 (Applying inverse block-encoded matrices, adapted from
[22]).

If $U$ is an $(\alpha,a,\varepsilon)$-block-encoding of an invertible matrix
$A\in\mathbb{C}^{2^{s}\times 2^{s}}$ and $U_{\mathbf{b}}$ is a unitary that
prepares the $s$-qubit quantum state $\ket{\mathbf{b}}$, then it takes
$O\left(\frac{\alpha\kappa}{\|A\|}\log(1/\varepsilon)\right)$ queries to $U$
and $U_{\mathbf{b}}$ to obtain a quantum state $\ket{\mathbf{y}}$ such that
$\|\ket{\mathbf{y}}-\ket{A^{-1}\mathbf{b}}\|<\varepsilon$.

###### Proof.

We apply the main theorem of [22], which assumed $\alpha=\|A\|=1$ to achieve a
query complexity of $O(\kappa\log(1/\varepsilon))$. In general, $\kappa$ needs
to be replaced by the inverse of the minimum singular value of the block-
encoded part. In our case, this singular value is
$\frac{\|A\|}{\alpha\kappa}$. ∎

As noted in Section 3, the query complexities of the above procedures, and
more generally quantum algorithms in the block-encoding framework [3], only
depend on the ratio $\alpha/\|A\|$. Hence, a block-encoding is optimal
(Definition 3.1) when $\alpha=\Theta(\|A\|)$. In the following, we apply such
optimal block-encodings in Section 4 to two specific applications of kernel
matrices.

###  5.1 Gravitational potential calculation: Quantum fast multipole method

We apply the block-encoding obtained in Section 4.1 to solve the gravitational
potential problem [23] on a quantum computer. Since we use a quantum version
of the hierarchical splitting, this is a quantum analogue to the classical
fast multipole method [39]. Consider a 1D system of $N$ particles which are
approximately uniformly distributed in the domain $[0,O(N)]$, where the mass
and location of particle $j$ are $m_{j}$ and $x_{j}$, respectively.
Furthermore, assume that $m_{j}$ are $\Theta(1)$ and the smallest pairwise
distance between the particles is $\Omega(1)$. We would like to compute the
accumulated potential at every particle $i$, stored in a vector $\Phi$ with
entries $\Phi_{i}=\sum_{j\neq i}m_{j}|x_{i}-x_{j}|^{-1}$.

######  Proposition 5.3 (Quantum fast multipole method).

Given an oracle $\mathcal{O}_{k}$ that performs
$\mathcal{O}_{k}:\ket{i}\ket{j}\ket{z}\rightarrow\ket{i}\ket{j}\ket{z\oplus\tilde{k}_{ij}},$
where $\tilde{k}_{ij}$ is a $b$-qubit string description of $k(x_{i},x_{j})$,
and a procedure $P_{m}$ that prepares the quantum state $\ket{\mathbf{m}}$
which stores the normalized masses of the particles, we can obtain a quantum
state $\ket{\phi}$ that $\varepsilon$-approximates the normalized vector
$\Phi/\|\Phi\|$, whose $j$-th entry is the accumulated potential at particle
$j$, using $O(1)$ queries to $P_{m}$,
$O\left(\operatorname{polylog}\frac{N}{\varepsilon}\right)$ one- and two-qubit
gates, and $O\left(b,\operatorname{polylog}\frac{N}{\varepsilon}\right)$
additional qubits.

###### Proof.

We first use $P_{m}$ to create the state
$\ket{\mathbf{m}}=\sum_{j=0}^{N-1}\frac{m_{j}}{\|\mathbf{m}\|}\ket{j}$ (assume
$\log N$ is an integer). Let $\mathbf{K}=(k(x_{i},x_{j}))_{i,j}$, observe that
$\ket{\Phi}=\mathbf{K}\ket{\mathbf{m}}=\frac{\Phi}{\|\mathbf{m}\|}$. Using
Lemma 4.3, we can construct a unitary $U$, which is an
$(\alpha,a,\varepsilon)$-block-encoding of $\mathbf{K}$, with the specified
gate complexity and wherein $\alpha=\Theta(\log N)$ and $a=O(\log N)$. Then,
we apply Lemma 5.1, with a more careful analysis on the post-selection step.
Here, the probability of succesfully post-selecting the subspace $\ket{0^{a}}$
is $\|\frac{\Phi}{\alpha\|\mathbf{m}\|}\|^{2}=\Theta(1)$ since
$\|\mathbf{m}\|=\Theta(\sqrt{N})$ and $\|\Phi\|=\Omega(\sqrt{N}\log N)$ (by
noting that the particle masses are $\Theta(1)$ and using similar calculations
to those in Remark 4.1). This gives us the specified query complexity. ∎

We remark that the (non-unitary) state preparation procedure $P_{m}$ can be
efficiently implemented when $m_{j}=\Theta(1)$ (see Theorem 1 of [40]). In
addition, if particles are not (nearly) uniformly distributed, we can
discretize the domain to form a finer mesh as described in Remark 4.4. Here we
have assumed these distances are $\Omega(1)$, hence the size of the mesh (and
the size of the kernel matrix) is $O(N)$. This method also enables computing
the potential at locations other than the sources themselves. This proposition
also applies to $N$-body calculations in higher-dimensional systems as the
hierarchical splitting and Lemma 4.3 can be generalized to these cases (see
Section D.3).

###  5.2 Integral equation: Quantum algorithm for linear systems of kernel
matrix

We provide an application of our hierarchical block-encoding procedure in
solving a quantum linear systems of kernel matrices. These problems arise in
many research areas, such as Gaussian processes [15] and integral equations
[14]. Consider the following integral equation over the circular strip of unit
radius and height and width $\lambda\ll 1$ as shown in Figure 2.

Figure 2: (A) A thin circular strip of unit radius and height $\lambda\ll 1$.
When using the collocation method to solve an integral equation on its
surface, one divides the strip into $N$ “panels” such that $\lambda=O(1/N)$.
(B) The “cyclic” hierarchical splitting of the kernel matrix that arises when
applying the collocation method to solve for an integral equation defined over
the strip (Equation 22).

$g(\mathbf{x})=f(\mathbf{x})+\int_{\Omega}k(\mathbf{x},\mathbf{x^{\prime}})f(\mathbf{x^{\prime}})d\mathbf{x^{\prime}},$
(22)

where
$k(\mathbf{x},\mathbf{x^{\prime}})=\|\mathbf{x}-\mathbf{x^{\prime}}\|^{-p}$
with $0<p\leq 2$, $g(\mathbf{x})$ is a known smooth function, and
$f(\mathbf{x})$ is the unknown function we would like to solve. Equations of
this form are commonly found in, e.g.,, Dirichlet problems [14], where
$k(\mathbf{x},\mathbf{x^{\prime}})$ is the Coulomb potential, $g(\mathbf{x})$
is the given potential on the strip, and $f(\mathbf{x})$ is the unknown charge
density. Here, we use the collocation method to numerically solve this
integral equation which represents the unknown function in a basis
$f(\mathbf{x})=\sum_{i=1}^{N}f_{i}\varphi_{i}(x)$ and solves for the
coefficients $f_{i}$ such that the integral equation is satisfied at specified
locations. Particularly, we use the piecewise-constant function as the basis
$\varphi_{i}(\theta)=1$ if $\theta\in[2\pi i/N,2\pi(i+1)/N)$ and 0 otherwise.
We refer to the segment $[2\pi i/N,2\pi(i+1)/N)$ as _panel_ $i$. Then, the
integral on the RHS of Equation 22 can be rewritten as

$\sum_{i=1}^{N}f_{i}\int_{\text{panel
}i}\frac{1}{\|\mathbf{x}-\mathbf{x^{\prime}}\|^{p}}d\mathbf{x^{\prime}},$ (23)

In the collocation method, we solve for the integral equation at the centroids
of the panels located at $c_{i}=2\pi(i+0.5)/N$. In other words, we solve the
following linear system

$\mathbf{g}=(I+\mathbf{K})\mathbf{f},$ (24)

where bold text indicates discretized values, e.g.,, $\mathbf{g}_{i}=g(c_{i})$
and $\mathbf{K}_{ij}=\int_{\text{panel
}j}\|c_{i}-\mathbf{x^{\prime}}\|^{-p}d\mathbf{x^{\prime}}$.

Observe that this collocation method can handle the singularity when $i=j$ as
$\mathbf{K}_{ii}=O(\lambda^{3-p})$. Furthermore, when $|i-j|$ is large,
$\mathbf{K}_{ij}\approx\text{distance}^{-p}\cdot\text{panel
volume}\approx\frac{\lambda^{2}}{N(2\sin(\pi|i-j|/N))^{p}}$. For panels $i$,
$j$ close to each other (e.g., $|i-j|<O(1)$), simple adaptive methods can be
used to approximate the integral. For instance, one can subdivide the panel
$j$ into eight smaller panels and approximate the integral as the sum of these
eight sub-panels; we neglect the details and refer the interested reader to
[41]. In summary, the entries of $\mathbf{K}$ are simple to calculate and they
approximately are

$\mathbf{K}_{ij}\approx\begin{cases}\frac{\lambda^{2}}{N(2\sin(\pi|i-j|/N))^{p}}&\text{
for }i\neq j\\\ O(\lambda^{3-p})&\text{ for }i=j\end{cases}.$ (25)

The error bound and the uniqueness of the solution in the collocation method
are well-studied in classical literature, we refer to [14, 41, 42] for more
details.

We now proceed to block-encode this kernel matrix. Observe that
$\frac{1}{N(2\sin(\pi|i-j|/N))^{p}}\leq\frac{1}{(\pi|i-j|)^{p}}$ for any
panels $i,j$ such that $|i-j|\leq N/2$, hence we can apply a “cyclic” version
of the hierarchical splitting in Section 2 to _optimally_ block-encode
$\mathbf{K}$ (see Figure 2B). As analyzed in Section 4.1, this block-encoding
only has a circuit depth of $O(\operatorname{polylog}\frac{N}{\varepsilon})$
and a normalization factor of $\alpha=\Theta(\|\mathbf{K}\|)$, as the operator
norm can be lowerbounded using the same technique as in Remark 4.1 by noticing
that $\frac{1}{N(2\sin(\pi|i-j|/N))^{p}}\geq\frac{1}{(2\pi|i-j|)^{p}}$ for any
$|i-j|\leq N/2$.

Finally, it is straightforward to use Lemma B.3 (linear combination of block-
encoded matrices) to obtain an optimal block-encoding of $I+\mathbf{K}$. Then,
we directly apply Lemma 5.2 to solve the linear system in Equation 24. Since
the block-encoding is optimal, the query complexity is
$O(\kappa\log(1/\varepsilon))$, for a total runtime of
$O(\kappa\operatorname{polylog}(N/\varepsilon))$ when including the circuit
complexity of the block-encoding (Lemma 4.3).

######  Proposition 5.4 (Solving discretized integral equations).

Consider the integral equation on a ring as described in Equation 22. We first
cast the equation into the $N\times N$ linear system in Equation 24 using the
collocation method with $N$ discretized points. Then, we can obtain its
numerical solution as a quantum state in time
$O(\kappa\operatorname{polylog}(N/\varepsilon))$, where $\kappa$ is the
condition number of the matrix $I+\mathbf{K}$ in Equation 24 and $\varepsilon$
is the error bound of the quantum operation.

We address two practical aspects when applying 5.4. First, we must prepare
$\mathbf{g}$ as a quantum state. In Appendix E, we present an efficient
procedure for preparing $\ket{\mathbf{g}}$ in the case where $\mathbf{g}$ is
sampled from a smooth function. At a high level, our procedure prepares the
quantum state in the Fourier regime, which is approximately sparse due to the
smoothness of $g$; then we apply the quantum Fourier transform to obtain the
desired real-regime state. To prove 4.5, we also present another state
preparation procedure using ideas from $\mathcal{H}$-matrix splitting in
Section C.4 which could be of independent interest. Second, the query
complexity of the block-encoding depends on the condition number $\kappa$ of
the matrix $I+\mathbf{K}$, which can be bounded as follows

$\kappa\leq\frac{1+\|\mathbf{K}\|}{1-\|\mathbf{K}\|}.$ (26)

Note that here, $\|\mathbf{K}\|$ scales with $\lambda$. Since we require
$\lambda\ll 1$ for the collocation method to work, the operator norm
$\|\mathbf{K}\|$, and thus $\kappa$, is bounded. Therefore, the runtime of 5.4
is $O(\operatorname{polylog}\frac{N}{\varepsilon})$.

###  5.3 Notes on sparsification approach

Since the entries of kernel matrices often decay along a row or column, one
may assume that they can be approximated by only considering a sparse set of
matrix entries and applying the naive block-encoding method for sparse
matrices (Lemma B.5). Here, we analyze this approach showing that it is only
favorable when kernel entries decay exponentially and not polynomially. In
fact, the runtime of this method is exponentially worse in the target error
bound $\varepsilon$ compared to that of our approach. Consider the
polynomially decaying kernel in Equation 13. One could attempt to approximate
this kernel matrix by a band matrix $\mathbf{K}_{b}$ of sparsity $d$. Defining
the error to be $\|\mathbf{K}-\mathbf{K}_{b}\|$, it can be seen that the error
is of magnitude $\int_{d}^{N}x^{-p}dx$. Since this integral diverges when
$p\leq 1$, one can only sparsify the matrix when $p>1$. In this case, to
achieve an error bound of $\varepsilon_{s}$, the sparsity $d$ needs to be
$\Omega(\varepsilon_{s}^{\frac{1}{1-p}})$. Then, one can apply Lemma B.5
(block-encoding of sparse matrices) to directly block-encode the sparsified
matrix with a normalization factor of
$\Theta(\varepsilon_{s}^{\frac{1}{1-p}})$ using $O(\log
N+\operatorname{polylog}(\frac{1}{\varepsilon}))$ extra gates and $O(1)$
oracle queries. This results in an exponentially worse runtime in terms of the
error bound compared to our method, which has a normalization factor of
$O(1)$.

For the exponential kernels, however, the sparsification method requires
keeping $\Theta(\log(1/\varepsilon_{s}))$ or more entries per row or column.
With Lemma B.5, the sparsified block-encoding has a normalization factor of
$\Theta(\log(1/\varepsilon_{s}))$ using the same gate and query complexities
as above. Thus, exponential kernels may admit a sparsification approach with
an extra factor of $\Theta(\log(1/\varepsilon))$ in the runtime.

###  5.4 Notes on condition number

The condition number $\kappa$ and its relation to the dimension $N$ of the
matrix play a central role in the runtime of our algorithms. Numerical
experiments shown in Section C.5 indicate that for purely polynomially
decaying kernels without any regularization and ignoring self-interacting
terms (diagonal entries set to zero), the condition number grows polynomially
with $N$, i.e., the minimum singular value decays with $N$. In practice,
however, the input vector often lives in the well-conditioned subspace where
these smaller singular values can be ignored as observed in the gravitation
example earlier. For force calculation problems, this is the case when
particles have comparable masses/charges. For numerical integration or
integral equations, this is equivalently the case when the input functions
change slowly.

We note that Lemma 5.2 for matrix inversion uses the discrete adiabatic
theorem, whose runtime depends directly on the condition number $\kappa$ of
the matrix operation. Hence, in the case where the matrix operation $A$ is
ill-conditioned but the input vector lives in the well-conditioned subspace,
we can use techniques that allow for filtering out small eigenvalues before
the inversion step such as the QSVT [3]. The query complexity of the block-
encoding in matrix inversion in this case is
$O(\frac{\log(1/\zeta,1/\varepsilon)}{\zeta\|A^{-1}\ket{b}\|})$, where $\zeta$
is a chosen threshold bounding the singular values in the well-conditioned
subspace.

Fortunately, for applications in numerical analysis or machine learning, the
kernel matrix is often regularized by an added matrix, e.g., the addition of
an identity matrix as in Fredholm integral equations of the second kind or the
noise matrices $\sigma_{n}^{2}I$ in Gaussian processes and kernel ridge
regression [14, 15]. In the integral equation setting, this regularization is
equivalent to scaling the magnitude of self-interacting terms (diagonal
entries). Furthermore, if the input vector is potentially in the ill-
conditioned subspace, preconditioning methods in [43] may be used to improve
the condition number.

##  6 Conclusion

Our results show that factoring certain kernel matrices into hierarchical or
$\mathcal{H}$-matrices leads to efficient algorithms for implementing such
kernel matrices on a quantum computer, thus expanding the scope of application
of quantum computers to a class of matrices which are neither low-rank nor
sparse. In fact, our approach to implementing such hierarchical matrices, in a
sense, blends prior quantum methods for implementing unitary block encodings
by hierarchically splitting a matrix into a sparse matrix component close to
the diagonal and progressively larger blocks farther away from the diagonals.
Looking forward, it is an interesting question whether any extensions to the
$\mathcal{H}$-matrix framework can lead to further improvements in the quantum
setting. Classically, the most significant extension of the
$\mathcal{H}$-matrix framework is perhaps the development of
$\mathcal{H}^{2}$-matrices, which incorporates a further level of hierarchical
nested bases to provide a logarithmic speedup in the matrix dimension for
performing matrix operations [24, 44]. Furthermore, recent classical work has
parameterized $\mathcal{H}$-matrices to integrate them into neural networks
and design learning algorithms with inductive biases that match the form of
the hierarchical splitting [37, 45]. Future work can study whether such
parameterizations can be performed as well in the quantum setting via
variational quantum algorithms or other quantum analogues to classical neural
networks [46].

## Acknowledgements

We thank the anonymous reviewers for helpful feedback and corrections on this
paper. This project was funded by AFOSR, ARO, and DARPA. BTK acknowledges
funding from the MIT Energy Initiative fellowship.

## References

  * [1] Aram W. Harrow, Avinatan Hassidim and Seth Lloyd  “Quantum Algorithm for Linear Systems of Equations”  In _Physical Review Letters_ 103.15 American Physical Society (APS), 2009  DOI: https://doi.org/10.1103/PhysRevLett.103.150502
  * [2] Andrew M. Childs, Robin Kothari and Rolando D. Somma  “Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision”  In _SIAM Journal on Computing_ 46.6 Society for Industrial & Applied Mathematics (SIAM), 2017, pp. 1920–1950  DOI: https://doi.org/10.1137/16M1087072
  * [3] András Gilyén, Yuan Su, Guang Hao Low and Nathan Wiebe  “Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics”  In _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ , 2019, pp. 193–204  DOI: https://doi.org/10.48550/arXiv.1806.01838
  * [4] Iordanis Kerenidis and Anupam Prakash  “Quantum Recommendation Systems”, 2016  DOI: https://doi.org/10.48550/arXiv.1603.08675
  * [5] Leonard Wossnig, Zhikuan Zhao and Anupam Prakash  “Quantum Linear System Algorithm for Dense Matrices”  In _Physical Review Letters_ 120.5 American Physical Society (APS), 2018  DOI: http://dx.doi.org/10.1103/PhysRevLett.120.050502
  * [6] András Gilyén, Zhao Song and Ewin Tang  “An improved quantum-inspired algorithm for linear regression”  In _Quantum_ 6 Verein zur Förderung des Open Access Publizierens in den Quantenwissenschaften, 2022, pp. 754  DOI: https://doi.org/10.22331/q-2022-06-30-754
  * [7] Changpeng Shao and Ashley Montanaro  “Faster quantum-inspired algorithms for solving linear systems”, 2021  DOI: https://doi.org/10.48550/arXiv.2103.10309
  * [8] David P. Woodruff  “Sketching as a tool for numerical linear algebra”  In _Foundations and Trends in Theoretical Computer Science_ 10.1–2 Now Publishers, 2014, pp. 1–157  DOI: 10.1561/0400000060
  * [9] Lin-Chun Wan et al.  “Asymptotic quantum algorithm for the Toeplitz systems”  In _Physical Review A_ 97.6 American Physical Society (APS), 2018  DOI: 10.1103/physreva.97.062322
  * [10] A Mahasinghe and J B Wang  “Efficient quantum circuits for Toeplitz and Hankel matrices”  In _Journal of Physics A: Mathematical and Theoretical_ 49.27 IOP Publishing, 2016, pp. 275301  DOI: 10.1088/1751-8113/49/27/275301
  * [11] Grecia Castelazo et al.  “Quantum algorithms for group convolution, cross-correlation, and equivariant transformations”  In _Physical Review A_ 106.3 APS, 2022, pp. 032402  DOI: 10.1103/PhysRevA.106.032402
  * [12] Andrew M. Childs and Wim Dam  “Quantum algorithms for algebraic problems”  In _Reviews of Modern Physics_ 82.1 American Physical Society (APS), 2010, pp. 1–52  DOI: 10.1103/revmodphys.82.1
  * [13] R Beatson and Leslie Greengard  “A short course on fast multipole methods”  In _Wavelets, multilevel methods, and elliptic PDEs_ , Numerical Mathematics and Scientific Computation  Oxford University Press, 1997, pp. 1–37  URL: https://nyuscholars.nyu.edu/en/publications/a-short-course-on-fast-multipole-methods
  * [14] Kendall Atkinson and Weimin Han  “Numerical solution of fredholm integral equations of the second kind”  In _Theoretical Numerical Analysis_ Springer, 2009, pp. 473–549  DOI: https://doi.org/10.1007/978-1-4419-0458-4_12
  * [15] Carl Edward Rasmussen  “Gaussian processes in machine learning”  In _Summer school on machine learning_ , 2003, pp. 63–71  Springer  DOI: https://doi.org/10.1007/978-3-540-28650-9_4
  * [16] W. Hackbusch  “A Sparse Matrix Arithmetic Based on $\mathcal{H}$-Matrices. Part I: Introduction to $\mathcal{H}$-Matrices”  In _Computing_ 62.2, 1999, pp. 89–108  DOI: 10.1007/s006070050015
  * [17] W. Hackbusch and B. N. Khoromskij  “A Sparse $\mathcal{H}$-Matrix Arithmetic. Part II: Application to Multi-Dimensional Problems”  In _Computing_ 64.1, 2000, pp. 21–47  DOI: 10.1007/PL00021408
  * [18] Guang Hao Low and Isaac L. Chuang  “Hamiltonian Simulation by Qubitization”  In _Quantum_ 3 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2019, pp. 163  DOI: 10.22331/q-2019-07-12-163
  * [19] Shantanav Chakraborty, András Gilyén and Stacey Jeffery  “The power of block-encoded matrix powers: improved regression techniques via faster Hamiltonian simulation”  In _ICALP_ , 2019  DOI: https://doi.org/10.4230/LIPIcs.ICALP.2019.33
  * [20] Guang Hao Low and Isaac L. Chuang  “Optimal Hamiltonian Simulation by Quantum Signal Processing”  In _Physical Review Letters_ 118.1 American Physical Society (APS), 2017  DOI: 10.1103/physrevlett.118.010501
  * [21] Joran Apeldoorn and András Gilyén  “Improvements in Quantum SDP-Solving with Applications”  In _ICALP_ , 2019  DOI: https://doi.org/10.4230/LIPIcs.ICALP.2019.99
  * [22] Pedro CS Costa et al.  “Optimal Scaling Quantum Linear-Systems Solver via Discrete Adiabatic Theorem”  In _PRX Quantum_ 3.4 APS, 2022, pp. 040303  DOI: 10.1103/PRXQuantum.3.040303
  * [23] Josh Barnes and Piet Hut  “A hierarchical O(N log N) force-calculation algorithm”  In _Nature_ 324.6096, 1986, pp. 446–449  DOI: 10.1038/324446a0
  * [24] Steffen Börm, Lars Grasedyck and Wolfgang Hackbusch  “Introduction to hierarchical matrices with applications”  In _Engineering Analysis with Boundary Elements_ 27.5, 2003, pp. 405–422  DOI: https://doi.org/10.1016/S0955-7997(02)00152-2
  * [25] Markus Fenn and Gabriele Steidl  “FMM and H-matrices: A Short Introduction to the Basic Idea”, 2002  URL: https://madoc.bib.uni-mannheim.de/744/1/TR-02-008.pdf
  * [26] W. Hackbusch, B. Khoromskij and S. A. Sauter  “On $\mathcal{H}^{2}$-matrices”  In _Lectures on Applied Mathematics_ Berlin, Heidelberg: Springer Berlin Heidelberg, 2000, pp. 9–29  DOI: https://doi.org/10.1007/978-3-642-59709-1_2
  * [27] J. Carrier, L. Greengard and V. Rokhlin  “A Fast Adaptive Multipole Algorithm for Particle Simulations”  In _SIAM J. Sci. Stat. Comput._ 9.4 USA: Society for IndustrialApplied Mathematics, 1988, pp. 669–686  DOI: 10.1137/0909044
  * [28] Jaswinder Pal Singh, Chris Holt, John L Hennessy and Anoop Gupta  “A parallel adaptive fast multipole method”  In _Proceedings of the 1993 ACM/IEEE Conference on Supercomputing_ , 1993, pp. 54–65  DOI: 10.1145/169627.169651
  * [29] E. Tyrtyshnikov  “Mosaic-Skeleton approximations”  In _Calcolo_ 33, 1996, pp. 47–57  DOI: 10.1007/BF02575706
  * [30] Achi Brandt  “Multilevel computations of integral transforms and particle interactions with oscillatory kernels”  In _Computer Physics Communications_ 65.1, 1991, pp. 24–38  DOI: https://doi.org/10.1016/0010-4655(91)90151-A
  * [31] Gregory Beylkin, Ronald R. Coifman and Vladimir Rokhlin  “Fast wavelet transforms and numerical algorithms I”  In _Communications on Pure and Applied Mathematics_ 44, 1991, pp. 141–183  DOI: https://doi.org/10.1002/cpa.3160440202
  * [32] Robin Kothari  “Efficient algorithms in quantum query complexity”, 2014  URL: https://uwspace.uwaterloo.ca/handle/10012/8625
  * [33] Yihui Quek and Patrick Rebentrost  “Fast algorithm for quantum polar decomposition, pretty-good measurements, and the Procrustes problem”, 2021  DOI: 10.48550/arXiv.2106.07634
  * [34] Lov Grover and Terry Rudolph  “Creating superpositions that correspond to efficiently integrable probability distributions”, 2002  DOI: 10.48550/arXiv.quant-ph/0208112
  * [35] Michael A. Nielsen and Isaac Chuang  “Quantum Computation and Quantum Information”  Cambridge University Press, 2011  DOI: 10.1017/CBO9780511976667
  * [36] George B. Arfken, Hans J. Weber and Frank E. Harris  “Chapter 10 - Green’s Functions”  In _Mathematical Methods for Physicists (Seventh Edition)_ Boston: Academic Press, 2013, pp. 447–467  DOI: https://doi.org/10.1016/B978-0-12-384654-9.00010-4
  * [37] Yuwei Fan, Lin Lin, Lexing Ying and Leonardo Zepeda-Nún͂ez  “A Multiscale Neural Network Based on Hierarchical Matrices”  In _Multiscale Modeling & Simulation_ 17.4, 2019, pp. 1189–1213  DOI: 10.1137/18M1203602
  * [38] Gilles Brassard, Peter Høyer, Michele Mosca and Alain Tapp  “Quantum amplitude amplification and estimation”  In _Quantum Computation and Information_ American Mathematical Society, 2002, pp. 53–74  DOI: 10.1090/conm/305/05215
  * [39] L. Greengard and V. Rokhlin  “A fast algorithm for particle simulations”  In _Journal of Computational Physics_ 73.2, 1987, pp. 325–348  DOI: 10.1016/0021-9991(87)90140-9
  * [40] Kosuke Mitarai, Masahiro Kitagawa and Keisuke Fujii  “Quantum analog-digital conversion”  In _Physical Review A_ 99.1 APS, 2019, pp. 012301  DOI: 10.1103/PhysRevA.99.012301
  * [41] P. K. Kythe  “An introduction to boundary element methods (1st ed.)”  CRC Press, 1995  URL: https://www.routledge.com/An-Introduction-to-Boundary-Element-Methods/Kythe/p/book/9780367449148
  * [42] W. Hackbusch  “The Panel Clustering Method for BEM”  In _Discretization Methods in Structural Mechanics_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1990, pp. 299–306  DOI: 10.1007/978-3-642-49373-7_28
  * [43] Yu Tong, Dong An, Nathan Wiebe and Lin Lin  “Fast inversion, preconditioned quantum linear system solvers, fast Green’s-function computation, and fast evaluation of matrix functions”  In _Physical Review A_ 104.3 American Physical Society (APS), 2021  DOI: 10.1103/physreva.104.032422
  * [44] Lin Lin, Jianfeng Lu and Lexing Ying  “Fast construction of hierarchical matrix representation from matrix–vector multiplication”  In _Journal of Computational Physics_ 230.10 Elsevier, 2011, pp. 4071–4087  DOI: 10.1016/j.jcp.2011.02.033
  * [45] Yuwei Fan et al.  “A multiscale neural network based on hierarchical nested bases”  In _Research in the Mathematical Sciences_ 6.2 Springer, 2019, pp. 1–28  DOI: 10.1007/s40687-019-0183-3
  * [46] Marco Cerezo et al.  “Variational quantum algorithms”  In _Nature Reviews Physics_ 3.9 Nature Publishing Group, 2021, pp. 625–644  DOI: 10.1038/s42254-021-00348-9
  * [47] Dominic W. Berry, Graeme Ahokas, Richard Cleve and Barry C. Sanders  “Efficient Quantum Algorithms for Simulating Sparse Hamiltonians”  In _Communications in Mathematical Physics_ 270.2, 2007, pp. 359–371  DOI: 10.1007/s00220-006-0150-x
  * [48] Andris Ambainis  “Variable time amplitude amplification and a faster quantum algorithm for solving systems of linear equations”, 2010  DOI: 10.48550/arXiv.1010.4458
  * [49] Yiğit Subaşı, Rolando D. Somma and Davide Orsucci  “Quantum Algorithms for Systems of Linear Equations Inspired by Adiabatic Quantum Computing”  In _Phys. Rev. Lett._ 122 American Physical Society, 2019, pp. 060504  DOI: 10.1103/PhysRevLett.122.060504
  * [50] Dong An and Lin Lin  “Quantum Linear System Solver Based on Time-optimal Adiabatic Quantum Computing and Quantum Approximate Optimization Algorithm”  In _ACM Transactions on Quantum Computing_ 3.2 Association for Computing Machinery (ACM), 2022, pp. 1–28  DOI: 10.1145/3498331
  * [51] Lin Lin and Yu Tong  “Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems”  In _Quantum_ 4 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2020, pp. 361  DOI: 10.22331/q-2020-11-11-361
  * [52] Jonathan M. Borwein and Peter B. Borwein  “Pi and the AGM: A Study in the Analytic Number Theory and Computational Complexity”  USA: Wiley-Interscience, 1987  URL: https://www.wiley.com/en-us/Pi+and+the+AGM%3A+A+Study+in+Analytic+Number+Theory+and+Computational+Complexity-p-9780471315155
  * [53] Sanjeev Arora and Boaz Barak  “Computational Complexity: A Modern Approach”  USA: Cambridge University Press, 2009  URL: https://doi.org/10.1017/CBO9780511804090
  * [54] Steffen Börm, Lars Grasedyck and Wolfgang Hackbusch  “Lecture notes on Hierarchical Matrices”, 2003  URL: https://www.mis.mpg.de/publications/other-series/ln/lecturenote-2103.html
  * [55] Wajih Halim Boukaram, George Turkiyyah and David E. Keyes  “Hierarchical Matrix Operations on GPUs: Matrix-Vector Multiplication and Compression”, 2019  DOI: https://doi.org/10.48550/arXiv.1902.01829

##  Appendix A Prior works on block-encoding and quantum linear systems
algorithms

In this section, we provide a brief overview and discuss the runtimes of
previous works in Table 1. As a reminder, the quantum linear system problem
(QLSP) seeks an approximate solution $\ket{x^{\prime}}$ to

$A\vec{x}=\vec{b}$ (27)

such that

$\left\|\ket{x^{\prime}}-\frac{\vec{x}}{\|\vec{x}\|}\right\|\leq\varepsilon.$
(28)

Let $N$ be the dimension; $d$ be the sparsity; $\kappa$ be the condition
number; and $k$ be the rank of the matrix $A$, which is assumed to be square,
normalized ($\|A\|=1$), and Hermitian 111If not, one can consider the
equivalent linear system $\begin{pmatrix}0&A\\\
A^{\dagger}&0\end{pmatrix}\vec{y}=\begin{pmatrix}\vec{b}\\\ 0\end{pmatrix}$,
whose solution is $\vec{y}=\begin{pmatrix}0\\\ \vec{x}\end{pmatrix}$.. The
original HHL algorithm [1] runs in time
$O(d^{2}\kappa^{2}\operatorname{polylog}(N)/\varepsilon)$ by performing a
Hamiltonian simulation of $A$ (i.e. $e^{-iAt}$) and quantum phase estimation
(QPE) subroutines. Here, the Hamiltonian simulation [47] and QPE steps
together contribute a factor of
$O(d^{2}\kappa\operatorname{polylog}(N)/\varepsilon)$, while the post-
selection step contributes the remaining factor of $O(\kappa)$.

Later, [2] noticed that the QPE step can be circumvented by utilizing the fact
that the inverse function $\frac{1}{x}$ can be written, to the desired
precision $\varepsilon/\kappa$, as a linear combination of
$O(\kappa\operatorname{polylog}(\kappa/\varepsilon))$ Chebyshev polynomials.
Combining this expansion with the variable-time amplitude amplification
technique in [48], they achieved a runtime of
$O(d\kappa\operatorname{polylog}(N\kappa/\varepsilon))$. This polynomial
expansion approach is generalized in the block-encoding framework by [3, 19]
to perform almost any polynomial matrix transformation with the same
complexity. All these works use a naive block-encoding in the oracle-access
model (Lemma 48 in [3], same as Lemma B.5 or Lemma 3.2).

In another approach, [5] used a QRAM data structure to implement a unitary $W$
whose eigenvalues are related to the singular values of $A$ (this algorithm
does not require $A$ to be square or Hermitian) via its Frobenius norm
$\|A\|_{F}$. Thus, performing QPE on $W$ yields a quantum singular value
estimation algorithm for $A$. Matrix multiplication and inversion can then be
implemented immediately with a rotation operator controlled on the register
storing the singular values, just like done in HHL. The runtime of both
operations, including the post-selection step, is
$O(\kappa^{2}\|A\|_{F}\operatorname{polylog}(N)/\varepsilon).$ For a full-rank
matrix whose singular values are $O(1)$, we have $\|A\|_{F}=\sqrt{N}$. This
QRAM model, however, is strictly stronger than the oracle-access model as it
requires the entries be loaded and pre-processed in the data structure.

Recent adiabatic algorithms for QLSP [49, 50, 51] also enjoy success and
culminate in [22], which achieved an optimal query complexity of
$O(\kappa\log(1/\varepsilon))$ by constructing a sophisticated schedule
function for the adiabatic evolution. This query complexity is better than the
algorithms in [2, 3, 5] and has been shown to be optimal [22]. However, all of
these assume that the matrix $A$ is perfectly block-encoded inside a unitary
$U=\begin{pmatrix}\frac{A}{\|A\|}&\cdot\\\ \cdot&\cdot\end{pmatrix}$ which is
generally hard to construct without prior knowledge of the structure of $A$.
The exact query complexity and total runtime in the general case is stated in
Lemma 5.2 and governed by the ratio $\frac{A}{\alpha}$, where $\alpha$ is the
normalization factor in the block-encoding.

The runtimes shown in Table 1 are obtained by combining the optimal adiabatic
solver in [22] (see Lemma 5.2) with the corresponding block-encoding used in
[2, 3, 19, 4, 5]. In particular, the naive oracle-access block-encoding in [2,
3, 19] yields a normalization factor of $\alpha=d=N$ for dense matrices. The
QRAM-based block-encoding in [4, 5] yields $\alpha=\|A\|_{F}=\sqrt{N}$ for
full-rank matrices. The runtime for matrix multiplication is simply obtained
by computing the probability of success in the post-selection step after
applying the block-encoding, as described in Section 3 and Lemma 5.1.

Quantum-inspired algorithms are a class of classical algorithms equipped with
sample-query access, a classical analogue to QRAM. The currently fastest
quantum-inspired algorithms for QLSP [6, 7] are based on stochastic gradient
descent (specifically randomized Kaczmarz method) for a runtime of
$\widetilde{O}(\|A\|_{F}^{6}\kappa^{8}/\varepsilon^{2})$ (where the notation
$\widetilde{O}$ suppresses the log factors). Assuming $\kappa=O(1)$ and the
matrix has rank $k$, this runtime is $\widetilde{O}(k^{6}/\varepsilon^{2})$,
which is probitive when the matrix is full-rank.

We note that all the algorithms above are “general-purpose” since they do not
assume any inherent structures in the matrix $A$. Whereas, our work builds on
these past works and focuses on optimizing the performance of QLSP algorithms
for dense and full-rank kernel matrices and potentially other classes of
matrices that suit the $\mathcal{H}$-matrix framework.

##  Appendix B Details of block-encodings

We list some previously established results in the block-encoding framework
and prove the helper lemmas presented in the main text. Most of the following
results are drawn from Section 4 of [3]. As a reminder, a block-encoding is
defined as follows.

######  Definition B.1 (Block-encoding [3]).

Suppose that $A$ is an $s$-qubit operator, $\alpha,\varepsilon>0$, then we say
that the $(s+a)$-qubit unitary $U$ is an $(\alpha,a,\varepsilon)$-block-
encoding of $A$, if

$\left\|A-\alpha(\bra{0}^{\otimes a}\otimes I_{s})U(\ket{0}^{\otimes a}\otimes
I_{s})\right\|\leq\varepsilon,$

where $\|\cdot\|$ denotes the operator norm of a matrix.

###  B.1 Linear combinations and multiplications of block-encodings

Combinations of block-encoded matrices can be efficiently combined linearly or
multiplied with each other. First we define state-preparation-pairs used to
prepare the coefficients in a linear combination.

######  Definition B.2 (State preparation pair [3]).

Let $y\in\mathbb{C}^{m}$ and $\|y\|_{1}\leq\beta$, the pair of unitaries
$(P_{L},P_{R})$ is called a $(\beta,n,\varepsilon)$-state-preparation-pair of
$y$ if $P_{L}\ket{0}^{\otimes n}=\sum_{j=0}^{2^{n}-1}c_{j}\ket{j}$ and
$P_{R}\ket{0}^{\otimes n}=\sum_{j=1}^{2^{n}-1}d_{j}\ket{j}$ such that
$\sum_{j=0}^{m-1}\left|\beta c_{j}^{*}d_{j}-y_{j}\right|\leq\varepsilon_{1}$
and $c_{j}^{*}d_{j}=0$ for any $j\in\\{m,\ldots,2^{n}-1\\}$.

We can then implement a linear combination of block-encoded matrices using an
above state preparation pair.

######  Lemma B.3 (Linear combination of block-encoded matrices, adapted from
[3]).

Let $A=\sum_{j=0}^{m-1}y_{j}A_{j}$ be an $s$-qubit operator where
$\|y\|_{1}\leq\beta$. Suppose $(P_{L},P_{R})$ is a
$(\beta,n,\varepsilon_{1})$-state-preparation-pair for $y$,
$W=\sum_{j=0}^{m-1}\ket{j}\bra{j}\otimes
U_{j}+\left(\left(I-\sum_{j=0}^{m-1}\ket{j}\bra{j}\right)\otimes I_{a}\otimes
I_{s}\right)$ is an $(n+a+s)$-qubit unitary such that for all $j\in[m]$ we
have that $U_{j}$ is an $\left(\alpha,a,\varepsilon_{2}\right)$-block-encoding
of $A_{j}$. Then the unitary $\widetilde{W}=(P_{L}^{\dagger}\otimes
I_{a}\otimes I_{s})W(P_{R}\otimes I_{a}\otimes I_{s})$ is an
$\left(\alpha\beta,a+n,\alpha\varepsilon_{1}+\beta\varepsilon_{2}\right)$-block-
encoding of $A$.

###### Proof.

We have:

$\displaystyle\left\|A-\alpha\beta\left(\bra{0}^{n}\bra{0}^{a}\otimes
I_{s}\right)\widetilde{W}\left(\ket{0}^{n}\ket{0}^{a}\otimes
I_{s}\right)\right\|$ $\displaystyle=$
$\displaystyle\left\|A-\alpha\sum_{j=0}^{m-1}\beta\left(c_{j}^{*}d_{j}\right)\left(\bra{0}^{a}\otimes
I_{s}\right)U_{j}\left(\ket{0}^{a}\otimes I_{s}\right)\right\|$
$\displaystyle\leq$
$\displaystyle\alpha\varepsilon_{1}+\left\|A-\alpha\sum_{j=0}^{m-1}y_{j}\left(\bra{0}^{a}\otimes
I_{s}\right)U_{j}\left(\ket{0}^{a}\otimes I_{s}\right)\right\|$
$\displaystyle\leq$
$\displaystyle\alpha\varepsilon_{1}+\sum_{j=0}^{m-1}|y_{j}|\left\|A_{j}-\alpha\left(\bra{0}^{a}\otimes
I_{s}\right)U_{j}\left(\ket{0}^{a}\otimes I_{s}\right)\right\|$
$\displaystyle\leq$
$\displaystyle\alpha\varepsilon_{1}+\beta\varepsilon_{2}.\qed$

We note that in the original work [3], Lemma 52 states that the error of the
block-encoding is $\alpha\varepsilon_{1}+\alpha\beta\varepsilon_{2}$, which is
a typo. Lemma B.3 also applies to cases where $P_{L}$ and $P_{R}$ are non-
unitary. As long as block-encodings of $P_{L}$, $P_{R}$ are provided, we can
still implement the unitary $\widetilde{W}$ in Lemma B.3 by using the
following lemma, which allows one to multiply block-encoded matrices.

######  Lemma B.4 (Product of block-encoded matrices [3]).

If $U$ is an $(\alpha,a,\delta)$-block-encoding of an $s$-qubit operator $A$,
and $V$ is an $(\beta,b,\varepsilon)$-block-encoding of an $s$-qubit operator
$B$ then $\left(I_{b}\otimes U\right)\left(I_{a}\otimes V\right)$ is an
$(\alpha\beta,a+b,\alpha\varepsilon+\beta\delta)$-block-encoding of $AB$. Note
that here $I_{a}$ ($I_{b}$) acts on the ancilla qubits of $U$ ($V$).

###### Proof.

$\displaystyle\|AB-\alpha\beta(\bra{0^{a+b}}\otimes I_{s})(I_{b}\otimes
U)(I_{a}\otimes V)(\ket{0^{a+b}}\otimes I_{s})\|$ $\displaystyle=$
$\displaystyle\|AB-\underbrace{\alpha(\bra{0^{a}}\otimes
I_{s})U(\ket{0^{a}}\otimes
I_{s})}_{\widetilde{A}}\underbrace{\beta(\bra{0^{b}}\otimes
I_{s})V(\ket{0^{b}}\otimes I_{s})}_{\widetilde{B}}\|$ $\displaystyle=$
$\displaystyle\|AB-\widetilde{A}B+\widetilde{A}B-\widetilde{A}\widetilde{B}\|$
$\displaystyle=$
$\displaystyle\|(A-\widetilde{A})B+\widetilde{A}(B-\widetilde{B})\|$
$\displaystyle\leq$
$\displaystyle\|A-\widetilde{A}\|\beta+\alpha\|B-\widetilde{B}\|$
$\displaystyle\leq$ $\displaystyle\alpha\varepsilon+\beta\delta.\qed$

###  B.2 Block-encoding of oracle-access matrices

The following lemma can be used to block-encode any matrix whose entries are
provided by an oracle.

######  Lemma B.5 (Block-encoding of oracle-access matrices, adapted from
Lemma 48 of [3]).

Let $A\in\mathbb{C}^{2^{s}\times 2^{s}}$ be a $d_{r}$-row-sparse and
$d_{c}$-column-sparse matrix. Suppose that we have the access to the following
$2(s+1)$-qubit oracles

$\displaystyle\mathcal{O}_{r}:\ket{i}\ket{k}\rightarrow\ket{i}\ket{r_{ik}}$
$\displaystyle 0\leq i<2^{s},0\leq k<d_{r},$
$\displaystyle\mathcal{O}_{c}:\ket{l}\ket{j}\rightarrow\ket{c_{lj}}\ket{j}$
$\displaystyle 0\leq j<2^{s},0\leq l<d_{c},$

where $r_{ik}$ is the index for the $k$-th non-zero entry of the $i$-th row of
A, or if there are less than $k$ non-zero entries, then $r_{ik}=k+2^{s}$, and
similarly $c_{lj}$ is the index for the $l$-th non-zero entry of the $j$-th
column of $A$, or if there are less than $l$ non-zero entries, then
$c_{lj}=l+2^{s}$. Additionally, let $\hat{a}\geq\max_{i,j}|a_{ij}|$ and
suppose the following oracle is provided

$\mathcal{O}_{A}:\ket{i}\ket{j}\ket{0}^{\otimes
b}\rightarrow\ket{i}\ket{j}\ket{\tilde{a}_{ij}},\hskip 14.22636pt0\leq
i,j<2^{s},$

where $\tilde{a}_{ij}$ is the (exact) $b$-qubit description of
$a_{ij}/\hat{a}$, and if $i$ or $j$ is out of range then $\tilde{a}_{ij}=0$.
Then one can implement an $(\hat{a}\sqrt{d_{r}d_{c}},s+3,\varepsilon)$-block-
encoding of $A$ with a single use of $\mathcal{O}_{r},\mathcal{O}_{c}$, two
uses of $\mathcal{O}_{A}$ and additionally
$O(s+\operatorname{polylog}(\frac{\hat{a}\sqrt{d_{r}d_{c}}}{\varepsilon}))$
one- and two-qubit gates while using
$O(b,\operatorname{polylog}(\frac{\hat{a}\sqrt{d_{r}d_{c}}}{\varepsilon}))$
ancilla qubits (which are discarded before the post-selection step).

###### Proof.

Let $G_{d}$ be the $(s+1)$-qubit unitary that performs:
$\ket{0}\rightarrow\frac{1}{\sqrt{d}}\sum_{i=0}^{d-1}\ket{i}$, which can be
implemented with $O(s)$ Hadamard gates (assume $d$ is a power of 2 and note
that $d\leq 2^{s}$). Additionally, define the following $2(s+1)$-qubit
operators: $V_{R}=\mathcal{O}_{c}(G_{d_{c}}\otimes I_{s+1})$ and
$V_{L}=\mathcal{O}_{r}(I_{s+1}\otimes G_{d_{r}})\operatorname{SWAP}_{s+1}$,
where $\operatorname{SWAP}_{s+1}$ swaps two $(s+1)$-qubit registers (which
uses $s$ two-qubit swap gates). For any $0\leq i,j<2^{s}$, we have that

$\displaystyle\bra{0^{s+1}}\bra{i}V_{L}^{\dagger}V_{R}\ket{0^{s+1}}\ket{j}$
$\displaystyle=$
$\displaystyle\left(\frac{1}{\sqrt{d_{r}}}\sum_{k=0}^{d_{r}-1}\bra{i}\bra{r_{ik}}\right)\left(\frac{1}{\sqrt{d_{c}}}\sum_{l=0}^{d_{c}-1}\ket{c_{lj}}\ket{j}\right)$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{d_{r}d_{c}}}\text{ if
}a_{ij}\neq 0\text{ and }0\text{ otherwise.}$

Observe that after applying $V_{R}$, instead of applying $V_{L}^{\dagger}$
right away, we can first append a $b$-qubit register and call
$\mathcal{O}_{A}$ to query the matrix entry, then append the register
$\ket{0}$ and perform the controlled rotation:
$\ket{0}\ket{\tilde{a}_{ij}}\rightarrow(a_{ij}/\hat{a}\ket{0}+\sqrt{1-|a_{ij}/\hat{a}|^{2}}\ket{1})\ket{\tilde{a}_{ij}}$,
and finally call $\mathcal{O}_{A}$ again to uncompute the $b$-qubit ancilla
register. Conditioning on the ancilla qubit being $\ket{0}$, we obtain the
desired block-encoding. We only count $s+3$ qubits as ancilla qubits in the
block-encoding notation because the others can be discarded immediately after
the uncomputation step.

We analyze the circuit complexity of the controlled rotation step above.
Observe that, assuming the $b$-qubit description $\ket{\tilde{a}_{ij}}$ is
exact, if the controlled rotation can be done with an accuracy of
$O(\frac{\varepsilon}{\hat{a}\sqrt{d_{r}d_{c}}})$ then we achieve the overall
accuracy of the block-encoding. This step implements circuits to compute
elementary functions such as the square root and trigonometric functions with
$O(\log(\frac{\hat{a}\sqrt{d_{r}d_{c}}}{\varepsilon}))$ additional qubits.
Implementing these elementary functions on these additional qubits requires
$O(b,\operatorname{polylog}(\frac{\hat{a}\sqrt{d_{r}d_{c}}}{\varepsilon}))$
elementary logic gates [52, 53]. After computing the rotation angles, it takes
$O(\log(\frac{\hat{a}\sqrt{d_{r}d_{c}}}{\varepsilon}))$ controlled quantum
gates to perform the controlled rotation. ∎

We can see that the normalization factor of the above block-encoding depends
on the sparsity of the matrix. Without loss of generality, assume $\hat{a}=1$.
Then, the above lemma when applied to dense matrices produces a block-encoding
with a normalization factor of $O(N)$. This is not optimal (see definition in
Section 3) if the operator norm of the matrix is significantly smaller than
$O(N)$. Hence, we refer to it as the “naive” block-encoding. We provide some
examples in which this lemma works optimally in section C.6.

######  Remark B.6 (“Naive” block-encoding of dense oracle-access matrices
(same as Lemma 3.2)).

In the dense case, the oracles $\mathcal{O}_{r},\mathcal{O}_{c}$ are not
needed; we can simply use $s$ Hadamard gates to query all matrix entries. It
can be easily verified that in this case we can obtain an
$(\hat{a}2^{s},s+1,\varepsilon)$-block-encoding with two uses of
$\mathcal{O}_{A}$ and additionally
$O(s+\operatorname{polylog}(\frac{\hat{a}2^{s}}{\varepsilon}))$ one- and two-
qubit gates while using
$O(b,\operatorname{polylog}(\frac{\hat{a}2^{s}}{\varepsilon}))$ ancilla
qubits.

If the matrix $A$ is rectangular, that is $A\in\mathbb{C}^{2^{s}\times
2^{t}}$, where we assume $s>t$ without loss of generality, we can block-encode
$A$ by applying the above lemmas on the $2^{s}$ by $2^{s}$ matrix
$\widetilde{A}$ defined as

$\widetilde{A}_{ij}=\begin{cases}A_{ij}&\text{if }i<2^{s}\text{ and
}j<2^{t}\\\ 0&\text{otherwise}\end{cases}.$ (29)

Then for an input state $\ket{\psi}$, we can obtain $A\ket{\psi}$ by applying
$\widetilde{A}$ on the state $\ket{0^{s-t}}\otimes\ket{\psi}$.

###  B.3 Block-encoding of low-rank matrices

Recall Lemma 3.4, copied below.

######  Lemma B.7 (Block-encoding of low-rank operators with state preparation
unitaries, inspired by Lemma 1 of [33]).

Let
$A=\sum_{i=0}^{p-1}\sigma_{i}\mathbf{u}_{i}\mathbf{v}_{i}^{\dagger}\in\mathbb{C}^{2^{s}\times
2^{s}}$. Let $r=\lceil\log p\rceil$ and
$\sum_{i=0}^{p-1}|\sigma_{i}|\leq\beta$. Suppose the following $(r+s)$-qubit
unitaries are provided:

$G_{L}:\ket{i}\ket{0^{s}}\rightarrow\ket{i}\ket{\mathbf{u}_{i}},G_{R}:\ket{i}\ket{0^{s}}\rightarrow\ket{i}\ket{\mathbf{v}_{i}},$

where $0\leq i<p$ and $\ket{\mathbf{u}_{i}}$ $(\ket{\mathbf{v}_{i}})$ is the
quantum state whose amplitudes are the entries of $\mathbf{u}_{i}$
$(\mathbf{v}_{i})$. Let $(P_{L},P_{R})$ be a $(\beta,r,\varepsilon)$-state-
preparation-pair for the vector $[\sigma_{0},\ldots,\sigma_{p-1}]$. Then, one
can construct a $(\beta,r+s,\varepsilon)$-block-encoding of $A$ with one use
of each of $G_{L},G_{R}^{\dagger},P_{L}^{\dagger},P_{R}$ and a
$\operatorname{SWAP}$ gate on two $s$-qubit registers.

###### Proof.

Observe that the $(2s+r)$-qubit unitary $U=(P_{L}^{\dagger}\otimes
I_{2s})(G_{R}^{\dagger}\otimes
I_{s})(I_{r}\otimes\operatorname{SWAP}_{s})(G_{L}\otimes I_{s})(P_{R}\otimes
I_{2s})$ is a $(\beta,r+s,\varepsilon)$-block-encoding of $A$. Indeed, for any
$0\leq m,n<2^{s}$ we have that

$\displaystyle(\bra{0^{r+s}}\otimes\bra{m})U(\ket{0^{r+s}}\otimes\ket{n})$
$\displaystyle=$
$\displaystyle\left(\sum_{j=0}^{p-1}c_{j}^{*}\bra{j}\bra{\mathbf{v}_{j}}\bra{m}\right)(I_{r}\operatorname{SWAP}_{s})\left(\sum_{i=0}^{p-1}d_{i}\ket{i}\ket{\mathbf{u}_{i}}\ket{n}\right)$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{p-1}c_{i}^{*}d_{i}\bra{\mathbf{v_{i}}}n\rangle\langle
m\ket{\mathbf{u}_{i}}=\langle
m|\left(\sum_{i=0}^{p-1}c_{i}^{*}d_{i}\ket{\mathbf{u}_{i}}\bra{\mathbf{v_{i}}}\right)|n\rangle,$

and

$\displaystyle\left\|\beta\left(\sum_{i=0}^{p-1}c_{i}^{*}d_{i}\ket{\mathbf{u}_{i}}\bra{\mathbf{v_{i}}}\right)-A\right\|$
$\displaystyle=$ $\displaystyle\left\|\left(\sum_{i=0}^{p-1}(\beta
c_{i}^{*}d_{i}-\sigma_{i})\ket{\mathbf{u}_{i}}\bra{\mathbf{v_{i}}}\right)\right\|$
$\displaystyle\leq$ $\displaystyle\sum_{i=0}^{p-1}|(\beta
c_{i}^{*}d_{i}-\sigma_{i})|\leq\varepsilon.\qed$

The unitaries $G_{R},G_{L}$ can be constructed using existing state
preparation methods based on the structure of the entries in the singular
vectors (e.g., [34] utilized efficiently integrable distributions, [4] used a
QRAM data structure). As an example, when given oracle-access to the entries
of the vectors $\mathbf{u},\mathbf{v}$, one can use the following lemma to
block-encode rank-1 operators.

######  Lemma B.8 (Block-encoding of rank-1 operators with oracle access
entries).

Let $A=\mathbf{u}\mathbf{v}^{\dagger}\in\mathbb{C}^{2^{s}\times 2^{s}}$ be a
rank-1 operator. Suppose the following oracles are provided:

$\mathcal{O}_{\mathbf{u}}:\ket{i}\ket{z}^{b}\rightarrow\ket{i}\ket{z\oplus\tilde{u}_{i}},\mathcal{O}_{\mathbf{v}}:\ket{j}\ket{z}^{b}\rightarrow\ket{j}\ket{z\oplus\tilde{v}_{j}},$

where $\tilde{u}_{i}(\tilde{v}_{j})$ is the $b$-qubit exact description of
$u_{i}(v_{j})$. Let $\max_{i}|u_{i}|\leq\hat{u}$ and
$\max_{j}|v_{j}|\leq\hat{v}$. Then one can obtain a
$(2^{s}\hat{u}\hat{v},s+2b+2,0)$-block-encoding of $A$ with two uses of
$\mathcal{O}_{\mathbf{u}}$, $\mathcal{O}_{\mathbf{v}}$ each.

###### Proof.

Let $\operatorname{CR}_{\hat{u}}$ denote the conditioned rotation operator:
$\operatorname{CR}_{\hat{u}}\ket{\tilde{u}_{i}}\ket{0}=\ket{\tilde{u}_{i}}\left(\frac{u_{i}}{\hat{u}}\ket{0}+\sqrt{1-\left|\frac{u_{i}}{\hat{u}}\right|^{2}}\ket{1}\right)$.
Let $G_{\mathbf{u}}=(\mathcal{O}_{\mathbf{u}}\otimes I)(I_{s}\otimes
CR_{\hat{u}})(\mathcal{O}_{\mathbf{u}}\otimes I)(H^{\otimes s}\otimes
I_{b+1})$. Observe that

$(I_{s}\otimes\bra{0}^{b+1})G_{\mathbf{u}}\ket{0}^{s+b+1}=2^{-s/2}\sum_{i=0}^{2^{s}-1}\frac{u_{i}}{\hat{u}}\ket{i}.$

Define $G_{\mathbf{v}}$ similarly as above. Then we have222Here,
$G_{\mathbf{u}}$ and $G_{\mathbf{v}}$ act on each other register
$\ket{0}^{b+1}$.

$\displaystyle\bra{0}^{s}\bra{m}\bra{0}^{b+1}\bra{0}^{b+1}(G_{\mathbf{v}}^{\dagger}\otimes
I)(I\otimes\operatorname{SWAP}_{s})(G_{\mathbf{u}}\otimes
I)\ket{0}^{s}\ket{n}\ket{0}^{b+1}\ket{0}^{b+1}$
$\displaystyle=\frac{1}{2^{s}}\sum_{i,j=0}^{2^{s}-1}\braket{j}{n}\braket{m}{i}\frac{u_{i}}{\hat{u}}\frac{v_{j}}{\hat{v}}=\frac{1}{2^{s}\hat{u}\hat{v}}u_{m}v_{n}.\qed$

This low-rank block-encoding can also be applied for rectangular matrices.
Suppose $\mathbf{u}\in\mathbb{C}^{2^{m}}$ and
$\mathbf{v}\in\mathbb{C}^{2^{n}}$, where we assume without loss of generality
$m>n$, then we can apply the lemmas on the matrix
$\widetilde{A}=\mathbf{u}(\bra{0}^{m-n}\otimes\mathbf{v}^{\dagger})$.

######  Remark B.9.

The normalization factor $2^{s}\hat{u}\hat{v}$ is to ensure that the block-
encoding is unitary. Thus, Lemma B.8 works best when the entries $u_{i}$ and
$v_{j}$ change slowly. If $\mathbf{u}$ or $\mathbf{v}$ is sparse, we can apply
the techniques in Lemma B.5 to improve the normalization factor of this low-
rank block-encoding. In general, however, the above lemma is rather redundant
since we can directly apply Lemma B.5 naively by constructing the matrix
entry-access oracle from $\mathcal{O}_{\mathbf{u}}$ and
$\mathcal{O}_{\mathbf{v}}$. We simply present this lemma here for
completeness.

###  B.4 Block-encoding of block-sparse matrices (Lemma 3.5)

In this section, we provide a block-encoding procedure for block-sparse
matrices. First we consider the block-diagonal case and later use these
techniques for constructing general block-sparse matrices.

######  Lemma B.10 (Block-encoding of block-diagonal matrices).

Let $\\{A^{j}:0\leq j\leq 2^{t}-1\\}$ be a set of $s$-qubit operators and each
$U^{j}$ be an $(\alpha_{j},a,\varepsilon)$-block-encoding of $A^{j}$. Let
$W=\sum_{j=0}^{2^{t}-1}\ket{j}\bra{j}\otimes U^{j}$ and
$\hat{\alpha}=\max_{j}\alpha_{j}$. Furthermore, suppose that the following
oracle is provided:
$\mathcal{O}_{\alpha}:\ket{j}\ket{z}\rightarrow\ket{j}\ket{z\oplus\tilde{\alpha}_{j}}$,
where $\tilde{\alpha}_{j}$ is the (exact) $b$-qubit description of
$\alpha_{j}$ and $z$ is a $b$-qubit string. Then one can obtain an
$(\hat{\alpha},a+1,2\varepsilon)$-block-encoding of
$A=\sum_{j=0}^{2^{t}-1}\ket{j}\bra{j}\otimes A^{j}$ with two uses of
$\mathcal{O}_{\alpha}$, one use of $W$, and additionally
$O(\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ one- and two-
qubit gates while using
$O(b,\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ ancilla qubits
(which are discarded before the post-selection step).

###### Proof.

Let $\operatorname{CR}$ denote the controlled rotation operator:
$\operatorname{CR}\ket{\tilde{\alpha}_{j}}\ket{0}=\ket{\tilde{\alpha}_{j}}\left(\bar{\alpha}_{j}\ket{0}+\sqrt{1-|\bar{\alpha}_{j}|^{2}}\ket{1}\right)$,
with
$|\bar{\alpha}_{j}-\frac{\alpha_{j}}{\hat{\alpha}}|\leq\frac{\varepsilon}{\hat{\alpha}}$.
Implementing this operator requires
$O(\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ one- and two-
qubit gates while using
$O(b,\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ ancilla qubits
(see explanation in the proof of Lemma B.5). Observe that
$\widetilde{W}=(\mathcal{O}_{\alpha}\otimes
I_{s+a+1})(\operatorname{CR}\otimes I_{s+t+a})(\mathcal{O}_{\alpha}\otimes
I_{s+a+1})(I_{b+1}\otimes W)$ is the desired block-encoding333The subscripts
indicate the registers the identity act on.. Indeed, we have that

$\widetilde{W}(\ket{0}\otimes\ket{0}^{b}\otimes
I_{s+t+a})=\sum_{j=0}^{2^{t}-1}\left(\bar{\alpha}_{j}\ket{0}+\sqrt{1-|\bar{\alpha}_{j}|^{2}}\ket{1}\right)\otimes\ket{0}^{b}\otimes\ket{j}\bra{j}\otimes
U^{j},$

Therefore,

$\displaystyle\left\|\hat{\alpha}\left(\bra{0}^{1+b+a}\otimes
I_{s+t}\right)\widetilde{W}\left(\ket{0}^{1+b+a}\otimes
I_{s+t}\right)-A\right\|$ $\displaystyle=$
$\displaystyle\left\|\bigoplus_{0\leq
j<2^{t}}\left(\hat{\alpha}\bar{\alpha}_{j}\bra{0}^{a}U^{j}\ket{0}^{a}-A^{j}\right)\right\|$
$\displaystyle\leq$ $\displaystyle\left\|\bigoplus_{0\leq
j<2^{t}}\left(\alpha_{j}\bra{0}^{a}U^{j}\ket{0}^{a}-A^{j}\right)\right\|+\left\|\bigoplus_{0\leq
j<2^{t}}(\hat{\alpha}\bar{\alpha}_{j}-\alpha_{j})\bra{0}^{a}U^{j}\ket{0}^{a}\right\|$
$\displaystyle\leq$ $\displaystyle\max_{0\leq
j<2^{t}}\left\|\alpha_{j}\bra{0}^{a}U^{j}\ket{0}^{a}-A^{j}\right\|+\max_{0\leq
j<2^{t}}\left\|(\hat{\alpha}\bar{\alpha}_{j}-\alpha_{j})\bra{0}^{a}U^{j}\ket{0}^{a}\right\|$
$\displaystyle\leq$ $\displaystyle 2\varepsilon.$

We do not include $b$ in the notation of the block-encoding because this
register can actually be taken out of the system right after the uncomputation
step. ∎

In the above lemma, if the normalization factors $\alpha_{j}$ are the same for
all blocks, then the operator $W$ is already an
$(\hat{\alpha},a,\varepsilon)$-block-encoding of the desired block-diagonal
matrix.

Next, we apply the proof techniques above to generalize to the block-sparse
case. The following lemma is the same as Lemma 3.5.

######  Lemma B.11 (Block-encoding of block-sparse matrices).

Let $A=\sum_{i,j=0}^{2^{t}-1}\ket{i}\bra{j}\otimes A^{ij}$ be a $d_{r}$-row-
block-sparse and $d_{c}$-column-block-sparse matrix, where each $A^{ij}$ is an
$s$-qubit operator. Let $U^{ij}$ be an $(\alpha_{ij},a,\varepsilon)$-block-
encoding of $A^{ij}$. Suppose that we have the access to the following
$2(t+1)$-qubit oracles

$\displaystyle\mathcal{O}_{r}:\ket{i}\ket{k}\rightarrow\ket{i}\ket{r_{ik}}\qquad\forall
i\in\\{0,\ldots,2^{t}-1\\},k\in\\{0,\ldots,d_{r}-1\\},$
$\displaystyle\mathcal{O}_{c}:\ket{l}\ket{j}\rightarrow\ket{c_{lj}}\ket{j}\qquad\forall
j\in\\{0,\ldots,2^{t}-1\\},l\in\\{0,\ldots,d_{c}-1\\},$

where $r_{ik}$ is the index for the $k$-th non-zero block of the $i$-th block-
row of A, or if there are less than $k$ non-zero blocks, then
$r_{ik}=k+2^{t}$, and similarly $c_{lj}$ is the index for the $l$-th non-zero
block of the $j$-th block-column of $A$, or if there are less than $l$ non-
zero blocks, then $c_{lj}=l+2^{t}$. Additionally, suppose the following oracle
is provided:

$\displaystyle\mathcal{O}_{\alpha}:\ket{i}\ket{j}\ket{z}\rightarrow\ket{i}\ket{j}\ket{z\oplus\tilde{\alpha}_{ij}},\forall
i,j\in\\{0,\ldots,2^{t}-1\\}$

where $\tilde{\alpha}_{ij}$ is the $b$-qubit description of $\alpha_{ij}$ (if
$i$ or $j$ is out of range then $\tilde{\alpha}_{ij}=0$), and let the
$(2t+2+s+a)$-qubit unitary $W=\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j})\otimes U^{ij}+\left(I-\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j})\right)\otimes I_{s+a}$. Let
$\hat{\alpha}=\max_{i,j}\alpha_{ij}$. Then one can implement an
$(\hat{\alpha}\sqrt{d_{r}d_{c}},t+a+3,2\sqrt{d_{r}d_{c}}\varepsilon)$-block-
encoding of $A$, with one use of each of $\mathcal{O}_{r},\mathcal{O}_{c},$
and $W$, two uses of $\mathcal{O}_{\alpha}$,
$O(t+\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ additional
one- and two-qubit gates, and
$O(b,\operatorname{polylog}(\frac{\hat{\alpha}}{\varepsilon}))$ extra ancilla
qubits (which are discarded before the post-selection step).

###### Proof.

Let $\operatorname{CR}$ be the controlled rotation operator:
$\operatorname{CR}\ket{\tilde{\alpha}_{ij}}\ket{0}=\ket{\tilde{\alpha}_{ij}}\left(\bar{\alpha}_{ij}\ket{0}+\sqrt{1-|\bar{\alpha}_{ij}|^{2}}\ket{1}\right)$,
where
$|\bar{\alpha}_{ij}-\frac{\alpha_{ij}}{\hat{\alpha}}|\leq\frac{\varepsilon}{\hat{\alpha}}$.
Let $G_{d}$ be a $(t+1)$-qubit unitary such that
$G_{d}\ket{0^{t+1}}=\sum_{k=0}^{d-1}\frac{\ket{k}}{\sqrt{d}}$, where $d\leq
2^{t}$. Additionally, let $V_{R}=\mathcal{O}_{c}(G_{d_{c}}\otimes I_{t+1})$
and $V_{L}=\mathcal{O}_{r}(I_{t+1}\otimes
G_{d_{r}})\operatorname{SWAP}_{t+1}$. For any $0\leq i,j<2^{t}$, observe that

$\displaystyle\bra{0^{t+1}}\bra{i}V_{L}^{\dagger}V_{R}\ket{0^{t+1}}\ket{j}$
$\displaystyle=\left(\sum_{k=0}^{d_{r}-1}\bra{i}\frac{\bra{r_{ik}}}{\sqrt{d_{r}}}\right)\left(\sum_{k=0}^{d_{c}-1}\ket{c_{lj}}\frac{\ket{j}}{\sqrt{d_{c}}}\right)$
(30) $\displaystyle=\frac{1}{\sqrt{d_{r}d_{c}}}\text{ if }A^{ij}\neq 0\text{
and }0\text{ otherwise.}$

Let the $(2t+2+s+a)$-qubit unitary $W=\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j})\otimes U^{ij}+(I-\sum_{i,j:A^{ij}\neq
0}(\ket{i}\ket{j}\bra{i}\bra{j}))\otimes I_{s+a}$ and, similarly to Lemma
B.10, $\widetilde{W}=(\mathcal{O}_{\alpha}\otimes
I_{s+a+1})(\operatorname{CR}\otimes I_{s+2t+2+a})(\mathcal{O}_{\alpha}\otimes
I_{s+a+1})(I_{b+1}\otimes W)$ be a $(2t+s+a+b+3)$-qubit unitary 444The
subscripts indicate the registers the identity act on.. Then
$(V_{L}^{\dagger}\otimes I_{s+a+b+1})\widetilde{W}(V_{R}\otimes I_{s+a+b+1})$
is the desired block-encoding. Indeed, this follows directly from Equation 30
and a similar proof to that of Lemma B.10, noting that $\widetilde{W}$ leaves
the $\ket{c_{lj}}$ and $\ket{j}$ registers unchanged, and

$\displaystyle\left\|\sum_{i,j:A^{ij}\neq
0}\ket{i}\bra{j}\otimes(\hat{\alpha}\bar{\alpha}_{ij}\bra{0}^{a}U^{ij}\ket{0}^{a}-A^{ij})\right\|$
$\displaystyle\leq$ $\displaystyle\left\|\sum_{i,j:A^{ij}\neq
0}\ket{i}\bra{j}\otimes(\alpha_{ij}\bra{0}^{a}U^{ij}\ket{0}^{a}-A^{ij})\right\|+\left\|\sum_{i,j:A^{ij}\neq
0}\ket{i}\bra{j}\otimes(\hat{\alpha}\bar{\alpha}_{ij}-\alpha_{ij})\bra{0}^{a}U^{ij}\ket{0}^{a})\right\|$
$\displaystyle\leq$
$\displaystyle\sqrt{d_{r}d_{c}}\cdot\left(\max_{i,j:A^{ij}\neq
0}\left\|\alpha_{ij}\bra{0}^{a}U^{ij}\ket{0}^{a}-A^{ij}\right\|+\max_{i,j:A^{ij}\neq
0}\left\|(\hat{\alpha}\bar{\alpha}_{ij}-\alpha_{ij})\bra{0}^{a}U^{ij}\ket{0}^{a}\right\|\right)$
$\displaystyle\leq$ $\displaystyle 2\sqrt{d_{r}d_{c}}\varepsilon.$

Implementing $G_{d}$ and the controlled rotation operator requires the
indicated additional gate and ancilla qubit complexities (see explanation in
the proof of Lemma B.5). ∎

This lemma does not require sparsity to work. It generalizes Lemma B.5 (also
Lemma 48 of [3]) for block-encoding matrices with oracle access to entries. If
$\alpha_{ij}$ are the same for all blocks, we can omit the oracle
$\mathcal{O}_{\alpha}$ and the controlled rotation step, hence only $t+a+2$
ancilla qubits are needed.

##  Appendix C Detailed analysis of block-encodings of kernel matrices

###  C.1 Circuit complexity analysis of Section 4.1

In this section, we provide detailed analysis on circuit and query
complexities of the block-encoding procedure presented in Section 4.1. Our
problem is to optimally block-encode the kernel matrix
$\mathbf{K}=(k(x_{i},x_{j}))_{i,j=0}^{N-1}$, where $x_{i}=i$, $N=2^{L}$, and
$k(x,x^{\prime})=|x-x^{\prime}|^{-p}$ is a polynomially decaying kernel, given
the oracle

$\mathcal{O}_{k}:\ket{i}\ket{j}\ket{0}^{\otimes
b}\rightarrow\ket{i}\ket{j}\ket{\tilde{k}(x_{i},x_{j})},$ (31)

where $\tilde{k}(x_{i},x_{j})$ is the (exact) $b$-qubit description of
$k(x_{i},x_{j})$.

As a reminder, our procedure uses the hierarchical splitting in Figure 1.
Specifically, at level $\ell$ of the hierarchy, each admissible block has size
$2^{L-\ell}$. In addition, the maximum entry of an admissible block in
$\mathbf{K}^{(\ell)}$ is bounded by $d_{\min}^{-p}=2^{-(L-\ell)p}$. Thus, an
admissible block at level $\ell$ can be
$(\alpha_{\ell},L-\ell+1,\frac{\varepsilon}{3})$-block-encoded via the naive
procedure for dense matrices (Lemma 3.2 or Remark B.6), where
$\alpha_{\ell}=2^{(L-\ell)(1-p)}$. From here, each $\mathbf{K}^{(\ell)}$,
which is a block-sparse matrix of sparsity 3, can be
$(3\alpha_{\ell},L+3,\varepsilon)$-block-encoded via the procedure for block-
sparse matrices (Lemma 3.5). The adjacent interaction part $\mathbf{K}_{ad}$,
which is a tridiagonal matrix, can be $(3,L+3,\varepsilon)$-block-encoded with
Lemma B.5. Finally, we use the procedure for linear combination of block-
encoded matrices (Lemma B.3) to obtain a block-encoding of
$\mathbf{K}=\sum_{\ell=2}^{L}\mathbf{K}^{(\ell)}+\mathbf{K}_{ad}$ as
$U=\sum_{\ell=2}^{L}3\alpha_{\ell}U^{(\ell)}+3U_{ad}$, where the $U$’s denote
the corresponding block-encodings. This step requires a $(\log L)$-qubit state
preparation pair (Definition 3.3) which prepares the coefficients
$[3,3\alpha_{2},\ldots,3\alpha_{L}]$ in the linear combination. This entire
procedure results in a $(\alpha,L+\log L+3,\alpha\varepsilon)$-block-encoding
of the _exact_ kernel matrix $\mathbf{K}$, where,

$\alpha=\begin{cases}3(1+\frac{2^{1-p}-N^{1-p}}{2^{1-p}-4^{1-p}})\leq
O(1)&\text{ if }p>1\\\ 3\log N&\text{ if }p=1\\\
3(1+\frac{N^{1-p}-2^{1-p}}{4^{1-p}-2^{1-p}})&\text{ if }p<1\\\ \end{cases}.$
(32)

In the main text, we have shown that this procedure results in an optimal
block-encoding. That is, in the obtained $(\alpha,\log N+\log\log
N+3,\varepsilon)$-block-encoding, the normalization factor $\alpha$ is exactly
$\Theta(\|\mathbf{K}\|)$. Here, we show why exactly two queries to
$\mathcal{O}_{k}$ and $O(\operatorname{polylog}(\frac{N}{\varepsilon}))$ extra
one- and two-qubit gates are needed by carefully combining Lemma 3.5, Lemma
B.3, and Lemma 3.2.

#### Resources

We use 4 registers $A$ ($\log L$ qubits), $B$ ($L+1$ qubits), $C$ ($L+1$
qubits), $D$ (single qubit). For simplicity, assume $\log L$ is an integer. We
want to construct a $(2L+\log L+3)$-qubit unitary $U$ such that, for any
$0\leq m,n\leq 2^{L}-1$, we have

$\bra{0}_{A}\bra{0}_{B}\bra{m}_{C}\bra{0}_{D}U\ket{0}_{A}\ket{0}_{B}\ket{n}_{C}\ket{0}_{D}=\frac{k_{mn}}{\alpha}.$
(33)

Here, although $0\leq m,n\leq 2^{L}-1$, we allocate $(L+1)$ qubits to the
registers $B,C$ for later use.

#### State preparation pair

Observe that Lemma B.3 requires a state preparation pair $(P_{L},P_{R})$
(Definition 3.3) that prepares the coefficient vector
$\bm{\beta}=[3,3\alpha_{2},\ldots,3\alpha_{L}]$ (the first entry being
$\alpha_{ad}$ of $\mathbf{K}_{ad}$). Note that $\|\bm{\beta}\|_{1}=\alpha$,
the overall normalization factor. These operators are $\log L=\log\log N$
qubits, hence we assume they are provided for now (later we will show how to
efficiently construct them). In particular, we choose $P_{L}=P_{R}$ and

$P_{R}\ket{0}_{A}=\frac{1}{\sqrt{\alpha}}\sum_{\ell=0}^{L-1}\sqrt{\beta_{\ell}}\ket{\ell}_{A}.$
(34)

#### Implementing the “right” part

Suppose the registers are in the state

$\ket{0}_{A}\ket{0}_{B}\ket{n}_{C}\ket{0}_{D}.$ (35)

We apply $P_{R}$ to get (omitting the register names for ease of notation)

$\frac{1}{\sqrt{\alpha}}\sum_{\ell=1}^{L}\sqrt{\beta_{\ell-1}}\ket{\ell-1}\ket{0}\ket{n}\ket{0}.$
(36)

Next, we take care of Lemma 3.5 for the block-sparse levels
$\mathbf{K}^{(\ell)}$. Conditioned on register A being $\ell-1$, where
$2\leq\ell\leq L$, we focus on a particular $\ell$. Define the operator
$D_{3}:\ket{0^{L+1}}\rightarrow\frac{1}{\sqrt{3}}\sum_{c=1}^{3}\ket{0}\ket{c}\ket{0^{L-2}}$.
Apply $D_{3}$ on the register $B$ we get

$\frac{1}{\sqrt{3}}\sum_{c=1}^{3}\ket{\ell-1}\underbrace{\ket{0}\ket{c}\ket{0^{L-2}}}_{\text{register
}B}\ket{n}\ket{0}.$ (37)

Next, we apply the controlled block-row index oracle $\mathcal{O}_{r}$, which,
based on $\ell,c,n$, computes the block-row index $m_{b}$ of the $c$-th non-
zero admissible block in $\mathbf{K}^{(\ell)}$ that intersects column $n$
(there are at most $3$ such blocks since $\mathbf{K}^{(\ell)}$ is 3-block-
sparse) and writes the result $m_{b}$ to the first $\ell+1$ qubits of register
$B$. Note that this computation is easy and reversible. We denote this set of
non-zero block-row index $m_{b}$ as $I(\ell,n)$. If there are less than $c$
such blocks, $\mathcal{O}_{r}$ writes $c+2^{\ell}$ instead (this choice
ensures the oracle is reversible and also the first qubit of $\ket{m_{b}}$
acts as a flag that will be useful later).

$\frac{1}{\sqrt{3}}\sum_{m_{b}\in
I(\ell,n)}\ket{\ell-1}\underbrace{\ket{m_{b}}\ket{0^{L-\ell}}}_{\text{register
}B}\ket{n}\ket{0}.$ (38)

Next, we take care of the steps in Lemma 3.2. Also conditioned on register $A$
being $\ell-1$, where $2\leq\ell\leq L$, we apply Hadamard gates on the last
$L-\ell$ qubits of register $B$

$\frac{1}{\sqrt{3}}\frac{1}{\sqrt{2^{L-\ell}}}\sum_{m_{b}\in
I(\ell,n)}\sum_{m_{i}=0}^{2^{L-\ell}-1}\ket{\ell-1}\underbrace{\ket{m_{b}}\ket{m_{i}}}_{\text{register
}B}\ket{n}\ket{0}.$ (39)

Observe that when the first qubit of $\ket{m_{b}}$ is zero,
$\ket{m_{b}}\ket{m_{i}}$ together form a valid row index $m<2^{L}$. If this
qubit is one, then $\ket{m_{b}}\ket{m_{i}}$ is not a valid row index and we
can simply take $k_{mn}=0$ when querying the matrix entry oracle
$\mathcal{O}_{k}$. In particular, we query the oracle $\mathcal{O}_{k}$ on
registers $B,C$ to obtain $\ket{\tilde{k}_{mn}}$ in another appended $b$-qubit
register. Then, we perform the controlled rotation to rotate register $D$ as:
$\ket{0}\rightarrow\left(\frac{k_{mn}}{k^{(\ell)}_{\max}}\ket{0}+\sqrt{1-\left|\frac{k_{mn}}{k^{(\ell)}_{\max}}\right|^{2}}\ket{1}\right)$,
where $k^{(\ell)}_{\max}$ is the maximum entry in the admissible blocks at
level $\ell$, which we have shown to be $2^{-(L-\ell)p}$. Finally, we call
$\mathcal{O}_{k}$ one more time to uncompute and discard the appended
register. This gives

$\frac{1}{\sqrt{3}}\frac{1}{\sqrt{2^{L-\ell}}}\sum_{m_{b}\in
I(\ell,n)}\sum_{m_{i}=0}^{2^{L-\ell}-1}\ket{\ell-1}\underbrace{\ket{m_{b}}\ket{m_{i}}}_{\text{register
}B}\ket{n}\left(\frac{k_{mn}}{k^{(\ell)}_{\max}}\ket{0}+\sqrt{1-\left|\frac{k_{mn}}{k^{(\ell)}_{\max}}\right|^{2}}\ket{1}\right)$
(40)

for $2\leq\ell\leq L$.

On the other hand, conditioned on $\ell=1$ (register $A$ being $\ket{0}$), we
take care of the adjacent part $\mathbf{K}_{ad}$ rather simply. Apply the
$(2L+2)$-qubit operator $\mathcal{O}_{ad}$ which performs the map
$\ket{0}\ket{n}\rightarrow\frac{1}{\sqrt{3}}\sum_{i=-1}^{1}\ket{n+i}\ket{n}$
for $0\leq n<2^{L}$ (if $n+i=-1$ then it writes to the first register any
value that is not a valid row index, e.g., $2^{L+1}-1$) on registers $B,C$.
Then, we perform the oracle query, controlled rotation, and uncomputation
similarly as done above to obtain

$\frac{1}{\sqrt{3}}\sum_{i=-1}^{1}\ket{0}_{A}\ket{n+i}_{B}\ket{n}_{C}\left(k_{n+i,n}\ket{0}+\sqrt{1-|k_{n+i,n}|^{2}}\ket{1}\right),$
(41)

where $k_{n+i,n}=0$ if $n+i$ is an invalid row index.

To summarize, the entire procedure from Equation 36 to Equation 40 (and
Equation 41), denoted by $U_{R}$, performs the following

$\displaystyle U_{R}P_{R}\ket{0}\ket{0}\ket{n}\ket{0}=$
$\displaystyle\frac{1}{\sqrt{3\alpha}}\left[\sqrt{\beta_{0}}\sum_{i=-1}^{1}\ket{0}\ket{n+i}\ket{n}\left(k_{n+i,n}\ket{0}+\sqrt{1-|k_{n+i,n}|^{2}}\ket{1}\right)\right.$
(42)
$\displaystyle\left.+\sum_{\ell=2}^{L}\frac{\sqrt{\beta_{\ell-1}}}{\sqrt{2^{L-\ell}}}\sum_{m_{b}\in
I(\ell,n)}\sum_{m_{i}=0}^{2^{L-\ell}-1}\ket{\ell-1}\underbrace{\ket{m_{b}}\ket{m_{i}}}_{\text{equiv.
to some
}m}\ket{n}\left(\frac{k_{mn}}{k^{(\ell)}_{\max}}\ket{0}+\sqrt{1-\left|\frac{k_{mn}}{k^{(\ell)}_{\max}}\right|^{2}}\ket{1}\right)\right],$

where $\beta_{0}=3$ is the normalization factor of the adjacent level
$\mathbf{K}_{ad}$ and $\beta_{\ell-1}=3\alpha_{\ell}$ is the normalization
factor of the hierarchical level $\mathbf{K}^{(\ell)}$. Importantly, notice
that register $B$ has full support on all $2^{L}$ possible values of row index
$m$ in this superposition.

#### Implementing the “left” part

We would like to transform the state

$\ket{0}\ket{0}\ket{m}\ket{0}$ (43)

to obtain similar result to that of the “right” procedure above. The procedure
follows the same steps, except

  * •

$P_{R}$ needs to be replaced by $P_{L}$.

  * •

$\mathcal{O}_{r}$ needs to be replaced by $\mathcal{O}_{c}$, the controlled
block-column index oracle which, based on $\ell,m$, computes the block-column
indices $n_{b}$ of the non-zero admissible blocks in $\mathbf{K}^{(\ell)}$
that intersect row $m$ (there are at most $3$ such blocks since
$\mathbf{K}^{(\ell)}$ is 3-block-sparse) and writes the result $n_{b}$ to the
first $\ell+1$ qubits of register $B$. We denote this set of non-zero block-
column indices $n_{b}$ as $J(\ell,m)$.

  * •

There is no need to query the entry oracle $\mathcal{O}_{k}$ (nor the
accompanied controlled rotation).

  * •

Registers $B$ and $C$ are swapped at the end.

Using similar notation to that of Equation 42, the above modifications give

$\displaystyle\operatorname{SWAP}_{B,C}U_{L}P_{L}\ket{0}\ket{0}\ket{m}\ket{0}=\frac{1}{\sqrt{3\alpha}}$
$\displaystyle\left[\sqrt{\beta_{0}}\sum_{i=-1}^{1}\ket{0}\ket{m}\ket{m+i}\ket{0}\right.$
(44)
$\displaystyle\left.+\sum_{\ell=2}^{L}\frac{\sqrt{\beta_{\ell-1}}}{\sqrt{2^{L-\ell}}}\sum_{n_{b}\in
J(\ell,m)}\sum_{n_{j}=0}^{2^{L-\ell}-1}\ket{\ell-1}\ket{m}\underbrace{\ket{n_{b}}\ket{n_{j}}}_{\text{equiv.
to some }n}\ket{0}\right].$

#### Combining both parts

Combining Equation 42 with Equation 44 and noting that
$\beta_{0}=3,\beta_{\ell-1}=3\cdot 2^{(L-\ell)(1-p)}$ and
$k^{(\ell)}_{\max}=2^{-(L-\ell)p}$, we get

$\bra{0}\bra{0}\bra{m}\bra{0}P_{L}^{\dagger}U_{L}^{\dagger}\operatorname{SWAP}_{B,C}U_{R}P_{R}\ket{0}\ket{0}\ket{n}\ket{0}=\frac{k_{mn}}{\alpha},$
(45)

which is exactly what we set out to prove (Equation 33).

#### Circuit complexity

As seen above, this block-encoding only uses two calls to $\mathcal{O}_{k}$,
and one call to each of $\mathcal{O}_{r}$, $\mathcal{O}_{c}$ and $P_{R}$,
$P_{L}$. To achieve an overall error bound of $\varepsilon$, it suffices for
the controlled rotation on $\ket{\tilde{k}_{mn}}\ket{0}$ to have error
$O(\frac{\varepsilon}{N})$. This can be done with
$O(b,\operatorname{polylog}(\frac{N}{\varepsilon}))$ extra qubits and one- and
two-qubit gates (see explanation in proof of Lemma B.5).

#### Constructing state preparation pair

$P_{L}$ and $P_{R}$ act on only $\log\log N$ qubits, hence were assumed given
as unitaries. If this is not the case, we can still construct them (Equation
34) as block-encoded operators with only an extra normalization factor of
$\log N$ as follows. Let

$\displaystyle\mathcal{P}:\ket{0^{\log L}}\ket{0}$
$\displaystyle\rightarrow\frac{1}{\sqrt{L}}\sum_{\ell=0}^{L-1}\ket{\ell}\ket{0}$
(46)
$\displaystyle\rightarrow\frac{1}{\sqrt{L}}\sum_{\ell=0}^{L-1}\ket{\ell}\ket{\tilde{\beta}_{\ell}}\ket{0}\qquad\text{(query
$\beta_{\ell}$ into an appended register)}$
$\displaystyle\rightarrow\frac{1}{\sqrt{L}}\sum_{\ell=0}^{L-1}\ket{\ell}\left(\sqrt{\frac{\beta_{\ell}}{\hat{\beta}}}\ket{0}+\sqrt{1-\left|\frac{\beta_{\ell}}{\hat{\beta}}\right|}\ket{1}\right),$

where $\beta_{0}=3$ and $\beta_{\ell-1}=3\alpha_{\ell}=3\cdot
2^{(L-\ell)(1-p)}$ for $2\leq\ell\leq L$ and
$\hat{\beta}=\max_{\ell}|\beta_{\ell}|$. As explained in Lemma B.5, the
controlled rotation can be implemented with error bound $\varepsilon$ using
$O(\operatorname{polylog}(\frac{L}{\varepsilon}))$ extra one- and two-qubit
gates. Comparing with Equation 34 we see that $\mathcal{P}$ is a
$(\sqrt{\frac{\hat{\beta}L}{\alpha}},1,\varepsilon)$-block-encoding of
$P_{R}$. Thus, using Lemma B.4 (multiplying block-encoded matrices) we can
achieve Equation 45 with just two more ancilla qubits. The RHS of Equation 45
then becomes
$\frac{k_{mn}}{\alpha}\left(\frac{\hat{\beta}L}{\alpha}\right)^{-1}$. In other
words, the normalization factor will become a factor of
$\frac{\hat{\beta}L}{\alpha}$ larger than the optimal block-encoding. We can
bound this extra factor rather easily. As a reminder, $\alpha$ is given in
Equation 32.

  * •

For $p>1$: $\hat{\beta}=3$ and $\alpha\geq 3$, thus
$\frac{\hat{\beta}L}{\alpha}\leq L=\log N$.

  * •

For $p=1$: $\hat{\beta}=3$ and $\alpha=3\log N$, thus
$\frac{\hat{\beta}L}{\alpha}=1$.

  * •

For $p<1$: $\hat{\beta}=3\cdot 2^{(L-2)(1-p)}=O(N^{1-p})$ and
$\alpha=\Omega(N^{1-p})$, thus $\frac{\hat{\beta}L}{\alpha}\leq O(L)=O(\log
N)$.

Therefore, this extra factor only increases the runtimes of Lemma 5.1
(applying block-encoded matrices) and Lemma 5.2 (applying inverse of block-
encoded matrices) by a factor of $O(\log N)$.

###  C.2 Variants of hierarchical splitting

The form of the hierarchical splitting can be slightly modified for other
types of decaying kernels. For example, consider kernels of the form
$k(x,x^{\prime})=\frac{1}{|(x-x^{\prime})+c|^{p}}$ on the domain
$\Omega=[0,N-1)$ and $-N<c<N$ is an integer. In this case, we can use a
“shifted” hierarchical splitting to obtain an optimal block-encoding of the
kernel matrix $\mathbf{K}=(|i-j+c|^{-p})_{i,j=0}^{N-1}$, in which the
hierarchy is shifted in the horizontal direction (along the rows, see Figure
3.) Similarly, for kernels of the form
$k(x,x^{\prime})=\frac{1}{(|x-x^{\prime}|-c)^{p}}$ ($0\leq c<N$), we can apply
a hierarchical splitting shifted along the skew-diagonal direction, such that
the upper and lower triangular halves of the matrix are shifted in opposite
directions.

Figure 3: Hierarchical splitting corresponding to the kernel
$k(x,x^{\prime})=\frac{1}{|x-x^{\prime}+4|^{p}}$ with uniform distribution of
particles, which is obtained by shifting the splitting in Figure 1C 4 units to
the right. Note that the fine splittings in the first 4 columns are actually
only needed when the sum $x-x^{\prime}+4$ is modulo $N$.

###  C.3 Adaptive algorithm for non-uniform particle distribution

In particle simulation tasks, if the distribution of particles is not
approximately uniform, an adaptive hierarchical splitting can be employed. Let
the smallest interparticle spacing be $N^{-c}$, where $N$ is the number of
particles. Then, we can scale the mesh such that the interparticle distance is
at least $1$ and the mesh size is $N^{c}$. The added mesh sites are treated as
zero-mass particles. The kernel matrix has size $N^{c}\times N^{c}$ and the
hierarchical splitting constructed on this mesh has $c\log N$ levels. The
entire block-encoding procedure presented in Section 4 applies, with only an
extra factor of $\operatorname{poly}(c)$ in the normalization factor and the
required resources of the block-encoding. The optimality of the normalization
factor of this adaptive block-encoding is harder to state precisely, since it
depends on the distribution of particles which is not determined a priori.

###  C.4 Generalization of quantum hierarchical block-encoding

Here, we provide a block-encoding procedures in which one has direct access to
the structure of the entry magnitudes, thus generalizing hierarchical block-
encoding procedure presented in Section 4.1 and Section C.1.

First, we present a procedure for preparing a quantum state
$\ket{x}\in\mathbb{C}^{N}$, where $N=2^{L}$, with high probability of success.
For any $0\leq\ell<L-1$, let $I_{\ell}=\\{i:2^{-(\ell+1)}<|x_{i}|\leq
2^{-\ell}\\}$ and $I_{L-1}=\\{i:|x_{i}|\leq 2^{-(L-1)}\\}$. Also, let
$n_{\ell}=|I_{\ell}|$. Apart from the oracle
$\mathcal{O}_{x}:\ket{i}\ket{0}\rightarrow\ket{i}\ket{\tilde{x}_{i}}$, our
procedure assumes the following oracle is provided:

$\mathcal{O}_{index}:\ket{\ell}\ket{j}\rightarrow\ket{0}\ket{f(\ell,j)},\qquad
0\leq j<n_{\ell},$ (47)

where the first register is $\log L$ qubits and the second register is $L=\log
N$ qubits, and $f(\ell,j)$ is an efficiently computable and _reversible_
function that computes the index of the $j$-th element in the set $I_{\ell}$
(according to some fixed ordering, e.g. top to bottom if $\ket{x}$ is treated
as a column vector) if $0\leq j<n_{\ell}$, or else the oracle writes
$N+(j-n_{\ell}+1)+\sum_{l=0}^{\ell-1}(N-n_{l})$ as a bit string on the
concatenation of the two registers (so that the first register will be non-
zero). This ensures the oracle is reversible. An example of such efficient and
reversible $f(\ell,j)$ is in regular $\mathcal{H}$-matrix splitting. Our
procedure is as follows:

$\displaystyle\small\ket{0^{\log L}}\ket{0}$
$\displaystyle\rightarrow\sum_{\ell=0}^{L-1}\frac{\sqrt{n_{\ell}}2^{-\ell}}{\beta}\ket{\ell}\ket{0},\qquad\text{where
}\beta=\left(\sum_{\ell}n_{\ell}2^{-2\ell}\right)^{1/2}$
$\displaystyle\rightarrow\sum_{\ell=0}^{L-1}\frac{\sqrt{n_{\ell}}2^{-\ell}}{\beta}\ket{\ell}\sum_{j=1}^{n_{\ell}}\frac{\ket{j}}{\sqrt{n_{\ell}}}\qquad\text{(Hadamards
controlled on $\ell$)}$
$\displaystyle\rightarrow\sum_{\ell=0}^{L-1}\sum_{j=1}^{n_{\ell}}\frac{2^{-\ell}}{\beta}\ket{0}\ket{f(\ell,j)}\qquad\text{(call
oracle $\mathcal{O}_{index}$)}$
$\displaystyle\rightarrow\sum_{\ell,j}\frac{2^{-\ell}}{\beta}\ket{0}\ket{f(\ell,j)}\left(\frac{x_{f(\ell,j)}}{2^{-\ell}}\ket{0}+\sqrt{1-\left|\frac{x_{f(\ell,j)}}{2^{-\ell}}\right|^{2}}\ket{1}\right)$
(query $x_{f(\ell,j)}$ and rotate ancilla qubit)

Above, we have ignored the terms when $j$ is invalid ($j\geq n_{\ell}$) in the
sum because we can simply set $x_{f(\ell,j)}=0$ for those cases. Furthermore,
observe that the superposition has support on all indices. Thus, conditioning
on the ancilla being zero, we obtain the desire state $\ket{x}$. The
probability of successfully post-selecting is $\frac{1}{\beta^{2}}$, where

$\beta^{2}=\sum_{\ell=0}^{L-1}n_{\ell}2^{-2\ell}=\sum_{\ell=0}^{L-2}n_{\ell}2^{-2\ell}+n_{L-1}2^{-2(L-1)}\leq
4\sum_{i}|x_{i}|^{2}+\frac{4}{N}\leq 8.$ (48)

We also assumed that the first operation in the procedure, which maps
$\ket{0^{\log L}}$ to a weighted superposition over $\ket{\ell}$, can be done
perfectly. In fact, it can only be block-encoded with an extra normalization
factor of $\sqrt{L}$: for a vector $\ket{y}\in\mathbb{C}^{L}$ one can simply
do

$\displaystyle\ket{0^{\log L}}\ket{0}$
$\displaystyle\rightarrow\sum_{\ell=0}^{L-1}\frac{1}{\sqrt{L}}\ket{\ell}\ket{0},$