# Recursive Quantum Eigenvalue/Singular-Value Transformation:
Analytic Construction of Matrix Sign Function by Newton Iteration

Kaoru Mizuta<EMAIL_ADDRESS>RIKEN Center for Quantum Computing
(RQC), Hirosawa 2-1, Wako, Saitama 351-0198, Japan  Department of Applied
Physics, The University of Tokyo, Hongo 7-3-1, Bunkyo, Tokyo 113-8656, Japan
Keisuke Fujii  Graduate School of Engineering Science, Osaka University, 1-3
Machikaneyama, Toyonaka, Osaka 560-8531, Japan.  Center for Quantum
Information and Quantum Biology, Osaka University, Japan.  RIKEN Center for
Quantum Computing (RQC), Hirosawa 2-1, Wako, Saitama 351-0198, Japan  Fujitsu
Quantum Computing Joint Research Division at QIQB, Osaka University, 1-2
Machikaneyama, Toyonaka 560-0043, Japan

###### Abstract

Quantum eigenvalue transformation (QET) and its generalization, quantum
singular value transformation (QSVT), are versatile quantum algorithms that
allow us to apply broad matrix functions to quantum states, which cover many
of significant quantum algorithms such as Hamiltonian simulation. However,
finding a parameter set which realizes preferable matrix functions in these
techniques is difficult for large-scale quantum systems: there is no
analytical result other than trivial cases as far as we know and we often
suffer also from numerical instability. In this Letter, we propose recursive
QET or QSVT (r-QET or r-QSVT), in which we can execute complicated matrix
functions by recursively organizing block-encoding by low-degree QET or QSVT.
Owing to the simplicity of recursive relations, it works only with a few
parameters with exactly determining the parameters, while its iteration
results in complicated matrix functions. In particular, by exploiting the
recursive relation of Newton iteration, we construct the matrix sign function,
which can be applied for eigenstate filtering for example, in a tractable way.
We show that an analytically-obtained parameter set composed of only $8$
different values is sufficient for executing QET of the matrix sign function
with an arbitrarily small error $\varepsilon$. Our protocol will serve as an
alternative protocol for constructing QET or QSVT for some useful matrix
functions without numerical instability.

Introduction.— Quantum Eigenvalue Transformation (QET) is a versatile quantum
algorithm which enables to apply various matrix functions [1]. For a hermitian
matrix of interest $A$, QET executes parallel processing of its eigenvalues
and thereby allows us to apply broad matrix polynomial functions
$\sum_{n}c_{n}A^{n}$ to arbitrary quantum states. With its generalization to
general matrices, called quantum singular value transformation (QSVT), it
covers various today’s important quantum algorithms such as Hamiltonian
simulation [2, 3] and search algorithms [4], by properly constructing
polynomial approximations. Not only it provides unified understanding of
quantum algorithms [5], but also it can serve more efficient alternative
algorithms for various purposes.

QET can be executed by repetition of parameterized unitary gates on ancilla
qubits and unitary gates embedding the target matrix, called block-encoding.
While tunability of the parameters ensures realization of broad functions by
quantum signal processing (QSP) [2], we must accurately determine a proper
parameter set for a desired function. Although finding the parameters for
degree-$q$ polynomials within an error $\varepsilon$ can be executed by
$\mathrm{poly}\left(q,\log(1/\varepsilon)\right)$-time classical computation,
its numerical instability has become one of the central problems for accurate
implementation of QET/QSVT. In fact, several numerical algorithms trying to
solve this instability, such as parameter determination by optimization, have
recently been developed [6, 7, 8, 9]. By contrast, there are only a few
results on analytical parameter determination. As far as we know, they are
limited to trivial cases for Chebyshev polynomials, which are useful for
Grover’s search algorithm [4, 10, 11].

In this Letter, we propose recursive QET/QSVT (r-QET/r-QSVT) that can
potentially determine all the parameters in an analytical or numerically-much-
cheaper way. In these protocols, we recursively organize block-encoding by
low-degree QET/QSVT so that it can reproduce recursive relations of matrix
functions, and then obtain complicated matrix functions by iteration. For
instance, we can exploit Newton iteration for matrix functions as the
recursive relation [12]. Then, with a sufficient number of iterations for its
convergence, r-QET/r-QSVT enables us to organize nontrivial matrix functions
only with a smaller number of parameters that can be easily determined. As a
prominent consequence, we obtain a tractable implementation of matrix sign
functions with arbitrarily small error. Using Newton-Schulz iteration of Padé
family [13], we analytically obtain a parameter set that accurately builds QET
of matrix sign functions. Furthermore, the parameter set has constant unique
values which do not depend on either an allowable error $\varepsilon$ or any
parameter of the matrix. Although our construction expenses the computational
cost compared to the optimal protocol [14, 15] due to the strong limitation on
the parameters, it suffers from no numerical instability and even can
overwhelm the optimal one when we take recovery of coherent errors into
account [16]. With various recursive constructions of matrix functions such as
Newton iteration [12] and logistic map [17, 18], r-QET/r-QSVT will give a
promising candidate for executing complicated operations on quantum computers
in accurate and stable ways.

Quantum Eigenvalue Transformation (QET).— Throughout the main text, we
concentrate on QET and thus r-QET for hermitian matrices for simplicity (See
Supplemental Materials S2 for QSVT, whose discussion is completely parallel).
We begin with briefly introducing QET here. Let a hermitian matrix $A$ have
spectral decomposition $A=\sum_{a}a\ket{a}\bra{a}$ ($a\in\mathbb{R}$) on a
finite-dimensional Hilbert space $\mathcal{H}$. Block-encoding of $A$ is
defined by a unitary gate $O_{A}$ satisfying

$\braket{0}{O_{A}}{0}_{b}=\frac{A}{\alpha},\quad\alpha>0.$ (1)

Here, $\ket{0}_{b}$ denotes a reference state in an ancillary Hilbert space
$\mathcal{H}_{b}$. We set $\alpha=1$ by the renormalization $A\to A/\alpha$
below. Construction of block-encoding is known for a linear combination of
unitaries, a sparse-access matrix, and so on [3].

Combining parametrized unitary operations on the ancilla system,

$R_{\phi}=e^{i\phi(2\ket{0}\bra{0}_{b}-I_{b})}\otimes I,$ (2)

we define a degree-$q$ QET operator by

$\displaystyle\mathrm{QET}[A,\vec{\phi}]=$
$\displaystyle\quad\begin{cases}R_{\phi_{1}}O_{A}\prod_{i=1}^{(q-1)/2}\left[R_{\phi_{2i}}O_{A}^{\dagger}R_{\phi_{2i+1}}O_{A}\right],&\text{$q$:
odd},\\\
\prod_{i=1}^{q/2}\left[R_{\phi_{2i-1}}O_{A}^{\dagger}R_{\phi_{2i}}O_{A}\right],&\text{$q$:
even}.\end{cases}$ (3)

By properly tuning the parameter set $\vec{\phi}\in\mathbb{R}^{q}$, it can
realize various polynomial functions $f(x)=\sum_{n}c_{n}x^{n}$
($c_{n}\in\mathbb{C}$) as

$\braket{0}{\mathrm{QET}[A,\vec{\phi}]}{0}_{b}=f(A)=\sum_{n}c_{n}A^{n}.$ (4)

It is proven that there exists a parameter set $\vec{\phi}$ if and only if
$f(x)$ satisfies all the following conditions [1]:

  1. (i)

$f(x)$ has a degree at most $q$ and a parity $(-1)^{q}$.

  2. (ii)

$|f(x)|\leq 1$ for any $x\in[-1,1]$ and $|f(x)|\geq 1$ for any
$x\in(-\infty,-1]\cup[1,\infty)$.

  3. (iii)

(If $q$ is even) $|f(ix)f^{\ast}(ix)|\geq 1$ for any $x\in\mathbb{R}$, where
$f^{\ast}(x)$ is defined by $f^{\ast}(x)=\sum_{n}c_{n}^{\ast}x^{n}$.

By $\order{1}$ controlled operations $\mathrm{QET}[A,\vec{\phi}]$, generic
renormalized matrix functions $f(x)$ with $|f(x)|\leq 1/4$
(${}^{\forall}x\in[-1,1]$) are also realizable.

For a desired function $f(x)$ satisfying (i)-(iii), how can we find a proper
parameter set $\vec{\phi}$ ? As far as we know, the analytical result is
restricted to the Chebyshev polynomials $f(x)=T_{n}(x)$ with the trivial
angles $\vec{\phi}=((q-1)\pi/2,-\pi/2,-\pi/2,\ldots,-\pi/2)$, which can be
utilized for Grover’s search and its family such as amplitude amplification
[4]. In general, it requires finding all the roots of $1-f(x)f^{\ast}(x)$,
which is a degree-$q$ polynomial in $x^{2}$, and iteratively decomposing
$\mathrm{QET}[A,\vec{\phi}]$ into lower-degree QET operators [1]. However, for
useful functions such as $e^{-itA}$ (Hamiltonian simulation [2, 3]),
$A^{-1}/2\kappa$ (quantum linear system problem, QLSP [19]), and
$\mathrm{sign}(A)$ (eigenstate filtering [15]), the typical degree $q$ is
quite large as $q\in\mathrm{poly}\left(N,\log(1/\varepsilon)\right)$ depending
on the system size $N$ and the allowable error $\varepsilon$. Thus, it is
difficult to accurately compute $\vec{\phi}$ generally having
$\mathrm{poly}\left(N,\log(1/\varepsilon)\right)$ different values in a
numerically-stable way, although it can be done by
$\mathrm{poly}\left(N,\log(1/\varepsilon)\right)$-time classical calculation.
This instability has been partially resolved by refining root-finding problems
and decomposition into lower degrees [6, 7], or employing optimization [8, 9].
They have numerically succeeded up to $q\sim 10^{4}$ within $10^{2}\sim
10^{4}$ seconds.

