# Quantum algorithms for matrix geometric means

Nana Liu<EMAIL_ADDRESS>Institute of Natural Sciences, School of
Mathematical Sciences, MOE-LSC, Shanghai Jiao Tong University, Shanghai,
200240, P. R. China  Shanghai Artificial Intelligence Laboratory, Shanghai,
China  University of Michigan-Shanghai Jiao Tong University Joint Institute,
Shanghai 200240, China  Qisheng Wang<EMAIL_ADDRESS>Graduate
School of Mathematics, Nagoya University, Nagoya 464-8602, Japan  Mark M.
Wilde<EMAIL_ADDRESS>School of Electrical and Computer Engineering,
Cornell University, Ithaca, New York 14850, USA  Zhicheng Zhang
<EMAIL_ADDRESS>Centre for Quantum Software and Information, University of
Technology Sydney, Ultimo, NSW 2007, Australia

###### Abstract

Matrix geometric means between two positive definite matrices can be defined
equivalently from distinct perspectives – as solutions to certain nonlinear
systems of equations, as points along geodesics in Riemannian geometry, and as
solutions to certain optimisation problems. This diversity already suggests
the potential for varied applications, as well as acting as a bridge between
different domains. Here we devise new quantum subroutines to efficiently
prepare quantum unitary operators that embed the standard matrix geometric
mean and its generalisations called the weighted matrix geometric mean. This
enables the construction of solutions to the algebraic Riccati equation, which
is an important class of nonlinear systems of equations that appears in
machine learning, optimal control, estimation, and filtering. Using these
subroutines, we present a new class of quantum learning algorithms called
quantum geometric mean metric learning. This has applications in efficiently
finding the best distance measure and solving classification problems in the
weakly supervised limit and for anomaly detection, for both classical and
quantum problems. We also show how our method can be generalised to a
particular $p^{\text{th}}$-order system of nonlinear equations. These quantum
subroutines for matrix geometric means are also useful in other areas of
quantum information. For example, we show how to use them in the estimation of
geometric Rényi relative entropies and the Uhlmann fidelity by means of the
Fuchs–Caves observable. In particular, our quantum algorithms for estimating
the Uhlmann and Matsumoto fidelities have optimal dependence on the precision.
Finally, we provide a $\mathsf{BQP}$-complete problem based on matrix
geometric means that can be solved by our subroutines, thus characterising
their computational capability.

###### Contents

  1. I Introduction
    1. I.1 Summary of our results
  2. II Background
    1. II.1 Matrix geometric means
    2. II.2 Algebraic Riccati equations
    3. II.3 Block-encoding
  3. III Quantum subroutines for matrix geometric means, algebraic Riccati equations, and higher-order nonlinear equations
    1. III.1 Quantum subroutine for matrix geometric means
    2. III.2 $B=0$ algebraic Riccati equation
    3. III.3 $B\neq 0$ algebraic Riccati equation
    4. III.4 Higher-order polynomial equations
  4. IV Applications
    1. IV.1 Quantum geometric mean metric learning
      1. IV.1.1 Learning Euclidean metric from data
      2. IV.1.2 1-class quantum learning
      3. IV.1.3 Extension to weighted geometric mean metric learning
    2. IV.2 Estimation of quantum fidelity and geometric Rényi relative entropies
      1. IV.2.1 Fidelity
      2. IV.2.2 Geometric fidelity and geometric Rényi relative entropy
  5. V BQP-hardness
  6. VI Discussion
  7. A Proof of Lemma 5
  8. B Preliminary lemmas of the block-encoding formalism and other useful results
  9. C Proof of Lemma 7
  10. D Proof of Lemma 8
  11. E Proof of Lemma 9
  12. F Proof of Lemma 11
  13. G Proof of Lemma 12
  14. H Proof of Theorems 18 and 19
  15. I Proof of Lemmas 17 and 20
  16. J Proof of Lemma 23

##  I Introduction

Quantum computation is considered a rapidly emerging technology that has
important implications for the development of algorithms. Many quantum
algorithms that have theoretically demonstrated potential quantum advantage,
however, have been chiefly directed towards linear problems – in part because
quantum mechanics is itself linear. These include simulating solutions of
linear systems of equations [1], known as quantum linear algebra, and linear
ordinary and partial differential equations [2, 3, 4, 5].

However, many problems of scientific interest are nonlinear. While most
nonlinear systems of equations of interest for applications only appear after
discretising nonlinear ordinary and partial differential equations, there is
an important class of nonlinear system of equations that is not only relevant
to partial differential equations but is also of independent interest. This
class consists of the algebraic Riccati equations, which are nonlinear matrix
equations with quadratic nonlinearity [6]. These are also the stationary
states of the Riccati matrix differential equations, which are essential for
many applications in applied mathematics, science, and engineering problems.
These nonlinear matrix equations are particularly relevant for optimal
control, stability theory, filtering (e.g., Kalman filter [7]), network
theory, differential games, and estimation problems [8].

It turns out that solutions to the algebraic Riccati equations are closely
connected with the concept of a matrix geometric mean. For example, the unique
solution to the simplest algebraic Riccati equation can be precisely expressed
as the standard matrix geometric mean, as we will recall later. The matrix
geometric means are matrix generalisations of the scalar geometric mean and
have a long history in mathematics [9, 10]; there are diverse approaches to
this same concept. For example, the standard matrix geometric mean can be
defined as the output of an optimisation problem. The matrix geometric mean
between two matrices also has an elegant geometric interpretation as a
midpoint along the geodesic joining these two matrices that live in Riemannian
space [6]. The Monge map between two Gaussian distributions, appearing in
optimal transport, can also be expressed in terms of the matrix geometric mean
[11]. The standard and weighted matrix geometric means appear in quantum
information in the form of quantum entropic [12, 13, 14] and fidelity [15, 16,
17] measures.

However, computing the matrix geometric mean involves matrix multiplication
and also nonlinear operations like taking inverses and square roots of
matrices. Here classical numerical schemes can be inefficient, with costs that
are polynomial in the size of the matrix [18]. The processing of several
matrix multiplications can, under certain conditions, be more efficient
through quantum processing. Our aim here is to construct quantum subroutines
that embed the standard and weighted matrix geometric means into unitary
operators and to determine the conditions under which these embeddings can
indeed be conducted efficiently. There are many such possible unitary
operators, and we choose a formalism called block-encoding [19, 20, 21].

The block-encoding of a non-unitary matrix $Y$ is a unitary matrix $U_{Y}$
whose upper left-hand corner is proportional to $Y$. The construction of this
unitary matrix allows realisation by means of a quantum circuit, which
describes unitary evolution. The matrix $Y$ can be subsequently recovered by
extracting only the top-left corner through measurement. This provides a
convenient building block for constructing sums and (integer and non-integer)
powers of matrices $Y$ by concatenating its block-encodings via unitary
circuits. This formalism allows us to form the block-encoding of the standard
and weighted matrix geometric means, which are products of matrices and their
roots. From these block-encodings, one can also recover their expectation
values with respect to certain states. These different expectation values are
then relevant for various applications, like in machine learning and quantum
fidelity estimation.

Under certain assumptions, we show how these can be efficiently implementable
on quantum devices. This efficiency arises from the fact that matrix
multiplications can be more efficient with quantum algorithms. This
observation has an important consequence. It means that a quantum device can
efficiently prepare solutions of the (nonlinear) algebraic Riccati equations.
The expectation values of these solutions can also be shown to be efficiently
recoverable for different applications. Our approach differs from many past
works in three key respects: (a) ours is the first quantum subroutine, to the
best of our knowledge, to prepare solutions of nonlinear matrix equations
without using iterative methods. The solutions themselves are matrices and not
vectors, which differs from other quantum algorithms for nonlinear systems of
equations, for example [22, 23, 24]; (b) the solutions are not embedded in a
pure quantum state, but rather an observable, thus introducing a novel
embedding of the solution. This is important when solutions themselves are in
matrix form (for matrix equations), which differs from the quantum embeddings
of solutions of discretised nonlinear ordinary and partial differential
equations (solutions not in matrix form) [22, 23, 24]; (c) we show the
efficient recovery of outputs for nonlinear systems of equations directly
relevant for applications.

One class of applications is in the area of machine learning. Machine learning
algorithms often require an assignment of a metric, or distance measure, in
order to compute distances between data points. The values of these distances
then become central to the outcome, for instance, in making a prediction for
classification. This means that the choice of the metric itself is important,
but the best metric can depend on the actual data. Learning the metric from
given data – called metric learning – can also be formulated as a learning
problem. While most of these metric learning algorithms require iterative
techniques like gradient descent to minimise the proposed loss function, a
class of metric learning algorithms called geometric mean metric learning [25]
admits closed-form solutions. It has also been shown to attain higher
classification accuracy with greater speed than previous methods. Here we
devise efficient quantum algorithms, using our quantum subroutine for the
matrix geometric mean, for geometric mean metric learning for both classical
and quantum data. For quantum data, we propose new algorithms that can be used
for the anomaly detection of quantum states, which differs from previous
algorithms [26]. The applicability extends also to asymmetric cases for which
there is a higher cost to be paid for false negatives or true positives. This
is in fact related to the weighted matrix geometric mean.

There is also an important connection between the solution of the geometric
mean metric learning problem and the Fuchs–Caves observable [15], which
appears in quantum fidelity estimation. This allows for a re-derivation of
quantum fidelity from the point of view of machine learning. We show that our
quantum subroutines for the matrix geometric mean can also be used in the
efficient estimation of geometric Rényi relative entropies and the quantum
fidelity by means of the Fuchs–Caves observable. This new way of estimating
the quantum fidelity has polynomially better performance in precision than
previously known fidelity estimation algorithms. It is also shown to be
optimal with respect to precision.

We can also extend our method to a more general class of nonlinear systems of
equations of $p^{\text{th}}$-degree. These are $p^{\text{th}}$-degree
polynomial generalisations of the simplest algebraic Riccati equations. We
show that the unique solutions of these equations are weighted matrix
geometric means. We similarly devise quantum subroutines to prepare their
block-encodings. The weighted matrix geometric mean for two quantum states has
an elegant geometric interpretation as the state at $(1/p)^{\text{th}}$ of the
length along the geodesic connecting two quantum states in Riemannian space.
We also show these are relevant to the weighted version of our new quantum
learning algorithm. Furthermore, preparing block-encodings of the weighted
matrix geometric means allows us to construct, to the best of our knowledge,
the first quantum algorithm for estimating the geometric Rényi relative
entropies.

###  I.1 Summary of our results

For convenience, we provide a brief summary of our results here. Our first
contribution consists of basic quantum subroutines in Section III for matrix
geometric means (see Definition 1) and their weighted generalisation (see
Definition 2).

Solving algebraic Riccati equations. We then consider the problem of solving
the algebraic matrix Riccati equation

$YAY-B^{\dagger}Y-Y^{\dagger}B-C=0,$ (1)

where $A$, $B$, and $C$ are $D\times D$ complex-valued matrices. We delineate
quantum algorithms with time complexity $O\lparen\operatorname{poly}\log
D\rparen$ for solving Eq. (1) for well-conditioned matrices, in Section III.2
and Section III.3. Here, we say a matrix $A$ is well-conditioned if $A\geq
I/(\operatorname{poly}\log D)$. The higher-order case $Y\lparen
AY\rparen^{p-1}=C$ is studied in Section III.4. In Section V, we show that it
is $\mathsf{BQP}$-complete to solve the equation $YAY=C$, a special case of
Eq. (1), in which case the solution is $Y=A^{-1}\\#C$ (see Definition 1 for
the meaning of this notation).

Geometric mean metric learning. We introduce quantum algorithms for learning
the metric in machine learning, by phrasing this as an optimisation problem
using a geometric perspective. Unlike other metric learning algorithms, this
optimisation problem has a closed-form solution. This follows the geometric
mean metric learning method [25]. The solution turns out to be expressible in
terms of the matrix geometric mean $Y=A^{-1}\\#C$. We design quantum
algorithms for the learning task for classical data (Section IV.1.1) as well
as for quantum data (Section IV.1.2). We present the conditions under which
the quantum algorithm is more efficient than the corresponding classical
algorithm. The learning algorithm for quantum data is uniquely quantum in
nature and has no classical counterpart.

(Uhlmann) fidelity estimation. Based on the Fuchs–Caves observable [15], we
design a new quantum algorithm for fidelity estimation in Section IV.2.1 via
the fidelity formula
$F\lparen\rho,\sigma\rparen=\Tr\lparen\lparen\sigma^{-1}\\#\rho\rparen\sigma\rparen$,
which involves the matrix geometric mean. We prove that our algorithm has
$\widetilde{O}\lparen 1/\epsilon\rparen$ dependence on the precision
$\epsilon$ and that it is optimal up to polylogarithmic factors.

Geometric Rényi relative entropy. In Section IV.2.2, we present the first
quantum algorithm for computing the geometric Rényi relative entropy, to the
best of our knowledge. In particular, we design a quantum algorithm for
computing the geometric fidelity
$\hat{F}_{1/2}\lparen\rho,\sigma\rparen\coloneqq\Tr\lparen\rho\\#\sigma\rparen$
(also known as the Matsumoto fidelity [16, 17]) with $\widetilde{O}\lparen
1/\epsilon\rparen$ dependence on the precision $\epsilon$, and we prove that
it is optimal up to polylogarithmic factors.

Organisation of this paper. In Section II, we begin with a review of the
standard matrix geometric mean, weighted matrix geometric mean, the algebraic
Riccati equation, and block-encoding. In Section III we compute the costs
required to prepare block-encodings of the solutions of algebraic Riccati
equations and their $p^{\text{th}}$-order generalisations. Applications are
presented Section IV. In Section V we show how our new quantum subroutines for
the matrix geometric mean can solve a $\mathsf{BQP}$-complete problem. We end
in Section VI with discussions.

##  II Background

In this section we give a brief overview of the standard and weighted matrix
geometric means and their role in solving algebraic Riccati equations (see
[27, Chapters 4 & 6] and [9, 10] for more details). We then provide a
definition of block-encoding. Throughout the paper, unless otherwise stated,
we deal with Hermitian matrices.

###  II.1 Matrix geometric means

######  Definition 1 (Matrix geometric mean).

Fix $D\in\mathbb{N}$. Given two $D\times D$ positive definite matrices $A$ and
$C$, the matrix geometric mean of $A$ and $C$ is defined as

$\displaystyle A\\#C\coloneqq A^{1/2}(A^{-1/2}CA^{-1/2})^{1/2}A^{1/2}>0.$ (2)

Note that the matrix geometric mean between $A^{-1}$ and $C$ is thus defined
by

$\displaystyle A^{-1}\\#C=A^{-1/2}(A^{1/2}CA^{1/2})^{1/2}A^{-1/2}>0.$ (3)

Alternatively, the matrix geometric mean $A\\#C$ can be equivalently be
written as

$\displaystyle A\\#C=\max\Biggl{\\{}Y\geq 0:\begin{pmatrix}A&Y\\\
Y&C\end{pmatrix}\geq 0\Biggr{\\}},$ (4)

where the ordering of Hermitian matrices is given by the Löwner partial order.

The matrix geometric mean appears in quantum information, for example, like
the Fuchs–Caves observable [15], in quantum fidelity and entropy operators
like the Tsallis relative operator entropy [28], and quantum fidelity measures
between states [16, 12, 13, 17] and channels [29]. This concept can also be
generalised to the weighted matrix geometric mean.

######  Definition 2 (Weighted matrix geometric mean).

Fix $p>0$. The weighted matrix geometric mean with weight $1/p$ is defined as

$\displaystyle A\\#_{1/p}\,C\coloneqq
A^{1/2}(A^{-1/2}CA^{-1/2})^{1/p}A^{1/2}.$ (5)

The weighted matrix geometric mean between $A^{-1}$ and $C$ is then equal to

$\displaystyle A^{-1}\\#_{1/p}\,C=A^{-1/2}(A^{1/2}CA^{1/2})^{1/p}A^{-1/2}.$
(6)

The canonical matrix geometric mean corresponds to the weighted geometric mean
with weight $1/p=1/2$.

We will use the definitions in Eqs. (3) and (6) here and throughout because,
as we will see later on, they are relevant to solutions of classes of
nonlinear matrix equations like the algebraic Riccati equations.

For positive definite matrices (which include full-rank density matrices), the
standard and weighted matrix geometric means have elegant geometric
interpretations. It is known that the inner product on the real vector space
formed by the set of Hermitian matrices gives rise to a Riemannian metric [27,
Chapter 6]. This Riemannian metric is defined on the manifold $M_{H}$ formed
by the set of positive definite matrices. Following [27, Eqs. (6.2) & (6.4)],
a trajectory $\gamma:[a,b]\rightarrow M_{H}$ on this manifold is a piecewise
differential path on $M_{H}$ whose length is defined by
$L(\gamma)\coloneqq\int_{a}^{b}\left\|\gamma^{-1/2}(t)\gamma^{\prime}(t)\gamma^{-1/2}(t)\right\|_{2}\,dt$.
Then the distance $\delta(A^{-1},C)=\inf_{\gamma}L(\gamma)$ between any two
positive definite matrices $A^{-1}$ and $C$ on this manifold is defined to be
shortest length joining these two points. Then we have the following result.

######  Lemma 3 ([27, Theorem 6.1.6]).

If $A^{-1}$ and $C$ are two positive definite matrices, then there exists a
unique geodesic joining $A^{-1}$ and $C$. This geodesic has the following
parameterisation with $t\in[0,1]$:

$\displaystyle\gamma_{\operatorname{geod}}(t)=A^{-1/2}(A^{1/2}CA^{1/2})^{1/t}A^{-1/2},\qquad
t\in[0,1].$ (7)

