††thanks: corresponding author

# Approximating Hamiltonian dynamics with the Nyström method

Alessandro Rudi<EMAIL_ADDRESS>INRIA - Sierra project team, Paris,
France  Leonard Wossnig<EMAIL_ADDRESS>Department of Computer
Science, University College London, London, United Kingdom and Rahko Ltd.,
London, United Kingdom  Carlo Ciliberto<EMAIL_ADDRESS>Department of
Computer Science, University College London, London, United Kingdom  Andrea
Rocchetto<EMAIL_ADDRESS>Department of Computer Science, University of
Texas at Austin, Austin, United States, Department of Computer Science,
University of Oxford, Oxford, United Kingdom, and Department of Computer
Science, University College London, London, United Kingdom  Massimiliano
Pontil<EMAIL_ADDRESS>Computational Statistics and Machine
Learning, IIT, Genoa, Italy  Simone Severini<EMAIL_ADDRESS>Department
of Computer Science, University College London, London, United Kingdom

###### Abstract

Simulating the time-evolution of quantum mechanical systems is BQP-hard and
expected to be one of the foremost applications of quantum computers. We
consider classical algorithms for the approximation of Hamiltonian dynamics
using subsampling methods from randomized numerical linear algebra. We derive
a simulation technique whose runtime scales polynomially in the number of
qubits and the Frobenius norm of the Hamiltonian. As an immediate application,
we show that sample based quantum simulation, a type of evolution where the
Hamiltonian is a density matrix, can be efficiently classically simulated
under specific structural conditions. Our main technical contribution is a
randomized algorithm for approximating Hermitian matrix exponentials. The
proof leverages a low-rank, symmetric approximation via the Nyström method.
Our results suggest that under strong sampling assumptions there exist
classical poly-logarithmic time simulations of quantum computations.

##  1 Introduction

Special purpose quantum simulators permit the efficient implementation of
unitary dynamics governed by physically meaningful families of Hamiltonians,
while the general task is BQP-hard – since we can implement any quantum
computation by a sequence of Hamiltonian evolutions. The properties
responsible for efficiency correspond to plausible structural restrictions
such as locality of interaction between subsystems [Llo96], or sparsity, which
permits simulation complexities sub-logarithmic in the inverse error [AT03].
Moreover, sparse operators provide examples of a class of quantum circuits
with efficient weak classical simulation [SN13], whereas general strong
simulation is known to be #P -hard [Van10]. The intuition supporting this
widely used empirical phenomenon is based on the fact that sparse matrices
have good storage requirements and an _easier_ combinatorial structure
mirroring the small number of physical interactions.

In this paper, we consider the problem of classical simulation of quantum
Hamiltonian dynamics. In particular, we will be interested in the case of
time-independent Hamiltonians. In this setting, the problem is mathematically
equivalent to the task of approximating the matrix exponential of an Hermitian
matrix, i.e. the Hamiltonian of the system.

The task is particularly relevant for quantum many-body physics where the
exponential scaling of the wave function limits to only a handful of qubits
the applicability of ordinary differential equation solvers. In order to
simulate bigger system it is thus necessary to leverage specific properties of
the system. An important class of methods exploits the entanglement structure
of the states to obtain efficient representations, known as tensor networks
[VMC08, Orú14], that can then be evolved either via Trotterisation [Vid04] or
with a time-dependent variational principle [VGC04, Hae+11].

We conclude this introduction by noting that, from a computational
perspective, the problem of simulating Hamiltonian dynamics is related, but
different, to the problem of simulating quantum circuits. In the latter case
the task no longer involves the computation of costly matrix exponentials and
reduces to approximating the action of a circuit, described as a sequence of
unitary operations, on a state vector. Recent techniques developed for this
problem involve the notion of stabilizer rank [Bra+18] or neural network
quantum states [JBC18].

###  1.1 Results

We give a classical, randomized algorithm for the strong simulation of quantum
Hamiltonian dynamics. Recall that in a strong quantum simulation one computes
the amplitude of a particular outcome, whereas in a weak simulation one only
samples from the output distribution of a quantum circuit. While the first is
known to have unconditional and exponential lower bounds [HNS18], and it is
therefore hard for both classical and quantum computers, the latter can be
performed efficiently by a quantum computer (for circuits of polynomial size).
In this sense, our approach is tackling a problem that is harder than the one
solved by quantum computers (which can only sample from the probability
distribution induced by measurements on the output state). Surprisingly, we
are still able to find cases in which the time evolution can be simulated
efficiently.

More specifically, given an efficient classical description of a Hamiltonian
$H$ and of an $n=\log_{2}N$-qubit quantum state $\psi$, we propose conditions
under which the evolution can be simulated efficiently. The algorithm takes as
input an index $i\in\\{1,\dots,2^{n}\\}$ and returns the complex amplitude
corresponding to the projection of the state vector in the $i$-th element of
the computational basis.

Our algorithm relies on an input model that is commonly used in quantum
simulation (condition 1 and 2) supplemented by a condition that is frequently
encountered in the randomised numerical linear algebra literature (condition
3). More specifically:

  1. 1.

We require $H$ to be row-computable. This condition states that every row of
$H$ must have at most a number $s=O(\polylog(N))$ of non-zero entries.
Furthermore, there exists a classical efficient algorithm that, given a row-
index $i$, outputs a list of the non-zero entries of the row.

  2. 2.

We require the initial state $\psi$ to have have at most $q=O(\polylog(N))$
non-zero entries. Furthermore, there exists a classical efficient algorithm
that outputs a list of the non-zero entries.

  3. 3.

We require $H$ to be efficiently row-searchable. This condition informally
states that we can efficiently sample randomly selected indices of the rows of
$H$ in a way proportional to the norm of the row (general case) or the
diagonal element of the row (positive semidefinite case).

It is important to remark that, although the notion of row-searchability is
commonly assumed to hold in the randomised numerical linear algebra literature
(see for example Section $4$ in the seminal paper by Frieze, Kannan, and
Vempala [FKV04]), in the context of quantum systems this is a much stronger
requirement and cannot be assumed to hold for every matrix. Indeed, whilst for
a polynomially sized matrix it is always possible to evaluate all the row-
norms efficiently in the number of non-zero entries of the matrix, this is
generally no longer the case for the exponentially sized matrices found in
quantum mechanics.

Under this input model we prove the following theorem that here we state
informally.

######  Theorem 1 (Informal statement of the main result).

If $H$ is a row-computable, row-searchable Hamiltonian on $n$ qubits, and if
$\psi$ is an $n$-qubit quantum state with an efficient classical description
then, there exists an algorithm that, with probability $1-\delta$,
approximates any chosen amplitude of the state $e^{iHt}\psi$ in time

$O\left(sq+\frac{t^{9}\lVert H\rVert^{4}_{F}\lVert
H\rVert^{7}}{\epsilon^{4}}\left(n+\log\frac{1}{\delta}\right)^{2}\right),$

where $\epsilon$ determines the quality of the approximation,
$\lVert\cdot\rVert_{F}$ is the Frobenius norm, $\lVert\cdot\rVert$ is the
spectral norm, $s$ is the maximum number of non-zero elements in the rows of
$H$, and $q$ is the number of non-zero elements in $\psi$.

By analyzing the dependency of the runtime on the Frobenius norm we can
determine under which conditions we can obtain efficient Hamiltonian
simulations. Informally, we obtain:

######  Corollary 2 (Informal).

If $H$ is a row-computable, row-searchable Hamiltonian on $n$ qubits such that
$\lVert H\rVert^{2}_{F}-\frac{1}{N}\mathrm{Tr}\left(H\right)^{2}\leq
O(\polylog(N))$, and if $\psi$ is an $n$-qubit quantum state with an efficient
classical description, then there exists an efficient algorithm that
approximates any chosen amplitude of the state $e^{iHt}\psi$.

The algorithm proceeds by performing a two steps approximation. First, the
Hamiltonian $H$ is approximated in terms of a low rank operator $\widehat{H}$
which is more amenable to computations by sampling the rows with a probability
proportional to the row norm. Second, the time evolution
$e^{i\widehat{H}t}\psi$ is approximated by a truncated Taylor expansion of the
matrix exponential. By combining the structure of $\widehat{H}$ and the
spectral properties of the truncated exponential we are able to guarantee that
the above procedure can be efficiently performed.

In a more specific way, we separately consider the case of a generic
Hamiltonian $H$ and the restricted case of positive semidefinite Hamiltonians,
for which a refined analysis is possible. The proposed algorithm leverages on
a low-rank approximation of the Hamiltonian $H$ to efficiently approximate the
matrix exponential $e^{iHt}$. Such approximation is performed by randomly
sampling $M=O(\polylog(N))$ rows according to the distribution determined by
the row-searchability condition and then collating them in a matrix
$A\in\mathbb{C}^{M\times N}$.

When $H$ is positive semidefinite, we consider the approximation
$\widehat{H}=AB^{+}A^{*}$, where $B^{+}$ is the pseudoinverse of the positive
semidefinite matrix $B\in\mathbb{C}^{M\times M}$ obtained by selecting the
rows of $A$ whose indices correspond to those originally sampled for the rows
of $H$. The approximation of the time evolution $e^{i\widehat{H}t}\psi$ is
then performed in terms of a Taylor expansion of the matrix exponential
function, truncated to the $K$-th order. Leveraging the structure of
$\widehat{H}$ it is possible to formulate the truncated expansion only in
terms of linear operations involving the matrices $A^{*}A$, $B^{+}$ and $B$
and the vector $A^{*}\psi$. Under the row-computability assumption, all these
operations can be performed efficiently.

In the general Hermitian setting, rows of $H$, that are sampled to form $A$,
are first rescaled according to their sampling probability. Differently from
the positive semidefinite case, we resort to a slightly more involved variant
of the approximation of the matrix exponential, in particular we use
$\widehat{H}^{2}:=AA^{*}$ to approximate $H^{2}$. This approach involves two
auxiliary functions in which the matrix exponential is decomposed and for each
of these two functions we evaluate its truncated Taylor expansion. We show
that this approximation can be formulated only in terms of linear operations
involving $A^{*}A$ and $A^{*}\psi$. Again, under the row-computability
assumption, these operations can be performed efficiently.

###  1.2 Applications and outlook

It remains an open question—and an interesting research direction—to determine
if there exist physically relevant families of Hamiltonians that respect the
requirements for efficient simulation outlined in this paper. In general
terms, our methods appear to function well in cases of “small variations” (for
example, if all the eigenvalues of the Hamiltonian are contained in an
exponentially small band) and in this sense share some similarities with
perturbation theory.

As a specific application of our theorem, we propose the case of sample based
Hamiltonian simulation, that is the simulation of quantum dynamics where the
Hamiltonian is a density matrix [LMR14]. This type of simulation has recently
found applications in various quantum algorithms for machine learning tasks
such as linear regression [SSP16] and least squares support vector machines
[RML14]. Note that when these algorithms are used to analyze classical data,
they assume that the data can be efficiently encoded into a density matrix.
Specifically, as a direct consequence of our main theorem, we obtain the
following result:

######  Theorem 3 (Informal).

If $\rho$ is a row-computable, row-searchable density matrix and if $\psi$ is
an $n$-qubit quantum state with an efficient classical description, then there
exists an efficient algorithm that approximates any chosen amplitude of the
state $e^{i\rho t}\psi$.

We remark that all our results do not require the Hamiltonian to be sparse
(that is, to have at most $\polylog(N)$ non-zero entries) but only to be row-
sparse (that is, to have at most $\polylog(N)$ non-zero entries _per row_).
The latter requirement is compatible with matrices that are non-sparse.