Recursive QET (r-QET) and Newton iteration.— Here, we propose the protocol
named recursive QET (r-QET), and provide the formulation combined with Newton
iteration for matrix functions [12]. It aims to implement complicated matrix
functions by QET with keeping tractability of the parameter set $\vec{\phi}$
in terms of analytical or numerically-cheap computation. Our strategy is to
employ recursive relations: while each step operation is executed in a simple
way, its repetition forms rather complicated functions.

Suppose that we want to execute complicated matrix functions of $A$ with its
block-encoding $O_{A}$. Based on the fact that QET generates block-encoding
from block-encoding as Eq. (4), we organize r-QET by recursively defining a
series of block-encodings $\\{O_{X_{n}}\\}_{n}$ by

$O_{X_{n+1}}=\mathrm{QET}[O_{X_{n}},\vec{\phi}^{g}],\quad n=0,1,2,\ldots.$ (5)

Here, the initial input $O_{X_{0}}$ is dependent on $O_{A}$ (e.g.
$O_{X_{0}}=O_{A}$), and we have options for the parameter set $\vec{\phi}_{g}$
including its values and dimension. With these options, the above construction
forms a recursive relation of matrices $X_{n}=\braket{0}{O_{X_{n}}}{0}_{b}$,

$X_{n+1}=g(X_{n}),\quad X_{0}=\braket{0}{O_{X_{0}}}{0}_{b},$ (6)

with a variety of polynomial functions $g$. Recursive relations of matrices
like Eq. (6) are known to be available for complicated matrix functions
exemplified by a matrix logistic map [17, 18].

One of the most promising candidates for the recursive relation is Newton
iteration, which was originally invented for solving nonlinear equations [12].
Newton iteration for matrices enables us to efficiently compute various matrix
functions $f(A)$ with properly choosing the function $g$ as a result of
iteration $\lim_{N\to\infty}X_{n}=f(A)$. For instance, $g(X)=(3X-X^{3})/2$ and
$X_{0}=A$ generates the matrix sign function $X_{n}\to\mathrm{sign}(A)$
(defined later) and the one by $g(X)=2X-XAX$ and $X_{0}=\theta A$
($0<\theta\ll 1$) generates the matrix inversion $X_{n}\to A^{-1}$ [12]. When
r-QET combined with Newton iteration, we organize the parameter set
$\vec{\phi}_{g}$ so that the corresponding polynomial function $g$ reproduces
Newton iteration. The iteration continues until it achieves an allowable error
$\varepsilon$ as $\norm{X_{n}-f(A)}\leq\varepsilon$ ($\norm{\cdot}$ denotes
the operator norm), and thus the iteration number $n$ depends on
$\varepsilon$. Then, the resulting unitary gate $O_{X_{n}}$ provides an
accurate block-encoding for the function $f(A)$.

By iterative substitution of the recursive relation Eq. (5), the block-
encoding $O_{X_{n}}$ is rewritten by