This geodesic has length given by

$\displaystyle\delta(A^{-1},C)=L(\gamma_{\operatorname{geod}})=\left\|\log(A^{1/2}CA^{1/2})\right\|_{2}.$
(8)

In the above, $\|X\|_{2}\coloneqq\sqrt{\operatorname{Tr}[X^{\dagger}X]}$
denotes the Schatten $2$-norm, whereas $\|\cdot\|$ refers to operator norm
throughout our paper.

From this viewpoint, the matrix geometric mean
$A^{-1}\\#C=\gamma_{\operatorname{geod}}(t=1/2)$ can clearly be interpreted as
the midpoint along the geodesic joining $A^{-1}$ and $C$. Similarly, the
weighted geometric mean with weight $1/p$ can be interpreted as the point
along the manifold when $t=1/p$.

###  II.2 Algebraic Riccati equations

Let us begin with a general form of the algebraic Riccati equation for the
unknown $D\times D$ matrix $Y$:

$\displaystyle Y^{\dagger}AY-B^{\dagger}Y-Y^{\dagger}B-C=0,$ (9)

where $A$, $B$, and $C$ are $D\times D$ matrices with complex-valued entries.
This can be understood as a matrix version of the famous (scalar) quadratic
equation $ay^{2}-2by-c=0$. Solutions of equations like (9) are not always
guaranteed to exist, and certain conditions are required to prove the
existence of, for instance, Hermitian solutions [30]. See [31] for conditions
on solvability. Even if existence can be shown, the solutions may not be
unique or could alternatively be uncountably many [32, 33, 34]. However, there
are unique solutions under certain conditions. For instance, if all the matrix
entries are real-valued, then for symmetric positive semidefinite $A,C$ and
symmetric positive $Y$, there is a unique positive definite solution if and
only if an associated matrix $H=\begin{pmatrix}-B&A\\\ C&B^{T}\end{pmatrix}$
has no imaginary eigenvalues [35].

In this paper, we confine our attention to simpler cases, for example in
Lemmas 4 and 5, when there are unique solutions.

######  Lemma 4 (Solution of simple algebraic Riccati equation).

Consider the following algebraic Riccati equation when $A$ and $C$ are
positive definite matrices and $Y$ is Hermitian:

$\displaystyle YAY=C.$ (10)

This equation has a unique positive definite solution given by the standard
matrix geometric mean:

$\displaystyle Y=A^{-1}\\#C\coloneqq
A^{-1/2}(A^{1/2}CA^{1/2})^{1/2}A^{-1/2}>0.$ (11)

###### Proof.

This lemma is well known from [36, 37], but we provide a brief proof for
completeness. Starting from the Riccati equation in (10) and by using the fact
that $A$ is positive definite with a unique square root, consider that

$\displaystyle YAY=C\quad$ $\displaystyle\Leftrightarrow\quad
YA^{1/2}A^{1/2}Y=C$ (12) $\displaystyle\Leftrightarrow\quad
A^{1/2}YA^{1/2}A^{1/2}YA^{1/2}=A^{1/2}CA^{1/2}$ (13)
$\displaystyle\Leftrightarrow\quad(A^{1/2}YA^{1/2})^{2}=A^{1/2}CA^{1/2}.$ (14)

Since the matrix $A^{1/2}CA^{1/2}$ is positive definite and the equality in
the last line above has been shown, both $A^{1/2}CA^{1/2}$ and
$(A^{1/2}YA^{1/2})^{2}$ have a unique positive definite square root, implying
that

$\displaystyle A^{1/2}YA^{1/2}=(A^{1/2}CA^{1/2})^{1/2}\quad$
$\displaystyle\Leftrightarrow\quad Y=A^{-1/2}(A^{1/2}CA^{1/2})^{1/2}A^{-1/2},$
(15)

thus justifying that $Y=A^{-1}\\#C$ is the unique positive definite solution
as claimed. ∎

See [38] for a discussion of (10) in the infinite-dimensional case.

If $A$ and $C$ are both positive definite with unit trace
$\operatorname{Tr}(A)=1=\operatorname{Tr}(C)$, then $A$ and $C$ can also be
interpreted as density matrices. Then the operator $A^{-1}\\#C$ is also known
as the Fuchs–Caves observable [39], which is of relevance in the study of
quantum fidelity. We will return to this point later. See also [40, Section V]
for an interpretation of (10) when $A$ and $C$ are density matrices.

We can also extend Lemma 4 to the $B\neq 0$ case, and the following holds.

######  Lemma 5.

If $A$ and $C$ are positive definite, $B$ is an arbitrary matrix, and
$(A^{-1}B)=(A^{-1}B)^{\dagger}$, then a Hermitian solution to Eq. (9) can be
expressed as

$\displaystyle Y=A^{-1}\\#(B^{\dagger}A^{-1}B+C)+A^{-1}B.$ (16)

###### Proof.

See Appendix A. ∎

Classical algorithms for solving algebraic Riccati equations are typically
inefficient [18] with respect to the size of the problem, i.e., polynomial in
$D$. We will be looking at conditions for which a quantum algorithm for
solving algebraic Riccati equations can be executed with less complexity.

###  II.3 Block-encoding

Classical information can be embedded in quantum systems in the form of
quantum states, either pure or mixed, or in the form of quantum processes. A
closed quantum system evolves under a unitary transformation, represented by a
unitary matrix. In this paper, we will be focusing on how a matrix solution to
a matrix equation can be embedded in a unitary matrix. Unlike other quantum
subroutines that prepare solutions of a linear system of equations embedded in
the amplitudes of a pure quantum state, here we first embed the solution $Y$
into a unitary matrix.

There are different ways of embedding an arbitrary matrix into a unitary
matrix. For instance, it is guaranteed by the Sz. Nagy dilation theorem (see,
e.g., [41, Theorem 1.1]) that such a unitary matrix should always exist. We
choose a flexible dilation known as block-encoding [19, 20, 21]. A unitary
matrix $U_{Y}$ is called a block-encoding of a matrix $Y$ if it satisfies the
following definition.

######  Definition 6 (Block-encoding).

Fix $n,a\in\mathbb{N}$ and $\epsilon,\alpha\geq 0$. Let $Y$ be an $n$-qubit
operator. An $(n+a)$-qubit unitary $U_{Y}$ is an $(\alpha,a,\epsilon)$-block-
encoding of an operator $Y$ if

$\left\|Y-\alpha\bra{0}_{a}U_{Y}\ket{0}_{a}\right\|\leq\epsilon.$ (17)

Here $|0\rangle_{a}$ are the $|0\rangle$ states in the computational basis of
the $a$-ancilla qubits. The block-encoding formalism allows one to construct,
for example, block-encodings of sums of matrices, linear combinations of
block-encoded matrices, and polynomial approximations of negative and positive
power functions of matrices [19]. We list several associated lemmas in
Appendix B for convenience.

##  III Quantum subroutines for matrix geometric means, algebraic Riccati
equations, and higher-order nonlinear equations

Let us focus on cases where the solutions to the algebraic Riccati equations
can be captured by the matrix geometric mean in Lemmas 4 and 5. The
computation of the matrix geometric mean involves the computation of the
square roots of matrices and several matrix multiplications. For $D\times D$
matrices, typically these costs will scale polynomially with $D$ for a
classical algorithm. However, quantum algorithms for matrix multiplications of
block-encoded matrices can be performed more efficiently when compared to the
number of classical numerical steps. These series of block-encoded matrix
multiplications can be achieved in the quantum case via the block-encoding
formalism.

Let us begin with the algebraic Riccati equation in Eq. (9):

$\displaystyle YAY-B^{\dagger}Y-Y^{\dagger}B-C=0.$ (18)

It is our goal below first to construct a block-encoding of the solution $Y$,
denoted $U_{Y}$, under the conditions obeyed in Lemmas 4 and 5. This we
consider as a subroutine that we can then employ in various applications.

Below we assume that we also have access to the block-encodings of $A$, $B$,
$C$ – denoted $U_{A}$, $U_{B}$, $U_{C}$, respectively – as well as their
inverses $U_{A}^{\dagger}$, $U_{B}^{\dagger}$, $U_{C}^{\dagger}$ and
controlled versions. For example, if $A$, $B$, and $C$ are positive semi-
definite matrices with unit trace, these can be considered as density
matrices. Then from Lemma 30 we can prepare block-encodings $U_{A}$, $U_{B}$,
$U_{C}$ by access to the unitaries that prepare purifications of $A$, $B$, and
$C$, with only a single query to each purification and $O(\log D)$ gates. In
more general scenarios, we can leave the preparation of these block-encodings
to a later stage, which also depends on the particular application. Below,
$\kappa_{A}$ and $\kappa_{C}$ denote the condition numbers for $A$ and $C$,
respectively. It is important to clarify that, as assumed in [1], all of our
quantum algorithms for matrix geometric means assume that

$\displaystyle I$ $\displaystyle\geq A\geq I/\kappa_{A},$ (19) $\displaystyle
I$ $\displaystyle\geq C\geq I/\kappa_{C},$ (20) $\displaystyle I$
$\displaystyle\geq B^{\dagger}B.$ (21)

This means that $\kappa_{A}$ and $\kappa_{C}$ are really equal to the inverses
of the minimum eigenvalues of $A$ and $C$, respectively, and
$\left\|A\right\|,\left\|B\right\|,\left\|C\right\|\leq 1$. The upper bounds
above are automatically satisfied whenever $A$ and $C$ are density matrices.

###  III.1 Quantum subroutine for matrix geometric means

As a warm-up, we present a quantum subroutine for implementing block-encodings
of the weighted matrix geometric means.

######  Lemma 7 (Block-encoding of weighted matrix geometric mean).

Suppose that $U_{A},U_{C}$ are $(1,a,0)$-block-encodings of matrices $A,C$,
respectively, where $A\geq I/\kappa_{A}$, $C\geq I/\kappa_{C}$ and $I$ is the
identity matrix. For $\epsilon\in(0,1/2)$, one can implement a
$(\kappa_{A}^{1/p}\gamma_{p},5a+20,\epsilon)$-block-encoding of $Y$ for every
fixed real $p\neq 0$, where

$\gamma_{p}=\begin{cases}1&p>0,\\\
\kappa_{A}^{-1/p}\kappa_{C}^{-1/p}&p<0,\end{cases}$ (22)

and

$Y=A\\#_{1/p}C=A^{1/2}\left(A^{-1/2}CA^{-1/2}\right)^{1/p}A^{1/2},$ (23)

using

  * •

$\widetilde{O}\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$
queries to $U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$;

  * •

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates; and

  * •

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

######  Remark 1.

In the above and in what follows, ‘queries to $U$’ refers to access not only
to $U$, but also to its inverse $U^{\dagger}$, controlled-$U$, and
controlled-$U^{\dagger}$. Here and in the following, $\widetilde{O}(\cdot)$
suppresses logarithmic factors of functions appearing in $(\cdot)$.

###### Proof sketch of Lemma 7.

See Appendix C for a detailed proof. As an illustration for the construction
of our quantum subroutines, we outline the basic idea. Other quantum
subroutines later presented in this section are obtained using similar ideas.
Our approach consists of three main steps:

  1. 1.

Implement a block-encoding of $A^{-1/2}$, using roughly
$\widetilde{O}\lparen\kappa_{A}\rparen$ queries to a block-encoding of $A$
(for simplicity, we ignore the $\epsilon$-dependence in our brief explanation
here). This is done by applying quantum singular value transformation (QSVT)
[19] with polynomial approximations of negative power functions (see Lemma
27).

  2. 2.

Implement a block-encoding of $\lparen A^{-1/2}CA^{-1/2}\rparen^{1/p}$, using
roughly $\widetilde{O}\lparen\kappa_{A}\kappa_{C}\rparen$ queries to a block-
encoding of $A^{-1/2}CA^{-1/2}$. This is done by applying QSVT with polynomial
approximations of positive power functions (see Lemma 28). Note that a block-
encoding of $A^{-1/2}CA^{-1/2}$ can be implemented using $O\lparen 1\rparen$
queries to block-encodings of $A^{-1/2}$ and $C$ by the method for realising
the product of block-encoded matrices (see Lemma 24).

  3. 3.

Similar to Step 2, implement a block-encoding of $A^{1/2}\lparen
A^{-1/2}CA^{-1/2}\rparen^{1/p}A^{1/2}$, using $O\lparen 1\rparen$ queries to
block-encodings of $A^{1/2}$ and $\lparen A^{-1/2}CA^{-1/2}\rparen^{1/p}$,
where a block-encoding of $A^{1/2}$ can be implemented using
$\widetilde{O}\lparen\kappa_{A}\rparen$ queries to a block-encoding of $A$.

To conclude, the overall query complexity is roughly
$\widetilde{O}\lparen\kappa_{A}\rparen\cdot\widetilde{O}\lparen\kappa_{A}\kappa_{C}\rparen+\widetilde{O}\lparen\kappa_{A}\rparen=\widetilde{O}\lparen\kappa_{A}^{2}\kappa_{C}\rparen$.
Note that the construction is mainly based on QSVT and thus is also time
efficient. So the overall time complexity is equal to the query complexity
only up to polylogarithmic factors. ∎

###  III.2 $B=0$ algebraic Riccati equation

Let us begin with the unique positive definite solution to the algebraic
Riccati equation with $B=0$, i.e., Eq. (9), which can be expressed as the
matrix geometric mean $Y=A^{-1}\\#C$, according to Lemma 4, where $A$ and $C$
are positive definite matrices. Then we have the following lemma, which
characterises a block-encoding of the solution in a quantum circuit.

######  Lemma 8.

Suppose that $U_{A},U_{C}$ are $(1,a,0)$-block-encodings of matrices $A,C$,
respectively, with $A\geq I/\kappa_{A}$ and $C\geq I/\kappa_{C}$. For
$\epsilon\in(0,1/2)$, one can implement a $(\kappa_{A},5a+20,\epsilon)$-block-
encoding of $Y$, where

$Y=A^{-1}\\#C=A^{-1/2}\left(A^{1/2}CA^{1/2}\right)^{1/2}A^{-1/2},$ (24)

using

  * •

$\widetilde{O}\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$
queries to $U_{C}$ and
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$;

  * •

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates; and

  * •

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

###### Proof.

See Appendix D. ∎

###  III.3 $B\neq 0$ algebraic Riccati equation

Here we want to construct a block-encoding of a Hermitian solution to the
algebraic Riccati equation via the standard matrix geometric mean, according
to Lemma 5. We then have the following lemma.

######  Lemma 9.

Suppose that $U_{A},U_{B},U_{C}$ are $(1,a,0)$-block-encodings of matrices
$A,B,C$, respectively, with $A\geq I/\kappa_{A}$, $C\geq I/\kappa_{C}$ and
$A^{-1}B=\left(A^{-1}B\right)^{\dagger}$. For $\epsilon\in(0,1/2)$, one can
implement a $(\kappa_{A}^{3/2},b,\epsilon)$-block-encoding of $Y$, where
$b=O\\!\left(a+\log\left(1/\epsilon\right)\right)$ and

$\displaystyle Y$
$\displaystyle=A^{-1}\\#\left(B^{\dagger}A^{-1}B+C\right)+A^{-1}B$ (25)
$\displaystyle=A^{-1/2}\left(A^{1/2}\left(B^{\dagger}A^{-1}B+C\right)A^{1/2}\right)^{1/2}A^{-1/2}+A^{-1}B,$
(26)

using

  * •

$\widetilde{O}\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$
queries to $U_{B}$ and $U_{C}$, and
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$;

  * •

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates; and

  * •

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

###### Proof.

See Appendix E. ∎

###  III.4 Higher-order polynomial equations

Algebraic Riccati equations are second-order nonlinear equations whose
solutions are given by the second-order matrix geometric mean, i.e., $p=2$. We
can also generalise our formalism to particular $p^{\text{th}}$-order
nonlinear matrix equations, whose solutions involve $p\in\\{3,4,\ldots\\}$
weighted matrix geometric mean. For example, consider the following
$p^{\text{th}}$-order nonlinear matrix equations, which we call the
$p^{\text{th}}$-order YAY algebraic equations

$\displaystyle Y(AY)^{p-1}=C,$ (27)

where $p$ is the highest order polynomial in $Y$. It is straightforwardly
checked that the solutions can be written in terms of the weighted geometric
mean from Definition 2:

$\displaystyle Y=A^{-1/2}(A^{1/2}CA^{1/2})^{1/p}A^{-1/2}=A^{-1}\\#_{1/p}C.$
(28)

See [42] for a discussion of this kind of equation in the infinite-dimensional
case.

######  Lemma 10 (Solution of simple $p^{\text{th}}$-order algebraic $YAY$
equation).

Fix $p\in\\{2,3,4,\ldots\\}$. Consider the $p^{\text{th}}$-order algebraic
$YAY$ equation when $A$ and $C$ are positive definite matrices:

$\displaystyle Y(AY)^{p-1}=C.$ (29)

This equation has a unique positive definite solution given by the following
weighted geometric mean:

$\displaystyle Y=A^{-1/2}(A^{1/2}CA^{1/2})^{1/p}A^{-1/2}=A^{-1}\\#_{1/p}C>0.$
(30)

###### Proof.

The proof is a generalisation of that for Lemma 4 and we provide it for
completeness. Starting from the equation in (29) and by using the fact that
$A$ is positive definite with a unique square root, consider that

$\displaystyle Y(AY)^{p-1}=C\quad$ $\displaystyle\Leftrightarrow\quad
Y(A^{1/2}A^{1/2}Y)^{p-1}=C$ (31)
$\displaystyle\Leftrightarrow\quad(A^{1/2}YA^{1/2})^{p}=A^{1/2}CA^{1/2},$ (32)