Going beyond immediate applications, we believe that the introduction of
randomised numerical algebra techniques in quantum mechanics may provide a new
direction from which to tackle the quantum many-body problem. In this spirit,
it is relevant to mention that shortly after our preprint was posted on arXiv,
Tang [Tan19] showed that, using a strong sampling data structure—and with a
significantly higher polynomial overhead—it is possible to derive a classical
poly-logarithmic time algorithm for recommendation problems that nearly
matches the asymptotic scaling of the quantum algorithm proposed by Kerenidis
and Prakash [KP16]. The work by Tang assumes the ability to sample from the
probability distribution defined by the squared entries of a matrix divided by
the $\ell_{2}$ norm of the matrix, the so called $\ell_{2}$-norm sampling. Our
row-searchable condition is fundamentally equivalent to this and both our
results provide evidence that a careful assessment of the state preparation
conditions is fundamental in order to determine whether a quantum algorithm
provides a true advantage over a classical variant.

More specifically, [Tan19] defines the sample access to the input in the
following way;

######  Definition 1 ([Tan19]).

We have sample access to data point $x\in\mathbb{C}^{N}$ if, given an index
$i\in[N]$, we can produce independent random samples $i\in[N]$ where $i$ is
sampled with probability $|x_{i}|^{2}/\lVert x\rVert^{2}$.

This definition can be easily extended to the columns (or rows) $H_{:,i}$
($H_{i,:}$) of a Hamiltonian $H$. In this case, we sample a column with
probability $p(i)=\lVert H_{:,i}\rVert^{2}/\lVert H\rVert_{F}^{2}$ (similarly
for a row). The sampling scheme we use for Hermitian matrices is equivalent to
this one (when the matrix is also PSD we use a more computationally efficient
variant). Additionally, we also provide a method based on binary-trees to
compute the sampling probabilities if these are not given by an oracle but
stored in a standard memory model.

Informally, we use the following sample access model;

######  Definition 2 (Informal statement of the sampling access model).

We have sample access to the Hermitian matrix $H\in\mathbb{C}^{N\times N}$ if,
given an index $i\in[N]$, we can randomly sample independent row indices
$i\in[N]$ according to probability $p(i)=\lVert H_{:,i}\rVert^{2}/\lVert
H\rVert_{F}^{2}$ in time $O(\mathrm{poly}(\log(N)))$ for a polynomially
sparse, or structured matrix $H$.

Definition 1 and Definition 2 are fundamentally the same. The main difference
is that, while we use a traditional memory structure and provide a fast way to
calculate the marginals using a binary tree (see Sec. 3), [Tan19] assumes a
memory structure which allows one to sample efficiently according to this
distribution.

Finally, it is relevant to mention that the techniques introduced by Tang have
been directly applied to the problem of classical Hamiltonian simulation in
[Chi+19]. While the time complexity of this algorithm has a $36$-th power
dependency on the Frobenius norm of the Hamiltonian, our algorithm scales with
a $4$-th power.

###  1.3 Related work

####  1.3.1 Classical and quantum approximation of matrix exponentials

The problem of matrix exponentiation has been extensively studied in the
linear algebra literature [Hig05, Hig09, AH09, AH11]. Presently, for the case
of general Hermitian matrices, there is no known algorithm with a runtime
logarithmic in the dimension of the input matrix. Such fast scaling is
required for the simulation of quantum dynamics where the dimensions of the
matrix that describes the evolution scale exponentially with the number of
qubits in the system.

The results described in this paper are based on randomized numerical linear
algebra techniques. These methods, along with results from spectral graph
theory, have lead to a variety of new classical algorithms to approximate
matrix exponentials [DKM06, Dri+11, Mah+11, Woo+14, RCR15]. For specific types
of matrices, these techniques give efficient runtimes. For example, Orrecchia
et al. [OSV12] have demonstrated that the spectral sparsifiers of Spielmann
and Teng [ST04, ST11] can be used to approximate exponentials of strictly
diagonally dominant matrices in time almost linear in the number of non-zero
entries of $H$.

Quantum computers on the other hand can approximate efficiently some kinds of
matrix exponentials. In particular, there exist time-efficient quantum
algorithms for simulating the dynamic of row-sparse Hamiltonians that have
only a linear dependency in the row-sparsity. For an important class of
algorithms the simulation exploits an efficient edge-coloring of the graph
associated with the Hamiltonian matrix $H$ [Chi+03, CK11]. Once this edge
coloring is found, the Hamiltonian can be decomposed into a sum of sparse
Hamiltonians assuming that $H$ is sparse. It is known that these terms can be
simulated separately using the Trotter-Suzuki formula [Tro59, Suz76]. Improved
methods have been developed which result in algorithms with a runtime of
$\tilde{O}(s\ \poly(\log(Ns/\epsilon)))$ where, $s$ is the sparsity, $N$ the
dimension and $\epsilon$ the maximum error in the solution [CW12, BCK15,
LC17].

####  1.3.2 Randomized numerical linear algebra and Nyström methods

Randomized numerical linear algebra (RandNLA) seeks to solve large-scale
linear algebra problems exploiting randomisation as a computational resource.
Central is the notion of a sketch. A sketch is a smaller or sparser
approximation of the original input instance, such as a matrix of data points,
that is used to compute quantities of interest. The review of Drineas and
Mahoney covers the main ideas and tools of RandNLA [DM18]. As a simple example
of the techniques used in RandNLA we present here the case of approximate
matrix multiplication via random sampling. In general, given two $n\times n$
square matrices $A$ and $B$, computing $AB$ takes $O(n^{3})$ time. Drineas,
Kannan, and Mahoney showed that, given an efficient sampling method for the
rows and columns of the matrices, it is possible to efficiently approximate
$AB$ in time $O(c\,n^{2})$, where $c$ is the number of rows and columns
randomly sampled from the matrices [DKM06]. The algorithm is structurally very
simple. Draw $c$ random samples of columns of $A$ and rows of $B$ according to
a probability distribution $p$. Group the samples into two, properly scaled,
smaller matrices $C$ and $R$. When the probability distribution $p$ is chosen
appropriately and the columns and rows accordingly re-scaled, it is possible
to show that that $\lVert AB-CR\rVert_{F}=O(\lVert A\rVert_{F}\lVert
B\rVert_{F}/\sqrt{c})$.

Similar ideas can be used to compute low-rank approximations of matrices that
can be then used in a wide range of applications [DKM06a]. In general, these
methods produce sketches that do not preserve given symmetries of the matrix,
a fundamental requirement for applications in quantum mechanics. A technique
where the symmetry of the sketched matrix is preserved is the Nyström method,
a RandNLA tool developed for the approximation of kernel matrices in
statistical learning theory. Roughly speaking, the method allows one to
construct a lower dimensional, symmetric, positive semidefinite approximation
of a given matrix given a sampling schemes for its columns. More specifically,
let $K\in\mathbb{R}^{n\times n}$ be a symmetric, rank $r$, PSD matrix,
$K_{:,j}$ the $j$-th column vector of $K$, and $K_{i,:}$ the $i$-th row vector
of $K$. The singular value decomposition of $K$ is $K=U\Sigma U^{\top}$, where
the columns of $U$ are orthogonal and
$\Sigma=\operatorname{diag}(\sigma_{1},\dots,\sigma_{r})$ is the matrix of the
singular values $\sigma_{1},\dots,\sigma_{r}$ of $K$. The Moore-Penrose
pseudoinverse of $K$ is
$K^{+}=\sum_{t=1}^{r}\sigma_{t}^{-1}U_{:,t}{U_{:,t}}^{\top}.$ The Nyström
method finds a low-rank approximation of $K$ that preserves the symmetry and
PSD property of the matrix. Let $C$ denote the $n\times l$ matrix formed by
(uniformly) sampling $l\ll n$ columns of $K$, $W$ denote the $l\times l$
matrix consisting of the intersection of these $l$ columns with the
corresponding $l$ rows of $K$, and $W_{k}$ denote the best rank-$k$
approximation of $W$

$W_{k}=\operatorname{argmin}_{V\in\mathbb{R}^{l\times
l},\operatorname{rank}(V)=k}\lVert V-W\rVert_{F}.$

The Nyström method generates a rank-$k$ approximation $\tilde{K}$ of $K$ for
$k<n$ defined by:

$\tilde{K}_{k}=CW_{k}^{+}C^{\top}\approx K$

The running time of the algorithm is $O(nkl)$ [KMT12].

The performance of the algorithm can be improved using non-uniform sampling
schemes. For a sampling scheme equivalent to the one we use in this paper we
have

######  Theorem 4 (Theorem $3$ in [DM05]).

Let $K$ be a $n\times n$ symmetric PSD matrix, let $k<l$ be a rank parameter,
and let $\tilde{K}$ be the approximation constructed using the Nyström method
by sampling $l$ columns with probabilities $p_{i}=|K_{:,i}|^{2}/\lVert
K\rVert_{F}^{2}$. Let $\epsilon>0$ and
$\eta=1+\sqrt{8\operatorname{log}(1/\delta)}$, if $l>64k\eta^{2}/\epsilon^{4}$
then, with probability at least $1-\delta$

$\left\lVert K-\tilde{K}_{k}\right\rVert_{F}\leq\left\lVert
K-K_{k}\right\rVert_{F}+\epsilon\sum_{i=1}^{n}K_{ii}^{2}.$

The Nyström method has proved to be a powerful tool in a range of applications
where the matrices are approximately low rank. The method in its present form
was developed by Williams and Seeger [WS01] as a sampling-based algorithm to
solve regression and classification problems involving Gaussian processes.
This problem requires the approximation of symmetric, positive semidefinite
matrices that can be well low-rank approximated [WS01, Wil+02]. The technique
proved to be closely related to a method for solving linear integral equations
developed by Nyström [Nys30] and hence the name Nyström method.

It is worth mentioning the Nyström extension, a further refinement of the
technique that has found numerous applications ranging from large-scale
machine learning problems, to applications in statistics and signal processing
[WS01, Wil+02, ZK10, TKR08, Fow+04, KMT12, BW07, BW08, KMT09, LKL10, MTJ11,
ZK10, ZTK08]. Typical extensions that substantially improve the performance,
e.g. lead to lower reconstruction error, introduce non-uniform importance
sampling distributions or random mixing of the input before sampling the
columns.

###  1.4 Organization

Section 2 introduces relevant notation and definitions. Section 3 discusses
the row-searchability condition and the algorithm to sample efficiently from
the rows of the Hamiltonian. Section 4 proves the theorem for the approximate
simulation of exponentials of PSD matrices. Section 5 proves the theorem for
the approximate simulation of exponentials of Hermitian matrices. Section 6
discusses applications to the simulation of the evolution of density matrices.

##  2 Preliminaries

We denote vectors with lower-case letters. For a vector $x\in\mathbb{C}^{n}$,
let $x_{i}$ denotes the $i$-th element of $x$. A vector is sparse if most of
its entries are $0$. For an integer $k$, let $[k]$ denotes the set
$\\{1,\dots,k\\}$.

For a matrix $A\in\mathbb{C}^{m\times n}$ let $a^{j}:=A_{:,j}$, $j\in[n]$
denote the $j$-th column vector of $A$, $a_{i}:=A_{i,:}$, $i\in[m]$ the $i$-th
row vector of $A$, and $a_{ij}=A(i,j)$ the $(i,j)$-th element. We denote by
$A_{i:j}$ the sub-matrix of $A$ that contains the rows from $i$ to $j$.