$O_{X_{n}}=\mathrm{QET}[\mathrm{QET}[\ldots[\mathrm{QET}[O_{X_{0}},\vec{\phi}^{g}],\vec{\phi}^{g}],\ldots,\vec{\phi}^{g}],$
(7)

which has the form of $O_{X_{n}}=\mathrm{QET}[O_{X_{0}},\vec{\phi}_{n}]$. The
$(\deg(g^{n}))$-dimensional parameter set $\vec{\phi}_{n}$ is determined
solely by the $(\deg(g))$-dimensional one $\vec{\phi}^{g}$. Since the
desirable error $\varepsilon$ solely affects the iteration number $n$, we can
obtain a parameter set of QET achieving arbitrarily small error in an
analytical or numerically-stable way.

While we concentrate on the usage of the standard QET for implementing the
function $g$ here, we note that other types of QET/QSVT are also available.
For instance, if the recursive relation $g(X)$ include some matrices other
than $X$ (e.g. $g(X)=2X-XAX$ for matrix inversion), we need at-least multi-
variate QET/QSVT [20, 21].

Figure 1: (a) Some members of Padé family. (b) Intuitive picture of Newton
iteration based on $p_{2}(x)$. Each eigenvalue of $X_{n}$, represented by
$x_{n}$, approaches $1$ or $-1$ depending on its sign.

Analytical r-QET for Matrix Sign Function.— In r-QET, we should carefully
choose the recursive relation so that the tractable function $g$ can be
realized by QET and the resulting function $X_{n}$ is meaningful. Here, we
show the power of r-QET combined with Newton iteration by analytically
constructing the parameter set of QET for matrix sign functions.

Suppose that a hermitian matrix $A=\sum_{a}a\ket{a}\bra{a}$ has a spectral gap
$2\Delta$ ($>0$) around zero as $a\in[-1,-\Delta]\cup[\Delta,1]$. Then, the
matrix sign function $\mathrm{sign}(A)$ is defined by

$\mathrm{sign}(A)=\sum_{a}\mathrm{sign}(a)\ket{a}\bra{a},\quad\mathrm{sign}(a)=\frac{a}{|a|},$
(8)

which is useful for various tasks such as eigenstate filtering [22]. The
matrix sign function $\mathrm{sign}(A)$ can be generated by Newton-Shultz
iteration using a series of rational functions called Padé family
$\\{p_{l}(x)\\}_{l\in\mathbb{N}}$ with the initial input $X_{0}=A$ [12, 13].
As the recursive relation $g(X)$, we adopt the second simplest case,

$p_{2}(X)=\frac{1}{8}(15X-10X^{3}+3X^{5}),\quad\deg(p_{2})=5,$ (9)

since the simplest one $p_{1}(X)=(3X-X^{3})/2$ violates the condition (ii)
required for QET [See Fig. 1 (a)]. The convergence to $\mathrm{sign}(A)$ is
confirmed by the fact that every eigenvalue of $X_{n}$, denoted by $x_{n}$,
moves to $+1$ or $-1$ based on the recursive relation $x_{n+1}=p_{2}(x_{n})$
[See Fig. 1 (b)].

It is sufficient to determine the parameter set
$\vec{\phi}^{p_{2}}\in\mathbb{R}^{5}$, which cosntructs the relation,

$O_{X_{n+1}}=\mathrm{QET}[O_{X_{n}},\vec{\phi}^{p_{2}}],\quad
O_{X_{0}}=O_{A}.$ (10)

Finding $\vec{\phi}^{q}$ can be attributed to finding all the roots of
$1-p_{2}(x)p_{2}^{\ast}(x)$. While it is generally hard in an analytical way
when the degree is larger than $5$, we find them by the factorization,

$1-p_{2}(x)p_{2}^{\ast}(x)=-\frac{9}{64}(1-x^{2})^{3}(x^{2}-s)(x^{2}-s^{\ast}),$
(11)

with $s=(11+3\sqrt{15}i)/6$. This leads to the following analytical parameter
set $\vec{\phi}^{p_{2}}$ (See Supplementary Materials S1 for its detailed
calculation):

$\displaystyle\phi_{1}^{p_{2}}=0,\quad\phi_{2}^{p_{2}}=\pi+\frac{1}{2}\arctan\frac{\sqrt{15}}{7},$
$\displaystyle\phi_{3}^{p_{2}}=\pi+\frac{1}{2}\arctan\sqrt{15},\quad\phi_{4}^{p_{2}}=-\frac{1}{2}\arctan\sqrt{15},$
$\displaystyle\phi_{5}^{p_{2}}=-\frac{1}{2}\arctan\frac{\sqrt{15}}{7}.$ (12)

Therefore, by repeated substitution Eq. (7), we obtain an analytical QET
operator $\mathrm{QET}[O_{A},\vec{\phi}_{n}]$ for the matrix sign function
$\mathrm{sign}(A)$. The parameter set $\vec{\phi}_{n}$ has only eight
different angles $\pm\phi_{i}^{p_{2}}$ for $i=2,3,4,5$, and they always appear
in the fixed orders of
$\phi_{2}^{p_{2}}\to\phi_{3}^{p_{2}}\to\phi_{4}^{p_{2}}\to\phi_{5}^{p_{2}}$ or
$-\phi_{5}^{p_{2}}\to-\phi_{4}^{p_{2}}\to-\phi_{3}^{p_{2}}\to-\phi_{2}^{p_{2}}$.

Cost for Matrix Sign Function.— Let us evaluate the cost for the matrix sign
function. We repeat the recursive relation until the desirable error
$\varepsilon\in[0,1]$ is achieved as
$\norm{X_{n}-\mathrm{sign}(A)}\leq\varepsilon$. The convergence rate to
$\mathrm{sign}(A)$ based on the Newton iteration, Eq. (9), is dominated by the
gap $\Delta$ as follows [13]:

$\displaystyle\norm{X_{n}-\mathrm{sign}(A)}\leq(1-\Delta^{2})^{3^{n}}.$ (13)

By using the relation $\log(1-\Delta^{2})^{-1}\geq\Delta^{2}$ for
$\Delta\in[0,1)$, it is sufficient to choose the iteration number by

$n=\left\lceil\log_{3}\left(\frac{1}{\Delta^{2}}\log(1/\varepsilon)\right)\right\rceil.$
(14)

The cost is measured by the query complexity, i,e., the number of the oracle
$O_{A}$ in the unitary gate $O_{X_{n}}$ giving
$X_{n}=\mathrm{sign}(A)+\order{\varepsilon}$. Let the query complexity of
$O_{X_{n}}$ denote $q_{n}$. Then, the recursive relation Eq. (10) immediately
implies $q_{n+1}=5q_{n}$ and $q_{0}=1$. Under the proper iteration number $n$
by Eq. (14), we need the query complexity given by

$\displaystyle q_{n}$ $\displaystyle=$ $\displaystyle 5^{n}$ (15)
$\displaystyle\leq$ $\displaystyle
5\times\left(\frac{1}{\Delta^{2}}\log(1/\varepsilon)\right)^{\log_{3}5}$
$\displaystyle\in$
$\displaystyle\Theta\left(\frac{1}{\Delta^{2\log_{3}5}}\log^{\log_{3}5}(1/\varepsilon)\right).$
(16)

Let us compare with the standard QET approach [14, 1, 5], which uses
polynomial approximation of the error function
$\mathrm{erf}(kx)\simeq\mathrm{sign}(x)$ with large $k>0$. This yields the
query complexity $\Theta\left(\Delta^{-1}\log(1/\varepsilon)\right)$ to
implement $\mathrm{sign}(A)+\order{\varepsilon}$, which is known to be optimal
both in $\Delta$ and $\varepsilon$. Considering the value $\log_{3}5\simeq
1.465$, the query complexity of r-QET, Eq. (16), is polynomially larger than
the optimal one. This difference comes from the flexibility of the parameter
set $\vec{\phi}$. The standard QET approach uses
$\order{\Delta^{-1}\log(1/\varepsilon)}$ different parameters $\vec{\phi}$
obtained by numerically solving an
$\order{\Delta^{-1}\log(1/\varepsilon)}$-degree equation. In contrast, our
approach employs only eight different values in a fixed order independent of
$\Delta$ and $\varepsilon$. At some expense of the query complexity, r-QET can
serve $\mathrm{sign}(A)$ up to arbitrarily small error $\varepsilon$, with
completely avoiding numerical instability and working only with a few kinds of
gates.

While r-QET fails to achieve the optimal query complexity, it can overwhelm
the optimal protocol when correcting coherent error is taken into account. In
QET, the multiplicative coherent error on the parameter set
$\phi_{i}\to\phi_{i}(1+\delta)$ ($i=1,2,\ldots,q$) is a possible obstacle to
accurate implementation of matrix functions. Recently, A. K. Tan, et al. [16]
have shown that the coherent error $\delta$ on degree-$q$ QET can be
suppressed up to $\order{\delta^{k+1}}$ ($k\in\mathbb{N}$) by additional query
complexity,

$q_{\mathrm{correct}}\in\order{2^{k}(c_{\vec{\phi}})^{k^{2}}q}.$ (17)

Here, $c_{\vec{\phi}}$ denotes the number of different values in
$\vec{\phi}\in\mathbb{R}^{q}$. The standard QET approach achieving the
optimality is expected to have $c_{\vec{\phi}}\in\order{q}$ different values
in $\vec{\phi}$, and thereby its total query complexity amounts to
$\mathcal{O}(\Delta^{-k^{2}-1}\log^{k^{2}+1}(1/\varepsilon))$. By contrast,
r-QET employs constant values with $c_{\vec{\phi}}=8$ regardless of any other
parameters. Therefore, the total query complexity including the recovering
remains Eq. (16) as long as $k\in\order{1}$, and it can overwhelm the
originally-optimal protocol for arbitrary $k\geq 2$.

Generalization of r-QET for matrix sign functions.— Our result for matrix sign
functions can be generalized to other members of the Padé family or
implementation of polar decomposition by r-QSVT.

In the first case, the recursive relation $X_{n+1}=p_{l}(X_{n})$ with
$X_{0}=A$ ($l\in\mathbb{N}$) by the Padé family,

$p_{l}(x)=x\sum_{k=0}^{l}\frac{(2k-1)!!}{2^{k}k!}(1-x^{2})^{k}.$ (18)

casts the matrix sign function as $\lim_{n\to\infty}X_{n}=\mathrm{sign}(A)$.
The conditions (i)-(iii) for QET can be satisfied only when $l$ is even [23].
r-QET requires only $2l$ different values for the parameter set for the
degree-$(2l+1)$ polynomial $p_{l}(X)$, where one of them can be zero. As well
as $l=2$, we can analytically determine the parameter set for $l=4$ by the
quartic formula. Even for larger $l\geq 6$, the numerical instability for the
parameter determination is much more suppressed than the standard QET. The
advantage of generalization to larger degrees is the computational cost. The
convergence rate to $\mathrm{sign}(A)$ becomes faster as
$\norm{X_{n}-\mathrm{sign}(A)}\leq(1-\Delta^{2})^{(l+1)^{n}}$ [13], and hence
the query complexity $q_{n}$ in $O_{A}$ becomes smaller as

$\displaystyle
q\in\Theta\left(\frac{1}{\Delta^{2(1+\nu_{l})}}\log^{1+\nu_{l}}(1/\varepsilon)\right),$
(19) $\displaystyle\nu_{l}=\frac{\log(2l+1)}{\log(l+1)}-1.$ (20)

r-QET achieves the scaling
$\Theta(\Delta^{-(2+o(1))}\log^{1+o(1)}(1/\varepsilon))$ for large $l$ [23].
While it is still not optimal in the gap $\Delta$, it can reach the optimal
one in the desirable error $\varepsilon$.

Let us discuss the generalization to r-QSVT. QSVT produces block-encoding
executing polynomial transformation of every singular value from block-
encoding of a generic matrix $A$, and hence r-QSVT can be composed by
recursive iteration of QSVT (See Supplemental Materials S2 in detail). When we
use the same parameter set $\vec{\phi}^{p_{2}}$ for the matrix sign function,
r-QSVT organizes block-encoding $\\{O_{X_{n}}\\}_{n}$ reproducing

$X_{n+1}=\frac{X_{n}}{8}\left\\{15-10X_{n}^{\dagger}X_{n}+3(X_{n}^{\dagger}X_{n})^{2}\right\\},\quad
X_{0}=A.$ (21)

Under this recursive iteration, the matrix converges to $X_{n}\to
A(\sqrt{A^{\dagger}A})^{-1}$, which is the unitary part of the polar
decomposition of the non-singular matrix $A$. With the same iteration number
Eq. (14), r-QSVT achieves an arbitrarily small error $\varepsilon$ and
analyticity (or numerical stability) of the parameter set for the polar
decomposition.

Discussion and Conclusion.— In this Letter, we propose recursive QET/QSVT that
executes recursive relations by QET/QSVT. All the parameters can be determined
by low-degree polynomials for recursive relations, which enables analytical or
numerically-stable calculation and also leads to feasibility of recovering
coherent noise. Particularly, the construction of matrix sign function when
combined with Newton iteration is the first analytical result on parameters
for useful functions other than the trivial Chebyshev polynomials. Indeed,
with the analytically-obtained parameters given by Eq. (12), we can execute
eigenstate filtering, and thereby solve quantum linear system problems [22]
(See also Supplemental Materials S3). It will not only contribute to
understanding the expressibility of QET with restricted parameters [8, 9], but
also open up an efficient path to construct QET for various purposes.

We conclude this Letter with some future directions of r-QET/r-QSVT. The first
direction is to explore accurate implementation of other matrix functions by
Newton iteration. We have concentrated on using the standard QET as a
subroutine for the matrix sign function for simplicity. However, Newton
iteration covers various matrix functions and can show quadratic or faster
convergence. Exploiting a series of QET/QSVT protocols, such as those for
Fourier series expansions [24, 25, 26] and multivariate polynomial functions
[20, 21], to reproduce their recursive relations will provide alternative ways
to execute various matrix functions with efficiently and accurately providing
the required parameters.

Second, it will be also important to seek for the usage of various recursive
relations. For instance, recursive relations based on the matrix logistic map
can generate either chaotic or non-chaotic behaviors of each eigenvalue [17,
18]. While its usefulness in quantum many-body systems is still controversial,
examining compatibility of such recursive relations with r-QET/r-QSVT will
open up a new way of quantum operations by QET/QSVT.

###### Acknowledgements.

K. M. is supported by RIKEN Special Postdoctoral Researcher Program. This work
is supported by MEXT Quantum Leap Flagship Program (MEXTQLEAP) Grant No.
JPMXS0118067394, JPMXS0120319794, and JST COI-NEXT program Grant No.
JPMJPF2014.

## References

  * 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_, STOC 2019 (Association for Computing Machinery, New York, NY, USA, 2019) pp. 193–204.
  * Low and Chuang [2017a] G. H. Low and I. L. Chuang, Optimal Hamiltonian Simulation by Quantum Signal Processing, Phys. Rev. Lett. 118, 010501 (2017a).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Grover [1997] L. K. Grover, Quantum Mechanics Helps in Searching for a Needle in a Haystack, Phys. Rev. Lett. 79, 325 (1997).
  * Martyn  _et al._ [2021] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, Grand Unification of Quantum Algorithms, PRX Quantum 2, 040203 (2021).
  * Haah [2019] J. Haah, Product decomposition of periodic functions in quantum signal processing, Quantum 3, 190 (2019).
  * Chao  _et al._ [2020] R. Chao, D. Ding, A. Gilyen, C. Huang, and M. Szegedy, Finding Angles for Quantum Signal Processing with Machine Precision, arXiv:2003.02831 [quant-ph] (2020).
  * Dong  _et al._ [2021] Y. Dong, X. Meng, K. B. Whaley, and L. Lin, Efficient phase-factor evaluation in quantum signal processing, Phys. Rev. A 103, 042419 (2021).
  * Wang  _et al._ [2022a] J. Wang, Y. Dong, and L. Lin, On the energy landscape of symmetric quantum signal processing, Quantum 6, 850 (2022a).
  * Høyer [2000] P. Høyer, Arbitrary phases in quantum amplitude amplification, Phys. Rev. A 62, 052304 (2000).
  * Long [2001] G. L. Long, Grover algorithm with zero theoretical failure rate, Phys. Rev. A 64, 022307 (2001).
  * Higham [2008] N. J. Higham,  _Functions of Matrices_ (Society for Industrial and Applied Mathematics, 2008).
  * Kenney and Laub [1991] C. Kenney and A. J. Laub, Rational Iterative Methods for the Matrix Sign Function, SIAM Journal on Matrix Analysis and Applications 12, 273 (1991), https://doi.org/10.1137/0612020 .
  * Low and Chuang [2017b] G. H. Low and I. L. Chuang, Hamiltonian simulation by uniform spectral amplification, arXiv:1707.05391 [quant-ph] (2017b).
  * Lin and Tong [2020a] L. Lin and Y. Tong, Near-optimal ground state preparation, Quantum 4, 372 (2020a).
  * Tan  _et al._ [2023] A. K. Tan, Y. Liu, M. C. Tran, and I. L. Chuang, Error Correction of Quantum Algorithms: Arbitrarily Accurate Recovery Of Noisy Quantum Signal Processing, arXiv:2301.08542 [quant-ph] (2023).
  * Navickas  _et al._ [2011] Z. Navickas, R. Smidtaite, A. Vainoras, and M. Ragulskis, The logistic map of matrices, Discrete Continuous Dyn. Syst. Ser. B 16, 927 (2011).
  * Łukasz Pawela and Życzkowski [2023] Łukasz Pawela and K. Życzkowski, Matrix logistic map: fractal spectral distributions and transfer of chaos (2023), arXiv:2303.06176 [nlin.CD] .
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum Algorithm for Systems of Linear Equations with Exponentially Improved Dependence on Precision, SIAM Journal on Computing 46, 1920 (2017), https://doi.org/10.1137/16M1087072 .
  * Rossi and Chuang [2022] Z. M. Rossi and I. L. Chuang, Multivariable quantum signal processing (m-QSP): prophecies of the two-headed oracle, Quantum 6, 811 (2022).
  * Borns-Weil  _et al._ [2023] Y. Borns-Weil, T. Saffat, and Z. Stier, A Quantum Algorithm for Functions of Multiple Commuting Hermitian Matrices, arXiv:2302.11139 [quant-ph] (2023).
  * Lin and Tong [2020b] L. Lin and Y. Tong, Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems, Quantum 4, 361 (2020b).
  * [23] Exactly speaking, we was not able to found the proof of the latter part of the condition (ii) for generic even $l$, $|p_{l}(x)|\geq 1$ for any $x\notin[-1,1]$, although it seems to be correct from their graphs. By computing the derivatives of $p_{l}(x)$ at $x=\pm 1$, we have proven that $p_{l}(x)$ is compatible with the standard QET at least up to $l=10^{3}$. While we believe that the exponent $\nu_{l}$ [See Eq. (19)] can be arbitrarily small under $l\to\infty$, it is ensured to reach $\nu_{10^{3}}\simeq 0.100$ at the present stage.
  * Dong  _et al._ [2022] Y. Dong, L. Lin, and Y. Tong, Ground-State Preparation and Energy Estimation on Early Fault-Tolerant Quantum Computers via Quantum Eigenvalue Transformation of Unitary Matrices, PRX Quantum 3, 040305 (2022).
  * de Lima Silva  _et al._ [2022] T. de Lima Silva, L. Borges, and L. Aolita, Fourier-based quantum signal processing, arXiv:2206.02826 [quant-ph] (2022).
  * Wang  _et al._ [2022b] X. Wang, Y. Wang, Z. Yu, and L. Zhang, Quantum Phase Processing: Transform and Extract Eigen-Information of Quantum Systems, arXiv:2209.14278 [quant-ph] (2022b).

Supplemental Materials for
“Recursive Quantum Eigenvalue/Singular-Value Transformation;
Analytic Construction of Matrix Sign Function by Newton iteration”

Kaoru Mizuta1,2 and Keisuke Fujii3,4,1,5

1RIKEN Center for Quantum Computing (RQC), Hirosawa 2-1, Wako, Saitama
351-0198, Japan
2Department of Applied Physics, The University of Tokyo, Hongo 7-3-1, Bunkyo,
Tokyo 113-8656, Japan
3Graduate School of Engineering Science, Osaka University,
1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan.
4Center for Quantum Information and Quantum Biology, Osaka University, Japan.
5Fujitsu Quantum Computing Joint Research Division at QIQB,
Osaka University, 1-2 Machikaneyama, Toyonaka 560-0043, Japan
(Dated: )

##  S1 Parameter determination of QET

In this section, we show how to determine the parameter set
$\vec{\phi}\in\mathbb{R}^{q}$ in QET, for this Letter to be self-contained.
The discussion here is mainly based on Ref. [1].

Suppose that we want to realize a degree-$q$ polynomial $f(x)$ satisfying
(i)-(iii) in the main text by degree-$q$ QET. This exploits a technique of
quantum signal processing (QSP) [2], since the achievable functions are
equivalent to QET/QSVT by so-called qubitization technique [3]. We begin with
determining the parameter set for QSP to realize the function $f(x)$. QSP is a
repetition of single-qubit rotation operators that embodies a polynomial
function of a certain signal $x\in[-1,1]$. A rotation-based degree-$q$ QSP
operator $\mathrm{QSP}[x,\vec{\phi}]$ on $\mathbb{C}^{2}$ is defined by

$\displaystyle\mathrm{QSP}[x,\vec{\phi}^{\prime}]$ $\displaystyle=$
$\displaystyle
e^{i\phi_{0}^{\prime}Z}\prod_{i=1}^{q}\left[W(x)e^{i\phi_{i}^{\prime}Z}\right],$
(S1) $\displaystyle W(x)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}x&i\sqrt{1-x^{2}}\\\
i\sqrt{1-x^{2}}&x\end{array}\right),$ (S4)