where the last line is obtained by left and right multiplying the previous
line by $A^{1/2}$. Since the matrix $A^{1/2}CA^{1/2}$ is positive definite and
the equality in the last line above has been shown, both $A^{1/2}CA^{1/2}$ and
$(A^{1/2}YA^{1/2})^{p}$ have a unique positive definite $p$th root, implying
that

$\displaystyle A^{1/2}YA^{1/2}=(A^{1/2}CA^{1/2})^{1/p}\quad$
$\displaystyle\Leftrightarrow\quad Y=A^{-1/2}(A^{1/2}CA^{1/2})^{1/p}A^{-1/2},$
(33)

thus justifying that $Y=A^{-1}\\#_{1/p}C$ is the unique positive definite
solution as claimed. ∎

To construct a block-encoding for the weighted geometric mean (and thus the
solution of (27)), we have proven the following lemma, which holds for every
non-zero real number $p$.

######  Lemma 11.

Suppose that $U_{A},U_{C}$ are $(1,a,0)$-block-encodings of matrices $A,C$,
respectively, with $A\geq I/\kappa_{A}$ and $C\geq I/\kappa_{C}$ and let
$p\neq 0$ be any fixed non-zero real number. For $\epsilon\in(0,1/2)$, one can
implement a $(\kappa_{A}\gamma_{p},5a+20,\epsilon)$-block-encoding of $Y$,
where

$Y=A^{-1}\\#_{1/p}C=A^{-1/2}\left(A^{1/2}CA^{1/2}\right)^{1/p}A^{-1/2},$ (34)

and

$\gamma_{p}=\begin{cases}1&p>0,\\\
\kappa_{A}^{-1/p}\kappa_{C}^{-1/p}&p<0,\end{cases}$ (35)

using

  * •

$\widetilde{O}\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$
queries to $U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$;

  * •

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates; and

  * •

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

###### Proof.

See Appendix F. ∎

##  IV Applications

Here we explore two classes of applications for preparing block-encodings of
the matrix geometric mean. The first class of applications is to learning
problems, in particular for metric learning from data, both quantum and
classical. Next we demonstrate how having access to the matrix geometric mean
also allows us to compute some fundamental quantities in quantum information,
like the quantum fidelity between two mixed states via the Fuchs–Caves
observable, as well as geometric Rényi relative entropies.

###  IV.1 Quantum geometric mean metric learning

In learning problems, there is typically a loss function $L$ that we want to
optimise. Suppose we have $D\times D$ positive definite matrices $Y$, $A$, and
$C$. We note that here the uniqueness result in Lemma 4 continues to hold.
Consider the following optimisation problem:

$\displaystyle\min_{Y\geq 0}L(Y),\qquad
L(Y)\coloneqq\operatorname{Tr}(YA)+\operatorname{Tr}(Y^{-1}C).$ (36)

It turns out that, for given $A$ and $C$, the unique $Y$ minimising $L(Y)$ is
$Y=A^{-1}\\#C$. In [25], this result was proven for real positive definite
matrices $A$ and $C$, and here we extend it to positive definite Hermitian
matrices. In [43], the same optimisation problem was considered in the context
of quantum fidelity, where it was shown that the optimal value of (36) is
equal to $\operatorname{Tr}[(A^{1/2}CA^{1/2})^{1/2}]$.

######  Lemma 12.

Fix $A$ and $C$ to be positive definite matrices. The unique solution to
$\min_{Y\geq 0}L(Y)$ where
$L(Y)=\operatorname{Tr}(YA)+\operatorname{Tr}(Y^{-1}C)$ is the matrix
geometric mean $Y=A^{-1}\\#C$.

###### Proof.

If $L(Y)$ is both strictly convex and strictly geodesically convex, then the
solution to $\nabla L(Y)=0$ will also be a global minimiser. For the proof of
strict convexity and strict geodesic convexity, see Appendix G. Now $\nabla
L(Y)=A-Y^{-1}CY^{-1}=0$ implies the algebraic Riccati equation $YAY=C$ or
$Y=A^{-1}\\#C$, which is the unique solution for positive definite matrices
$A$ and $C$. ∎

We will use this property and map two learning problems – one for classical
data and another for quantum data – onto this optimisation problem. Using the
block-encoding for the matrix geometric mean in Lemma 8, we then devise
quantum algorithms for learning a Euclidean metric from data as well as a
$1$-class classification algorithm for quantum states. We also extend to the
case of weighted learning, where there are unequal contributions to the loss
function in Eq. (36) from $\text{Tr}(YA)$ and $\text{Tr}(Y^{-1}C)$.

####  IV.1.1 Learning Euclidean metric from data

Machine learning algorithms rely on distance measures to quantify how similar
one set of data is to another. Naturally different distance measures can give
rise to different results, and so choosing the right metric is crucial for the
success of an algorithm. The distance measure itself can in fact be learned
for example in a weakly supervised scenario, and this is called metric
learning [44]. Here we are provided with the following two sets $\mathcal{S}$
(similar) and $\mathcal{D}$ (dissimilar) of pairs (training data)

$\displaystyle\mathcal{S}$
$\displaystyle\coloneqq\\{(\boldsymbol{x},\boldsymbol{x}^{\prime})\,|\,\boldsymbol{x},\boldsymbol{x}^{\prime}\text{
are in the same class}\\},$ (37) $\displaystyle\mathcal{D}$
$\displaystyle\coloneqq\\{(\boldsymbol{x},\boldsymbol{x}^{\prime})\,|\,\boldsymbol{x},\boldsymbol{x}^{\prime}\text{
are in different classes}\\},$ (38)

and $\\{(\boldsymbol{x}^{(k)},\boldsymbol{x}^{\prime(k)})\\}_{k}$ are the data
points, where $k$ labels all the pairs that either belong to $\mathcal{S}$ or
$\mathcal{D}$. An important example in metric learning is learning the
Euclidean metric from data, which can be reformulated as a simple optimisation
problem with a closed-form solution. Learning a Euclidean metric is a common
form of metric learning, where we can learn a Mahalanobis distance $d_{Y}$

$\displaystyle
d_{Y}(\boldsymbol{x},\boldsymbol{x}^{\prime})=(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}Y(\boldsymbol{x}-\boldsymbol{x}^{\prime})=\operatorname{Tr}(Y(\boldsymbol{x}-\boldsymbol{x}^{\prime})(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}),$
(39)

with $Y$ a real $D\times D$ symmetric positive definite matrix. To identify a
suitable $Y$, one requires a suitable cost function.

In geometric mean metric learning [25], we want $d_{Y}$ to be minimal between
data in the same class, i.e., $\mathcal{S}$. At the same time, when the data
are in different classes, i.e., $\mathcal{D}$, we want $d_{Y^{-1}}$ to be
minimal instead. Thus we want to minimise the sum
$\sum_{\mathcal{S}}d_{Y}+\sum_{\mathcal{D}}d_{Y^{-1}}$. This leads to an
optimisation problem of the form in Eq. (36)

$\displaystyle\min_{Y\geq 0}L(Y),\qquad
L(Y)\coloneqq\operatorname{Tr}(YA)+\operatorname{Tr}(Y^{-1}C),$ (40)
$\displaystyle
A=\sum_{(\boldsymbol{x},\boldsymbol{x}^{\prime})\in\mathcal{S}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T},$
(41) $\displaystyle
C=\sum_{(\boldsymbol{x},\boldsymbol{x}^{\prime})\in\mathcal{D}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T},$
(42)

where we assume $A$ and $C$ are positive definite. From Lemma 12, we see that
the optimal solution to Eq. (40) is the matrix geometric mean $Y=A^{-1}\\#C$.
We see that this is also in fact the solution of the $B=0$ algebraic Riccati
equation $YAY=C$. In Lemma 8, we saw that, given access to the block-encodings
of $U_{A}$, $U_{C}$, their inverses, and controlled versions, we can construct
the block-encoding of $Y$, denoted $U_{Y}$. Lemma 8 also shows that the query
and gate complexity costs are efficient in $D$, i.e., $O(\text{poly}\log D)$,
when the condition numbers for $A$ and $C$ are also polynomial in $\log D$.

While we see that it is possible to efficiently recover $U_{Y}$, it is not
sufficient for a direct application to machine learning. Learning $Y$ is part
of the learning stage, but reading off the classical components of $Y$
directly from $U_{Y}$ is inefficient. However, if we consider the testing
stage in machine learning, then we need to compute the actual distance $d_{Y}$
if we are given a new data pair $(\boldsymbol{y},\boldsymbol{y}^{\prime})$,
known as testing data. For this testing data, we do not know its
classification into $\mathcal{S}$ or $\mathcal{D}$ a priori. Thus the task is
to show that, having access to $U_{Y}$, it is then sufficient to compute
$d_{Y}$ without needing to read out the elements of $Y$. For example, a large
value of $d_{Y}(\boldsymbol{y},\boldsymbol{y}^{\prime})$ means that we should
classify $(\boldsymbol{y},\boldsymbol{y}^{\prime})\in\mathcal{D}$, whereas a
small value of $d_{Y}(\boldsymbol{y},\boldsymbol{y}^{\prime})$ means that we
should classify $(\boldsymbol{y},\boldsymbol{y}^{\prime})\in\mathcal{S}$. We
discuss later the preparation of the block-encodings of $U_{A}$, $U_{C}$.

Before proceeding, we first discuss the preparation of a quantum state that we
later require. Given a new data pair (testing data)
$(\boldsymbol{y},\boldsymbol{y}^{\prime})$ for which we want to compute
$d_{Y}$, where we use the optimal $Y$, we can define a corresponding pure
quantum state with $m=O(\log D)$ qubits
$|\psi\rangle_{y,y^{\prime}}=(1/\mathcal{N}_{\psi})\sum_{i=1}^{D}(\boldsymbol{y}-\boldsymbol{y}^{\prime})_{i}|i\rangle$,
with normalisation constant
$\mathcal{N}_{\psi}^{2}=\sum_{i=1}^{D}(\boldsymbol{y}-\boldsymbol{y}^{\prime})_{i}^{2}$.
Its amplitudes are proportional to $\boldsymbol{y}-\boldsymbol{y}^{\prime}$
for any pair $(\boldsymbol{y},\boldsymbol{y}^{\prime})$. We say that the state
$|\psi\rangle_{y,y^{\prime}}$ has sparsity $\sigma$ if $\sigma$ is the number
of non-zero entries in the amplitude. We can use optimal state preparation
schemes [45, 46, 47] to prepare $|\psi\rangle_{y,y^{\prime}}$.

######  Lemma 13.

A circuit [45] producing an $m$-qubit state
$|z\rangle=\sum_{i=1}^{2^{m}}z_{i}|i\rangle$ from $|0\rangle$ with given
classical entries $\\{z_{i}\\}_{i=1}^{2^{m}}$ can be implemented using
$O(m\sigma)$ CNOT gates and $O(\sigma(\log\sigma+m))$ one-qubit gates, where
the specification of the circuit can be found with a classical algorithm with
time complexity $O(m\sigma^{2}\log\sigma)$. It also uses a constant number of
ancilla qubits. Alternative preparation strategies [46, Remark 1],[47, Theorem
2] have a reduced gate (depth) complexity $\Theta(\log(m\sigma))$, but require
$O(m\sigma\log\sigma)$ ancilla qubits.

Since here $m=O(\log D)$, we see that as long as $\sigma$ is small, e.g.,
$\sigma=O(\operatorname{poly}\log D)$, then the total initial state
preparation cost, in either gate complexity and number of ancillas is
$O(\operatorname{poly}\log D)$. Next we compute $d_{Y}$.

######  Theorem 14.

Suppose we are given $U_{A}$, $U_{C}$, which are $(1,\log D,0)$-block-
encodings
$A=\sum_{(\boldsymbol{x},\boldsymbol{x}^{\prime})\in\mathcal{S}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}$
and
$C=\sum_{(\boldsymbol{x},\boldsymbol{x}^{\prime})\in\mathcal{D}}(\boldsymbol{x}-\boldsymbol{x}^{\prime})(\boldsymbol{x}-\boldsymbol{x}^{\prime})^{T}$,
respectively. We also assume access to their inverses and controlled versions.
We assume that the data obeys $\kappa_{A},\kappa_{C}=O(\operatorname{poly}\log
D)$. Given a testing data pair $(\boldsymbol{y},\boldsymbol{y}^{\prime})$, we
assume that the corresponding state $|\psi\rangle_{y,y^{\prime}}$ has sparsity
$O(\operatorname{poly}\log D)$ and
$\mathcal{N}_{\psi}=O(\operatorname{poly}\log D)$. Then computing
$d_{Y}(\boldsymbol{y},\boldsymbol{y}^{\prime})$ to precision $\epsilon$ has a
query and gate complexity $O(\operatorname{poly}(\log D,1/\epsilon))$.

###### Proof.

We first make the observation that, for optimal $Y$,

$\displaystyle
d_{Y}(\boldsymbol{y},\boldsymbol{y^{\prime}})=\mathcal{N}_{\psi}^{2}\langle\psi|_{y,y^{\prime}}Y|\psi\rangle_{y,y^{\prime}}\approx\mathcal{N}_{\psi}^{2}\kappa_{A}\operatorname{Tr}(\langle
0|U_{Y}|0\rangle|\psi\rangle_{y,y^{\prime}}\langle\psi|_{y,y^{\prime}}),$ (43)

where $Y=A^{-1}\\#C$ and $U_{Y}$ is a $(\kappa_{A},5a+20,\epsilon)$ block-
encoding of $Y$. The proportionality constant of $\kappa_{A}$ comes from Lemma
8, which shows that $\|Y-\kappa_{A}\langle
0|_{5a+20}U_{Y}|0\rangle_{5a+20}\|\leq\epsilon$, where $a=\log D$. For
simplicity, we neglect the subscript on the $|0\rangle$ states. This trace can
be interpreted as an expectation value of $|\psi\rangle_{y,y^{\prime}}$ with
$Y$ as the observable, and comes from the definition of $d_{Y}$ and
$|\psi\rangle_{y,y^{\prime}}$.

To compute this trace given $U_{Y}$, we observe $\operatorname{Tr}((S\otimes
T)X)=\operatorname{Tr}(X_{T}S)$, where if
$T=\sum_{n}\lambda_{n}|u_{n}\rangle\\!\langle v_{n}|$, then
$X_{T}=\sum_{n}\lambda_{n}\langle v_{n}|X|u_{n}\rangle$. So we can rewrite

$\displaystyle d_{Y}(\boldsymbol{y},\boldsymbol{y}^{\prime})$
$\displaystyle\approx\mathcal{N}^{2}_{\psi}\kappa_{A}\operatorname{Tr}(\langle
0|U_{Y}|0\rangle|\psi\rangle_{y,y^{\prime}}\langle\psi|_{y,y^{\prime}})$ (44)
$\displaystyle=\mathcal{N}^{2}_{\psi}\kappa_{A}\operatorname{Tr}((|\psi\rangle_{y,y^{\prime}}\langle\psi|_{y,y^{\prime}}\otimes|0\rangle\\!\langle
0|)U_{Y})$ (45)
$\displaystyle=\mathcal{N}^{2}_{\psi}\kappa_{A}\langle\Psi|_{y,y^{\prime}}U_{Y}|\Psi\rangle_{y,y^{\prime}},$
(46)

where
$|\Psi\rangle_{y,y^{\prime}}=|\psi\rangle_{y,y^{\prime}}\otimes|0\rangle$. The
last expectation value can be realised with a conventional swap test [48, 49]
between the states $|\Psi\rangle_{y,y^{\prime}}$ and
$U_{Y}|\Psi\rangle_{y,y^{\prime}}$. One can also use the destructive SWAP test
(i.e., Bell measurements and classical post-processing) [50]. Alternatively,
$\operatorname{Tr}((|\psi\rangle_{y,y^{\prime}}\langle\psi|_{y,y^{\prime}}\otimes|0\rangle\\!\langle
0|)U_{Y})$ can also be computed through a Hadamard test (Lemma 31), where one
is given the controlled-$U_{Y}$ and state
$|\psi\rangle_{y,y^{\prime}}\otimes|0\rangle$. For example, applying the
unitary $U_{Y}$ to $|\psi\rangle_{y,y^{\prime}}\otimes|0\rangle$ and using the
swap test with $|\psi\rangle_{y,y^{\prime}}\otimes|0\rangle$, we recover
$\langle\psi|_{x,x^{\prime}}Y|\psi\rangle_{x,x^{\prime}}$ to precision
$\epsilon$ with query and gate complexity $O(\operatorname{poly}(\log
D,1/\epsilon))$, when $\kappa_{A}=O(\operatorname{poly}\log D)$. Now,
$d_{Y}(\boldsymbol{y},\boldsymbol{y}^{\prime})=\mathcal{N}^{2}_{\psi}\langle\psi|_{y,y^{\prime}}Y|\psi\rangle_{y,y^{\prime}}$.
Since we only have $\sigma=O(\text{poly}\log D)$ non-zero entries in
$\boldsymbol{y}-\boldsymbol{y}^{\prime}$, the cost in the classical
computation of the normalisation constant is also of order
$O(\operatorname{poly}\log D)$. Assuming that
$\mathcal{N}^{2}_{\psi}=O(\operatorname{poly}\log D)$, then we recover $d_{Y}$
efficiently when given access to $U_{Y}$ and $|\psi\rangle_{y,y^{\prime}}$.