The supremum is denoted as $\mathrm{sup}$ and the infimum as $\mathrm{inf}$.
For a measure space $(X,\Sigma,\mu)$, and a measurable function $f$ an
essential upper bound of $f$ is defined as
$U_{f}^{\text{ess}}:=\\{l\in\mathbb{R}:\mu(f^{-1}(l,\infty)=0\\}$ if the
measurable set $f^{-1}(l,\infty)$ is a set of measure zero, i.e., if $f(x)\leq
l$ for almost all $x\in X$. Then the essential supremum is defined as
$\text{ess sup}f:=\mathrm{inf}\ U_{f}^{\text{ess}}$. We let the
$\mathrm{span}$ of a set $S=\\{v_{i}\\}_{1}^{k}\subseteq\mathbb{C}^{n}$ be
defined by
$\mathrm{span}\\{S\\}:=\left\\{x\in\mathbb{C}^{n}\,|\,\exists\,\\{\alpha_{i}\\}_{1}^{k}\subseteq\mathbb{C}\text{
with }x=\sum_{i=1}^{k}\alpha_{i}v_{i}\right\\}$. The set is linearly
independent if $\sum_{i}\alpha_{i}v_{i}=0$ if and only if $\alpha_{i}=0$ for
all $i$. The range of $A\in\mathbb{C}^{m\times n}$ is defined by
$\text{range}(A)=\\{y\in\mathbb{R}^{m}:y=Ax\text{ for some
}x\in\mathbb{C}^{n}\\}=\mathrm{span}(A^{1},\ldots,A^{n})$. Equivalently the
range of $A$ is the set of all linear combinations of the columns of $A$. The
nullspace $\mathrm{null}(A)$ (or kernel $\mathrm{ker}(A)$) is the set of
vectors such that $Av=0$. Given a set
$S=\\{v_{i}\\}_{1}^{k}\subseteq\mathbb{C}^{n}$. The null space of $A$ is
$\mathrm{null}(A)=\\{x\in\mathbb{R}^{b}:Ax=0\\}$.

The rank of a matrix $A\in\mathbb{C}^{m\times n}$, $\mathrm{rank}(A)$ is the
dimension of $\mathrm{range}(A)$ and is equal to the number of linearly
independent columns of $A$; Since this is equal to $\mathrm{rank}(A^{T})$ it
also equals the number of linearly independent rows of $A$, and satisfies
$\mathrm{rank}(A)\leq\min\\{m,n\\}$. The trace of a matrix is the sum of its
diagonal elements $\mathrm{Tr}\left(A\right)=\sum_{i}a_{ii}$. The support of a
vector $\mathrm{supp}(v)$ is the set of indices $i$ such that $v_{i}=0$ and we
call it sparsity of the vector. For a matrix we denote the sparsity as the
number of non zero entries, while row or column sparsity refers to the number
of non-zero entries per row or column. A symmetric matrix $A$ is positive
semidefinite (PSD) if all its eigenvalues are non-negative. For a PSD matrix
$A$ we write $A\succeq 0$. Similarly $A\succeq B$ is the partial ordering
which is equivalent to $A-B\succeq 0$.

We use the following standard norms. The Frobenius norm $\lVert
A\rVert_{F}=\sqrt{\sum_{i=1}^{m}\sum_{j=1}^{n}A_{ij}^{2}}$, and the spectral
norm $\lVert A\rVert=\sup_{x\in\mathbb{C}^{n},\;x\neq 0}\frac{|Ax|}{|x|}$.
Note that that $\lVert
A\rVert_{F}^{2}=\mathrm{Tr}\left(A^{T}A\right)=\mathrm{Tr}\left(AA^{T}\right)$.
Both norms are submultiplicative and unitarily invariant and they are related
to each other as $\lVert A\rVert\leq\lVert A\rVert_{F}\leq\sqrt{n}\lVert
A\rVert$.

The singular value decomposition of $A$ is $A=U\Sigma V^{*}$ where $U,V$ are
unitary matrices and $U^{*}$ defines the complex conjugate transpose, also
called Hermitian conjugate, of $U$. We denote the pseudo-inverse of a matrix
$A$ with singular value decomposition $U\Sigma V^{*}$ as
$A^{+}:=V\Sigma^{+}U^{*}$.

##  3 From row-searchability to efficient row-sampling

All our algorithms require row-searchable Hamiltonians. In this section we
define the row-searchability condition and describe an efficient algorithm to
sample from the rows of row-searchable Hamiltonians.

Let $n\in\mathbb{N}$. We first introduce a binary tree of subsets spanning
$\\{0,1\\}^{n}$. In the following, with abuse of notation, we identify binary
tuples with the associated binary number. Let $L$ be a binary string with
$|L|\leq n$, where $|L|$ denotes the length of the string. We denote with
$S(L)$ the set

$\displaystyle
S(L)=\\{L\\}\times\\{0,1\\}^{n-|L|}=\\{(L_{1},\dots,L_{|L|},v_{1},\dots,v_{n-|L|})\leavevmode\nobreak\
|\leavevmode\nobreak\ v_{1},\dots,v_{n-|L|}\in\\{0,1\\}\\}.$ (1)

We are now ready to state the row-searchability property for a matrix $H$.

######  Definition 3 (Row-searchability).

Let $H$ be a Hermitian matrix of dimension $2^{n}$, for $n\in\mathbb{N}$. $H$
is row-searchable if, for any binary string $L$ with $|L|\leq n$, it is
possible to compute the following quantity in $O(\poly(n))$

$w(S(L))=\sum_{i\in S(L)}h(i,H_{:,i}),$ (2)

where $h$ is the function computing the weight associated to the $i$-th column
$H_{:,i}$. For positive semidefinite $H$ we use $h(i,H_{:,i})=H_{i,i}$. i.e.
the diagonal element $i$ while for general Hermitian $H$ we use
$h(i,H_{:,i})=\|H_{:,i}\|^{2}$.

Row-searchability intuitively works as follows. If we are given a binary tree,
where the leaves contain the individual probabilities and the parents at each
level contain the marginal over their children nodes, then we can, for a
randomly sampled number in $[0,1]$ traverse this tree in $\log(N)$ time to
find the leave node that is sampled, i.e. the indices of the column of $H$.
More specifically, row-searchability requires the evaluation of $w(S(L)))$ as
defined in Eq. (2) which computes marginals of the diagonal of $H$, where the
co-elements, i.e. the elements where we are not summing over, are defined by
the tuple $L$. Hence, for empty $L$, $w(S(L))=\mathrm{Tr}\left(H\right)$.

Algorithm 1 MATLAB code for the sampling algorithm

Input: wS(L) corresponds to the function $w(S(L))$ defined in Eq. 2.

Output: $L$ is the sampled row index

    L = [];
    q = rand()*wS(L);

    for i=1:n
      if q >= wS([L 0])
          L = [L 0];
      else
          L = [L 1];
          q = q - wS([L 0]);
      end
    end

Note that the function $h$ considered is related to leverage score sampling,
which is a common approach in randomized algorithm for linear algebra [Mah+11,
Woo+14].

In Alg. 1 we provide an algorithm, that, given a row-searchable $H$, is able
to sample an index with probability $p(j)=h(j,H_{:,j})/w(\\{0,1\\}^{n})$. Let
$q$ be a random number uniformly sampled in $[0,T]$, where
$T=w(\\{0,1\\}^{n})$ is the sum of the weights associated to all the rows. The
algorithm uses logarithmic search, starting with $L$ empty and adding
iteratively $1$ or $0$, to find the index $L$ such that
$w(\\{0,\dots,L-1\\})\leq q\leq w(\\{0,\dots,L\\})$. The total time required
to compute one index, is $O(nQ(n))$ where $Q(n)$ is the maximum time required
to compute a $w(S(L))$ for $L\in\\{0,1\\}^{n}$. Note that if $w(S(L))$ can be
computed efficiently for any $L\in\\{0,1\\}^{n}$ then $Q(n)$ is polynomial and
the cost of the sampling procedure will be polynomial.

######  Remark 1 (Row-searchability more general than sparsity).

Note that if $H$ has a polynomial number of non-zero elements, then $w(S(L))$
can be always computed in polynomial time. Indeed given $L$, we go through the
list of elements describing $H$ and select only the ones whose row-index
starts with $L$ and then compute $w(S(L))$, both step requiring polynomial
time. However $w(S(L))$ can be computed efficiently even for Hamiltonians that
are not polynomially sparse. For example, take the diagonal Hamiltonian
defined by $H_{ii}=1/i$ for $i\in[2^{n}]$. This $H$ is not polynomially sparse
and in particular it has an exponential number of non-zero elements, but still
$w(S(L))$ can be computed in polynomial time, here in particular in $O(1)$.

##  4 Algorithm for PSD row-searchable Hermitian matrices

Given a $2^{n}\times 2^{n}$ matrix $H\succeq 0$, our goal is to produce an
approximation of the state

$\psi(t)=\exp{(iHt)}\psi.$ (3)

In particular, we will provide an expression of the form
$\widehat{\exp}(i\widehat{H}t)\psi$, where $\widehat{\exp}$ and $\widehat{H}$
are a suitable approximation respectively of the exponential function and $H$.
We give here an algorithm for $H\succeq 0$ which we then generalize in the
following section to arbitrary Hermitian $H$, if the row-searchability
condition, i.e. condition 3, is fulfilled.

Let $h$ be the diagonal of the positive semidefinite $H$ and let
$t_{1},\dots,t_{M}$, with $M\in\mathbb{N}$ be indices independently sampled
with repetition from $\\{1,\dots,2^{n}\\}$, with probabilities

$p(q)=h_{q}/\sum_{i}h_{i},$ (4)

e.g. via Alg. 1. Then, define the matrix $B\in\mathbb{C}^{M\times M}$ such
that $B_{i,j}=H_{t_{i},t_{j}}$, for $1\leq i,j\leq M$. Finally, denote by
$A\in\mathbb{C}^{2^{n}\times M}$ the matrix $A_{i,j}=H_{i,t_{j}}$ for $1\leq
i\leq 2^{n}$ and $1\leq j\leq M$. The approximated matrix is defined as
$\widehat{H}=AB^{+}A^{*}$, where $(\cdot)^{+}$ is the pseudoinverse. Let us
also define a function $g(x)=(e^{itx}-1)/x$. Then we have $e^{itx}=1+g(x)x$.
Note that $g$ is an analytic function, in particular,

$g(x)=\sum_{k\geq 1}\frac{(it)^{k}}{k!}x^{k-1}.$

Then

$e^{i\widehat{H}t}=I+g(\widehat{H})\widehat{H}=I+g(AB^{+}A^{*})AB^{+}A^{*}=I+Ag(B^{+}A^{*}A)B^{+}A^{*},$
(5)

where the last step is due to the fact that, given an analytic function
$q(x)=\sum_{k\geq 0}\alpha_{k}x^{k}$, we have

$\displaystyle q(AB^{+}A^{*})AB^{+}A^{*}$ $\displaystyle=\sum_{k\geq
1}\alpha_{k}(AB^{+}A^{*})^{k}AB^{+}A^{*}$ $\displaystyle=A\sum_{k\geq
1}\alpha_{k}(B^{+}A^{*}A)^{k}B^{+}A^{*}=Aq(B^{+}A^{*}A)B^{+}A^{*}.$

By writing $D=B^{+}A^{*}A$, the algorithm is now

$\widehat{\psi}_{M}(t)=\psi+Ag(D)B^{+}A^{*}\psi.$

Now we approximate $g$ with $g_{K}(x)$, which limits the series defining $g$
to the first $K$ terms, for $K\in\mathbb{N}$. Moreover, by rewriting
$g_{K}(D)B^{+}(A^{*}\psi)$ in an iterative fashion, we have

$b_{j}=\frac{(it)^{K-j}}{(K-j)!}v+Db_{j-1},\quad v=B^{+}(A^{*}\psi),$

with $b_{0}=\frac{(it)^{K}}{K!}v$ and so $b_{K-1}=g(D)B^{+}A^{*}\psi$.
Finally, the new approximate state is given by

$\widehat{\psi}_{K,M}(t)=\psi+Ab_{K-1}.$ (6)

Algorithm 2 MATLAB code for approximating Hamiltonian dynamics when $H$ is PSD

Input: $M$, $T=t_{1},\dots,t_{M}$ list of indices computed via Alg. 1. The
function compute_H_subMatrix, given two lists of indices, computes the
associated submatrix of $H$. compute_psi_subVector, given a list of indices,
computes the associated subvector of $\psi$.