where $\vec{\phi}^{\prime}\in\mathbb{R}^{q+1}$ is a tunable parameter set and
$Z$ is the Pauli-$Z$ operator. The QSP operator can be written in the form of

$\mathrm{QSP}[x,\vec{\phi}^{\prime}]=\left(\begin{array}[]{cc}f(x)&ih(x)\sqrt{1-x^{2}}\\\
ih^{\ast}(x)\sqrt{1-x^{2}}&f^{\ast}(x),\end{array}\right)$ (S5)

and there exists a corresponding parameter set
$\vec{\phi}^{\prime}\in\mathbb{R}^{q+1}$ if and only if $f(x)$ satisfies
(i)-(iii) in the main text. Then, a polynomial function $h(x)$ satisfies

  1. (i′)

$h(x)$ has a degree at most $q-1$

  2. (ii′)

$h(x)$ has a parity $(-1)^{q-1}$.

  3. (iii′)

$f(x)f^{\ast}(x)+(1-x^{2})h(x)h^{\ast}(x)=1$ for any $x$.

To obtain a parameter set $\vec{\phi}^{\prime}$ for the rotation-based QSP, we
should obtain the polynomial function $h(x)$. This can be executed by finding
all the roots $\\{s_{i}\\}_{i=1}^{q}$ of $1-f^{\ast}(x)f(x)$, which is a
degree-$q$ polynomial in $x^{2}$. From the conditions (i)-(iii) in the main
text, the polynomial $1-f^{\ast}(x)f(x)$ is always factorized in the form of

$1-f^{\ast}(x)f(x)=\alpha(1-x^{2})\prod_{s_{i};\mathrm{Im}(s_{i})\geq
0}(x^{2}-s_{i})(x^{2}-s_{i}^{\ast})$ (S6)

with $\alpha>0$, and thus

$h(x)=\sqrt{\alpha}\prod_{s_{i}\in S;\mathrm{Im}(s_{i})\geq 0}(x^{2}-s_{i}),$
(S7)

gives a proper choice satisfying (i′)-(iii′).

Once one obtains $h(x)=\sum_{n=0}^{q-1}h_{n}x^{n}$ for the desired polynomial
$f(x)=\sum_{n=0}^{q}f_{n}x^{n}$, the parameter set $\vec{\phi}^{\prime}$ is
determined in a recursive way. We define $\phi_{q}^{\prime}$ by

$e^{2i\phi_{q}^{\prime}}=\frac{f_{q}}{h_{q-1}},$ (S8)

and then we update the polynomials $f(x)$ and $h(x)$ by

$\displaystyle\tilde{f}(x)$ $\displaystyle=$ $\displaystyle
e^{-i\phi_{q}^{\prime}}xf(x)+e^{i\phi_{q}^{\prime}}(1-x^{2})h(x),$ (S9)
$\displaystyle\tilde{h}(x)$ $\displaystyle=$ $\displaystyle
e^{i\phi_{q}^{\prime}}xh(x)-e^{-i\phi_{q}^{\prime}}f(x).$ (S10)

The pair of $\tilde{f}(x)$ and $\tilde{h}(x)$ satisfies (i)-(iii) and
(i′)-(iii′) for the degree $q-1$, and hence it can be realized by a
degree-$(q-1)$ QSP. The parameter $\phi_{q}^{\prime}$ is determined by Eq.
(S8) with using $\tilde{f}(x)$ and $\tilde{h}(x)$ instead of $f(x)$ and
$h(x)$. We repeat this procedure until the obtained polynomials become
trivial.

Finally, we transform the angles $\vec{\phi}^{\prime}$ into those for QET.
With a proper choice of two-dimensional basis, QET can be identified with a
reflection-based QSP operation $\mathbb{C}^{2}$, which is defined by

$\displaystyle\overline{\mathrm{QSP}}[x,\vec{\phi}]$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{q}\left[e^{i\phi_{i}Z}R(x)\right],\quad\vec{\phi}\in\mathbb{R}^{q},$
(S11) $\displaystyle R(x)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}x&\sqrt{1-x^{2}}\\\
\sqrt{1-x^{2}}&-x\end{array}\right).$ (S14)

Namely, QET and reflection-based QSP realize the same function by
$\braket{0}{\mathrm{QET}[A,\vec{\phi}]}{0}$ and
$\braket{0}{\overline{\mathrm{QSP}}[x,\vec{\phi}]}{0}$ with the same angles
$\vec{\phi}$. Achievable functions by them are equivalent to those by
rotation-based QSP, where we choose the angles by