We saw that preparing the state $|\psi\rangle_{y,y^{\prime}}$ according to
Lemma 13 incurs a cost $O(\operatorname{poly}\log D)$. From Lemma 8 we can
construct a $(\kappa_{A},O(\log D),\epsilon)$-block-encoding of $Y$ with gate
and query complexity
$O(\operatorname{poly}(\kappa_{A},\kappa_{C},\log(1/\epsilon)))$. Since
$\kappa_{A},\kappa_{C}=O(\text{poly}\log D)$, then the theorem is proved. ∎

Thus, if our assumptions are obeyed, the quantum cost for computation of
$d_{Y}$ can be $O(\text{poly}\log D)$, whereas classical numerical algorithms
for computing the matrix geometric mean alone has cost
$O(\operatorname{poly}D)$ for $D\times D$ matrices [51, 52].

In Theorem 14, we also assumed access to $U_{A}$, $U_{C}$. We show below the
preparation of a block-encoding of density matrices which are proportional to
$A$ and $C$ and how this can be used to compute $d_{Y}$ efficiently. First
consider Lemma 30, which shows how to create a block-encoding of a density
matrix. We first observe that we can define density matrices $\rho_{A}$ and
$\rho_{C}$ where rewrite

$\displaystyle\rho_{A}=\frac{A}{\text{Tr}(A)},\qquad
A=\sum_{k\in\mathcal{S}}\mathcal{N}^{2}_{\psi_{k}}|\psi_{k}\rangle\\!\langle\psi_{k}|,\qquad\text{Tr}(A)=\sum_{k\in\mathcal{S}}\mathcal{N}^{2}_{\psi_{k}},$
$\displaystyle\rho_{C}=\frac{C}{\text{Tr}(C)},\qquad
C=\sum_{k\in\mathcal{D}}\mathcal{N}^{2}_{\psi_{k}}|\psi_{k}\rangle\\!\langle\psi_{k}|,\qquad\text{Tr}(C)=\sum_{k\in\mathcal{D}}\mathcal{N}^{2}_{\psi_{k}},$
(47)

where
$|\psi_{k}\rangle=(1/\mathcal{N}_{\psi_{k}})\sum_{i=1}^{D}(\boldsymbol{x}^{(k)}-\boldsymbol{x}^{\prime(k)})_{i}|i\rangle$
and
$\mathcal{N}^{2}_{\psi_{k}}=\sum_{i=1}^{D}(\boldsymbol{x}^{(k)}-\boldsymbol{x}^{\prime(k)})_{i}^{2}$
is the corresponding normalisation. Then from Lemma 30, if we are given
unitaries $V_{A}$ and $V_{C}$ that prepare purifications of $\rho_{A}$ and
$\rho_{C}$, respectively, it is possible to create $U_{\rho_{A}}$ and
$U_{\rho_{C}}$ using one query to $V_{A}$ and $V_{C}$ respectively and $O(\log
D)$ gates. One such class of states purifying $\rho_{A}$ and $\rho_{C}$ are

$\displaystyle|\Sigma_{A}\rangle$
$\displaystyle\coloneqq\sum_{k\in\mathcal{S}}\sqrt{p^{(A)}_{k}}|k\rangle|\psi_{k}\rangle,$
(48) $\displaystyle|\Sigma_{C}\rangle$
$\displaystyle\coloneqq\sum_{k\in\mathcal{D}}\sqrt{p^{(C)}_{k}}|k\rangle|\psi_{k}\rangle,$
(49)

where

$\displaystyle p_{k}^{(A)}$
$\displaystyle\coloneqq\mathcal{N}^{2}_{\psi_{k}}/\sum_{l\in\mathcal{S}}\mathcal{N}_{\psi_{l}}^{2},$
(50) $\displaystyle p_{k}^{(C)}$
$\displaystyle\coloneqq\mathcal{N}^{2}_{\psi_{k}}/\sum_{l\in\mathcal{D}}\mathcal{N}_{\psi_{l}}^{2}.$
(51)

To prepare $|\Sigma_{A}\rangle$ and $|\Sigma_{C}\rangle$ we require the
controlled unitaries $V_{A}=\sum_{k\in\mathcal{S}}|k\rangle\\!\langle
k|\otimes W_{k}^{(A)}$ and $V_{C}=\sum_{k\in\mathcal{D}}|k\rangle\\!\langle
k|\otimes W_{k}^{(C)}$ acting on states
$\sum_{k\in\mathcal{S}}\sqrt{p_{k}^{(A)}}|k\rangle|0\rangle$ and
$\sum_{k\in\mathcal{D}}\sqrt{p_{k}^{(A)}}|k\rangle|0\rangle$ respectively.
Here $W_{k}^{(A)}$ and $W_{k}^{(C)}$ are the state preparation circuits from
Lemma 13 that create $|\psi_{k}\rangle$ where $k\in\mathcal{S}$ and
$k\in\mathcal{D}$, respectively. Since $W_{k}^{(A)}$ and $W_{k}^{(C)}$ are
known circuits and assuming $\sigma=O(\operatorname{poly}\log D)$, it is
similarly efficient and also straightforward to realise $V_{A}$ and $V_{C}$.
Then from Lemma 30 it is possible to create $(1,O(\log D),0)$-block-encodings
of $\rho_{A}$ and $\rho_{C}$ with gate and query complexity
$O(\operatorname{poly}\log D)$, denoted $U_{\rho_{A}}$ and $U_{\rho_{C}}$,
respectively. In the case where $\text{Tr}(A)=1=\text{Tr}(C)$, then this
automatically gives us the unitaries $U_{A}$ and $U_{C}$ required in Theorem
14.

For general classical data, $\text{Tr}(A)=1=\text{Tr}(C)$ would not hold in
general. However, since $A=\text{Tr}(A)\rho_{A}$ and $C=\text{Tr}(C)\rho_{C}$,
the proof in Theorem 14 holds in the same way if we began with $U_{\rho_{A}}$
and $U_{\rho_{C}}$, from which we can create $U_{Y^{\prime}}$ where
$Y^{\prime}\equiv\rho_{A}^{-1}\\#\rho_{C}=(\text{Tr}(C)/\text{Tr}(A))^{1/2}Y$.
This implies

$\langle 0|U_{Y}|0\rangle\approx
Y=(\text{Tr}(A)/\text{Tr}(C))^{1/2}Y^{\prime}\approx(\text{Tr}(A)/\text{Tr}(C))^{1/2}\langle
0|U_{Y^{\prime}}|0\rangle.$ (52)

Following through the same proof idea as in Theorem 14 allows us to extract
$d_{Y}^{\prime}$. To recover $d_{Y}$, we just use
$d_{Y}\approx(\text{Tr}(A)/\text{Tr}(C))^{1/2}d_{Y}^{\prime}$. These
normalisations can be efficiently recovered by assuming the states
$|\psi_{k}\rangle$ have low sparsity $\sigma_{k}=O(\operatorname{poly}\log D)$
for each $k$. This also means that the normalisations $\text{Tr}(A)$ and
$\text{Tr}(C)$ are efficient to compute. So long as
$(\text{Tr}(A)/\text{Tr}(C))^{1/2}$ is $O(\operatorname{poly}\log D)$, then
$d_{Y}$ is efficiently estimable.

####  IV.1.2 1-class quantum learning

Here we propose a new quantum classification problem that is a 1-class
problem. This means that given a quantum state, we only want to know whether
this state belongs to a class $\mathcal{A}$ or not. This problem occurs in
many areas in machine learning, in particular in anomaly detection, where
$\mathcal{A}$ is the class of states that are considered anomalous. Here we
can be provided with the following training data:

$\displaystyle\rho=\frac{1}{N}\sum_{i=1}^{N}\rho_{i},\qquad\rho_{i}\in\mathcal{A},$
(53)
$\displaystyle\sigma=\frac{1}{M}\sum_{i=1}^{M}\sigma_{i},\qquad\sigma_{i}\notin\mathcal{A},$
(54)

where $\\{\rho_{i}\\}_{i}$ and $\\{\sigma_{i}\\}_{i}$ are sets of
$D$-dimensional states. In anomaly detection scenarios, there are usually much
fewer examples of anomalous states than ‘normal’ states, so that $N\ll M$.
However, we will not focus on subtleties associated with imbalanced training
data here.

Suppose that we have an incoming quantum state $\xi$ and we want to flag this
as belonging to the class $\mathcal{A}$ or not. Then it is useful to learn an
‘observable’ or a ‘witness’ $Y$ such that its expectation value
$\operatorname{Tr}(Y\xi)$ is large when $\xi$ is flagged as anomalous,
belonging to $\mathcal{A}$, but this value is small when $\xi$ is ‘normal’.
Thus we can set up an optimisation problem of the form

$\displaystyle\min_{Y\geq 0}L(Y),\qquad
L(Y)\coloneqq\operatorname{Tr}(Y\sigma)+\operatorname{Tr}(Y^{-1}\rho).$ (55)

It is sensible in the above to minimise $\operatorname{Tr}(Y^{-1}\rho)$ in
$L(Y)$ above since it is simple to show that a small value of
$\operatorname{Tr}(Y^{-1}\rho)$ implies a large value of
$\operatorname{Tr}(Y\rho)$. Since $\rho$ is a density matrix, it can be shown
that $\operatorname{Tr}(Y^{-1}\rho)\geq\operatorname{Tr}(Y\rho)^{-1}$, which
is a consequence of the operator Jensen inequality (see [53, Eqs. (29)–(35)]).
Thus $\operatorname{Tr}(Y^{-1}\rho)\leq\lambda$ implies
$\operatorname{Tr}(Y\rho)\geq 1/\lambda$.

When $\rho$ and $\sigma$ are (positive definite) density matrices, the unique
solution to Eq. (55) is given by the matrix geometric mean
$Y=\sigma^{-1}\\#\rho$. We can therefore proceed as before to compute
$\operatorname{Tr}(Y\xi)$, except now we do not need to be concerned with
state preparation of $\rho$ and $\sigma$, and we can assume that we are given
copies of $\rho$ and $\sigma$. Thus, given access to $U_{Y}$, we can estimate
the following expectation:

$\displaystyle\operatorname{Tr}(Y\xi)\approx\kappa_{\sigma}\operatorname{Tr}((\xi\otimes|0\rangle\\!\langle
0|)U_{Y}),$ (56)

where the $\kappa_{\sigma}$ constant follows from Lemma 8 and the error in the
above estimate is upper-bounded by $\epsilon$. We then have the following
result.

######  Theorem 15.

Suppose that we are given the block-encodings $U_{\rho}$ and $U_{\sigma}$,
where $\rho\in\mathcal{A}$, $\sigma\notin\mathcal{A}$ and that we are also
given access to multiple copies of $\xi$. Suppose further that
$\kappa_{\rho},\kappa_{\sigma}=O(\operatorname{poly}\log D)$. Then computing
$\operatorname{Tr}(Y\xi)$ for the optimal $Y$ in Eq. (55) to precision
$\epsilon>0$ has a query and gate complexity $O(\operatorname{poly}(\log
D,1/\epsilon))$.

###### Proof.

From Lemma 8 we can construct a $(\kappa_{\sigma},O(\log D),\epsilon)$-block-
encoding of $Y$ with gate and query complexity
$O(\operatorname{poly}(\kappa_{\rho},\kappa_{\sigma},\log(1/\epsilon)))$.
Considering that $\kappa_{\rho},\kappa_{\sigma}=O(\operatorname{poly}\log D)$,
applying the unitary $U_{Y}$ to $\xi\otimes|0\rangle\\!\langle 0|$, and using
the Hadamard test (Lemma 31) with $\xi\otimes|0\rangle\\!\langle 0|$, we
recover $\operatorname{Tr}(Y\xi)$ to precision $\epsilon$ with query and gate
complexity $O(\operatorname{poly}(\log D,1/\epsilon))$. ∎

We emphasise that this problem is entirely quantum in nature as we are given
directly only quantum data.

######  Remark 2.

The assumption of $U_{\rho}$ and $U_{\sigma}$ as block-encodings of $\rho$ and
$\sigma$, respectively, is without loss of generality in practice. There are
two quantum input models for quantum states that are commonly employed in
quantum algorithms:

  * •

Quantum query access model. In this model, quantum unitary oracles
$\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ are given such that they
prepare purifications of $\rho$ and $\sigma$, respectively. By the technique
of purified density matrix in [21] (see Lemma 30), we can implement $U_{\rho}$
and $U_{\sigma}$ from $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ with
query and gate complexity $\widetilde{O}\lparen 1\rparen$. Therefore, Theorem
15 can be adapted to the quantum query access model with query and gate
complexity $O(\operatorname{poly}(\log D,1/\epsilon))$.

  * •

Quantum sample access model. In this model, independent and identical copies
of $\rho$ and $\sigma$ are given. By the technique of density matrix
exponentiation [54, 55], we can implement unitary operators that are block-
encodings of $\rho$ and $\sigma$ using their copies (which was first noted in
[56] and later investigated in [57, 58, 59]). In this way, Theorem 15 can be
adapted to the quantum sample access model with sample and gate complexity
$O(\operatorname{poly}(\log D,1/\epsilon))$.

A very interesting observation to note here is that the matrix geometric mean
solution $Y$ to Eq. (55) is precisely the Fuchs–Caves observable [15], which
is important for distinguishing two states $\rho$ and $\sigma$. From this
observation, we can motivate the Fuchs–Caves observable as the observable that
gives rise to a kind of ‘optimal witness’ that distinguishes $\rho$ and
$\sigma$ and the value of this ‘witness’ is precisely quantum fidelity, as
shown in the next section. This provides an alternative motivation for the
form of quantum fidelity between two mixed states from a metric learning
viewpoint. In fact, a protocol involving a measurement of the Fuchs–Caves
observable also achieves an upper bound on sample complexity for the quantum
hypothesis testing problem in distinguishing $\rho$ and $\sigma$ [60, Appendix
F]. Thus, up to constant factors, the strategy also minimises the number of
copies of each state used for a given tolerated precision in distinguishing
the states. We note that the loss function also appears in Eq. (6) in [43],
but this is motivated from a different perspective.

####  IV.1.3 Extension to weighted geometric mean metric learning

The two terms in the loss function in Eq. (53), involving $\sigma$ and $\rho$
respectively, have equal weights. This means the learning algorithm deems
closeness to $\rho$ and farness to $\sigma$ of equal ‘importance’. However,
there are scenarios, especially in anomaly detection, where asymmetry is
preferable. For example, this occurs when there is a higher cost in getting
false negatives.

Modifying Eq. (53) by simply multiplying each of the two terms by different
constants $\alpha,\beta$ leads to
$L(Y)=\alpha\text{Tr}(Y\sigma)+\beta\text{Tr}(Y^{-1}\rho)$. However, this only
rescales the optimal solution $Y\rightarrow(\beta/\alpha)^{1/2}Y$ by a
constant factor, as observed in [25]. A new loss function is therefore
necessary for the asymmetric case.

Following [25], one can first observe that the solution $Y=\sigma^{-1}\\#\rho$
is in fact also a solution to the following optimisation problem when $t=1/2$:

$\displaystyle\min_{\tilde{Y}\geq 0}L_{t}(\tilde{Y}),\qquad
L_{t}(\tilde{Y})\coloneqq(1-t)\delta(\tilde{Y},\sigma^{-1})+t\delta(\tilde{Y},\rho),\qquad
t\in[0,1],$ (57)

where $\delta$ is the geodesic distance defined in Eq. (8). While the
mathematical proof is more involved, this fact can easily be understood from
the geometric viewpoint. Here $Y=\sigma^{-1}\\#\rho$ can be understood as the
midpoint along the unique geodesic in Riemannian space joining $\sigma^{-1}$
and $\rho$. When $t=1/2$, the optimal $\tilde{Y}$ is then the point along this
geodesic that simultaneously minimises the distance between $\tilde{Y}$ and
$\sigma^{-1}$, as well as $\tilde{Y}$ and $\rho$. This clearly must be the
midpoint. Similar geometric reasoning leads one to generalise to $t\neq 1/2$
where the solution to Eq. (57) is the weighted matrix geometric mean
$\tilde{Y}=\sigma^{-1}\\#_{t}\rho$. That this is the unique solution to Eq.
(57) is a special case (the $n=2$ case) in [61] and proofs can also be found
in [27, Chapter 6]. Also see [25] for a discussion in the context of geometric
mean metric learning.

We can proceed similarly to 1-class quantum learning algorithm with equal
weights as described in the previous section. The goal is also to output
$\text{Tr}(\tilde{Y}\xi)$ for some input test state $\xi$. Here we require
instead the construction of block-encodings of the weighted matrix geometric
mean, as given in Lemma 11. However, for any $t>0$ ($p>0$ in Lemma 11), we see
that there is no scaling difference for constructing the block-encoding for
the weighted version. Thus the cost, up to constant and logarithmic factors,
is identical for the quantum weighted geometric mean metric learning algorithm
as for the unweighted version in Theorem 15.

###  IV.2 Estimation of quantum fidelity and geometric Rényi relative
entropies

Here we describe our quantum algorithms for estimating quantum fidelity and
geometric Rényi relative entropies using our quantum subroutines for preparing
block-encodings of the standard and weighted matrix geometric means.

####  IV.2.1 Fidelity

The fidelity between two mixed quantum states is defined by [62]

$F\lparen\rho,\sigma\rparen\coloneqq\Tr\left\lparen\left\lparen\sigma^{1/2}\rho\sigma^{1/2}\right\rparen^{1/2}\right\rparen,$
(58)

which is a commonly considered measure of the closeness of or similarity
between two quantum states. Estimating the value of fidelity is a fundamental
task in quantum information theory. When given matrix descriptions of the
states $\rho$ and $\sigma$, it can be calculated directly using the formula
above or as the solution to a semidefinite optimisation problem [43].
Recently, several time-efficient quantum algorithms for fidelity estimation
have been developed when one has access to state-preparation circuits of
$\rho$ and $\sigma$ [63, 57, 64].