Output: vector $b$, s.t. $b=e^{iHt}\psi$.

    B = compute_H_subMatrix(T, T);

    D = zeros(M,M);
    v = zeros(M,1);
    for i=1:(2^n/M)
        E = compute_H_subMatrix((i-1)*M+1:M*i, T);
        D = D + E’*E;
        v = v + E’*compute_Psi_subVector((i-1)*M+1:M*i);
    end
    u
    D = B\D;

    b = zeros(M,1);
    for j=1:K
        b = (1i*t)^(K-j)/factorial(K-j) * v + D*b;
    end

A MATLAB implementation of this procedure is presented in Alg. 2.
Let the row sparsity $s$ of $H$ be of order $\poly(n)$. Then the total cost of
applying this operator is given by $O(M^{2}\poly(n)+KM^{2}+M^{3}))$ in time,
where the terms $M^{3}$ and $M^{2}\poly(n)$ are resulting from the calculation
of $D$. To compute the total cost in space, note that we do not have to save
$H$ or $A$ in memory, but only $B,D$ and the vectors $v,b_{j}$, for a total
cost of $O(M^{2})$. Indeed $D$ can be computed in the following way. Assuming,
without loss of generality, to have $2^{n}/M\in\mathbb{N}$, then

$D=B^{-1}\sum_{i=1}^{2^{n}/M}A_{M(i-1)+1:Mi}^{*}A_{M(i-1)+1:Mi},$

where $A_{a:b}$ is the submatrix of $A$ containing the rows from $a$ to $b$. A
similar reasoning holds for the computation of the vector $v$. In this
computation we have assumed that the sample probabilities are give to us and
that we can efficiently sample from the matrix $H$ according to these
probabilities. In order to make our algorithm practical we hence need to give
an algorithm for performing the sampling. This properties of the SPDS case
algorithm are summarized in the following theorem.

######  Theorem 5 (Algorithm for simulating PSD row-searchable Hermitian
matrices).

Let $\epsilon,\delta\in(0,1]$, let $K,M\in\mathbb{N}$ and $t>0$. Let $H$ be
positive semidefinite, where $K$ is the number of terms in the truncated
series expansions of $g(\hat{H})$ and $M$ the number of samples we take for
the approximation. Let $\psi(t)$ be the true evolution (Eq. 3) and let
$\widehat{\psi}_{K,M}(t)$ be the output of our Alg. 2 (Eq. 6). When

$K\geq e\,t\lVert H\rVert+\log\frac{2}{\epsilon},\qquad\quad
M\geq\max\left(405\mathrm{Tr}\left(H\right),\leavevmode\nobreak\
\frac{72\mathrm{Tr}\left(H\right)t}{\epsilon}\log\frac{36\mathrm{Tr}\left(H\right)t}{\epsilon\delta}\right),$
(7)

then the following holds with probability $1-\delta$,

$\lVert\psi(t)-\widehat{\psi}_{K,M}(t)\rVert\leq\epsilon.$