$\phi_{1}=\phi_{0}^{\prime}+\phi_{q}^{\prime}+(q-1)\frac{\pi}{2},\quad\phi_{i}=\phi_{i-1}^{\prime}-\frac{\pi}{2}\quad(2\leq
i\leq q).$ (S15)

Numerical instability arises in the root-finding problem of the degree-$q$
polynomial, Eq. (S6), and the error is accumulated via iterative computation
by Eqs. (S8)-(S10). In Newton-iteration-based QET, we do not suffer from them
since we employ low-degree QET. When we consider
$f(x)=p_{2}(x)=(15x-10x^{3}+3x^{3})/8$ to construct a matrix sign function as
the main text, the factorization Eq. (11) analytically gives the pairwise
function $h(x)$ by

$h(x)=\frac{3}{8}x^{4}-\frac{1}{16}(17+3\sqrt{15}i)x^{2}+\frac{1}{16}(11+3\sqrt{15}i).$
(S16)

The iterative determination of the angles by Eqs. (S8)-(S10) can be
analytically completed, which results in Eq. (12). We ensure the analyticity
of the angles $\vec{\phi}$ also when we use the higher-order Padé family
$f(x)=p_{4}(x)$, where $h(x)$ can be analytically computed by the
factorization,

$\displaystyle 1-(p_{4}(x))^{2}=-(1-x^{2})^{5}\sum_{i=0}^{4}a_{i}(x^{2})^{i},$
(S17) $\displaystyle a_{0}=1,\quad a_{1}=-\frac{17305}{16384},\quad
a_{2}=\frac{14235}{16384},$ (S18) $\displaystyle
a_{3}=-\frac{6475}{16384},\quad a_{4}=\frac{1225}{16384},$ (S19)

and the quartic formula.

##  S2 Extension to QSVT

In this section, we briefly discuss the extension to r-QSVT, especially for
the matrix sign function. Construction similar to the main text results in
QSVT with tractable parameters for polar decomposition of generic matrices.

We begin with introducing QSVT [1], with restricting to odd degrees for
simplicity. For a generic finite-dimensional matrix
$A\in\mathbb{C}^{d_{1}\times d_{2}}$, its block-encoding is defined by a
unitary operator $O_{A}$ such that

$\Pi_{1}O_{A}\Pi_{2}=A$ (S20)

is satisfied with certain $d_{1}$\- and $d_{2}$-dimensional projections
$\Pi_{1}$ and $\Pi_{2}$. With parametrized rotations designated by

$R_{\phi}^{\Pi_{1}}=e^{i\phi(2\Pi_{1}-I)},\quad
R_{\phi}^{\Pi_{2}}=e^{i\phi(2\Pi_{2}-I)},$ (S21)

a degree-$q$ QSVT operator for odd $q\in\mathbb{N}$ with a parameter set
$\vec{\phi}\in\mathbb{R}^{q}$ is defined by

$\mathrm{QSVT}[O_{A},\vec{\phi}]=R_{\phi_{1}}^{\Pi_{1}}O_{A}\prod_{i=1}^{(q-1)/2}\left[R_{\phi_{2i}}^{\Pi_{2}}O_{A}^{\dagger}R_{\phi_{2i+1}}^{\Pi_{1}}O_{A}\right].$
(S22)

The matrix $A$ is always described by singular value decomposition (SVD) as
$A=\sum_{\sigma}\sigma\ket{u_{\sigma}}\bra{v_{\sigma}}$, where $\sigma$
($>0$), $\ket{u_{\sigma}}\in\mathbb{C}^{d_{1}}$, and
$\ket{v_{\sigma}}\in\mathbb{C}^{d_{2}}$ are respectively called a singular
value, a left singular vector, and a right singular vector. QSVT executes
polynomial processing of every singular value given by

$\Pi_{1}\mathrm{QSVT}[O_{A},\vec{\phi}]\Pi_{2}=\sum_{\sigma}f(\sigma)\ket{u_{\sigma}}\bra{v_{\sigma}},$
(S23)

where the scalar function $f(x)$ satisfies all the conditions (i)-(iii) in the
main text. We note that this is not equal to a polynomial function $f(A)$ in
general, unlike QET.

Recursive QSVT (r-QSVT) is organized completely in a similar manner to r-QET;
We construct a series of unitaries $\\{O_{X_{n}}\\}_{n}$ by

$O_{X_{n+1}}=\mathrm{QSVT}[O_{X_{n}},\vec{\phi}^{g}],$ (S24)

which begins from a unitary $O_{X_{0}}$ such that
$\Pi_{1}O_{X_{0}}\Pi_{2}=X_{0}$. It casts a recursive relation as

$\displaystyle X_{n+1}$ $\displaystyle=$
$\displaystyle\Pi_{1}O_{X_{n+1}}\Pi_{2}$ (S25) $\displaystyle=$
$\displaystyle\sum_{\sigma_{n}}g(\sigma_{n})\ket{u_{\sigma_{n}}}\bra{v_{\sigma_{n}}},$

where SVD of $X_{n}$ is given by
$X_{n}=\sum_{\sigma_{n}}\sigma_{n}\ket{u_{\sigma_{n}}}\bra{v_{\sigma_{n}}}$.
If we choose the function $g$ based on Newton iteration, some complicated
matrix functions of $A$ will be realized. For instance, when the parameter set
$\vec{\phi}^{g}$ is chosen as Eq. (9), which is given for the matrix sign
function, it forms the recursive relation,

$\displaystyle X_{n+1}$ $\displaystyle=$
$\displaystyle\sum_{\sigma_{n}}\frac{1}{8}(15\sigma_{n}-10\sigma_{n}^{3}+3\sigma_{n}^{5})\ket{u_{\sigma_{n}}}\bra{v_{\sigma_{n}}}$
$\displaystyle=$
$\displaystyle\frac{1}{8}\left\\{15-10(X_{n}X_{n}^{\dagger})+3(X_{n}X_{n}^{\dagger})^{2}\right\\}X_{n},$

with $X_{0}=A$. This is nothing but the one for computing the polar
decomposition $U_{A}\in\mathbb{C}^{d_{1}\times d_{2}}$, which is a unitary
such that $A=U_{A}P_{A}$ with a positive semidefinite matrix
$P_{A}\in\mathbb{C}^{d_{2}\times d_{2}}$. Namely, we obtain
$\Pi_{1}O_{X_{n}}\Pi_{2}\to U_{A}$ by the recursive construction Eq. (S24).
Since the convergence rate is the same as the one of matrix sign functions
under the assumption ${}^{\forall}\sigma\in[\Delta,1]$, the query complexity
to $O_{A}$ to achieve an allowable error $\varepsilon$ is given by Eq. (16).

We can see the realization of polar decomposition from another aspect. The
components of the polar decomposition are respectively represented by

$U_{A}=\sum_{\sigma}\ket{u_{\sigma}}\bra{v_{\sigma}},\quad
P_{A}=\sum_{\sigma}\sigma\ket{v_{\sigma}}\bra{v_{\sigma}}.$ (S27)

Let us consider a QSVT operator $\mathrm{QSVT}[O_{A},\vec{\phi}]$ with
organizing $\vec{\phi}$ from $\vec{\phi}^{g}$ as well as the matrix sign
function. When the iteration number $n$ is large enough to satisfy Eq. (14),
QSVT processes every singular value by
$f(\sigma)=\mathrm{sign}(\sigma)+\order{\varepsilon}=1+\order{\varepsilon}$.
As a result, the QSVT operator $\mathrm{QSVT}[O_{A},\vec{\phi}]$ can
accurately reproduce the polar decomposition $U_{A}$ given by the formula Eq.
(S27). In a similar way, r-QSVT with other even-order members of Padé family
$g_{2l^{\prime}}(x)$ outputs the polar decomposition of generic matrices, with
yielding the query complexity Eq. (19).

##  S3 Quantum algorithms with
matrix sign function of r-QET

In the main text, we have shown that r-QET/r-QSVT provides analytical
parameter sets for matrix sign functions. Here, we provide that the same
parameter sets are available for eigenstate filtering, eigenstate preparation,
and quantum linear system problems (QLSP) as a complementary result by
combining the results of Ref. [22].

###  S3.1 Eigenstate filtering and eigenstate preparation

Eigenstate filtering is to filter out eigenstates having eigenenergies larger
than a certain cutoff. Assuming that a hermitian matrix $A=\sum_{a\neq
0,|a|\leq 1}a\ket{a}\bra{a}$ has a spectral gap at $a=0$, it is equivalent to
apply the projection,

$P_{+}(A)=\frac{1+\mathrm{sign}(A)}{2},$ (S28)

We use the matrix sign function generated by r-QET, such that
$\braket{0}{O_{X_{n}}}{0}_{b}=\mathrm{sign}(A)+\order{\varepsilon}$. Then, the
controlled operation defined by

$\displaystyle O_{\Delta}(A)=$ $\displaystyle\quad
e^{i(\pi/4)Y_{b^{\prime}}}\left(\ket{0}\bra{0}_{b^{\prime}}\otimes
I+\ket{1}\bra{1}_{b^{\prime}}\otimes
O_{X_{n}}\right)e^{i(\pi/4)Y_{b^{\prime}}}.$ (S29)

reproduces the filtering with additional one qubit $b^{\prime}$ by
$\braket{00}{O_{\Delta}(A)}{00}_{bb^{\prime}}=P_{+}(A)+\order{\varepsilon}$.
The query complexity in $O_{A}$ for the eigenstate filtering is

$\order{\Delta^{-2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\varepsilon)},$ (S30)