Here, we introduce a new approach for fidelity estimation that is based on the
Fuchs–Caves observable [15]. For two quantum states $\rho$ and $\sigma$, this
observable is given by $M=\sigma^{-1}\\#\rho$. Then, the fidelity between
$\rho$ and $\sigma$ can be represented as the expectation of $M$ with respect
to $\sigma$ (cf. [39, Eq. (9.159)]):

$F\lparen\rho,\sigma\rparen=\Tr\lparen M\sigma\rparen.$ (59)

######  Theorem 16 (Fidelity estimation via Fuchs–Caves observable).

Suppose that $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ prepare
purifications of mixed quantum states $\rho$ and $\sigma$, respectively. Then,
we can estimate $F\lparen\rho,\sigma\rparen$ within additive error $\epsilon$
using
$\widetilde{O}\lparen\min\\{\kappa_{\rho}^{2},\kappa_{\sigma}^{2}\\}\cdot\kappa_{\sigma}\kappa_{\rho}/\epsilon\rparen$
queries to $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$, where
$\kappa_{\rho},\kappa_{\sigma}>0$ are such that $\rho\geq I/\kappa_{\rho}$ and
$\sigma\geq I/\kappa_{\sigma}$.

###### Proof.

Suppose that $\rho$ and $\sigma$ are $n$-qubit mixed quantum states and
$\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ are $\lparen n+a\rparen$-qubit
unitary operators. By Lemma 30, we can implement two unitary operators
$U_{\rho}$ and $U_{\sigma}$ that are $\lparen 1,n+a,0\rparen$-block-encodings
of $\rho$ and $\sigma$ using $O\lparen 1\rparen$ queries to
$\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$, respectively. Then, by
applying Lemma 8, we can implement a
$\lparen\kappa_{\sigma},b,\delta\rparen$-block-encoding $U_{M}$ of
$M=\sigma^{-1}\\#\rho$, using
$\widetilde{O}\lparen\kappa_{\sigma}\kappa_{\rho}\log\lparen
1/\delta\rparen\rparen$ queries to $U_{\rho}$ and
$\widetilde{O}\lparen\kappa_{\sigma}^{2}\kappa_{\rho}\log^{2}\lparen
1/\delta\rparen\rparen$ queries to $U_{\sigma}$, where $b=O\lparen
n+a\rparen$, and $\kappa_{\rho}$ and $\kappa_{\sigma}$ satisfy $\rho\geq
I/\kappa_{\rho}$ and $\sigma\geq I/\kappa_{\sigma}$.

By the Hadamard test (given in Lemma 31), there is a quantum circuit $C$ that
outputs $0$ with probability
$\frac{1}{2}\left(1+\Re{\Tr\lparen\bra{0}_{b}U_{M}\ket{0}_{b}\sigma\rparen}\right)$,
using one query to $U_{M}$ and one sample of $\sigma$. By noting that

$\absolutevalue*{\kappa_{\sigma}\Re{\Tr\lparen\bra{0}_{b}U_{M}\ket{0}_{b}\sigma\rparen}-\Tr\lparen
M\sigma\rparen}\leq\Theta\lparen\delta\rparen,$ (60)

we conclude that an $O\lparen\epsilon/\kappa_{\sigma}\rparen$-estimate of
$\Re{\Tr\lparen\bra{0}_{b}U_{M}\ket{0}_{b}\sigma\rparen}$ with
$\delta=\Theta\lparen\epsilon/\kappa_{\sigma}\rparen$ suffices to obtain an
$\epsilon$-estimate of $\Tr\lparen M\sigma\rparen$ (which is the fidelity
according to Eq. (59)). By quantum amplitude estimation (given in Lemma 33),
this can be done using $O\lparen\kappa_{\sigma}/\epsilon\rparen$ queries to
$C$.

In summary, an $\epsilon$-estimate of $F\lparen\rho,\sigma\rparen$ can be
obtained by using
$\widetilde{O}\lparen\kappa_{\sigma}^{3}\kappa_{\rho}/\epsilon\rparen$ queries
to $\mathcal{O}_{\sigma}$ and
$\widetilde{O}\lparen\kappa_{\sigma}^{2}\kappa_{\rho}/\epsilon\rparen$ queries
to $\mathcal{O}_{\rho}$. The proof is completed by taking the minimum over
symmetric cases (i.e., simply flipping the role of $\rho$ and $\sigma$ since
the fidelity formula is symmetric under this exchange). ∎

The current best quantum query complexity of fidelity estimation is
$\widetilde{O}\lparen r^{2.5}/\epsilon^{5}\rparen$, due to [57], where $r$ is
the lower rank of the two input mixed quantum states. In comparison, our
quantum algorithm for fidelity estimation based on the Fuchs–Caves observable,
as given in Theorem 16, has a better dependence on the additive error
$\epsilon$, if $\kappa_{\rho}$ and $\kappa_{\sigma}$ are known in advance.

Moreover, we note that the $\epsilon$-dependence of the quantum algorithm
given in Theorem 16 is optimal (up to polylogarithmic factors), as stated in
Lemma 17 below.

######  Lemma 17 (Optimal $\epsilon$-dependence of fidelity estimation).

Suppose that $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ prepare
purifications of mixed quantum states $\rho$ and $\sigma$, respectively,
satisfying $\rho\geq I/\kappa_{\rho}$ and $\sigma\geq I/\kappa_{\sigma}$ for
$\kappa_{\rho},\kappa_{\sigma}>0$. Then, every quantum algorithm that
estimates $F\lparen\rho,\sigma\rparen$ within additive error $\epsilon$
requires query complexity $\Omega\lparen 1/\epsilon\rparen$ even if
$\kappa_{\rho}=\kappa_{\sigma}=\Theta\lparen 1\rparen$.

###### Proof.

See Appendix I. ∎

####  IV.2.2 Geometric fidelity and geometric Rényi relative entropy

Here we present, to the best of our knowledge, the first quantum algorithm for
computing the geometric $\alpha$-Rényi relative entropy, as introduced in [12,
13]. For $\alpha\in(0,1)\cup(1,2]$, the geometric $\alpha$-Rényi relative
entropy is defined as (see, e.g., [14, Eq. (9)] and [65, Eq. (7.6.1)])

$\hat{D}_{\alpha}\lparen\rho\|\sigma\rparen\coloneqq\frac{1}{\alpha-1}\log\hat{F}_{\alpha}\lparen\rho,\sigma\rparen,$
(61)

where

$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen\coloneqq\Tr\left\lparen\rho\\#_{1-\alpha}\sigma\right\rparen=\Tr\left\lparen\sigma\\#_{\alpha}\rho\right\rparen$
(62)

is known as the geometric $\alpha$-Rényi relative quasi-entropy. When
$\alpha\in(0,1)$, we also refer to $\hat{F}_{\alpha}\lparen\rho,\sigma\rparen$
as the geometric $\alpha$-fidelity. In particular, for the case of
$\alpha=1/2$, the quantity $\hat{F}_{1/2}\lparen\rho,\sigma\rparen$ is the
geometric fidelity (also known as the Matsumoto fidelity) [16, 17]. The
$\alpha$-geometric Rényi relative entropy has several uses in quantum
information theory, especially in analysing protocols involving feedback [14,
29, 66].

Here we present quantum algorithms in Theorems 18 and 19 for computing the
geometric Rényi relative (quasi-)entropy.

######  Theorem 18.

Suppose that $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ prepare
purifications of mixed quantum states $\rho$ and $\sigma$, respectively. Then,
for $\alpha\in\lparen 0,1\rparen\cup(1,2]$, we can estimate
$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen$ to within additive error
$\epsilon$ using

  * •

$\widetilde{O}\lparen\min\\{\kappa_{\rho},\kappa_{\sigma}\\}^{\min\\{1+\alpha,2-\alpha\\}}\cdot\kappa_{\rho}\kappa_{\sigma}/\epsilon\rparen$
queries for $\alpha\in(0,1)$, and

  * •

$\widetilde{O}\lparen\min\\{\kappa_{\rho}\kappa_{\sigma}^{\alpha-1},\kappa_{\rho}^{\alpha-1}\kappa_{\sigma},\kappa_{\rho}^{1+\alpha},\kappa_{\sigma}^{1+\alpha}\\}\cdot\kappa_{\rho}\kappa_{\sigma}/\epsilon\rparen$
queries for $\alpha\in(1,2]$

to $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$, where
$\kappa_{\rho},\kappa_{\sigma}>0$ satisfy $\rho\geq I/\kappa_{\rho}$ and
$\sigma\geq I/\kappa_{\sigma}$.

In particular, when $\alpha=1/2$, $\hat{F}_{1/2}\lparen\rho,\sigma\rparen$ is
the geometric fidelity (also known as the Matsumoto fidelity), which can be
estimated using
$\widetilde{O}\lparen\min\\{\kappa_{\rho},\kappa_{\sigma}\\}^{3/2}\cdot\kappa_{\rho}\kappa_{\sigma}/\epsilon\rparen$
queries to $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$.

###### Proof.

See Appendix H. ∎

######  Theorem 19.

Suppose that $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ prepare
purifications of mixed quantum states $\rho$ and $\sigma$, respectively. Then,
for $\alpha\in\lparen 0,1\rparen\cup(1,2]$, we can estimate
$\hat{D}_{\alpha}\lparen\rho\|\sigma\rparen$ within additive error $\epsilon$
using

  * •

$\widetilde{O}\lparen\min\\{\kappa_{\rho},\kappa_{\sigma}\\}^{\min\\{1+\alpha,2-\alpha\\}}\cdot\kappa_{\rho}\kappa_{\sigma}^{2-\alpha}/\epsilon\rparen$
queries for $\alpha\in(0,1)$, and

  * •

$\widetilde{O}\lparen\min\\{\kappa_{\rho}\kappa_{\sigma}^{\alpha-1},\kappa_{\rho}^{\alpha-1}\kappa_{\sigma},\kappa_{\rho}^{1+\alpha},\kappa_{\sigma}^{1+\alpha}\\}\cdot\kappa_{\rho}^{\alpha}\kappa_{\sigma}/\epsilon\rparen$
queries for $\alpha\in(1,2]$

to $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$, where
$\kappa_{\rho},\kappa_{\sigma}>0$ satisfy $\rho\geq I/\kappa_{\rho}$ and
$\sigma\geq I/\kappa_{\sigma}$.

###### Proof.

See Appendix H. ∎

Notably, we show that our quantum algorithm for estimating the geometric
fidelity $\hat{F}_{1/2}\lparen\rho,\sigma\rparen$ achieves an optimal
$\epsilon$-dependence. The optimality also holds for
$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen$ with $\alpha\in(0,1)$.

######  Lemma 20 (Optimal $\epsilon$-dependence of geometric $\alpha$-fidelity
estimation).

Suppose that $\mathcal{O}_{\rho}$ and $\mathcal{O}_{\sigma}$ prepare
purifications of mixed quantum states $\rho$ and $\sigma$, respectively, with
$\rho\geq I/\kappa_{\rho}$ and $\sigma\geq I/\kappa_{\sigma}$, where
$\kappa_{\rho},\kappa_{\sigma}>0$. Then, for any constant $\alpha\in(0,1)$,
any quantum algorithm that estimates
$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen$ within additive error $\epsilon$
requires query complexity $\Omega\lparen 1/\epsilon\rparen$ even if
$\kappa_{\rho}=\kappa_{\sigma}=\Theta\lparen 1\rparen$, where
$\Omega\lparen\cdot\rparen$ hides a constant factor that depends only on
$\alpha$.

###### Proof.

See Appendix I. ∎

It remains an open problem for optimality still holds if $\alpha\in(1,2]$.
However, note that when $\alpha\in(1,2]$, the inequality
$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen\geq 1$ holds, and so
$\hat{F}_{\alpha}\lparen\rho,\sigma\rparen$ cannot be interpreted as a
fidelity for these values of $\alpha$; thus, different techniques are required
in order to establish optimality.

##  V BQP-hardness

In this section, we consider the hardness of computing the matrix geometric
mean. Precisely, we show that our quantum algorithm for matrix geometric means
(given in Lemma 8) can be used to solve a $\mathsf{BQP}$-complete problem
(defined in Problem 1). Roughly speaking, this problem pertains to testing a
certain property of the matrix geometric mean of two well-conditioned sparse
matrices.

######  Problem 1 (Matrix geometric mean).

For functions $\kappa_{A}\colon\mathbb{N}\to\mathbb{N}$ and
$\kappa_{C}\colon\mathbb{N}\to\mathbb{N}$, let
$\textup{MGM}\lparen\kappa_{A},\kappa_{C}\rparen$ be a decision problem
defined as follows. For a size-$n$ instance of
$\textup{MGM}\lparen\kappa_{A},\kappa_{C}\rparen$, let $N=2^{n}$ and let
$A,C\in\mathbb{C}^{N\times N}$ be $O\lparen 1\rparen$-sparse positive definite
matrices with $I/\kappa_{A}\lparen n\rparen\leq A\leq I$ and
$I/\kappa_{C}\lparen n\rparen\leq C\leq I$, given by a
$\operatorname{poly}\lparen n\rparen$-size uniform classical circuit
$\mathcal{C}_{n}$ such that, for every $1\leq j\leq N$, the circuit
$\mathcal{C}_{n}\lparen j\rparen$ computes the positions and values of the
non-zero entries in the $j$-th row of $A$ and $C$. Let
$Y\in\mathbb{C}^{N\times N}$ be the matrix geometric mean of $A$ and $C$ such
that $YAY=C$. The task is to decide which of the following is the case,
promised that one of the two holds:

  * •

_Yes_ : $\langle\psi|M|\psi\rangle\geq 2/3$;

  * •

_No_ : $\langle\psi|M|\psi\rangle\leq 1/3$,

where $\ket{\psi}\coloneqq\frac{Y^{2}\ket{0}}{\lVert Y^{2}\ket{0}\rVert}$ and
$M=|0\rangle\\!\langle 0|\otimes I_{N/2}$ measures the first qubit.

######  Theorem 21.

$\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is
$\mathsf{BQP}$-complete.

###### Proof.

The proof consists of two parts: Lemma 22 and Lemma 23.

  1. 1.

In Lemma 22, we state that $\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is $\mathsf{BQP}$-hard;
the proof employs a reduction of the Quantum Linear Systems Problem (QLSP).

  2. 2.

In Lemma 23, we state that $\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is in $\mathsf{BQP}$; the
proof employs the quantum algorithm for the matrix geometric mean given in
Lemma 8.

∎

######  Lemma 22.

$\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is $\mathsf{BQP}$-hard.

###### Proof.

We consider the Quantum Linear Systems Problem (QLSP) defined as follows.

######  Problem 2 (QLSP).

For functions $\kappa\colon\mathbb{N}\to\mathbb{N}$, let
$\textup{QLSP}\lparen\kappa\rparen$ be a decision problem defined as follows.
For a size-$n$ instance of $\textup{QLSP}\lparen\kappa\rparen$, let $N=2^{n}$
and $A\in\mathbb{C}^{N\times N}$ be an $O\lparen 1\rparen$-sparse Hermitian
matrix such that $I/\kappa\lparen n\rparen\leq A\leq I$, given by a
$\operatorname{poly}\lparen n\rparen$-size uniform classical circuit
$\mathcal{C}_{n}$ such that for every $1\leq j\leq N$, $\mathcal{C}_{n}\lparen
j\rparen$ computes the positions and values of the non-zero entries in the
$j$-th row of $A$. The task is to decide which of the following is the case,
promised that one of the two holds:

  * •

_Yes_ item: $\langle\psi|M|\psi\rangle\geq 2/3$;

  * •

_No_ item: $\langle\psi|M|\psi\rangle\leq 1/3$,

where $\ket{\psi}\coloneqq\frac{A^{-1}\ket{0}}{\lVert A^{-1}\ket{0}\rVert}$
and $M=|0\rangle\\!\langle 0|\otimes I_{N/2}$ measures the first qubit.

It was shown in [1] that $\textup{QLSP}\lparen\operatorname{poly}\lparen
n\rparen\rparen$ is $\mathsf{BQP}$-complete. Here, we reduce
$\textup{QLSP}\lparen\operatorname{poly}\lparen n\rparen\rparen$ to
$\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$, and therefore show the
$\mathsf{BQP}$-hardness of $\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$.

Consider any instance (matrix) $A\in\mathbb{C}^{N\times N}$ of
$\textup{QLSP}\lparen\kappa\rparen$, where $N=2^{n}$ and
$\kappa=\operatorname{poly}\lparen n\rparen$. We choose
$C=I\in\mathbb{C}^{N\times N}$ to be the identity matrix, which is a
$1$-sparse Hermitian matrix and each of whose rows can be easily computed.
Note that the matrix geometric mean $Y$ of $A^{-1}$ and $C$ is
$Y=A^{-1}\\#C=A^{-1/2}$. Then, it can be seen that $Y^{2}=A^{-1}$ and thus
$\ket{\psi_{Y}}=Y^{2}\ket{0}/\lVert Y^{2}\ket{0}\rVert=A^{-1}\ket{0}/\lVert
A^{-1}\ket{0}\rVert=\ket{\psi_{A}}$. Consequently, any quantum algorithm that
determines whether $\bra{\psi_{Y}}M\ket{\psi_{Y}}\geq 2/3$ or
$\bra{\psi_{Y}}M\ket{\psi_{Y}}\leq 1/3$ with success probability at least
$2/3$ can be used to determine whether $\bra{\psi_{A}}M\ket{\psi_{A}}\geq 2/3$
or $\bra{\psi_{A}}M\ket{\psi_{A}}\leq 1/3$. In summary,
$\textup{QLSP}\lparen\kappa\rparen$ can be reduced to
$\textup{MGM}\lparen\kappa,1\rparen$ through the above encoding. Therefore,
$\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is $\mathsf{BQP}$-hard. ∎