Note that with the result above, we have that $\widehat{\psi}_{K,M}(t)$ in Eq.
(6) (Alg. 2) approximates $\psi(t)$, with error at most $\epsilon$ and with
probability at least $1-\delta$, requiring a computational cost that is
$O(\frac{st^{2}\mathrm{Tr}\left(H\right)^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$
in time and
$O(\frac{t^{2}\mathrm{Tr}\left(H\right)^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$
in memory.

In the following we now prove the first main result of this work. To prove
Theorem 5 we need the following lemmas. Lemma 1 performs a basic decomposition
of the error in terms of the distance between $H$ and the approximation
$\widehat{H}$ as well as in terms of the approximation $g_{K}$ with respect to
$g$. Lemma 2 provides an analytic bound on the distance between $H$ and
$\widehat{H}$, expressed in terms of the expectation of eigenvalues or related
matrices which are then concentrated in Lemma 3.

######  Lemma 1.

Let $K,M\in\mathbb{N}$ and $t>0$, then

$\lVert\psi(t)-\widehat{\psi}_{K,M}(t)\rVert\leq t\lVert
H-\widehat{H}\rVert+\frac{(t\lVert\widehat{H}\rVert)^{K+1}}{(K+1)!}.$

###### Proof.

By definition we have that $e^{ixt}=1+g(x)x$ with $g(x)=\sum_{k\geq
1}x^{k-1}(it)^{k}/k!$ and $g_{K}$ is the truncated version of $g$. By adding
and subtracting $e^{i\widehat{H}t}$, we have

$\lVert
e^{iHt}\psi-(I+g_{K}(\widehat{H})\widehat{H})\psi\rVert\leq\lVert\psi\rVert\
(\lVert e^{iHt}-e^{i\widehat{H}t}\rVert+\lVert
e^{i\widehat{H}t}-(I+g_{K}(\widehat{H})\widehat{H})\rVert).$ (8)

By [Nak03],

$\lVert e^{iHt}-e^{i\widehat{H}t}\rVert\leq t\lVert H-\widehat{H}\rVert,$ (9)

moreover, by [Mat93], and since $\widehat{H}$ is Hermitian and hence all the
eigenvalues are real, we have

$\lVert
e^{i\widehat{H}t}-(I+g_{K}(\widehat{H})\widehat{H})\rVert\leq\frac{(t\lVert\widehat{H}\rVert)^{K+1}}{(K+1)!}\sup_{l\in[0,1]}\lVert
i^{K+1}e^{il\widehat{H}t}\rVert\leq\frac{(t\lVert\widehat{H}\rVert)^{K+1}}{(K+1)!}.$
(10)

Finally note that $\lVert\psi\rVert=1$. ∎

To study the norm $\lVert H-\widehat{H}\rVert$ note that, since $H$ is
positive semidefinite, there exists an operator $S$ such that $H=SS^{*}$, so
$H_{i,j}=s_{i}^{*}s_{j}$ with $s_{i},s_{j}$ the $i$-th and $j$-th row of $S$.
Denote with $C$ and $\widetilde{C}$ the operators

$C=S^{*}S,\quad\widetilde{C}=\frac{1}{M}\sum_{j=1}^{M}\frac{\mathrm{Tr}\left(H\right)}{h_{t_{j}}}s_{t_{j}}s_{t_{j}}^{*}.$

We then obtain the following result.

######  Lemma 2.

The following holds with probability $1$. For any $\tau>0$,

$\lVert
H-\widehat{H}\rVert\leq\frac{\tau}{1-\beta(\tau)},\quad\quad\beta(\tau)=\lambda_{\max}((C+\tau
I)^{-1/2}(C-\widetilde{C})(C+\tau I)^{-1/2}),$ (11)

moreover $\lVert\widehat{H}\rVert\leq\lVert H\rVert$.

###### Proof.

Define the selection matrix $V\in\mathbb{C}^{M\times 2^{n}}$, that is always
zero except for one element in each row which is $V_{j,t_{j}}=1$ for $1\leq
j\leq M$. Then we have that

$A=HV^{*},\quad B=VHV^{*},$

i.e., $A$ is again given by the rows according to the sampled indiced
$t_{1},\ldots,t_{M}$ and $B$ is the submatrix obtained from taking the rows
and columns according to the same indices. In particular by denoting with
$\widehat{P}$ the operator $\widehat{P}=S^{*}V^{*}(VSS^{*}V^{*})^{+}VS$, and
recalling that $H=SS^{*}$ and $C=S^{*}S$, we have

$\widehat{H}=AB^{+}A^{*}=SS^{*}V^{*}(VSS^{*}V^{*})^{+}VSS^{*}=S\widehat{P}S^{*}.$

By definition $\widehat{P}$ is an orthogonal projection operator, indeed it is
symmetric and, by definition $Q^{+}QQ^{+}=Q^{+}$, for any matrix $Q$, then

$\widehat{P}^{2}=S^{*}V^{*}[(VSS^{*}V^{*})^{+}(VSS^{*}V^{*})(VSS^{*}V^{*})^{+}]VS=S^{*}V^{*}(VSS^{*}V^{*})^{+}VS=\widehat{P}.$
(12)

Indeed this is a projection in the row space of the matrix $R:=S^{*}V^{*}$,
since with the singular value decomposition $R:=U_{R}\Sigma_{R}V_{R}^{*}$ we
have

$\hat{P}=R^{*}(R^{*}R)^{+}R=V_{R}\Sigma_{R}U_{R}^{*}U_{R}\Sigma^{-2}_{R}U_{R}^{*}U_{R}\Sigma_{R}V_{R}^{*}=V_{R}V_{R}^{*},$
(13)

which spans the same space as $R$. Finally, since
$(I-\widehat{P})=(I-\widehat{P})^{2}$, and $\lVert Z^{*}Z\rVert=\lVert
Z\rVert^{2}$, we have

$\lVert H-\widehat{H}\rVert=\lVert S(I-\widehat{P})S^{*}\rVert=\lVert
S(I-\widehat{P})^{2}S^{*}\rVert=\lVert(I-\widehat{P})S^{*}\rVert^{2}.$ (14)

Note that $\widetilde{C}$ can be rewritten as $\widetilde{C}=S^{*}V^{*}LVS$,
with $L$ a diagonal matrix, with
$L_{jj}=\frac{\mathrm{Tr}\left(H\right)}{Mh_{t_{j}}}$. Moreover $t_{j}$ is
sampled from the probability $p(q)=h_{q}/\mathrm{Tr}\left(H\right)$, so
$h_{t_{j}}>0$ with probability 1, then $L$ has a finite and strictly positive
diagonal, so $\widetilde{C}$ has the same range of $\widehat{P}$. Now, with
$C=S^{*}S$, we are able to apply Proposition 3 and Proposition 7 of [RCR15a],
and obtain

$\lVert(I-\widehat{P})S^{*}\rVert^{2}\leq\frac{\tau}{1-\beta(\tau)},\quad\quad\beta(\tau)=\lambda_{\max}((C+\tau
I)^{-1/2}(C-\widetilde{C})(C+\tau I)^{-1/2}).$ (15)

Finally, note that, since $\widehat{P}$ is a projection operator we have that
$\lVert\widehat{P}\rVert=1$, so

$\lVert\widehat{H}\rVert=\lVert
S\widehat{P}S^{*}\rVert\leq\lVert\widehat{P}\rVert\lVert S\rVert^{2}\leq\lVert
S\rVert^{2}=\lVert H\rVert,$

where the last step is due to the fact that $H=SS^{*}$. ∎

######  Lemma 3.

Let $\delta\in(0,1]$ and $\tau>0$. When

$M\geq\max\left(405\mathrm{Tr}\left(H\right),\leavevmode\nobreak\
67\mathrm{Tr}\left(H\right)\log\frac{\mathrm{Tr}\left(H\right)}{2\delta}\right),\quad\quad\tau=\frac{9\mathrm{Tr}\left(H\right)}{M}\log\frac{M}{2\delta},$
(16)

then with probability $1-\delta$ it holds that

$\lambda_{\max}((C+\tau I)^{-1/2}(C-\widetilde{C})(C+\tau
I)^{-1/2})\leq\frac{1}{2}.$

###### Proof.

Define the random variable
$\zeta_{j}=\sqrt{\frac{\mathrm{Tr}\left(H\right)}{h_{t_{j}}}}s_{t_{j}},$ for
$1\leq j\leq M$. Note that

$\lVert\zeta_{j}\rVert\leq\sqrt{\frac{\mathrm{Tr}\left(H\right)}{h_{t_{j}}}}\lVert
s_{t_{j}}\rVert\leq\sqrt{\mathrm{Tr}\left(H\right)},$

almost surely. Moreover,

$\mathbb{E}\zeta_{j}\zeta_{j}^{*}=\sum_{q=1}^{2^{n}}p(q)\frac{\mathrm{Tr}\left(H\right)}{h_{q}}s_{q}s_{q}^{*}=\sum_{q=1}^{2^{n}}s_{q}s_{q}^{*}=S^{*}S=C.$
(17)

By definition of $\zeta_{j}$, we have

$\widetilde{C}=\frac{1}{M}\sum_{j=1}^{M}\zeta_{j}\zeta_{j}^{*}.$

Since $\zeta_{j}$ are independent for $1\leq j\leq M$, uniformly bounded, with
expectation equal to $C$, and with $\zeta_{j}^{*}(C+\tau
I)^{-1}\zeta_{j}\leq\lVert\zeta_{j}\rVert^{2}\tau^{-1}\leq\mathrm{Tr}\left(H\right)\tau^{-1}$,
we can apply Proposition 8 of [RCR15a], that uses non-commutative Bernstein
inequality for linear operators [Tro12], and obtain

$\lambda_{\max}((C+\tau I)^{-1/2}(C-\widetilde{C})(C+\tau
I)^{-1/2})\leq\frac{2\alpha}{3M}+\sqrt{\frac{2\alpha}{Mt}},$ (18)

with probability at least $1-\delta$, with
$\alpha=\log\frac{4\mathrm{Tr}\left(C\right)}{\tau\delta}$. Since

$\mathrm{Tr}\left(C\right)=\mathrm{Tr}\left(S^{*}S\right)=\mathrm{Tr}\left(SS^{*}\right)=\mathrm{Tr}\left(H\right),$

by Remark 1 of [RCR15a], we have that

$\lambda_{\max}((C+\tau I)^{-1/2}(C-\widetilde{C})(C+\tau
I)^{-1/2})\leq\frac{1}{2},$

with probability $1-\delta$, when
$M\geq\max(405\kappa^{2},67\kappa^{2}\log\frac{\kappa^{2}}{2\delta})$ and
$\tau$ satisfies $\frac{9\kappa^{2}}{M}\log\frac{M}{2\delta}\leq\tau\leq\lVert
C\rVert$ (note that $\lVert C\rVert=\lVert H\rVert$), where $\kappa^{2}$ is a
bound for the following quantity

$\begin{split}\inf_{\tau>0}[(\lVert
C\rVert+\tau)(\mathrm{ess}\sup\zeta_{j}^{*}(C+\tau I)^{-1}\zeta_{j})]\\\
\leq\mathrm{Tr}\left(H\right)\inf_{\tau>0}\frac{\lVert
H\rVert+\tau}{\tau}\leq\mathrm{Tr}\left(H\right):=\kappa^{2},\end{split}$ (19)

where $\mathrm{ess}\sup$ here denotes the essential supremum. ∎

Now we are ready to prove the Theorem for PSD matrices.

###### Proof of Theorem 5.

By Lemma 1, we have

$\lVert e^{iHt}\psi-(I+g_{K}(\widehat{H})\widehat{H})\psi\rVert\leq t\lVert
H-\widehat{H}\rVert+\frac{(t\lVert\widehat{H}\rVert)^{K+1}}{(K+1)!}.$

Let $\tau>0$. By Lemma 2, we know that $\lVert\widehat{H}\rVert\leq\lVert
H\rVert$ and that

$\lVert H-\widehat{H}\rVert\leq\frac{\tau}{1-\beta(\tau)},$
$\quad\beta(\tau)=\lambda_{\max}((C+\tau I)^{-1/2}(C-\widetilde{C})(C+\tau
I)^{-1/2}),$

with probability $1$. Finally by Lemma 3, we have that the following holds
with probability $1-\delta$,

$\lambda_{\max}((C+\tau I)^{-1/2}(C-\widetilde{C})(C+\tau
I)^{-1/2})\leq\frac{1}{2},$

when

$M\geq\max\left(405\mathrm{Tr}\left(H\right),\leavevmode\nobreak\
67\mathrm{Tr}\left(H\right)\log\frac{\mathrm{Tr}\left(H\right)}{2\delta}\right)$

and

$\quad\tau=\frac{9\mathrm{Tr}\left(H\right)}{M}\log\frac{\mathrm{Tr}\left(H\right)}{2\delta}.$

So we have

$\lVert
e^{iHt}\psi-(I+g_{K}(\widehat{H})\widehat{H})\psi\rVert\leq\frac{18\mathrm{Tr}\left(H\right)t}{M}\log\frac{M}{2\delta}+\frac{(t\lVert
H\rVert)^{K+1}}{(K+1)!},$

with probability $1-\delta$.
Now we select $K$ such that $\frac{(t\lVert
H\rVert)^{K+1}}{(K+1)!}\leq\frac{\epsilon}{2}$. Since, by the Stirling
approximation, we have

$(K+1)!\geq\sqrt{2\pi}(K+1)^{K+3/2}e^{-K-1}\geq(K+1)^{K+1}e^{-K-1}.$

Since

$(1+x)\log(1/(1+x))\leq-x,$

for $x>0$ we can select $K=et\lVert H\rVert+\log\frac{2}{\epsilon}-1$, such
that we have

$\displaystyle\log\left(\frac{(t\lVert H\rVert)^{K+1}}{(K+1)!}\right)$
$\displaystyle\leq(K+1)\log\frac{et\lVert H\rVert}{K+1}$ $\displaystyle\leq
et\lVert H\rVert\left(1+\frac{\log\frac{2}{\epsilon}}{et\lVert
H\rVert}\right)\log\frac{1}{1+\frac{\log\frac{2}{\epsilon}}{et\lVert
H\rVert}}$ $\displaystyle\leq\log\frac{\epsilon}{2}.$

Finally we require $M$, such that

$\frac{18\mathrm{Tr}\left(H\right)t}{M}\log\frac{M}{2\delta}\leq\frac{\epsilon}{2},$

and select

$M=\frac{72\mathrm{Tr}\left(H\right)t}{\epsilon}\log\frac{36\mathrm{Tr}\left(H\right)t}{\epsilon\delta}.$

Then we have that

$\frac{18\mathrm{Tr}\left(H\right)t}{M}\log\frac{M}{2\delta}\leq\frac{\epsilon}{2}\frac{\log\frac{36\mathrm{Tr}\left(H\right)t}{\epsilon\delta}+\log\log\frac{36\mathrm{Tr}\left(H\right)t}{\epsilon\delta}}{2\log\frac{36\mathrm{Tr}\left(H\right)t}{\epsilon\delta}}\leq\frac{\epsilon}{2}.$

∎

We now extend this result to the more general case of arbitrary Hermitian
matrices that fulfill the row-searchability condition, i.e. the ability to
sample according to some leverage of the rows. This will lead to the second
main result of this work.

##  5 Algorithm for row-searchable Hermitian matrices

In this section we provide the algorithm for simulating Hermitian (possibly
non-psd) matrices and we provide guarantees on the efficiency when $H$ is row-
searchable. Let in the following $s$ be the maximum number of non-zero
elements of the rows of $H$, $\epsilon$ the error in the approximation of the
output states of the algorithm w.r.t. the ideal $\psi(t)$, and $t$ the
evolution time of the simulation. Let further $K$ be the order of the
truncated series expansions and $M$ the number of samples we take for the
approximation. As before we first outline the algorithm and then prove its
properties.

For arbitrary matrices $H$ we will use the following algorithm. Sample
$M\in\mathbb{N}$ independent indices $t_{1},\dots t_{M}$, with probability
$p(i)=\frac{\|h_{i}\|^{2}}{\|H\|_{F}^{2}}$, $1\leq i\leq 2^{n}$, where $h_{i}$
is the $i$-th row of $H$ (sample via Alg. 1). Then denote with $A$, the matrix
$2^{n}\times M$ defined as

$A=\left[\frac{1}{\sqrt{Mp(t_{1})}}h_{t_{1}},\dots,\frac{1}{\sqrt{Mp(t_{1})}}h_{t_{M}}\right].$

Then we will use $H^{2}=AA^{*}$ as the approximation for the Hamiltonian.
Define two functions that we will use to approximate $e^{ix}$,

$f(x)=\frac{\cos(\sqrt{x})-1}{x},\quad
g(x)=\frac{\sin(\sqrt{x})-\sqrt{x}}{x\sqrt{x}},$

moreover denote with $f_{K}$ and $g_{K}$ the $K$-truncated Taylor expansions
of $f$ and $g$, for $K\in\mathbb{N}$

$f_{K}(x)=\sum_{j=0}^{K}\frac{(-1)^{j+1}x^{j}}{(2j+2)!},\quad
g_{K}(x)=\sum_{j=0}^{K}\frac{(-1)^{j+1}x^{j}}{(2j+3)!}.$

In particular note that

$e^{ix}=1+ix+f(x^{2})x^{2}+ig(x^{2})x^{3}.$

Analogously to the previous algorithm, we hence estimate $e^{ix}$ via $f_{K}$
and $g_{K}$. The final algorithm will be

$\widehat{\psi}_{K,M}(t)=\psi+itu+t^{2}Af_{K}(t^{2}A^{*}A)v+it^{3}Ag_{K}(t^{2}A^{*}A)z,$
(20)

with $u=\hat{H}\psi$, $v=A^{*}\psi$, $z=A^{*}u$. Note that the product
$f_{k}(A^{*}A)v$ and $Ag_{K}(A^{*}A)z$ are done by exploiting the Taylor
series form of the two functions and performing only matrix vector products in
the same way as in Alg. 2. Denote with $s$ the maximum number of non-zero
elements in the rows of $H$, with $q$ the number of non-zero elements in
$\psi$. The final algorithm requires $O(sq)$ in space and time to compute $u$,
$O(M\min(s,q))$ in time and $O(M)$ in space to compute $v$ and $O(Ms)$ in time
and space to compute $z$. We therefore obtain a total computational complexity
of

$\displaystyle{\rm time:}\leavevmode\nobreak\ \leavevmode\nobreak\
O\left(sq+M\min(s,q)+sMK\right),$ (21) $\displaystyle\quad{\rm
space:}\leavevmode\nobreak\ \leavevmode\nobreak\ O\left(s(q+M)\right).$ (22)

Note that if $s>M$ is it possible to further reduce the memory requirements at
the cost of more computational time, by computing $B=A^{*}A$ that can be done
in blocks and require $O(sM^{2})$ in time and $O(M^{2})$ in memory, and then
compute

$\widehat{\psi}_{K,M}(t)=\psi+itu+t^{2}Af_{K}(t^{2}B)v+it^{3}Ag_{K}(t^{2}B)z.$

In that case the computational cost would be

$\displaystyle{\rm time:}\leavevmode\nobreak\ \leavevmode\nobreak\
O\left(sq+M\min(s,q)+M^{2}(s+K)\right),$ (23) $\displaystyle\quad{\rm
space:}\leavevmode\nobreak\ \leavevmode\nobreak\ O\left(sq+M^{2}\right).$ (24)

The properties of the this algorithm are summarized in the following theorem
(this is a formal statement of Theorem 1):

######  Theorem 6 (Algorithm for simulating row-samplable Hermitian matrices).

Let $\delta,\epsilon\in(0,1]$. Let $t>0$ and $K,M\in\mathbb{N}$, where $K$ is
the number of terms in the truncated series expansions of $g(\widehat{H})$ and
$M$ the number of samples we take for the approximation, and let $t>0$. Let
$\psi(t)$ be the true evolution (Eq. 3) and let $\widehat{\psi}_{K,M}(t)$ be
computed as in Eq. 20. When

$\displaystyle M\geq\frac{256t^{4}(1+t^{2}\lVert H\rVert^{2})\lVert
H\rVert^{2}_{F}\lVert H\rVert^{2}}{\epsilon^{2}}\log\frac{4\lVert
H\rVert^{2}_{F}}{\delta\lVert H\rVert^{2}},$ (25) $\displaystyle\quad K\geq
4t\sqrt{\lVert H\rVert^{2}+\epsilon}+\log\frac{4(1+t\lVert
H\rVert)}{\epsilon},$ (26)

then

$\lVert\widehat{\psi}_{K,M}(t)-\psi(t)\rVert\leq\epsilon,$

with probability at least $1-\delta$.

Note that with the result above, we have that $\widehat{\psi}_{K,M}(t)$ in Eq.
(20) approximates $\psi(t)$, with error at most $\epsilon$ and with
probability at least $1-\delta$, requiring a computational cost that is
$O\left(sq+M\min(s,q)+M^{2}(s+K)\right)$ in time and $O\left(sq+M^{2}\right)$
is memory.

Combining Eq. 23 and 24 with Eq. 25 and 26, the whole computational complexity
of the algorithm described in this section, is

$\displaystyle{\rm time:}\leavevmode\nobreak\ \leavevmode\nobreak\
O\left(sq+\frac{t^{9}\lVert H\rVert^{4}_{F}\lVert
H\rVert^{7}}{\epsilon^{4}}\left(n+\log\frac{1}{\delta}\right)^{2}\right),$
(27) $\displaystyle{\rm space:}\leavevmode\nobreak\ \leavevmode\nobreak\
O\left(sq+\frac{t^{8}\lVert H\rVert^{4}_{F}\lVert
H\rVert^{6}}{\epsilon^{4}}\left(n+\log\frac{1}{\delta}\right)^{2}\right),$
(28)

where the quantity $\log\frac{4\lVert H\rVert^{2}_{F}}{\delta\lVert
H\rVert^{2}}$ in Eq. 25 was bounded using the following inequality

$\log\frac{\lVert H\rVert^{2}_{F}}{\lVert
H\rVert^{2}}\leq\log\frac{2^{n}\lambda_{MAX}^{2}}{\lambda_{MAX}^{2}}=n,$

where $\lambda_{MAX}$ is the biggest eigenvalue of $H$.

Observe now that simulation of the time evolution of $\alpha I$ does only
change the phase of the time evolution, where $I\in\mathbb{C}^{N\times N}$ is
the identity matrix and $\alpha$ some real parameter. We can hence perform the
time evolution of $\tilde{H}:=H-\alpha I$, since for any efficient classical
description of the input state we can apply the time evolution of the diagonal
matrix $e^{-i\alpha It}$. We can then optimize the parameter $\alpha$ such
that the Frobenius norm of the operator $\tilde{H}$ is minimized, i.e.

$\displaystyle\alpha=\underset{\alpha}{\text{argmin}}\lVert\tilde{H}\rVert^{2}_{F}=\underset{\alpha}{\text{argmin}}\lVert
H-\alpha I\rVert^{2}_{F},$ (29)

from which we obtain the condition
$\alpha=\frac{\mathrm{Tr}\left(H\right)}{2^{n}}$. Since our algorithm requires
that $\lVert\tilde{H}\rVert_{F}$ is bounded by $\polylog N$. Using the
spectral theorem, and the fact that the Frobenius norm is unitarily invariant,
this in turn gives us after a bit of algebra the condition

$\displaystyle\lVert
H\rVert_{F}^{2}-\frac{1}{N}\mathrm{Tr}\left(H\right)^{2}\leq O(\polylog(N)),$
(30)

for which we can simulate the Hamiltonian $H$ efficiently.
We now prove the second main result of this work and establish the correctness
of the above results.

###### Proof of Theorem 6.

Denote with

$\displaystyle\widehat{Z}_{K}(Ht,At)$
$\displaystyle=I+itH+t^{2}Af_{K}(t^{2}A^{*}A)A^{*}+it^{3}Ag_{K}(t^{2}A^{*}A)A^{*}H,$
$\displaystyle\widehat{Z}(Ht,At)$
$\displaystyle=I+itH+t^{2}Af(t^{2}A^{*}A)A^{*}+it^{3}Ag(t^{2}A^{*}A)A^{*}H.$

By definition of $\widehat{\psi}_{K,M}(t)$ and the fact that
$\lVert\psi\rVert=1$, we have

$\lVert\widehat{\psi}_{K,M}(t)-\psi(t)\rVert\leq\lVert\widehat{Z}_{K}(At,Ht)-e^{iHt}\rVert\lVert\psi\rVert\\\
\leq\lVert\widehat{Z}_{K}(At,Ht)-\widehat{Z}(Ht,At)\rVert+\lVert\widehat{Z}(Ht,At)-e^{iHt}\rVert.$

We first study $\lVert\widehat{Z}(Ht,At)-e^{iHt}\rVert$. Define $l(x)=f(x)x$
and $m(x)=g(x)x$. Note that, by the spectral theorem, we have

$\widehat{Z}(Ht,At)=I+itH+t^{2}Af(t^{2}A^{*}A)A^{*}+it^{3}Ag(t^{2}A^{*}A)A^{*}H=I+itH+t^{2}f(t^{2}AA^{*})AA^{*}+it^{3}g(t^{2}AA^{*})AA^{*}H=I+itH+l(t^{2}AA^{*})+itm(t^{2}AA^{*})H.$

Since

$e^{ixt}=1+ixt+l(t^{2}x^{2})+itm(t^{2}x^{2})x,$

we have

$\lVert\widehat{Z}(Ht,At)-e^{iHt}\rVert=\lVert
l(t^{2}AA^{*})-l(t^{2}H^{2})+itm(t^{2}AA^{*})H-itm(t^{2}H^{2})H\rVert\leq\lVert
l(t^{2}AA^{*})-l(t^{2}H^{2})\rVert+t\lVert
m(t^{2}AA^{*})-m(t^{2}H^{2})\rVert\lVert H\rVert.$

To bound the norms in $l,m$ we will apply Thm. 1.4.1 of [AP16]. The theorem
state that if a function $f\in L^{\infty}(\mathbb{R})$, i.e. $f$ is in the
function space which elements are the essentially bounded measurable
functions, it is entirely on $\mathbb{C}$ and satisfies $|f(z)|\leq
e^{\sigma|z|}$ for any $z\in\mathbb{C}$. Then $\lVert
f(A)-f(B)\rVert\leq\sigma\lVert f\rVert_{L^{\infty}(\mathbb{R})}\lVert
A-B\rVert$. Note that

$\displaystyle|l(z)|$
$\displaystyle=\left|\sum_{j=1}^{\infty}\frac{(-1)^{j}z^{j}}{(2j)!}\right|\leq\sum_{j=1}^{\infty}\frac{|z|^{j}}{(2j)!}\leq\sum_{j=1}^{\infty}\frac{|z|^{j}}{j!}\leq
e^{|z|},$ $\displaystyle|l(z)|$
$\displaystyle=\left|\sum_{j=1}^{\infty}\frac{(-1)^{j}z^{j}}{(2j+1)!}\right|\leq\sum_{j=1}^{\infty}\frac{|z|^{j}}{(2j+1)!}\leq\sum_{j=1}^{\infty}\frac{|z|^{j}}{j!}\leq
e^{|z|}.$

Moreover it is easy to see that $\lVert
l\rVert_{L^{\infty}(\mathbb{R})},\lVert m\rVert_{L^{\infty}(\mathbb{R})}\leq
2.$ So

$\lVert\widehat{Z}(Ht,At)-e^{iHt}\rVert\leq 2(1+t\lVert H\rVert)\lVert
t^{2}AA^{*}-t^{2}H^{2}\rVert=2t^{2}(1+t\lVert H\rVert)\lVert
AA^{*}-H^{2}\rVert.$

Now note that, by defining the random variable
$\zeta_{i}=\frac{1}{p(t_{i})}h_{t_{i}}h_{t_{i}}^{*},$ we have that

$\displaystyle AA^{*}=\frac{1}{M}\sum_{i=1}^{M}\zeta_{i},$
$\displaystyle\mathbb{E}[\zeta_{i}]=\sum_{q=1}^{2^{n}}p(q)\frac{1}{p(q)}h_{t_{i}}h_{t_{i}}^{*}=H^{2},\leavevmode\nobreak\
\leavevmode\nobreak\ \forall i.$

Let $\tau>0$. By applying Thm. 1 of [Hsu14] (or Prop. 9 in [RCR15a]), for
which

$\lVert AA^{*}-H^{2}\rVert\leq\sqrt{\frac{\lVert H\rVert^{2}_{F}\lVert
H\rVert^{2}\tau}{M}}+\frac{2\lVert H\rVert^{2}_{F}\lVert H\rVert^{2}\tau}{M},$

with probability at least $1-4\frac{\lVert H\rVert^{2}_{F}}{\lVert
H\rVert^{2}}\tau/(e^{\tau}-\tau-1)$. Now since

$1-4\frac{\lVert H\rVert^{2}_{F}}{\lVert
H\rVert^{2}}\tau/(e^{\tau}-\tau-1)\geq 1-e^{\tau},$

when $\tau\geq e$, by selecting

$\tau=2\log\frac{4\lVert H\rVert^{2}_{F}}{\lVert H\rVert\delta},$

we have that the equation above holds with probability at least $1-\delta$.
Let $\eta>0$, by selecting

$M=4\eta^{-2}\lVert H\rVert^{2}_{F}\lVert H\rVert^{2}\tau,$

we then obtain

$\lVert AA^{*}-H^{2}\rVert\leq\eta,$

with probability at least $1-\delta$.

Now we study $\lVert\widehat{Z}_{K}(At,Ht)-\widehat{Z}(Ht,At)\rVert$. Denote
with $a$, $b$ the functions $a(x)=l(x^{2})$, $b(x)=m(x^{2})$ and with
$a_{K},b_{K}$ the associated $K$-truncated Taylor expansions. Note that
$a(x)=\cos(x)-1$, while $b(x)=(\sin(x)-x)/x$. Now by definition of
$\widehat{Z}_{K}$ and $\widehat{Z}$, we have

$\displaystyle\lVert\widehat{Z}_{K}(At,Ht)-$
$\displaystyle\widehat{Z}(Ht,At)\rVert$ $\displaystyle\leq\lVert
a_{K}(t\sqrt{AA^{*}})-a(t\sqrt{AA^{*}})\rVert+t\lVert H\rVert\lVert
b_{K}(t\sqrt{AA^{*}})-b(t\sqrt{AA^{*}})\rVert.$

Note that, since

$\sum_{j}x^{2j}/(2j)!=\cosh(x)\leq 2e^{|x|},$

and

$\displaystyle|(a_{K}-a)(x)|$
$\displaystyle=\left|\sum_{j=K+2}(-1)^{j}\frac{x^{2j}}{(2j)!}\right|$
$\displaystyle\leq\frac{|x|^{2K+4}}{(2K+4)!}\sum_{j=0}\frac{|x|^{2j}}{(2j)!}\frac{(2K+4)!(2j)!}{(2j+2K+4)!}$
$\displaystyle\leq\frac{2|x|^{2K+4}e^{|x|}}{(2K+4)!},$
$\displaystyle|(b_{K}-b)(x)|$
$\displaystyle=\left|\sum_{j=K+2}(-1)^{j}\frac{x^{2j}}{(2j+1)!}\right|$
$\displaystyle\leq\frac{|x|^{2K+4}}{(2K+4)!}\sum_{j=0}\frac{|x|^{2j}}{(2j)!}\frac{(2K+4)!(2j)!}{(2j+2K+5)!}$
$\displaystyle\leq\frac{2|x|^{2K+4}e^{|x|}}{(2K+4)!}.$

Let $R>0,\beta\in(0,1]$. Now note that, by Stirling approximation, $c!\geq
e^{c\log\frac{c}{e}}$, so by selecting
$K=\frac{e^{2}}{2}R+\log(\frac{1}{\beta})$, we have for any $|x|\leq R$,

$\displaystyle\log\left|\frac{2|x|^{2K+4}e^{|x|}}{(2K+4)!}\right|$
$\displaystyle\leq|x|+(2K+4)\log\frac{e|x|}{2K+4}$ $\displaystyle\leq
R+(2K+4)\log\frac{eR}{2K+4}$ $\displaystyle\leq
R-\left(e^{2}R+4+\log\frac{1}{\beta}\right)\log\left(e+\frac{4+\log{1}{\beta}}{eR}\right)$
$\displaystyle\leq R-\left(e^{2}R+4+\log\frac{1}{\beta}\right)$
$\displaystyle\leq-(e^{2}-1)R-4-\log\frac{1}{\beta}$
$\displaystyle\leq-\log\frac{1}{\beta}.$

So, by choosing $K\geq\log\frac{1}{\beta}+e^{2}R/2$, we have
$|a_{K}(x)-a(x)|,|b_{K}(x)-b(x)|\leq\beta$. With this we finally obtain

$\lVert\widehat{Z}_{K}(At,Ht)-\widehat{Z}(Ht,At)\rVert\leq 2(1+t\lVert
H\rVert)\beta,$

when $K\geq e^{2}t\lVert A\rVert/2+\log\frac{1}{\beta}$, and therefore we have

$\lVert\widehat{\psi}_{K,M}(t)-\psi(t)\rVert\leq 2t^{2}(1+t\lVert
H\rVert)\eta+2(1+t\lVert H\rVert)\beta,$

with probability at least $1-\delta$, when

$M\geq 8\eta^{-2}\lVert H\rVert^{2}_{F}\lVert H\rVert^{2}\log\frac{4\lVert
H\rVert^{2}_{F}}{\lVert H\rVert^{2}\delta},\quad K\geq 4t\lVert
A\rVert+\log\frac{1}{\beta}.$

In particular, by choosing $\eta=\epsilon/(4t^{2}(1+t\lVert H\rVert)$ and
$\beta=\epsilon/(4(1+t\lVert H\rVert))$, we have

$\lVert\widehat{\psi}_{K,M}(t)-\psi(t)\rVert\leq\epsilon,$

with probability at least $1-\delta$.

With this result in mind note that, in the event where $\lVert
AA^{*}-H^{2}\rVert\leq\epsilon$, we have that

$|\lVert AA^{*}\rVert-\lVert H\rVert^{2}|\leq\lVert
AA^{*}-H^{2}\rVert\leq\epsilon,$

and therefore, $\lVert A\rVert\leq\sqrt{\lVert H\rVert^{2}+\epsilon}$. ∎

##  6 Application to density matrix simulation

Sample-based Hamiltonian simulation is a method for simulating Hamiltonians
which are density matrices [LMR14]. Such method has been used in several
recent quantum machine learning algorithms, including quantum support vector
machines [RML14], quantum gradient descent / Newton’s method [Reb+19], and
quantum linear regression [SSP16]. These techniques all make use of a quantum
algorithm that implements the time evolution $e^{i\rho t}\ket{\psi}$ governed
by a Hamiltonian corresponding to the density matrix $\rho$ of a pure
$n$-qubit state $\ket{\psi}$ of dimension $2^{n}$. Specifically, given an
oracle that returns superpositions of the entries of a density matrix
$\rho\in\mathbb{C}^{2^{n}\times 2^{n}}$, the quantum algorithm for density
matrix exponentiation implements a unitary matrix $U$, such that
$\|U\ket{\psi}-e^{i\rho t}\ket{\psi}\|\leq\epsilon$ and requires
$O(t^{2}/\epsilon)$ copies of $\rho$. This has been shown to be optimal
[Kim+17].

The total runtime of the density matrix simulation algorithm is
$O(t^{2}T(U_{\rho})/\epsilon)$, where $T(U_{\rho})$ is the time required to
prepare a quantum superposition of the entries of the matrix $\rho$. If the
state preparation can be achieved in $O(\poly(n))$ time, the total runtime of
the algorithm reduces to $O(\poly(n)t^{2}/\epsilon)$, which is polylogarithmic
in the dimension of $\rho$. A data structure that grants this sort of fast
access is common in the quantum machine learning literature [Bia+17, Cil+18].
As discussed in the introduction, it was recently noted by Tang [Tan19] that
some common implementations of this data structure, such as [KP16], are
equivalent to the ability to perform $\ell_{2}$-norm sampling on the entries
of $\rho$. This requirement is fundamentally equivalent to the row-searchable
condition used in this paper.

More formally, the sample-based Hamiltonian simulation Theorem states that

######  Theorem 7 (Sample-based Hamiltonian simulation [LMR14]).

Given access to multiple copies of an $n$-qubit state $\rho$, there is an
efficient quantum circuit implementing a unitary $U$, such that
$\|U\psi-e^{i\rho t}\psi\|\leq\epsilon$, where $0\leq\epsilon\leq 1$, $t>0$,
and $\psi$ is an arbitrary pure state. The algorithm requires
$O(t^{2}/\epsilon)$ copies of $\rho$.

Note that the full running time of the algorithm depends on the memory access
mode and the time to prepare $\rho$.

Theorem 5, i.e. the algorithm for positive semidefinite matrices, suggests a
straightforward classical analogue of sample-based Hamiltonian simulation
which outputs a classical description of the evolved quantum state, under the
condition that we can efficiently compute marginals of the diagonal entries of
$\rho$. Note that this is the case for any density matrix with structured
diagonal entries (_e.g._ , non-decreasing order). More precisely, we obtain a
classical efficient algorithm for density matrix simulation with the following
properties (the next result is a formal statement of Theorem 3):

######  Theorem 8 (Classical sample-based Hamiltonian simulation).

Let $\epsilon,\delta\in(0,1]$, and let $\rho$ be a row-searchable $n$-qubit
density matrix with at most $s$ non-zero entries in each row. Given an
efficient classical description of an arbitrary input state $\psi$ then, with
probability $1-\delta$, Algorithm 2 returns an approximation of any chosen
amplitude of the state $\widehat{\psi}$ such that
$\lVert\widehat{\psi}-e^{i\rho t}\psi\rVert\leq\epsilon$. The algorithm runs
in time $O(\frac{st^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$ and requires
$O(\frac{t^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$ of memory.

###### Proof of Theorem 8.

The Theorem is an immediate consequence of the unitarity of the trace of
density matrices and Theorem 5 from which we know that it is possible to
compute an approximation $\widehat{\psi}(t)$ of $\psi(t)$, with error at most
$\epsilon$ and with probability at least $1-\delta$, in
$O(\frac{st^{2}\mathrm{Tr}\left(H\right)^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$
time and
$O(\frac{t^{2}\mathrm{Tr}\left(H\right)^{2}}{\epsilon^{2}}\log^{2}\frac{1}{\delta})$
memory. ∎

#### Acknowledgements

We acknowledge Dorit Aharonov, Leonardo Banchi, Sougato Bose, Fernando
Brandão, Giuseppe Carleo, Vladimir Korepin, Laura Mančinska, Tzu-Chieh Wei for
helpful conversations. A.R. is supported by EPSRC and by QinetiQ. C.C. and
M.P. are supported by EPSRC. S.S. is supported by the Royal Society, EPSRC,
the National Natural Science Foundation of China, and the grant ARO-MURI
W911NF-17-1-0304 (US DOD, UK MOD and UK EPSRC under the Multidisciplinary
University Research Initiative). L.W. is supported by the Royal Society.

## References

  * [AT03] Dorit Aharonov and Amnon Ta-Shma  “Adiabatic Quantum State Generation and Statistical Zero Knowledge”  In _Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing_ San Diego, CA, USA: Association for Computing Machinery, 2003, pp. 20–29  DOI: 10.1145/780542.780546
  * [AP16] Alexei Borisovich Aleksandrov and Vladimir Vsevolodovich Peller  “Operator Lipschitz functions”  In _Russian Mathematical Surveys_ 71.4 IOP Publishing, 2016, pp. 605  DOI: 10.1070/RM9729
  * [BW07] Mohamed-Ali Belabbas and Patrick J Wolfe  “Fast low-rank approximation for covariance matrices”  In _2nd IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2007._ , 2007, pp. 293–296  IEEE  DOI: 10.1109/CAMSAP.2007.4498023
  * [BW08] Mohamed-Ali Belabbas and Patrick J Wolfe  “On sparse representations of linear operators and the approximation of matrix products”  In _42nd Annual Conference on Information Sciences and Systems_ , 2008, pp. 258–263  IEEE  DOI: 10.1109/CISS.2008.4558532
  * [BCK15] Dominic W Berry, Andrew M Childs and Robin Kothari  “Hamiltonian simulation with nearly optimal dependence on all parameters”  In _IEEE 56th Annual Symposium on Foundations of Computer Science (FOCS)_ , 2015, pp. 792–809  IEEE  DOI: 10.1109/FOCS.2015.54
  * [Bia+17] Jacob Biamonte et al.  “Quantum machine learning”  In _Nature_ 549.7671 Nature Publishing Group, 2017, pp. 195  DOI: 10.1038/nature23474
  * [Bra+18] Sergey Bravyi et al.  “Simulation of quantum circuits by low-rank stabilizer decompositions”  In _arXiv preprint arXiv:1808.00128_ , 2018
  * [Chi+19] Nai-Hui Chia et al.  “Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning”  In _arXiv preprint arXiv:1910.06151_ , 2019
  * [CK11] Andrew M. Childs and Robin Kothari  “Simulating Sparse Hamiltonians with Star Decompositions”  In _Proceedings of the 5th Conference on Theory of Quantum Computation, Communication, and Cryptography_ Leeds, UK: Springer-Verlag, 2011, pp. 94–103  DOI: 10.1007/978-3-642-18073-6_8
  * [CW12] Andrew M. Childs and Nathan Wiebe  “Hamiltonian Simulation Using Linear Combinations of Unitary Operations”  In _Quantum Info. Comput._ 12.11-12 Paramus, NJ: Rinton Press, Incorporated, 2012, pp. 901–924  DOI: 10.5555/2481569.2481570
  * [Chi+03] Andrew M. Childs et al.  “Exponential Algorithmic Speedup by a Quantum Walk”  In _Proceedings of the Thirty-Fifth Annual ACM Symposium on Theory of Computing_ San Diego, CA, USA: Association for Computing Machinery, 2003, pp. 59–68  DOI: 10.1145/780542.780552
  * [Cil+18] Carlo Ciliberto et al.  “Quantum machine learning: a classical perspective”  In _Proc. R. Soc. A_ 474.2209 The Royal Society, 2018, pp. 20170551  DOI: 10.1098/rspa.2017.0551
  * [DKM06] Petros Drineas, Ravi Kannan and Michael W Mahoney  “Fast Monte Carlo algorithms for matrices I: Approximating matrix multiplication”  In _SIAM Journal on Computing_ 36.1 SIAM, 2006, pp. 132–157  DOI: 10.1137/S0097539704442684
  * [DKM06a] Petros Drineas, Ravi Kannan and Michael W Mahoney  “Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix”  In _SIAM Journal on computing_ 36.1 SIAM, 2006, pp. 158–183  DOI: 10.1137/S0097539704442696
  * [DM05] Petros Drineas and Michael W. Mahoney  “On the Nyström Method for Approximating a Gram Matrix for Improved Kernel-Based Learning”  In _J. Mach. Learn. Res._ 6 JMLR.org, 2005, pp. 2153–2175  DOI: 10.5555/1046920.1194916
  * [DM18] Petros Drineas and Michael W Mahoney  “Lectures on randomized numerical linear algebra”  In _The Mathematics of Data_ 25 American Mathematical Soc., 2018, pp. 1
  * [Dri+11] Petros Drineas, Michael W. Mahoney, S. Muthukrishnan and Tamás Sarlós  “Faster Least Squares Approximation”  In _Numer. Math._ 117.2 Berlin, Heidelberg: Springer-Verlag, 2011, pp. 219–249  DOI: 10.1007/s00211-010-0331-6
  * [Fow+04] Charless Fowlkes, Serge Belongie, Fan Chung and Jitendra Malik  “Spectral Grouping Using the Nyström Method”  In _IEEE Trans. Pattern Anal. Mach. Intell._ 26.2 USA: IEEE Computer Society, 2004, pp. 214–225  DOI: 10.1109/TPAMI.2004.1262185
  * [FKV04] Alan Frieze, Ravi Kannan and Santosh Vempala  “Fast Monte-Carlo Algorithms for Finding Low-Rank Approximations”  In _J. ACM_ 51.6 New York, NY, USA: Association for Computing Machinery, 2004, pp. 1025–1041  DOI: 10.1145/1039488.1039494
  * [Hae+11] Jutho Haegeman et al.  “Time-dependent variational principle for quantum lattices”  In _Physical review letters_ 107.7 APS, 2011, pp. 070601  DOI: 10.1103/PhysRevLett.107.070601
  * [Hig05] Nicholas J. Higham  “The Scaling and Squaring Method for the Matrix Exponential Revisited”  In _SIAM J. Matrix Anal. Appl._ 26.4 USA: Society for IndustrialApplied Mathematics, 2005, pp. 1179–1193  DOI: 10.1137/04061101X
  * [Hig09] Nicholas J. Higham  “The Scaling and Squaring Method for the Matrix Exponential Revisited”  In _SIAM Rev._ 51.4 USA: Society for IndustrialApplied Mathematics, 2009, pp. 747–764  DOI: 10.1137/090768539
  * [Hsu14] Daniel Hsu  “Weighted sampling of outer products”  In _arXiv preprint arXiv: 1410.4429_ , 2014
  * [HNS18] Cupjin Huang, Michael Newman and Mario Szegedy  “Explicit lower bounds on strong quantum simulation”  In _arXiv preprint arXiv:1804.10368_ , 2018
  * [JBC18] Bjarni Jónsson, Bela Bauer and Giuseppe Carleo  “Neural-network states for the classical simulation of quantum computing”  In _arXiv preprint arXiv:1808.05232_ , 2018
  * [KP16] Iordanis Kerenidis and Anupam Prakash  “Quantum recommendation systems”  In _arXiv preprint arXiv:1603.08675_ , 2016
  * [Kim+17] Shelby Kimmel et al.  “Hamiltonian simulation with optimal sample complexity”  In _npj Quantum Information_ 3.1 Nature Publishing Group, 2017, pp. 13  DOI: 10.1038/s41534-017-0013-7
  * [KMT09] Sanjiv Kumar, Mehryar Mohri and Ameet Talwalkar  “On Sampling-Based Approximate Spectral Decomposition”  In _Proceedings of the 26th Annual International Conference on Machine Learning_ Montreal, Quebec, Canada: Association for Computing Machinery, 2009, pp. 553–560  DOI: 10.1145/1553374.1553446
  * [KMT12] Sanjiv Kumar, Mehryar Mohri and Ameet Talwalkar  “Sampling methods for the Nyström method”  In _Journal of Machine Learning Research_ 13.Apr, 2012, pp. 981–1006  DOI: 10.5555/2188385.2343678
  * [LKL10] Mu Li, James T. Kwok and Bao-Liang Lu  “Making Large-Scale Nyström Approximation Possible”  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ Haifa, Israel: Omnipress, 2010, pp. 631–638  DOI: 10.5555/3104322.3104403
  * [Llo96] Seth Lloyd  “Universal Quantum Simulators”  In _Science_ 273.5278, 1996, pp. 1073–1078  DOI: 10.1126/science.273.5278.1073
  * [LMR14] Seth Lloyd, Masoud Mohseni and Patrick Rebentrost  “Quantum principal component analysis”  In _Nature Physics_ 10.9 Nature Publishing Group, 2014, pp. 631  DOI: 10.1038/nphys3029
  * [LC17] Guang Hao Low and Isaac L Chuang  “Hamiltonian Simulation by Uniform Spectral Amplification”  In _arXiv preprint arXiv:1707.05391_ , 2017
  * [MTJ11] Lester Mackey, Ameet Talwalkar and Michael I. Jordan  “Divide-and-Conquer Matrix Factorization”  In _Proceedings of the 24th International Conference on Neural Information Processing Systems_ Granada, Spain: Curran Associates Inc., 2011, pp. 1134–1142  DOI: 10.5555/2986459.2986586
  * [Mah+11] Michael W Mahoney  “Randomized algorithms for matrices and data”  In _Foundations and Trends® in Machine Learning_ 3.2 Now Publishers, Inc., 2011, pp. 123–224  DOI: 10.1561/2200000035
  * [Mat93] Roy Mathias  “Approximation of matrix-valued functions”  In _SIAM journal on matrix analysis and applications_ 14.4 SIAM, 1993, pp. 1061–1063  DOI: 10.1137/0614070
  * [AH09] Awad H Al-Mohy and Nicholas J Higham  “A new scaling and squaring algorithm for the matrix exponential”  In _SIAM Journal on Matrix Analysis and Applications_ 31.3 SIAM, 2009, pp. 970–989  DOI: 10.1137/09074721X
  * [AH11] Awad H Al-Mohy and Nicholas J Higham  “Computing the action of the matrix exponential, with an application to exponential integrators”  In _SIAM journal on scientific computing_ 33.2 SIAM, 2011, pp. 488–511  DOI: 10.1137/100788860
  * [Nak03] Ritsuo Nakamoto  “A norm inequality for Hermitian operators”  In _The American mathematical monthly_ 110.3 Mathematical Association of America, 2003, pp. 238  DOI: 10.1080/00029890.2003.11919961
  * [Nys30] Evert J Nyström  “Über die praktische Auflösung von Integralgleichungen mit Anwendungen auf Randwertaufgaben”  In _Acta Mathematica_ 54.1 Springer, 1930, pp. 185–204  DOI: 10.1007/BF02547521
  * [OSV12] Lorenzo Orecchia, Sushant Sachdeva and Nisheeth K. Vishnoi  “Approximating the Exponential, the Lanczos Method and an Õ(m)-Time Spectral Algorithm for Balanced Separator”  In _Proceedings of the Forty-Fourth Annual ACM Symposium on Theory of Computing_ New York, New York, USA: Association for Computing Machinery, 2012, pp. 1141–1160  DOI: 10.1145/2213977.2214080
  * [Orú14] Román Orús  “A practical introduction to tensor networks: Matrix product states and projected entangled pair states”  In _Annals of Physics_ 349 Elsevier, 2014, pp. 117–158  DOI: 10.1016/j.aop.2014.06.013
  * [RML14] Patrick Rebentrost, Masoud Mohseni and Seth Lloyd  “Quantum support vector machine for big data classification”  In _Physical review letters_ 113.13 APS, 2014, pp. 130503  DOI: 10.1103/PhysRevLett.113.130503
  * [Reb+19] Patrick Rebentrost et al.  “Quantum gradient descent and Newton’s method for constrained polynomial optimization”  In _New Journal of Physics_ 21.7 IOP Publishing, 2019, pp. 073023  DOI: 10.1088/1367-2630/ab2a9e
  * [RCR15] Alessandro Rudi, Raffaello Camoriano and Lorenzo Rosasco  “Less is More: Nyström Computational Regularization”  In _Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume 1_ Montreal, Canada: MIT Press, 2015, pp. 1657–1665  DOI: 10.5555/2969239.2969424
  * [RCR15a] Alessandro Rudi, Raffaello Camoriano and Lorenzo Rosasco  “Less is More: Nyström Computational Regularization”  In _arXiv preprint arXiv:1507.04717_ , 2015
  * [SSP16] Maria Schuld, Ilya Sinayskiy and Francesco Petruccione  “Prediction by linear regression on a quantum computer”  In _Physical Review A_ 94.2 APS, 2016, pp. 022342  DOI: 10.1103/PhysRevA.94.022342
  * [SN13] Martin Schwarz and Maarten Van den Nest  “Simulating quantum circuits with sparse output distributions”  In _arXiv preprint arXiv:1310.6749_ , 2013
  * [ST04] Daniel A. Spielman and Shang-Hua Teng  “Nearly-Linear Time Algorithms for Graph Partitioning, Graph Sparsification, and Solving Linear Systems”  In _Proceedings of the Thirty-Sixth Annual ACM Symposium on Theory of Computing_ Chicago, IL, USA: Association for Computing Machinery, 2004, pp. 81–90  DOI: 10.1145/1007352.1007372
  * [ST11] Daniel A Spielman and Shang-Hua Teng  “Spectral sparsification of graphs”  In _SIAM Journal on Computing_ 40.4 SIAM, 2011, pp. 981–1025  DOI: 10.1137/08074489X
  * [Suz76] Masuo Suzuki  “Generalized Trotter’s formula and systematic approximants of exponential operators and inner derivations with applications to many-body problems”  In _Communications in Mathematical Physics_ 51.2 Springer, 1976, pp. 183–190  DOI: 10.1007/BF01609348
  * [TKR08] Ameet Talwalkar, Sanjiv Kumar and Henry Rowley  “Large-scale manifold learning”  In _IEEE Conference on Computer Vision and Pattern Recognition_ , 2008, pp. 1–8  IEEE  DOI: 10.1109/CVPR.2008.4587670
  * [Tan19] Ewin Tang  “A Quantum-Inspired Classical Algorithm for Recommendation Systems”  In _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ Phoenix, AZ, USA: Association for Computing Machinery, 2019, pp. 217–228  DOI: 10.1145/3313276.3316310
  * [Tro12] Joel A. Tropp  “User-Friendly Tail Bounds for Sums of Random Matrices”  In _Found. Comput. Math._ 12.4 Berlin, Heidelberg: Springer-Verlag, 2012, pp. 389–434  DOI: 10.1007/s10208-011-9099-z
  * [Tro59] Hale F Trotter  “On the product of semi-groups of operators”  In _Proceedings of the American Mathematical Society_ 10.4, 1959, pp. 545–551  DOI: 10.1090/S0002-9939-1959-0108732-6
  * [Van10] Maarten Van Den Nest  “Classical simulation of quantum computation, the Gottesman-Knill theorem, and slightly beyond”  In _Quantum Information & Computation_ 10.3 Rinton Press, Incorporated, 2010, pp. 258–271  DOI: 10.5555/2011350.2011356
  * [VGC04] Frank Verstraete, Juan J Garcia-Ripoll and Juan Ignacio Cirac  “Matrix product density operators: Simulation of finite-temperature and dissipative systems”  In _Physical review letters_ 93.20 APS, 2004, pp. 207204  DOI: 10.1103/PhysRevLett.93.207204
  * [VMC08] Frank Verstraete, Valentin Murg and J Ignacio Cirac  “Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems”  In _Advances in Physics_ 57.2 Taylor & Francis, 2008, pp. 143–224  DOI: 10.1063/1.5026985
  * [Vid04] Guifré Vidal  “Efficient Simulation of One-Dimensional Quantum Many-Body Systems”  In _Phys. Rev. Lett._ 93 American Physical Society, 2004, pp. 040502  DOI: 10.1103/PhysRevLett.93.040502
  * [WS01] Christopher K.. Williams and Matthias Seeger  “Using the Nyström Method to Speed Up Kernel Machines”  In _Advances in Neural Information Processing Systems 13_ MIT Press, 2001, pp. 682–688  DOI: 10.5555/3008751.3008847
  * [Wil+02] CKI. Williams, CE. Rasmussen, A. Schwaighofer and V. Tresp  “Observations on the Nyström Method for Gaussian Process Prediction”, 2002  Max-Planck-Gesellschaft
  * [Woo+14] David P Woodruff  “Sketching as a tool for numerical linear algebra”  In _Foundations and Trends® in Theoretical Computer Science_ 10.1–2 Now Publishers, Inc., 2014, pp. 1–157  DOI: 10.1561/0400000060
  * [ZK10] Kai Zhang and James T Kwok  “Clustered Nyström method for large scale manifold learning and dimension reduction”  In _IEEE Transactions on Neural Networks_ 21.10 IEEE, 2010, pp. 1576–1587  DOI: 10.1109/TNN.2010.2064786
  * [ZTK08] Kai Zhang, Ivor W. Tsang and James T. Kwok  “Improved Nyström Low-Rank Approximation and Error Analysis”  In _Proceedings of the 25th International Conference on Machine Learning_ Helsinki, Finland: Association for Computing Machinery, 2008, pp. 1232–1239  DOI: 10.1145/1390156.1390311