when we use the matrix sign function by r-QET with the function $p_{l}(x)$
[See Eq. (19)]. While it is worse than the optimal cost
$\order{\Delta^{-1}\log(1/\varepsilon)}$, achieved by the standard QET [22],
r-QET determines all the parameters composed of $4l+1$ different rotation
angles in numerically-cheaper ways. For instance, it requires only $9$
different rotational angles, which are analytically given by
$\\{\pm\phi_{i}\\}_{i=2,3,4,5}$ in Eqs. (12) and $\pi/4$ in Eq. (S29), in the
simplest case $l=2$.

Eigenstate preparation is to prepare a certain isolated eigenstate of $A$
(here, we set the eigenvalue at zero). Assuming that $A$ has no eigenvalue in
$[-\Delta,\Delta]$ except for the zero eigenvalue, it is completed by applying
the projection,

$P_{0}(A)=\ket{0}\bra{0}=P_{+}\left(\frac{A+\Delta/2}{1+\Delta/2}\right)P_{+}\left(-\frac{A-\Delta/2}{1+\Delta/2}\right),$
(S31)

to a proper initial state $\ket{\psi_{0}}$. In the above filtering, we use the
block-encoding of $(A\pm\Delta/2)/(1+\Delta/2)$ for normalization, which
requires one query to a controlled-$O_{A}$ operation. We also need single-
qubit unitary gates on an additional qubit, whose rotation angle is
$\pm\arctan\sqrt{\Delta/2}$. When the initial state $\ket{\psi_{0}}$, prepared
by a unitary $U_{0}$, has the overlap $\gamma=|\braket{0}{\psi_{0}}|$, r-QET
with $p_{l}(x)$ requires

$\order{\frac{1}{\gamma}\Delta^{-2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\gamma\varepsilon)}$
(S32)

queries to $O_{A}$ and $\order{\gamma^{-1}}$ queries to $U_{0}$. Here, we use
the amplitude amplification [5] which allows the quadratic speedup from the
time $\order{\gamma^{-2}}$, expected by the success probability. At the
expense of the cost to some extent compared to the optimal one
$\order{\gamma^{-1}\Delta^{-1}\log(1/\varepsilon)}$ [22], r-QET works only
with $4l+3$ different rotation angles. In addition, they can be analytically
obtained in the simplest cases $l=2,4$, or avoid numerical instability for
larger $l$.

###  S3.2 Quantum linear system problem (QLSP)

Quantum linear system problem (QLSP) is to obtain $\ket{x}=A^{-1}\ket{b}$ from
a given non-singular hermitian matrix $A$ and a given vector $\ket{b}$. Here,
we define the condition number of $A$ by $\kappa=\Delta^{-1}$, and the vector
$\ket{b}$ is prepared by a unitary gate $U_{b}$. Recently, it has been shown
that QLSP can be solved by using eigenstate filtering and eigenstate
preparation as subroutines [22]. This relies on the fact that $\ket{x}$ is
obtained by an isolated eigenstate as

$\displaystyle\tilde{A}\left(\begin{array}[]{c}\ket{x}\\\
0\end{array}\right)=0,$ (S35)
$\displaystyle\tilde{A}=\left(\begin{array}[]{cc}0&A(I-\ket{b}\bra{b})\\\
(I-\ket{b}\bra{b})A&0\end{array}\right).$ (S38)

The block-encoding of $\tilde{A}$ is organize by $\order{1}$ queries to
$O_{A}$ and $U_{b}$. The algorithm prepares a state having an $\order{1}$
overlap with ${}^{\mathrm{t}}(\ket{x},0)$ by Zeno effect reproduced by
eigenstate filtering, and then make a projection onto
${}^{\mathrm{t}}(\ket{x},0)$ within $\order{\varepsilon}$ by eigenstate
preparation in $\tilde{A}$.

In the first step, the algorithm begins with preparing the zero-energy
eigenstate ${}^{\mathrm{t}}(\ket{b},0)$ of the trivial Hamiltonian
$\tilde{A}_{0}=X\otimes(I-\ket{b}\bra{b})$ Then, to approximately obtain the
eigenstate of the target Hamiltonian $\tilde{A}$, it uses a path

$\tilde{A}(s)=f_{s}\tilde{A}+(1-f_{s})\tilde{A}_{0},$ (S39)

with a specific scheduling function $f_{s}$ satisfying $f_{s=0}=0$ and
$f_{s=1}=1$. The algorithm mimics Zeno effect along this path by a series of
eigenstate filtering, which is represented by

$\prod_{t=0}^{M}P_{0}\left(\tilde{A}(t/M)\right)\left(\begin{array}[]{c}\ket{b}\\\
0\end{array}\right)=c\left(\begin{array}[]{c}\ket{x}\\\
0\end{array}\right)+\dots.$ (S40)

To achieve the overlap $c\in\order{1}$ with ${}^{\mathrm{t}}(\ket{x},0)$, we
need $M\in\order{(\log\kappa)^{2}}$-times repetition and the error smaller
than $\varepsilon_{\mathrm{EF}}\in\order{M^{-2}}$ of the projections $P_{0}$
in eigenstate filtering. As a result, this step requires at-most

$\displaystyle\order{M\Delta^{-2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\varepsilon_{\mathrm{EF}})}$
$\displaystyle\quad=\order{\kappa^{2(1+\nu_{l})}(\log\kappa)^{2}(\log\log\kappa)^{1+\nu_{l}}}$
(S41)

queries to the unitary gates $O_{A}$ and $U_{b}$. In order to organize the
block-encoding of $\tilde{A}(s)$ for every $s=0,1/M,2/M,\ldots,1$, we also
need $M+1\in\order{(\log\kappa)^{2}}$ different angles, which are analytically
determined by $\arctan\sqrt{f_{s}/(1-f_{s})}$.

After the reproduction of Zeno effect by eigenstate filtering, the algorithm
applies the projection $P_{0}(\tilde{A})$ within the allowable error
$\varepsilon$. This step yields
$\order{\Delta^{-2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\varepsilon)}=\order{\kappa^{2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\varepsilon)}$
query complexity. Due to the $\order{1}$ overlap after Zeno effect, this step
outputs the accurate solution of QSLP, $\mathrm{t}(\ket{x},0)$, with
$\order{1}$ success probability probability. In total, the QLSP algorithm
exploiting eigenstate filtering with r-QET is executed by

$\tilde{\mathcal{O}}\left(\kappa^{2(1+\nu_{l})}\log^{1+\nu_{l}}(1/\varepsilon)\right)$
(S42)

queries to the unitaries $O_{A}$ and $U_{b}$. This scaling is worse than that
of the original one using the standard optimal QET for eigenstate filtering
[22], whose scaling is $\tilde{\mathcal{O}}(\kappa\log(1/\varepsilon))$.
However, the algorithm based on r-QET has a strong advantage in the parameter
determination also in QLSP. While the standard one requires
$\tilde{O}(\kappa\log(1/\varepsilon))$ different parameters found by
numerically-unstable computation, r-QET works only with
$4l+3+M+1\in\order{(\log\kappa)^{2}}$ different parameters, and all of them
can be found analytically or in much numerically-cheaper ways.


be obtained with the true design parameters. Results on 200 transmissions from
the test set are shown in figure 6.7.

For the transmissions in the first bin, the difference in MAE compared to what
we would have for the true parameters was lower than 0.01. Most of the
transmissions end up in this bin, so we can conclude that all three loss
functions lead to good results. Of the three, the combined loss is superior to
the others, but only by a small margin. We observe that there are some
outliers with a high MAE. When looking at these transmissions, we saw that
they barely oscillate. It is to be expected that the Fourier transform of
these transmissions does not lead to much improvement.

At last, we compare the method when we make a guess for $\delta_{0}$ every 100
epochs or when we only use it to make an initial guess. We used the MSE
Combined loss function to make the comparison. It is clear that making a new
guess for $\delta_{0}$ leads to more transmission with a low MAE.

###  6.2 Designing optimal $\theta,n,l$

We now perform inverse design on $\theta,n$ and $l$. We need to make some
adaptations compared to our method in the earlier section. We still perform
gradient descent on the MSE of the predicted and desired transmissions
$T(\lambda)$. To get the initial values right, we look at the Fourier
transform of the transmission. In this case, there is no design parameter that
clearly determines the Fourier transform. Therefore, we need to compute the
minimum for every design parameters.

We compute this minimum in three steps. In the first step, we pick $n$ and $l$
randomly and choose the value of $\theta$ that minimizes the loss for these
fixed parameters. Then we use the value obtained for $\theta$ to compute a
minimum value of $n$, using a random $l$. At last, we use our guess for
$\theta$ and $n$ to compute $l$. In this 3-step process, we end up with
initial values for $\theta$, $n$ and $l$. This is illustrated in figure 6.8.
Again, there are three reasonable loss functions we can use to guess the
design parameters: the MSE on the transmissions, the MSE on the Fourier power
spectrum or a combination of both.

(a)

(b)

(c)

Figure 6.8: Before the gradient descent, we first minimize the loss for all
three parameters. The three images show the three consecutive minimizations we
perform. a) Step 1: We minimize $\theta$ for random values of $n$ and $l$. b)
Step 2: We minimize $n$ with the initial value obtained for $\theta$ and a
random $l$. c) Step 3: We minimize $l$ using the initial values we obtained
for $\theta$ and $n$.

We test this method on 200 transmissions from the test set. Again, we make a
new guess every 100 epochs. Results are shown in figure 6.9. The transmissions
in the first bin have to a MAE of less than 0.01. Of the three secondary
losses, the combined MSE gives the best results. We also see that using the
MSE on the transmission to guess the parameters works quite well. Looking at
the MSE on the Fourier transform alone gives the worst results.

We also check if we really need to make a guess every 100 epochs. We use the
best secondary loss function, the combined MSE, to make the comparison. We
compare to the case where the combined loss function is only used to determine
the initial guess. Results are shown in figure 6.9. The results are not as
significant as for $F$ and $\delta_{0}$. Making a new guess every 100 epochs
does lead to a lower MAE, but only by a small margin.