######  Lemma 23.

$\textup{MGM}\lparen\operatorname{poly}\lparen
n\rparen,\operatorname{poly}\lparen n\rparen\rparen$ is in $\mathsf{BQP}$.

###### Proof.

See Appendix J. ∎

##  VI Discussion

We constructed efficient block-encodings of the matrix geometric mean (and
weighted matrix geometric mean). These are unique solutions to the simplest
algebraic Riccati equations – quadratically nonlinear system of matrix
equations. Unlike the output of most quantum algorithms for linear systems of
equations, these solutions of the nonlinear matrix equations are not embedded
in pure quantum states, but rather in terms of observables from which we can
extract expectation values.

This allows us to introduce a new class of algorithms for quantum learning,
called quantum geometric mean metric learning. For example, this can be
applied in a purely quantum setting for picking out anomalous quantum states.
This can also be adapted to the case of flexible weights on the cost of
flagging an anomaly. The new quantum subroutines can also be used for the
first quantum algorithm, to the best of our knowledge, to compute the
geometric Rényi relative entropies and new quantum algorithms to compute
quantum fidelity by means of the Fuchs–Caves observable. In the latter case,
we demonstrate optimal scaling $\Omega(1/\epsilon)$ in precision.

While most of the applications introduced above are for quantum problems for
which there is no direct classical equivalent (although the quantum learning
algorithm can also be applied to learning Euclidean distances for classical
data), there are potential benefits that the new quantum subroutine can have
over purely classical methods. This could be exploited for future
applications. For example, classical numerical algorithms to compute the
matrix geometric mean have cost $O(\operatorname{poly}D)$ for $D\times D$
matrices [51, 52]. The same is also true for solving the differential matrix
Riccati equation and algebraic matrix Riccati equation [18] through iterative
methods and other methods based on finding the eigendecomposition of a larger
matrix [67]. For quantum processing on the other hand, we showed conditions
under which the block-encodings of some of these solutions can be obtained
with cost $O(\operatorname{poly}\log D)$.

For example, there are many classical problems for which it is important to
compute the matrix geometric mean between two matrices. They appear in imaging
[68, 69] and in the analysis of multiport electrical networks [70]. The
algebraic Riccati equation of the form in Eq. (9) also appears in optimal
control and Kalman filters. Under the assumptions in Lemma 5 when uniqueness
of its solution is also satisfied, it can be possible to construct its block-
encoding in Lemma 9. Although these assumptions are not generally satisfied,
this still gives an idea of the extent and reach of the matrix geometric mean.
Extensions of our algorithms to the matrix geometric mean consisting of more
than two matrices can also be explored, which already find applications in
areas like elasticity and radars [71, 72]. It is also intriguing to consider
purely quantum extensions of these problems. The main difficulty associated
with constructing block-encodings of multivariate geometric means is that they
are not known to have an analytical form as they do in the bivariate case;
rather, they are constructed as the solutions of nonlinear equations
generalizing the simple algebraic Riccati equation [10].

In addition to usefulness in applications, the standard and weighted matrix
geometric means also have an elegant interpretation in terms of geodesics in
Riemannian space. Despite the importance and beauty of Riemannian geometry in
mathematics and other areas in physics, sensing, and machine learning, it has
not appeared too much in quantum computation yet, apart from very notable
exceptions like [73]. This geometric perspective is useful in understanding
the weighted quantum learning algorithm, and we showed how it provided an
alternative motivation for the form of quantum fidelity via the Fuchs–Caves
observable. There is more potential here for the matrix geometric mean to
bring the ideas of geometry closer to quantum information and computation.

###### Acknowledgements.

NL acknowledges support from the Science and Technology Program of Shanghai,
China (21JC1402900), NSFC grant No. 12341104, the Shanghai Jiao Tong
University 2030 Initiative, and the Fundamental Research Funds for the Central
Universities. QW acknowledges support from the MEXT Quantum Leap Flagship
Program (MEXT Q-LEAP) Grant No. JPMXS0120319794. MMW acknowledges support from
the NSF under grants 2329662 and 2315398. ZZ acknowledges support from the
Sydney Quantum Academy, NSW, Australia.

## References

  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical Review Letters 103, 150502 (2009).
  * Childs  _et al._ [2021] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Jin  _et al._ [2022] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations via Schrödingerisation (2022), arXiv:2212.13969 .
  * Jin  _et al._ [2023a] S. Jin, N. Liu, and Y. Yu, Quantum simulation of partial differential equations: applications and detailed analysis, Physical Review A 108, 032603 (2023a).
  * An  _et al._ [2022] D. An, J.-P. Liu, D. Wang, and Q. Zhao, A theory of quantum differential equation solvers: limitations and fast-forwarding (2022), arXiv:2211.05246 .
  * Lancaster and Rodman [1995] P. Lancaster and L. Rodman,  _Algebraic Riccati Equations_ (Clarendon Press, 1995).
  * Salgado  _et al._ [1988] M. Salgado, R. Middleton, and G. C. Goodwin, Connection between continuous and discrete Riccati equations with applications to Kalman filtering, in  _IEE Proceedings D (Control Theory and Applications)_, Vol. 135 (IET, 1988) pp. 28–34.
  * Coppel [1974] W. Coppel, Matrix quadratic equations, Bulletin of the Australian Mathematical Society 10, 377 (1974).
  * Lawson and Lim [2001] J. D. Lawson and Y. Lim, The geometric mean, matrices, metrics, and more, The American Mathematical Monthly 108, 797 (2001).
  * Lawson and Lim [2021] J. D. Lawson and Y. Lim, The expanding universe of the geometric mean (2021), arXiv:2110.13371 [math.OA] .
  * Janati  _et al._ [2020] H. Janati, B. Muzellec, G. Peyré, and M. Cuturi, Entropic optimal transport between unbalanced Gaussian measures has a closed form, in  _Advances in Neural Information Processing Systems 33_ (2020) pp. 10468–10479.
  * Matsumoto [2013] K. Matsumoto, A new quantum version of $f$-divergence (2013), arXiv:1311.4722 .
  * Matsumoto [2018] K. Matsumoto, A new quantum version of $f$-divergence, in  _Reality and Measurement in Algebraic Quantum Theory_, Vol. 261, edited by M. Ozawa, J. Butterfield, H. Halvorson, M. Rédei, Y. Kitajima, and F. Buscemi (Springer Singapore, Singapore, 2018) pp. 229–273, series Title: Springer Proceedings in Mathematics & Statistics.
  * Fang and Fawzi [2021] K. Fang and H. Fawzi, Geometric Rényi divergence and its applications in quantum channel capacities, Communications in Mathematical Physics 384, 1615 (2021).
  * Fuchs and Caves [1995] C. A. Fuchs and C. M. Caves, Mathematical techniques for quantum communication theory, Open Systems & Information Dynamics 3, 345 (1995).
  * Matsumoto [2010] K. Matsumoto, Reverse test and quantum analogue of classical fidelity and generalized fidelity (2010), arXiv:1006.0302 .
  * Cree and Sikora [2020] S. Cree and J. Sikora, A fidelity measure for quantum states based on the matrix geometric mean (2020), arXiv:2006.06918 .
  * Ramesh  _et al._ [1989] A. V. Ramesh, S. Utku, and J. A. Garba, Computational complexities and storage requirements of some Riccati equation solvers, Journal of Guidance, Control, and Dynamics 12, 469 (1989).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. 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.
  * Chakraborty  _et al._ [2019] S. Chakraborty, A. Gilyén, and S. Jeffery, The power of block-encoded matrix powers: improved regression techniques via faster hamiltonian simulation, in  _46th International Colloquium on Automata, Languages, and Programming (ICALP 2019)_, Vol. 132 (2019) pp. 33:1–33:14.
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Jin  _et al._ [2023b] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum algorithms via linear representations for nonlinear ordinary and partial differential equations, Journal of Computational Physics 487, 112149 (2023b).
  * Joseph [2020] I. Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics, Physical Review Research 2, 043102 (2020).
  * Liu  _et al._ [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118, e2026805118 (2021).
  * Zadeh  _et al._ [2016] P. Zadeh, R. Hosseini, and S. Sra, Geometric mean metric learning, in  _Proceedings of the 33rd International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 48 (PMLR, New York, New York, USA, 2016) pp. 2464–2471.
  * Liu and Rebentrost [2018] N. Liu and P. Rebentrost, Quantum machine learning for quantum anomaly detection, Physical Review A 97, 042315 (2018).
  * Bhatia [2007] R. Bhatia,  _Positive Definite Matrices_ (Princeton University Press, 2007).
  * Furuichi  _et al._ [2004] S. Furuichi, K. Yanagi, and K. Kuriyama, Fundamental properties of Tsallis relative entropy, Journal of Mathematical Physics 45, 4868 (2004).
  * Katariya and Wilde [2021a] V. Katariya and M. M. Wilde, Geometric distinguishability measures limit quantum channel estimation and discrimination, Quantum Information Processing 20, 78 (2021a).
  * Ran and Vreugdenhil [1988] A. C. M. Ran and R. Vreugdenhil, Existence and comparison theorems for algebraic Riccati equations for continuous- and discrete-time systems, Linear Algebra and its Applications 99, 63 (1988).
  * Shurbet  _et al._ [1974] G. Shurbet, T. Lewis, and T. Boullion, Quadratic matrix equations, Ohio Journal of Science 74, 273 (1974).
  * Richardson and Kwong [1986] T. J. Richardson and R. H. Kwong, On positive definite solutions to the algebraic Riccati equation, Systems & Control Letters 7, 99 (1986).
  * Lancaster and Rodman [1980] P. Lancaster and L. Rodman, Existence and uniqueness theorems for the algebraic Riccati equation, International Journal of Control 32, 285 (1980).
  * Wimmer [1984] H. K. Wimmer, The algebraic Riccati equation: conditions for the existence and uniqueness of solutions, Linear Algebra and its Applications 58, 441 (1984).
  * Boyd and Barratt [1991] S. P. Boyd and C. H. Barratt,  _Linear Controller Design: Limits of Performance_ , Vol. 7 (Prentice Hall, 1991).
  * Kubo and Ando [1980] F. Kubo and T. Ando, Means of positive linear operators, Mathematische Annalen 246, 205 (1980).
  * Nakamura [2007] N. Nakamura, Geometric operator mean induced from the Riccati equation, Scientiae Mathematicae Japonicae 66, 83 (2007).
  * Pedersen and Takesaki [1972] G. K. Pedersen and M. Takesaki, The operator equation $tht=k$, Proceedings of the American Mathematical Society 36, 311 (1972).
  * Wilde [2017] M. M. Wilde,  _Quantum Information Theory_, 2nd ed. (Cambridge University Press, 2017).
  * Alsing  _et al._ [2024] P. M. Alsing, C. Cafaro, and S. Ray, Geodesics for mixed quantum states via their geometric mean operator (2024), arXiv:2404.04136 [quant-ph] .
  * Paulsen [2003] V. Paulsen,  _Completely Bounded Maps and Operator Algebras_ (Cambridge University Press, 2003).
  * Furuta [1988] T. Furuta, The operator equation $t(h^{1n}t)^{n}=k$, Linear Algebra and its Applications 109, 149 (1988).
  * Watrous [2013] J. Watrous, Simpler semidefinite programs for completely bounded norms, Chicago Journal of Theoretical Computer Science 10.4086/cjtcs.2013.008 (2013).
  * Kulis [2013] B. Kulis, Metric learning: A survey, Foundations and Trends in Machine Learning 5, 287 (2013).
  * Gleinig and Hoefler [2021] N. Gleinig and T. Hoefler, An efficient algorithm for sparse quantum state preparation, in  _2021 58th ACM/IEEE Design Automation Conference (DAC)_ (IEEE, 2021) pp. 433–438.
  * Sun  _et al._ [2023] X. Sun, G. Tian, S. Yang, P. Yuan, and S. Zhang, Asymptotically optimal circuit depth for quantum state preparation and general unitary synthesis, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 42, 3301 (2023).
  * Zhang  _et al._ [2022] X.-M. Zhang, T. Li, and X. Yuan, Quantum state preparation with optimal circuit depth: Implementations and applications, Physical Review Letters 129, 230504 (2022).
  * Barenco  _et al._ [1997] A. Barenco, A. Berthiaume, D. Deutsch, A. Ekert, R. Jozsa, and C. Macchiavello, Stabilization of quantum computations by symmetrization, SIAM Journal on Computing 26, 1541 (1997).
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, Quantum fingerprinting, Physical Review Letters 87, 167902 (2001).
  * Garcia-Escartin and Chamorro-Posada [2013] J. C. Garcia-Escartin and P. Chamorro-Posada, Swap test and Hong-Ou-Mandel effect are equivalent, Physical Review A 87, 052330 (2013).
  * Bini and Iannazzo [2011] D. A. Bini and B. Iannazzo, A note on computing matrix geometric means, Advances in Computational Mathematics 35, 175 (2011).
  * Iannazzo [2016] B. Iannazzo, The geometric mean of two matrices from a computational viewpoint, Numerical Linear Algebra with Applications 23, 208 (2016).
  * Katariya and Wilde [2021b] V. Katariya and M. M. Wilde, RLD Fisher information bound for multiparameter estimation of quantum channels, New Journal of Physics 23, 073040 (2021b).
  * Lloyd  _et al._ [2014] S. Lloyd, M. Mohseni, and P. Rebentrost, Quantum principal component analysis, Nature Physics 10, 631 (2014).
  * Kimmel  _et al._ [2017] S. Kimmel, C. Y.-Y. Lin, G. H. Low, M. Ozols, and T. J. Yoder, Hamiltonian simulation with optimal sample complexity, npj Quantum Information 3, 1 (2017).
  * Gilyén  _et al._ [2022] A. Gilyén, S. Lloyd, I. Marvian, Y. Quek, and M. M. Wilde, Quantum algorithm for Petz recovery channels and pretty good measurements, Physical Review Letters 128, 220502 (2022).
  * Gilyén and Poremba [2022] A. Gilyén and A. Poremba, Improved quantum algorithms for fidelity estimation (2022), arXiv:2203.15993 .
  * Wang and Zhang [2023] Q. Wang and Z. Zhang, Quantum lower bounds by sample-to-query lifting (2023), arXiv:2308.01794 .
  * Wang and Zhang [2024] Q. Wang and Z. Zhang, Time-efficient quantum entropy estimator via samplizer (2024), arXiv:2401.09947 .
  * Cheng  _et al._ [2024] H.-C. Cheng, N. Datta, N. Liu, T. Nuradha, R. Salzmann, and M. M. Wilde, Sample complexity of quantum hypothesis testing (2024), arXiv:2403.17868 .
  * Lawson and Lim [2011] J. Lawson and Y. Lim, Monotonic properties of the least squares mean, Mathematische Annalen 351, 267 (2011).
  * Uhlmann [1976] A. Uhlmann, The “transition probability” in the state space of a $*$-algebra, Reports on Mathematical Physics 9, 273 (1976).
  * Wang  _et al._ [2023] Q. Wang, Z. Zhang, K. Chen, J. Guan, W. Fang, J. Liu, and M. Ying, Quantum algorithm for fidelity estimation, IEEE Transactions on Information Theory 69, 273 (2023).
  * Wang  _et al._ [2024] Q. Wang, J. Guan, J. Liu, Z. Zhang, and M. Ying, New quantum algorithms for computing quantum entropies and distances, IEEE Transactions on Information Theory (2024), 2203.13522 .
  * Khatri and Wilde [2024] S. Khatri and M. M. Wilde, Principles of quantum communication theory: A modern approach (2024), arXiv:2011.04672 [quant-ph] .
  * Ding  _et al._ [2023] D. Ding, S. Khatri, Y. Quek, P. W. Shor, X. Wang, and M. M. Wilde, Bounding the forward classical capacity of bipartite quantum channels, IEEE Transactions on Information Theory 69, 3034 (2023).
  * Arnold and Laub [1984] W. F. Arnold and A. J. Laub, Generalized eigenproblem algorithms and software for algebraic Riccati equations, Proceedings of the IEEE 72, 1746 (1984).
  * Estatico and Di Benedetto [2013] C. Estatico and F. Di Benedetto, Shift-invariant approximations of structured shift-variant blurring matrices, Numerical Algorithms 62, 615 (2013).
  * Arsigny  _et al._ [2007] V. Arsigny, P. Fillard, X. Pennec, and N. Ayache, Geometric means in a novel vector space structure on symmetric positive-definite matrices, SIAM Journal on Matrix Analysis and Applications 29, 328 (2007).
  * Chansangiam [2012] P. Chansangiam, Operator means and applications, in  _Linear Algebra_, edited by H. A. Yasser (IntechOpen, Rijeka, 2012) Chap. 8.
  * Moakher [2006] M. Moakher, On the averaging of symmetric positive-definite tensors, Journal of Elasticity 82, 273 (2006).
  * Barbaresco [2009] F. Barbaresco, New foundation of radar Doppler signal processing based on advanced differential geometry of symmetric spaces: Doppler matrix CFAR and radar application, in  _International Radar Conference_ , Vol. 82 (Bordeaux France, 2009).
  * Nielsen  _et al._ [2006] M. A. Nielsen, M. R. Dowling, M. Gu, and A. C. Doherty, Quantum computation as geometry, Science 311, 1133 (2006).
  * Fujii  _et al._ [2009] J. I. Fujii, M. Fujii, and R. Nakamoto, Riccati equation and positivity of operator matrices, Kyungpook Mathematical Journal 49, 595 (2009).
  * Brassard  _et al._ [2002] G. Brassard, P. Høyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, Quantum Computation and Information 305, 53 (2002).
  * Ando [1979] T. Ando, Concavity of certain maps on positive definite matrices and applications to Hadamard products, Linear Algebra and its Applications 26, 203 (1979).
  * Belovs [2019] A. Belovs, Quantum algorithms for classical probability distributions, in  _Proceedings of the 27th Annual European Symposium on Algorithms_ (2019) pp. 16:1–16:11.

##  Appendix A Proof of Lemma 5

This follows by observing that (9) is a matrix version of the quadratic
equation and by following an argument similar to what is well known as
completing the square. Consider that

$\displaystyle\left(Y-A^{-1}B\right)^{\dagger}A\left(Y-A^{-1}B\right)$
$\displaystyle=\left(Y^{\dagger}-\left(A^{-1}B\right)^{{\dagger}}\right)A\left(Y-A^{-1}B\right)$
(63)
$\displaystyle=\left(Y^{\dagger}-B^{{\dagger}}\left(A^{-1}\right)^{{\dagger}}\right)A\left(Y-A^{-1}B\right)$
(64)
$\displaystyle=\left(Y^{\dagger}-B^{{\dagger}}A^{-1}\right)A\left(Y-A^{-1}B\right)$
(65)
$\displaystyle=Y^{\dagger}AY-Y^{\dagger}AA^{-1}B-B^{{\dagger}}A^{-1}AY+B^{{\dagger}}A^{-1}AA^{-1}B$
(66)
$\displaystyle=Y^{\dagger}AY-Y^{\dagger}B-B^{{\dagger}}Y+B^{{\dagger}}A^{-1}B.$
(67)

Then

$Y^{\dagger}AY-B^{{\dagger}}Y-Y^{{\dagger}}B-C=\left(Y-A^{-1}B\right)^{\dagger}A\left(Y-A^{-1}B\right)-B^{{\dagger}}A^{-1}B-C,$
(68)

and so (9) is equivalent to

$X^{\dagger}AX=D,$ (69)

where

$\displaystyle X$ $\displaystyle=Y-A^{-1}B,$ (70) $\displaystyle D$
$\displaystyle=B^{{\dagger}}A^{-1}B+C.$ (71)

Observe that $D$ is positive definite because $B^{{\dagger}}A^{-1}B$ is
positive semi-definite and $C$ is positive definite. So this is a reduction to
the original simplified form of the algebraic Riccati equation in (10), which
we know from Lemma 4 has the following unique positive definite solution:

$\displaystyle X$ $\displaystyle=A^{-1}\\#D$ (72)
$\displaystyle=A^{-1}\\#\left(B^{{\dagger}}A^{-1}B+C\right).$ (73)

This implies that

$\displaystyle Y$ $\displaystyle=X+A^{-1}B$ (74)
$\displaystyle=A^{-1}\\#\left(B^{{\dagger}}A^{-1}B+C\right)+A^{-1}B$ (75)

is a solution of (9).

######  Remark 3.

Contrary to what is stated in the proof of [74, Corollary 4], the solution of
(9), under the assumptions stated in Lemma 5, is not unique. Indeed,

$Y=-\left(A^{-1}\\#\left(B^{{\dagger}}A^{-1}B+C\right)\right)+A^{-1}B$ (76)

is also a legitimate solution. In fact, the following is a matrix version of
the famous quadratic formula:

$A^{-1}B\pm\left(A^{-1}\\#\left(B^{{\dagger}}A^{-1}B+C\right)\right),$ (77)

for which the scalar version is
$\frac{b}{a}\pm\sqrt{\frac{1}{a}\left(\frac{b^{2}}{a}+c\right)}$ corresponding
to a solution of $ay^{2}-2by-c=0$ (stated after (9)), under the assumption
that $a,c>0$.

##  Appendix B Preliminary lemmas of the block-encoding formalism and other
useful results

Let us introduce several preliminary lemmas of the block-encoding formalism,
which enable us to implement various arithmetic operations on the block-
encoded matrices. The first lemma states that, given block-encodings of two
matrices, we can obtain a block-encoding of their product.

######  Lemma 24 (Product of block-encoded matrices [19, Lemma 30]).

If $U$ is an $(\alpha,a,\delta)$-block-encoding of $A$ and $V$ is a
$(\beta,b,\epsilon)$-block-encoding of $B$, then there is a unitary $W$ that
is an $(\alpha\beta,a+b,\alpha\epsilon+\beta\delta)$-block-encoding of $AB$,
and can be implemented by one query to $U$ and $V$.

Taking the linear combination of several block-encoded matrices is also useful
and is stated in the following lemma.

######  Lemma 25 (Linear combination of block-encoded matrices [19, Lemma
29]).

Let $m\in\mathbb{N}$ and $\beta>0$ be constant, and let
$\boldsymbol{x}=(x_{1},\ldots,x_{m})\in\mathbb{R}^{m}$ be a vector such that
$\lVert\boldsymbol{x}\rVert_{1}\leq\beta$. Suppose that each $U_{j}$ is a
$(1,a,\epsilon)$-block-encoding of $A_{j}$ for $j=1$ to $m$. Then there is a
unitary $U$ that is a $(1,a+\eta\log(1/\epsilon),2\beta^{-1}\epsilon)$-block-
encoding of $\beta^{-1}\sum_{j=1}^{m}x_{j}A_{j}$, where $\eta$ is some
constant, and $U$ can be implemented by one query to each $U_{j}$ and
$\operatorname{polylog}(1/\epsilon)$ gates.

To construct our quantum algorithms for matrix geometric means, we need to
deal with the non-linear terms in the matrix geometric means. The tool to be
used is quantum singular value transformation [19], which in our case is an
(approximate) polynomial transformation of the block-encoded matrix, as stated
in the following lemma.

######  Lemma 26 (Polynomial eigenvalue transformation [19, Theorem 31]).

Let $U$ be a $(1,a,\epsilon)$-block-encoding of a Hermitian matrix $A$. If
$\delta\geq 0$ and $q(x)\in\mathbb{R}[x]$ is a polynomial of degree $d$ such
that $\absolutevalue*{q(x)}\leq 1$ for $x\in[-1,1]$, then there is a unitary
$\widetilde{U}$ that is a $(1,a+2,4d\sqrt{\epsilon}+\delta)$-block-encoding of
$q(A)/2$, and can be implemented by $d$ queries to $U$ and $O((a+1)d)$ gates.
A description of such an implementation can be computed classically in time
$O\\!\left(\operatorname{poly}\\!\left(d,\log(1/\delta)\right)\right)$.

We also need two polynomial approximation results for applying Lemma 26 in our
scenario. The following two lemmas show low-degree polynomials for
approximating the negative and positive power functions, respectively.

######  Lemma 27 (Polynomial approximations of negative power functions [19,
Corollary 67 in the full version]).