Figure 6.9: We performed inverse design on 200 transmissions in the test set.
We look at the MAE between the desired transmission and the final predicted
transmission. We then take difference between this MAE and the MAE we would
obtain with the true parameters. Transmissions in the first bin have a MAE of
less than $1\%$. In the first plot, we compared three loss functions to
initialize the design parameters. After 100 epochs, the loss was used to re-
initialize these parameters. The second plot compares our method when we re-
initialize the design parameters to the case where we only initialize once. We
used the combined MSE in the second plot.

As already discussed, the inverse design converges to a solution that is not
unique. To see this, we look at what happens for a transmission in the test
set shown in figure 6.10. The design parameters converge to a global minimum.
The optimal parameters are however not the same as the true parameters used to
generate the spectrum. This means that we have found a different design for
the Fabry-Pérot resonator with the same transmission $T(\lambda)$.

Figure 6.10: Inverse design for a transmission in the test set. The lower
panels show the desired transmission in blue and the transmission predicted
for the obtained design parameters in red. The upper panels show the evolution
of the design parameters. The panel in the lower right shows the evolution of
the MSE loss on the transmissions during training. We observe that $\theta$
makes a discontinuous jump when the parameter is re-initialized at 100 epochs.

###  6.3 Application: the block pulse and the elephant

Let us now conclude with what we have promised at the very beginning, inverse
design on an arbitrary transmission $T(\lambda)$. We start out with the case
of a block pulse. We used the inverse design method where we made a guess
every 100 epochs based on the combined MSE on the transmissions and the
Fourier transform. We did this for both $F$ and $\delta_{0}$ and the original
design parameters $\theta$, $n$ and $l$. The result is shown in figure 6.11.
The results are as good as we can hope for. The transmission adapts itself to
be near the block pulse, while still remaining a Fabry-Pérot transmission.

To show that we can really approximate arbitrary transmissions, we use an
example from the novel ”Le petit prince” by Antoine de Saint-Exupéry. The book
opens with a picture of an elephant that is being eaten by a snake. When this
picture is shown to adults in the book, they see an ordinary hat. If only they
knew.

We made an approximation of the elephant and performed inverse design on this
’transmission’. It was somewhat surprising to see that the algorithm was able
to handle this and did converge to a solution. The result is shown in figure
6.12. This shows that our inverse design method is capable of approximating
really any transmission.

(a)

(b)

Figure 6.11: Inverse design of a block pulse. The lower panels show the
desired block pulse in blue and the transmission predicted for the obtained
design parameters in red. The upper panels show the evolution of the design
parameters. The MSE on the transmissions during training is also plotted. a)
Inverse design on $F$ and $\delta_{0}$. We converge to a minimum. b) Inverse
design on $\theta$, $n$ and $l$. The design parameters converge to a minimum.

(a)

(b)

Figure 6.12: Inverse design of an elephant. The lower panels show the elephant
with its approximation in blue and the transmission predicted for the obtained
design parameters in red. The upper panels show the evolution of the design
parameters. The MSE on the transmissions during training is also plotted. a)
Inverse design on $F$ and $\delta_{0}$. We converge to a minimum. b) Inverse
design on $\theta$, $n$ and $l$. The design parameters converge to a minimum.

##  Chapter 7 Closing remarks

###  7.1 Main results

We successfully trained different neural networks to learn the transmission of
Fabry-Pérot. We showed that we can predict the wavelength-dependent
transmission. We also succeeded in making a neural network predicting the
angle-dependent transmission. This could be of use in the field of optical
computing. Finally, we showed that the parameters $F$ and $\delta_{0}$ that
appear in the analytical expression for the transmission can improve the
predictions on the wave-length transmissions. The best networks we created for
these three cases had a mean absolute error (MAE) of about 1% or less.

Based on the observation that it is very useful to know $F$ and $\delta_{0}$,
we wanted to recover them in a latent representation using a $\beta$-VAE. We
created a $\beta$-VAE that was able to make good reconstructions of the
transmissions $T(\lambda)$. Upon closer inspection, we saw that we could
retrieve a transformation of $F$ in one of these dimensions. The parameter
$\delta_{0}$ was however spread across multiple latent dimensions. The
correlation between $\delta_{0}$ and the latent dimensions also showed a non-
trivial pattern.

At last, we used the neural networks of chapter 4 to perform inverse design.
We tried to recover the design parameters from which a transmission
$T(\lambda)$ was made. We first tried gradient-descent with the MSE loss on
the transmission. We got stuck in local minima because the MSE loss on the
transmission is not a convex function of $\delta_{0}$. We then used the
Fourier transform of the transmissions to get us to the global minimum. This
allowed us to perform inverse design on all transmissions from the test set.
We were also able to design both $F$ and $\delta_{0}$ and the original design
parameters $\theta$, $n$ and $l$ for arbitrary transmissions $T(\lambda)$.

###  7.2 Limitations of current research

In chapter 4, we already experimented with different hyperparameters. There
are however a lot of hyperparameters that we did not look at. We expect that
bigger networks with more layers and more nodes could perform even better than
the networks we trained. Nonetheless, we did not believe it was very
interesting to further improve our results. We were already quite satisfied
with an MAE of about 1%.

The greatest shortcoming for the $\beta$-VAE was that $\delta_{0}$ was not
stored in a single latent dimension.

For the inverse design, we came up with the loss function of the Fourier
transform since we knew analytically that $\delta_{0}$ existed and was linked
to the oscillation of the transmission. We might not have thought of this if
we did not already know the properties of $\delta_{0}$.

###  7.3 Future research

The prediction of the transmission as well as the inverse design in this
thesis were performed on Fabry-Pérot. This is a resonator with only one layer.
It would be interesting to expand this to multilayer nanophotonic structures.
These are more challenging and new results to design them could be of
practical use. Moreover, it would be interesting to see if we can perform
inverse design without knowing the underlying parameters.

In future work on the $\beta$-VAE, we would like to understand what is
happening to $\delta_{0}$. We see two interesting options to go forward. The
first is to change the dimension of the latent representation and $\beta$.
Changing these parameters has a similar effect, they both control how much
information is stored in the latent representation. We could investigate over
how much dimensions $\delta_{0}$ is spread and if this information is spread
evenly among the dimensions.

The other option we see is to change the loss function. An MSE loss on the
Fourier transform proved to be very useful for inverse design. Since we are
again considering $\delta_{0}$ in the $\beta$-VAE, the Fourier transform might
be useful again.

We see that this thesis only began to explore inverse design in Nanophotonics.
There are still a lot of interesting directions to pursue that continue on
this work.

##  Appendix A Analytical expression Fabry-Pérot

The Fabry-Pérot system consists of a dielectric medium with index of
refraction $n$ in which a light ray resonates several times. The different
outgoing waves interfere, and this leads to a non-trivial total reflection $R$
and transmission $T$ of the system. If we assume that there is no loss in the
system, conservation of energy implies $R+T$ = 1. We can calculate the
resulting transmission $T$ by computing the sum of the different partial waves
$T_{i}$. The set-up of our calculation is shown in figure A.1. The first
outgoing wave is given by

$\displaystyle t_{0}$ $\displaystyle=Te^{i\frac{kl}{\text{cos}\theta_{mat}}},$
(A.1) $\displaystyle=Te^{i\frac{2\pi
n}{\lambda}\frac{l}{\text{cos}\theta_{mat}}}.$

In this expression, $k$ is the wavenumber of the light inside the resonator.
It can be written in function of the wavelength of the light and the angle of
refraction inside the resonator. The angle $\theta_{mat}$ in the equation
above is the angle inside the resonator. This angle is related to the angle
$\theta$ of the incoming light by Snell’s law

$n\text{sin}\theta_{mat}=n_{0}\text{sin}\theta.$ (A.2)

We assume that the medium outside the resonator is air so that we have
$n_{0}=1$. After a round trip in the resonator, a second outgoing wave
emerges, given by

$\displaystyle t_{1}^{\prime}$
$\displaystyle=TRe^{i\frac{3kl}{\text{cos}\theta_{mat}}}.$ (A.3)

We can not simply add the two waves $t_{0}$ and $t_{1}^{\prime}$ together,
since the wave $t_{0}$ has traveled over a distance in the time the light made
a round trip in the resonator. In order to compensate for this effect, we need
to retard the second wave. We get

$\displaystyle t_{1}$
$\displaystyle=TRe^{i\frac{3kl}{\text{cos}\theta_{mat}}-ik_{0}l_{0}}.$ (A.4)

Looking at figure A.1, we see that we have

$\displaystyle k_{0}=\frac{2\pi n_{0}}{\lambda},\hskip
14.22636ptl_{0}=2l\text{tan}\theta_{mat}\ \text{sin}\theta.$ (A.5)

Figure A.1: Fabry-Pérot resonator

In these equations, $k_{0}$ is the wavenumber of the light traveling in air
with $n_{0}=1$. The length $l_{0}$ is the distance the first wave $t_{0}$ will
have travelled in the time the light made a round trip in the resonator. To be
able to add up the different outgoing partial waves, we calculate the phase
difference between them

$\displaystyle\delta$
$\displaystyle=\frac{2kl}{\text{cos}\theta_{mat}}-k_{0}l_{0},$ (A.6)
$\displaystyle=\frac{2kl}{\text{cos}\theta_{mat}}-\frac{2\pi
n_{0}}{\lambda}2l\text{tan}\theta_{mat}\ \text{sin}\theta,$
$\displaystyle=\frac{2kl}{\text{cos}\theta_{mat}}-\frac{2\pi
n}{\lambda}2l\text{tan}\theta_{mat}\ \text{sin}\theta_{mat},$
$\displaystyle=\frac{2kl}{\text{cos}\theta_{mat}}-2kl\frac{\text{sin}^{2}\theta_{mat}}{\text{cos}\theta_{mat}},$
$\displaystyle=\frac{2kl\text{cos}^{2}\theta_{mat}}{\text{cos}\theta_{mat}},$
$\displaystyle=2kl\text{cos}\theta_{mat}.$