Let $f(x)=\left(x/\delta\right)^{-c}/2$. For $\delta,\epsilon\in(0,1/2)$ and
$c>0$, there is a polynomial $q(x)$ of degree
$O\\!\left((c+1)\delta^{-1}\log(1/\epsilon)\right)$ such that

  * •

$\absolutevalue*{q(x)-f(x)}\leq\epsilon$ for $x\in[\delta,1]$;

  * •

$\absolutevalue*{q(x)}\leq 1$ for $x\in[-1,\delta)$.

######  Lemma 28 (Polynomial approximations of positive power functions [20,
Lemma 10]).

Let $f(x)=x^{c}/2$. For $\delta,\epsilon\in(0,1/2)$ and $c\in(0,1)$, there is
a polynomial $q(x)$ of degree $O\\!\left(\delta^{-1}\log(1/\epsilon)\right)$
such that

  * •

$\absolutevalue*{q(x)-f(x)}\leq\epsilon$ for $x\in[\delta,1]$;

  * •

$\absolutevalue*{q(x)}\leq 1$ for $x\in[-1,\delta)$.

In practice, how to encode the desired matrices into block-encodings and how
to extract useful (classical) information from the block-encodings are of
great concern. For the encoding, a typical scenario is that we are given
sparse oracle access to a sparse matrix, and we can construct a block-encoding
of the matrix, as stated in the following lemma.

######  Lemma 29 (Block-encoding of sparse matrices, [19, Lemma 48 in the full
version]).

Suppose $A\in\mathbb{C}^{N\times N}$ is an $s$-sparse matrix such that every
entry $A_{j,k}$ satisfies $\absolutevalue{A_{j,k}}\leq 1$. Suppose sparse
oracles $\mathcal{O}_{s}$ and $\mathcal{O}_{A}$ are given such that

$\displaystyle\mathcal{O}_{s}\ket{j}\ket{k}$
$\displaystyle=\ket{j}\ket{l_{j,k}},$ (78)
$\displaystyle\mathcal{O}_{A}\ket{j}\ket{k}\ket{0}$
$\displaystyle=\ket{j}\ket{k}\ket{A_{j,k}},$ (79)

where $l_{j,k}$ denotes the column index of the $k$-th non-zero entry in the
$j$-th row. Here, we assume that the exact value of the entry $A_{j,k}$ is
given in a binary representation. Then, we can implement a quantum circuit
that is an $(s,\log_{2}N+3,\epsilon)$-block-encoding of $A$, using two queries
to $\mathcal{O}_{s}$, two queries to $\mathcal{O}_{A}$, and $O(\log
N+\log^{2.5}(s/\epsilon))$ one- and two-qubit quantum gates.

Another useful case of encoding commonly considered is that we are given
purified access to a density operator, and we can construct a block-encoding
of the density operator.

######  Lemma 30 (Block-encoding of density operators [21, Lemma 7], [19,
Lemma 25]).

Let $\rho$ be an $n$-qubit density operator, and let $V_{\rho}$ be an
$(n+a)$-qubit unitary that prepares a purification of $\rho$ such that
$\tr_{a}(V_{\rho}\ket{0}_{n+a}\\!\bra{0}V_{\rho}^{\dagger})=\rho$. Then there
is a $(2n+a)$-qubit unitary $\widetilde{V}$ that is a $(1,n+a,0)$-block-
encoding of $\rho$, and it can be implemented by one query to $V$ and $O(n)$
gates.

To extract classical information from the block-encodings, one needs to
perform quantum measurements. The Hadamard test is a useful and efficient way
to estimate the expectation value of a quantum observable on a given quantum
state. The following lemma shows a Hadamard test for block-encodings.

######  Lemma 31 (Hadamard test for block-encodings, [57, Lemma 9]).

Suppose that $U$ is a unitary operator that is a $\lparen 1,a,0\rparen$-block-
encoding of an $n$-qubit operator $A$. Then, there is a quantum circuit that
outputs $0$ with probability $\frac{1+\Re{\Tr\lparen A\rho\rparen}}{2}$, using
one query to $U$ and one sample of the mixed quantum state $\rho$.

The success probability of extracting classical information from a block-
encoding often depends on the scaling factor of the block-encoding. Robust
oblivious quantum amplitude amplification is a technique for amplifying the
scaling factor of block-encoding, as stated in the following lemma.

######  Lemma 32 (Robust oblivious amplitude amplification [58, Corollary
2.8], adapted from [19, Theorem 15]).

Let $U$ be a $(1,a,\epsilon)$-block-encoding of $A/\alpha$ such that
$\alpha>1$ and $\lVert A\rVert\leq 1$. Then there is a unitary $\widetilde{U}$
that is a $(1,a+2,8\alpha\epsilon)$-block-encoding of $A$ and can be
implemented by $O(\alpha)$ queries to $U$ and $O(a\alpha)$ gates.

Apart from the block-encoding formalism, let us introduce other useful
results. Quantum amplitude estimation allows one to estimate the amplitude of
a specific component of a quantum state, stated as follows.

######  Lemma 33 (Quantum amplitude estimation [75, Theorem 12]).

Suppose that unitary operator $U$ is given by

$U\ket{0}\ket{0}=\sqrt{p}\ket{0}\ket{\phi_{0}}+\sqrt{1-p}\ket{1}\ket{\phi_{1}},$
(80)

where $\ket{\phi_{0}}$ and $\ket{\phi_{1}}$ are normalized pure quantum
states, and $p\in\left[0,1\right]$. Then, we can obtain an estimate
$\widetilde{p}$ of $p$ such that

$\absolutevalue*{\widetilde{p}-p}\leq\frac{2\pi\sqrt{p(1-p)}}{M}+\frac{\pi^{2}}{M^{2}}$
(81)

with probability $\geq 8/\pi^{2}$ using $O\left\lparen M\right\rparen$ queries
to $U$. In particular, if we take $M=\Theta\lparen 1/\delta\rparen$, then
$\widetilde{p}$ is a $\delta$-estimate of $p$ with high probability.

Finally, the following two lemmas give relevant bounds on the condition number
of matrices, which will be useful in the complexity analysis of our
algorithms.

######  Lemma 34.

Let $A,B>0$ be two positive definite matrices such that
$\left\|A\right\|=\left\|B\right\|=1$. Then
$\kappa_{A+B}^{-1}\geq\kappa_{A}^{-1}+\kappa_{B}^{-1}$.

######  Lemma 35.

Let $A>0$ be a positive definite matrix, and let $B$ be a matrix of full rank,
such that $\left\|A\right\|=\left\|B\right\|=1$. Then
$\kappa_{B^{\dagger}AB}\leq\kappa_{A}\kappa_{B^{\dagger}B}$.

###### Proof.

First note that $B^{\dagger}B>0$ and $B^{\dagger}CB>0$ for every $C>0$. Since
$A\geq I/\kappa_{A}$, it follows that

$B^{\dagger}AB\geq\kappa_{A}^{-1}B^{\dagger}IB\geq\kappa_{A}^{-1}\kappa_{B^{\dagger}B}^{-1}I.$
(82)

It immediately follows that
$\kappa_{B^{\dagger}AB}^{-1}\geq\kappa_{A}^{-1}\kappa_{B^{\dagger}B}^{-1}$. ∎

##  Appendix C Proof of Lemma 7

In this appendix we prove Lemma 7. Let us first prove the following lemma,
which shows that we can implement a block-encoding of the matrix $\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}$.

######  Lemma 36.

Suppose that $U_{A},U_{C}$ are $(1,a,0)$-block-encodings of matrices $A,C$,
respectively, such that $A\geq I/\kappa_{A}$ and $C\geq I/\kappa_{C}$. For
$\epsilon\in(0,1/2)$, one can implement a $(1,3a+12,\epsilon)$-block-encoding
of $\kappa_{A}^{-1/p}\gamma_{p}^{-1}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}$ for any fixed real $p\neq 0$, where

$\gamma_{p}=\begin{cases}1&p>0,\\\
\kappa_{A}^{-1/p}\kappa_{C}^{-1/p}&p<0,\end{cases}$ (83)

using

  * •

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$ queries to
$U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$;

  * •

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates; and

  * •

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

###### Proof.

We first consider the case $p>0$. Let us construct
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$, a
block-encoding of $\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}$, step by step as follows. Along the way,
we also analyse the resources for each step. In the remainder of the paper, we
use the notation $\widetilde{O}_{a_{1},\ldots,a_{n}}(f)$ to denote
$O(f\operatorname{polylog}(b_{1},\ldots,b_{n}))$, where $a_{i},b_{i}$ are
parameters, $f$ is a function, and $b_{i}=a_{i}+a_{i}^{-1}$. Similarly, we use
$\widetilde{\Omega}_{a_{1},\ldots,a_{n}}(f)$ to denote
$\Omega(f/\operatorname{polylog}(b_{1},\ldots,b_{n}))$. In context without
ambiguity, we just omit the subscripts as usual.

  1. 1.

$U_{A}\rightarrow U_{(\kappa_{A}A)^{-1/2}}$:

     * •

Construction:

       1. (a)

Taking $c=1/2$, $\delta=\kappa_{A}^{-1}$, and $\epsilon=\epsilon_{1}$ in Lemma
27, we have a polynomial $q_{1}(x)$ of degree
$d_{1}=O\\!\left(\kappa_{A}\log(1/\epsilon_{1})\right)$ that approximates
$(\kappa_{A}x)^{-1/2}/2$.

       2. (b)

Taking $U=U_{A}$, $q=q_{1}(x)$, $\epsilon=0$, and $\delta=\epsilon_{1}$ in
Lemma 26, we have $U_{\left(\kappa_{A}A\right)^{-1/2}/4}$, a
$(1,a+2,\epsilon_{1})$-block-encoding of $q_{1}(A)/2$, which is therefore a
$(1,a+2,2\epsilon_{1})$-block-encoding of $\left(\kappa_{A}A\right)^{-1/2}/4$.

       3. (c)

Taking $U=U_{\left(\kappa_{A}A\right)^{-1/2}/4}$, $\alpha=4$ in Lemma 32, we
obtain $U_{\left(\kappa_{A}A\right)^{-1/2}}$, a
$(1,a+4,2^{6}\epsilon_{1})$-block-encoding of
$\left(\kappa_{A}A\right)^{-1/2}$.

     * •

Resources: $O(d_{1})$ queries to $U_{A}$, and $O(ad_{1})$ gates, and
$\operatorname{poly}(d_{1},\log(1/\epsilon_{1}))$ classical time.

  2. 2.

$U_{C},U_{(\kappa_{A}A)^{-1/2}}\rightarrow
U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}$:

     * •

Construction: by Lemma 24, given $U_{C}$ and $U_{(\kappa_{A}A)^{-1/2}}$, we
have $U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}$, a
$(1,3a+8,2^{7}\epsilon_{1})$-block-encoding of
$\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}$.

     * •

Resources: $O(1)$ queries to $U_{C}$ and $U_{(\kappa_{A}A)^{-1/2}}$.

  3. 3.

$U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}\rightarrow
U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$:

     * •

Construction:

       1. (a)

Taking $c=1/p$,
$\delta=\kappa_{A}^{-1}\kappa_{C}^{-1}\leq\kappa_{A}^{-1}\kappa_{A^{-1/2}CA^{-1/2}}^{-1}$
(by Lemma 35 and noting $\kappa_{A^{-1}}\leq 1$), and $\epsilon=\epsilon_{2}$
in Lemma 28, we have a polynomial $q_{2}(x)$ of degree
$d_{2}=O\\!\left(\kappa_{A}\kappa_{C}\log(1/\epsilon_{2})\right)$ that
approximates $x^{1/p}/2$.

       2. (b)

Taking $U=U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}$, $q=q_{2}(x)$,
$\epsilon=2^{7}\epsilon_{1}$, and $\delta=\epsilon_{2}$ in Lemma 26, we have
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}/4}$, a
$(1,3a+10,2^{5.5}d_{2}\sqrt{\epsilon_{1}}+\epsilon_{2})$-block-encoding of
$q_{2}\left(\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}\right)/2$, which is therefore a
$(1,3a+10,2^{5.5}d_{2}\sqrt{\epsilon_{1}}+2\epsilon_{2})$-block-encoding of
$\left(\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}\right)^{1/p}/4$.

       3. (c)

Taking $U=U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}/4}$, $\alpha=4$ in Lemma 32, we obtain
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$, a
$(1,3a+12,2^{10.5}d_{2}\sqrt{\epsilon_{1}}+2^{6}\epsilon_{2})$-block-encoding
of $\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}$.

     * •

Resources: $O(d_{2})$ queries to $U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}$, and
$O(ad_{2})$ gates, and $\operatorname{poly}(d_{2},\log(1/\epsilon_{2}))$
classical time.

Let us take

  * •

$\epsilon_{2}=2^{-7}\epsilon=\Omega(\epsilon)$.

  * •

$\epsilon_{1}=2^{-23}\epsilon
d_{2}^{-2}=\Omega\left(\kappa_{A}^{-2}\kappa_{C}^{-2}\epsilon\log^{-2}\left(1/\epsilon_{2}\right)\right)=\widetilde{\Omega}\left(\kappa_{A}^{-2}\kappa_{C}^{-2}\epsilon\right)$.

Then, we can calculate the final approximation error in
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$ as a
block-encoding of $\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}$:

$2^{10.5}d_{2}\sqrt{\epsilon_{1}}+2^{6}\epsilon_{2}\leq\epsilon.$ (84)

The number of ancilla qubits for constructing
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$ is
$3a+12$.

Finally, let us calculate the complexities of each step.

  1. 1.

$U_{A}\rightarrow U_{(\kappa_{A}A)^{-1/2}}$:

$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}\log\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}\log\left(1/\epsilon\right)\right)$
gates,

and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\log\left(1/\epsilon\right)\right)$
classical time.

  2. 2.

$U_{C},U_{(\kappa_{A}A)^{-1/2}}\rightarrow
U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}$:

$O(1)$ queries to $U_{C}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}\log\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}\log\left(1/\epsilon\right)\right)$
gates, and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\log\left(1/\epsilon\right)\right)$
classical time.

  3. 3.

$U_{\kappa_{A}^{-1}A^{-1/2}CA^{-1/2}}\rightarrow
U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$:

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$ queries to
$U_{C}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates, and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon\right)\right)$
classical time.

For the case $p<0$, the analysis is the same except that in Step 3a, we can
use Lemma 27 instead of Lemma 28. This only incurs an additional scaling
factor $\kappa_{A}^{1/p}\kappa_{C}^{1/p}$ into the final block-encoded matrix,
without significantly changing the complexity. ∎

Now we are ready to prove Lemma 7, which gives an implementation of a block-
encoding of the weighted matrix geometric mean in Eq. (2).

###### Proof of Lemma 7.

We first consider the case $p>0$. Let us construct $U_{\kappa_{A}^{-1/p}Y}$, a
$(1,5a+20,\epsilon)$-block-encoding of $\kappa_{A}^{-1/p}Y$, step by step as
follows, where

$Y=A\\#_{1/p}C=A^{1/2}\left(A^{-1/2}CA^{-1/2}\right)^{1/p}A^{1/2}.$ (85)

Along the way, we also analyse the resources for each step.

  1. 1.

$U_{A}\rightarrow U_{A^{1/2}}$:

     * •

Construction:

       1. (a)

Taking $c=1/2$, $\delta=\kappa_{A}^{-1}$, and $\epsilon=\epsilon_{1}$ in Lemma
28, we have a polynomial $q_{1}(x)$ of degree
$d_{1}=O\\!\left(\kappa_{A}\log(1/\epsilon_{1})\right)$ that approximates
$x^{1/2}/2$.

       2. (b)

Taking $U=U_{A}$, $q=q_{1}(x)$, $\epsilon=0$, and $\delta=\epsilon_{1}$ in
Lemma 26, we have $U_{A^{1/2}/4}$, a $(1,a+2,\epsilon_{1})$-block-encoding of
$q_{1}(A)/2$, which is therefore a $(1,a+2,2\epsilon_{1})$-block-encoding of
$A^{1/2}/4$.

       3. (c)

Taking $U=U_{A^{1/2}/4}$, $\alpha=4$ in Lemma 32, we obtain $U_{A^{1/2}}$, a
$(1,a+4,2^{6}\epsilon_{1})$-block-encoding of $A^{1/2}$.

     * •

Resources: $O(d_{1})$ queries to $U_{A}$, and $O(ad_{1})$ gates, and
$\operatorname{poly}(d_{1},\log(1/\epsilon_{1}))$ classical time.

  2. 2.

$U_{A},U_{C}\rightarrow U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$:

     * •

Construction:

Taking $\epsilon=\epsilon_{2}$ in Lemma 36, we can construct
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$, a
$(1,3a+12,\epsilon_{2})$-block-encoding of $U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$.

     * •

Resources:

According to Lemma 36, the resources for the above construction are:

       * –

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon_{2}\right)\right)$ queries
to $U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
queries to $U_{A}$;

       * –

$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
gates; and

       * –

$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon_{2}\right)\right)$
classical time.

  3. 3.

$U_{A^{1/2}},U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}\rightarrow U_{\kappa_{A}^{-1/p}Y}$:

     * •

Construction: by Lemma 24, given $U_{A^{1/2}}$ and
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$, we
have $U_{\kappa_{A}^{-1/2}Y}$, a
$(1,5a+20,2^{7}\epsilon_{1}+\epsilon_{2})$-block-encoding of
$\kappa_{A}^{-1/2}Y$.

     * •

Resources: $O(1)$ queries to $U_{A^{1/2}}$ and
$U_{\kappa_{A}^{-1/p}\left\lparen A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$.

Let us take $\epsilon_{1}=2^{-8}\epsilon$ and $\epsilon_{2}=2^{-1}\epsilon$.
Then the final approximation error in $U_{\kappa_{A}^{-1/p}Y}$ as a block-
encoding of $\kappa_{A}^{-1/p}Y$ is

$2^{7}\epsilon_{1}+\epsilon_{2}\leq\epsilon.$ (86)

The number of ancilla qubits for constructing $\kappa_{A}^{-1/p}Y$ is $5a+20$.
Finally, let us calculate the complexities of each step.

  1. 1.

$U_{A}\rightarrow U_{A^{1/2}}$:

$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}\log\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}\log\left(1/\epsilon\right)\right)$
gates,

and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\log\left(1/\epsilon\right)\right)$
classical time.

  2. 2.

$U_{A},U_{C}\rightarrow U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}$:

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon_{2}\right)\right)$ queries
to $U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
queries to $U_{A}$,
$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
gates, and
$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon_{2}\right)\right)$
classical time.

  3. 3.

$U_{A^{1/2}},U_{\kappa_{A}^{-1/p}\left\lparen
A^{-1/2}CA^{-1/2}\right\rparen^{1/p}}\rightarrow U_{\kappa_{A}^{-1/p}Y}$:

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon_{2}\right)\right)$ queries
to $U_{C}$,
$\widetilde{O}\\!\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
queries to $U_{A}$,
$\widetilde{O}\\!\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon_{2}\right)\right)$
gates, and
$\operatorname{poly}\\!\left(\kappa_{A},\kappa_{C},\log\left(1/\epsilon_{2}\right)\right)$
classical time.

By Definition 6, $U_{\kappa_{A}^{-1/p}Y}$ is also a
$(\kappa_{A}^{1/p},5a+20,\kappa_{A}^{1/p}\epsilon)$-block-encoding of $Y$.
Replacing the precision parameter immediately yields the results for the case
$p>0$ in Lemma 7. For the case $p<0$, the analysis is similar and is omitted.
∎

##  Appendix D Proof of Lemma 8

In this appendix we prove Lemma 8.

###### Proof of Lemma 8.

Let us construct $U_{\kappa_{A}^{-1}Y}$, a $(1,5a+20,\epsilon)$-block-encoding
of $\kappa_{A}^{-1}Y$, step by step as follows, where

$Y=A^{-1}\\#C=A^{-1/2}\left(A^{1/2}CA^{1/2}\right)^{1/2}A^{-1/2}.$ (87)

Along the way, we also analyse the resources for each step.

  1. 1.

$U_{A}\rightarrow U_{A^{1/2}}$:

     * •

Construction:

       1. (a)

Taking $c=1/2$, $\delta=\kappa_{A}^{-1}$, and $\epsilon=\epsilon_{1}$ in Lemma
28, we have a polynomial $q_{1}(x)$ of degree
$d_{1}=O\\!\left(\kappa_{A}\log(1/\epsilon_{1})\right)$ that approximates
$x^{1/2}/2$.

       2. (b)

Taking $U=U_{A}$, $q=q_{1}(x)$, $\epsilon=0$, and $\delta=\epsilon_{1}$ in
Lemma 26, we have $U_{A^{1/2}/4}$, a $(1,a+2,\epsilon_{1})$-block-encoding of
$q_{1}(A)/2$, which is therefore a $(1,a+2,2\epsilon_{1})$-block-encoding of
$A^{1/2}/4$.

       3. (c)

Taking $U=U_{A^{1/2}/4}$, $\alpha=4$ in Lemma 32, we obtain $U_{A^{1/2}}$, a
$(1,a+4,2^{6}\epsilon_{1})$-block-encoding of $A^{1/2}$.

     * •

Resources: $O(d_{1})$ queries to $U_{A}$, and $O(ad_{1})$ gates, and
$\operatorname{poly}(d_{1},\log(1/\epsilon_{1}))$ classical time.

  2. 2.

$U_{C},U_{A^{1/2}}\rightarrow U_{A^{1/2}CA^{1/2}}$:

     * •

Construction: by Lemma 24, given $U_{C}$ and $U_{A^{1/2}}$, we have
$U_{A^{1/2}CA^{1/2}}$, a $(1,3a+8,2^{7}\epsilon_{1})$-block-encoding of
$A^{1/2}CA^{1/2}$.

     * •

Resources: $O(1)$ queries to $U_{C}$ and $U_{A^{1/2}}$.

  3. 3.

$U_{A^{1/2}CA^{1/2}}\rightarrow U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}$:

     * •

Construction:

       1. (a)

Taking $c=1/2$,
$\delta=\kappa_{A}^{-1}\kappa_{C}^{-1}\leq\kappa_{A^{1/2}CA^{1/2}}^{-1}$ (by
Lemma 35), and $\epsilon=\epsilon_{2}$ in Lemma 28, we have a polynomial
$q_{2}(x)$ of degree
$d_{2}=O\\!\left(\kappa_{A}\kappa_{C}\log(1/\epsilon_{2})\right)$ that
approximates $x^{1/2}/2$.

       2. (b)

Taking $U=U_{A^{1/2}CA^{1/2}}$, $q=q_{2}(x)$, $\epsilon=2^{7}\epsilon_{1}$,
and $\delta=\epsilon_{2}$ in Lemma 26, we have
$U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}/4}$, a
$(1,3a+10,2^{5.5}d_{2}\sqrt{\epsilon_{1}}+\epsilon_{2})$-block-encoding of
$q_{2}\left(A^{1/2}CA^{1/2}\right)/2$, which is therefore a
$(1,3a+10,2^{5.5}d_{2}\sqrt{\epsilon_{1}}+2\epsilon_{2})$-block-encoding of
$\left(A^{1/2}CA^{1/2}\right)^{1/2}/4$.

       3. (c)

Taking $U=U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}/4}$, $\alpha=4$ in Lemma 32,
we obtain $U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}$, a
$(1,3a+12,2^{10.5}d_{2}\sqrt{\epsilon_{1}}+2^{6}\epsilon_{2})$-block-encoding
of $\left(A^{1/2}CA^{1/2}\right)^{1/2}$.

     * •

Resources: $O(d_{2})$ queries to $U_{A^{1/2}CA^{1/2}}$, and $O(ad_{2})$ gates,
and $\operatorname{poly}(d_{2},\log(1/\epsilon_{2}))$ classical time.

  4. 4.

$U_{A}\rightarrow U_{(\kappa_{A}A)^{-1/2}}$:

     * •

Construction:

       1. (a)

Taking $c=1/2$, $\delta=\kappa_{A}^{-1}$, and $\epsilon=\epsilon_{3}$ in Lemma
27, we have a polynomial $q_{3}(x)$ of degree
$d_{3}=O\\!\left(\kappa_{A}\log(1/\epsilon_{3})\right)$ that approximates
$(\kappa_{A}x)^{-1/2}/2$.

       2. (b)

Taking $U=U_{A}$, $q=q_{3}(x)$, $\epsilon=0$, and $\delta=\epsilon_{3}$ in
Lemma 26, we have $U_{\left(\kappa_{A}A\right)^{-1/2}/4}$, a
$(1,a+2,\epsilon_{3})$-block-encoding of $q_{3}(A)/2$, which is therefore a
$(1,a+2,2\epsilon_{3})$-block-encoding of $\left(\kappa_{A}A\right)^{-1/2}/4$.

       3. (c)

Taking $U=U_{\left(\kappa_{A}A\right)^{-1/2}/4}$, $\alpha=4$ in Lemma 32, we
obtain $U_{\left(\kappa_{A}A\right)^{-1/2}}$, a
$(1,a+4,2^{6}\epsilon_{3})$-block-encoding of
$\left(\kappa_{A}A\right)^{-1/2}$.

     * •

Resources: $O(d_{3})$ queries to $U_{A}$, and $O(ad_{3})$ gates, and
$\operatorname{poly}(d_{3},\log(1/\epsilon_{3}))$ classical time.

  5. 5.

$U_{(\kappa_{A}A)^{-1/2}},U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}\rightarrow
U_{\kappa_{A}^{-1}Y}$:

     * •

Construction: by Lemma 24, given $U_{(\kappa_{A}A)^{-1/2}}$ and
$U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}$, we have $U_{\kappa_{A}^{-1}Y}$, a
$(1,5a+20,2^{10.5}d_{2}\sqrt{\epsilon_{1}}+2^{6}\epsilon_{2}+2^{7}\epsilon_{3})$-block-
encoding of $\kappa_{A}^{-1}Y$.

     * •

Resources: $O(1)$ queries to $U_{(\kappa_{A}A)^{-1/2}}$ and
$U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}$.

Let us take

  * •

$\epsilon_{3}=2^{-9}\epsilon=\Omega(\epsilon)$.

  * •

$\epsilon_{2}=2^{-8}\epsilon=\Omega(\epsilon)$.

  * •

$\epsilon_{1}=2^{-23}\epsilon
d_{2}^{-2}=\Omega\left(\kappa_{A}^{-2}\kappa_{C}^{-2}\epsilon\log^{-2}\left(1/\epsilon_{2}\right)\right)=\widetilde{\Omega}\left(\kappa_{A}^{-2}\kappa_{C}^{-2}\epsilon\right)$.

Then, the final approximation error in $U_{\kappa_{A}^{-1}Y}$ as a block-
encoding of $\kappa_{A}^{-1}Y$ is

$2^{10.5}d_{2}\sqrt{\epsilon_{1}}+2^{6}\epsilon_{2}+2^{7}\epsilon_{3}\leq\epsilon.$
(88)

The number of ancilla qubits for constructing $U_{\kappa_{A}^{-1}Y}$ is
$5a+20$.

Finally, let us calculate the complexities of each step.

  1. 1.

$U_{A}\rightarrow U_{A^{1/2}}$:

$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}\log\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}\log\left(1/\epsilon\right)\right)$
gates,

and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\log\left(1/\epsilon\right)\right)$
classical time.

  2. 2.

$U_{C},U_{A^{1/2}}\rightarrow U_{A^{1/2}CA^{1/2}}$:

$O(1)$ queries to $U_{C}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}\log\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}\log\left(1/\epsilon\right)\right)$
gates, and
$\operatorname{poly}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A},\log\left(1/\epsilon\right)\right)$
classical time.

  3. 3.

$U_{A^{1/2}CA^{1/2}}\rightarrow U_{\left(A^{1/2}CA^{1/2}\right)^{1/2}}$:

$O\\!\left(\kappa_{A}\kappa_{C}\log\left(1/\epsilon\right)\right)$ queries to
$U_{C}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
queries to $U_{A}$,
$\widetilde{O}_{\kappa_{A},\kappa_{C},\epsilon}\left(a\kappa_{A}^{2}\kappa_{C}\log^{2}\\!\left(1/\epsilon\right)\right)$
gates, and