We used Snell’s law A.2 to go from the second to the third line. Then we used
trigonometric identities to further simplify the expression. Using the phase
difference between two outgoing waves, we can sum over all of them, we get

$\displaystyle t$ $\displaystyle=\sum_{n=0}^{\infty}TR^{n}e^{in\delta},$ (A.7)
$\displaystyle=\frac{T}{1-Re^{i\delta}}.$

Now that we have an expression for the sum of the partial waves, we can
calculate the total transmitted intensity

$\displaystyle T$
$\displaystyle=t^{\ast}t=\frac{T^{2}}{(1-Re^{i\delta})(1-Re^{-i\delta})},$
(A.8) $\displaystyle=\frac{T^{2}}{1-2Rcos\delta+R^{2}},$
$\displaystyle=\frac{(1-R)^{2}}{1-2R+4Rsin^{2}\frac{\delta}{2}+R^{2}},$
$\displaystyle=\frac{1}{1+\frac{4R}{(1-R)^{2}}sin^{2}\frac{\delta}{2}}.$

We can introduce a new parameter $F$ to simplify this expression. We end up
with

$T=\frac{1}{1+Fsin^{2}(\frac{\delta}{2})},$ (A.9)

where the two parameters are given by

$F=\frac{4R}{(1-R)^{2}},\hskip 28.45274pt\delta=\frac{2\pi
n}{\lambda}2lcos\theta_{mat}.$ (A.10)

The parameter $F$ depends on the reflectance $R$ on the sides of the
resonator. This can be calculated from the indices of refraction on both sides
and the angle under which the light comes in. For TE polarized light we obtain

$R=\Bigg{[}\frac{ncos\theta_{mat}-cos\theta}{ncos\theta_{mat}+cos\theta}\Bigg{]}^{2}.$
(A.11)

We always work with TE polarized light in this work. Equations A.9, A.10 and
A.11 are what we need to analytically compute the transmission of the Fabry-
Pérot resonator.

## References

  * [1] Giuseppe Carleo et al.  “Machine learning and the physical sciences”  In _Reviews of Modern Physics_ 91.4 American Physical Society (APS), 2019  arXiv:1903.10563
  * [2] OpenAI et al.  “Dota 2 with Large Scale Deep Reinforcement Learning”, 2019  arXiv:1912.06680 [cs.LG]
  * [3] Juan Carrasquilla and Roger G. Melko  “Machine learning phases of matter”  In _Nature Physics_ 13.5 Springer ScienceBusiness Media LLC, 2017, pp. 431–434  arXiv:1605.01735
  * [4] Nicholas Choma et al.  “Graph Neural Networks for IceCube Signal Classification”  In _CoRR_ abs/1809.06166, 2018  arXiv:1809.06166
  * [5] Dalya Baron  “Machine Learning in Astronomy: a practical overview”, 2019  arXiv:1904.07248 [astro-ph.IM]
  * [6] Dan Guest, Kyle Cranmer and Daniel Whiteson  “Deep Learning and Its Application to LHC Physics”  In _Annual Review of Nuclear and Particle Science_ 68.1 Annual Reviews, 2018, pp. 161–181  arXiv:1806.11484
  * [7] M.. al. Spühler  “A very short planar silica spot-size converter using a nonperiodic segmented waveguide.”  In _Light. Technol._ , 1998, pp. 1680–1685
  * [8] D.. Dobson and S.. Cox  “Maximizing band gaps in two-dimensional photonic crystals”  In _Appl. Math._ , 1999, pp. 2108–2120
  * [9] Alexander Y Piggott et al.  “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer”  In _Nature Photonics_ 9.6, 2015, pp. 374–377  arXiv:1504.00095
  * [10] Haejun Chung and Owen D. Miller  “Tunable metasurface inverse design for 80% switching efficiencies and 144∘ angular steering”, 2019  arXiv:1910.03132 [physics.optics]
  * [11] Logan Su et al.  “Inverse Design and Demonstration of a Compact on-Chip Narrowband Three-Channel Wavelength Demultiplexer”  In _ACS Photonics_ 5.2, 2018, pp. 301–305  arXiv:1709.08809
  * [12] Sean Molesky et al.  “Inverse design in nanophotonics”  In _Nature Photonics_ 12, 2018, pp. 659–670  arXiv:1801.06715
  * [13] John Peurifoy et al.  “Nanophotonic Particle Simulation and Inverse Design Using Artificial Neural Networks”  In _Science Advances_ 4.6 American Association for the Advancement of Science, 2017  arXiv:1712.03222
  * [14] Sunae So, Jungho Mun and Junsuk Rho  “Simultaneous inverse design of materials and parameters of core-shell nanoparticle via deep-learning: Demonstration of dipole resonance engineering”, 2019  arXiv:1904.02848 [physics.optics]
  * [15] Dianjing Liu, Yixuan Tan, Erfan Khoram and Zongfu Yu  “Training deep neural networks for the inverse design of nanophotonic structures”  In _ACS Photonics_ 5.4 American Chemical Society, 2017, pp. 1365–1369  arXiv:1710.04724
  * [16] Yashar Kiarashinejad, Sajjad Abdollahramezani and Ali Adibi  “Deep learning approach based on dimensionality reduction for designing electromagnetic nanostructures”  In _npj Computational Materials_ 6.1 Nature Research, 2020  arXiv:1902.03865
  * [17] Jiaqi Jiang et al.  “Free-form diffractive metagrating design based on generative adversarial networks”  In _ACS Nano_ 13.8, 2019, pp. 8872–8878  arXiv:1811.12436
  * [18] Frank F. Rosenblatt  “The perceptron: a probabilistic model for information storage and organization in the brain.”  In _Psychological review_ 65 6, 1958, pp. 386–408
  * [19] David E Rumelhart, Geoffrey E Hinton and Ronald J Williams  “Learning representations by back-propagating errors”  In _Nature_ 323, 1986, pp. 533–536
  * [20] Alex Krizhevsky, Ilya Sutskever and Geoffrey Hinton  “ImageNet Classification with Deep Convolutional Neural Networks”  In _Neural Information Processing Systems_ 25, 2012  DOI: 10.1145/3065386
  * [21] Pankaj Mehta et al.  “A high-bias, low-variance introduction to Machine Learning for physicists”  In _Physics Reports_ 810 Elsevier BV, 2019, pp. 1–124  arXiv:1803.08823
  * [22] Christopher M. Bishop  “Pattern Recognition and Machine Learning (Information Science and Statistics)”  Berlin, Heidelberg: Springer-Verlag, 2006
  * [23] Ian J. Goodfellow, Yoshua Bengio and Aaron Courville  “Deep Learning” http://www.deeplearningbook.org Cambridge, MA, USA: MIT Press, 2006
  * [24] Michael A. Nielsen  “Neural networks and deep learning”  Determination Press, 2015
  * [25] Y. Cun et al.  “Handwritten Digit Recognition with a Back-Propagation Network”  In _Advances in Neural Information Processing Systems 2_ San Francisco, CA, USA: Morgan Kaufmann Publishers Inc., 1990, pp. 396–404
  * [26] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun  “Deep Residual Learning for Image Recognition”  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ IEEE, 2016  arXiv:1512.03385
  * [27] Chigozie Nwankpa, Winifred Ijomah, Anthony Gachagan and Stephen Marshall  “Activation Functions: Comparison of trends in Practice and Research for Deep Learning”, 2018  arXiv:1811.03378 [cs.LG]
  * [28] Prajit Ramachandran, Barret Zoph and Quoc V. Le  “Searching for Activation Functions”, 2017  arXiv:1710.05941 [cs.NE]
  * [29] Sebastian Ruder  “An overview of gradient descent optimization algorithms”, 2016  arXiv: http://arxiv.org/abs/1609.04747
  * [30] Diederik P. Kingma and Jimmy Lei Ba  “Adam: A method for stochastic optimization”  In _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ International Conference on Learning Representations, ICLR, 2015  arXiv:1412.6980
  * [31] Geoffry.. Hinton and Ruslan R. Salakhutdinov  “Reducing the dimensionality of data with neural networks”  In _Science_ American Association for the Advancement of Science, 2006, pp. 504–507
  * [32] Lilian Weng  “From Autoencoder to Beta-VAE”  In _lilianweng.github.io_ , 2018  URL: http://lilianweng.github.io/lil-log/2018/08/12/from-autoencoder-to-beta-vae.html
  * [33] Dmytro Bondarenko and Polina Feldmann  “Quantum Autoencoders to Denoise Quantum Data”  In _Physical Review Letters_ 124 American Physical Society (APS), 2020  arXiv:1910.09169
  * [34] Diederik P Kingma and Max Welling  “Auto-Encoding Variational Bayes”, 2013  arXiv:1312.6114 [stat.ML]
  * [35] Irhum Shafkat  “Intuitively Understanding Variational Autoencoders”  medium, 2018  URL: https://towardsdatascience.com/intuitively-understanding-variational-autoencoders-1bfe67eb5daf
  * [36] Raban Iten et al.  “Discovering Physical Concepts with Neural Networks”  In _Physical Review Letters_ 124.1 American Physical Society, 2020  arXiv:1807.10300
  * [37] François Chollet  “Keras”, https://keras.io, 2015
  * [38] Pierre Ambs  “Optical Computing: A 60-Year Adventure”  In _Advances in Optical Technologies_ 2010 Hindawi Publishing Corporation, 2010, pp. 372652  URL: https://doi.org/10.1155/2010/372652