# Mean-squared-error-based adaptive estimation of pure quantum states and
unitary transformations

A. Rojas<EMAIL_ADDRESS>Instituto Milenio de Investigación en Óptica y
Departamento de Física, Universidad Concepción, casilla 160-C, Concepción,
Chile  L. Pereira  Instituto de Física Fundamental IFF-CSIC, Calle Serrano
113b, Madrid 28006, Spain.
S. Niklitschek  Facultad de Ciencias Físicas y Matemáticas, Departamento de
Estadística, Universidad de Concepción, casilla 160-C, Concepción, Chile  A.
Delgado  Instituto Milenio de Investigación en Óptica y Departamento de
Física, Universidad Concepción, casilla 160-C, Concepción, Chile

###### Abstract

In this article we propose a method to estimate with high accuracy pure
quantum states of a single qudit. Our method is based on the minimization of
the squared error between the complex probability amplitudes of the unknown
state and its estimate. We show by means of numerical experiments that the
estimation accuracy of the present method, which is given by the expectation
of the squared error on the sample space of estimates, is state independent.
Furthermore, the estimation accuracy delivered by our method is close to twice
the Gill-Massar lower bound, which represents the best achievable accuracy,
for all inspected dimensions. The minimization problem is solved via the
concatenation of the Complex simultaneous perturbation approximation, an
iterative stochastic optimization method that works within the field of the
complex numbers, and Maximum likelihood estimation, a well-known statistical
inference method. This can be carried out with the help of a multi-arm
interferometric array. In the case of a single qubit, a Mach-Zehnder
interferometer suffices. We also show that our estimation procedure can be
easily extended to estimate unknown unitary transformations acting on a single
qudit. Thereby, the estimation of unitary transformations achieves a higher
accuracy than that achieved by processes based on tomographic methods for
mixed states.

###### pacs:

Valid PACS appear here

##  I Introduction

Recently, the estimation of unknown quantum states has been studied from the
point of view of the estimation accuracy achievable by means of an ensemble of
$N_{T}$ identically prepared copies of the state to be estimated [1, 2, 3, 4,
5, 6, 7]. The ultime mixed-state estimation accuracy is given by the Gill-
Massar lower bound [8], which establishes the highest possible accuracy
achievable by means of separable measurements on the members of the ensemble.
For instance, the mean value $\bar{I}$ of the Uhlmann-Josza infidelity
$I=(Tr\sqrt{\sqrt{\rho}\tilde{\rho}\sqrt{\rho}})^{2}$ [9, 10] on the set of
estimates $\tilde{\rho}$ of the unknown state $\rho$ can be employed as a
metric for the estimation accuracy. In this case the Gill-Massar lower bound
becomes $\bar{I}\geq I_{GM}^{(mixed)}=(d^{2}-1)(d+1)/4N_{T}$ [11].

Adaptive two-stage standard quantum tomography saturates $I_{GM}^{(mixed)}$
for a single qubit [3]. In a first stage, standard quantum tomography [12, 13]
is carried out on an ensemble of size $N_{0}$ and a first estimate
$\tilde{\rho}_{0}$ is inferred. The eigenbasis of $\tilde{\rho}_{0}$ is used
to construct three new mutually unbiased bases, which are then employed to
perform a second stage of standard quantum tomography on an ensemble of size
$N_{T}-N_{0}$. Thereafter, the acquired data is post-processed via maximum
likelihood estimation [14, 15] to obtain a final estimate $\tilde{\rho}$. This
procedure clearly requires the capacity to adapt the measurement bases and
doubles the total number of measurement outcomes. Unfortunately, adaptive two-
stage standard quantum tomography departs from the Gill-Massar lower bound in
the case of a single qudit with $d>2$ [4].

In the particular case of pure states it has been shown that a much better
accuracy can be obtained. In this case, the Gill-Massar lower bound for the
infidelity is $I_{GM}^{(pure)}=(d-1)/N_{T}$. The 5-bases based quantum
tomographic method [6] produces an infidelity that lays in between
$I_{GM}^{(mixed)}$ and $I_{GM}^{(pure)}$ [7]. This method employs an adaptive
scheme where measurements on the canonical base are employed to define four
new measurement bases. The five bases determine univocally any pure state of a
single qudit and allow to certify the purity assumption.

An estimation accuracy closer to $I_{GM}^{(pure)}$ can be achieved by means of
formulating the problem of quantum state determination as an optimization
problem [16] and solving it by means of a combination of stochastic
optimization on the field of the complex numbers and maximum likelihood
estimation (MLE) [17, 18]. In this approach the infidelity is considered a
real function of complex arguments where the unknown state plays the role of a
set of fixed and unknown complex parameters. This function is optimized by
means of the Complex simultaneous stochastic approximation (CSPSA) method,
which allows to handle non-holomorphic functions with unknown parameters. This
optimization method requires the measurement of the infidelity at each
iteration. The information provided by the sequence of measurements can be
employed to enhance the rate of convergence of the optimization method when
processed via maximum likelihood estimation.

Here, we study the estimation of pure quantum states using the mean-squared
error (MSE) as figure of merit for the accuracy. It has been theoretically
proven and experimentally demonstrated that states with a small infidelity
might lead to very different physical properties [19, 20, 21]. Consequently,
the infidelity might turn to be inadequate to assess the estimation of high-
dimensional quantum systems. Therefore, it is advisable to explore other
accuracy metrics. We resort to the mean-squared error mainly because it can be
inferred from experimentally acquired data, it is inexpensive to compute, it
is an excellent metric in the context of optimization, and it is a desirable
measure in statistics and estimation theory [22]. We first show that the
squared error (SE) between the probability amplitudes of two pure quantum
states of a single qudit can be measured by means of a multi-arm
interferometric array. This allows us to employ the CSPSA method to optimize
SE. This iterative method and MLE are then combined to drive a sequence of
measurements in such a way that the SE rapidly decreases at each iteration.
Due to the intrinsic stochasticity of CSPSA, the estimation procedure
generates, for a fixed unknown state, a sample of estimates. Via Monte Carlo
numerical experiments we show that the accuracy of the estimation procedure,
that is, the mean of SE on the sample of estimates (or MSE), is nearly
independent of the state to be estimated. Moreover, mean and median of SE
agree on the sample of estimates, which indicates a symmetric distribution of
estimates without outliers. Numerical simulations show that after a few
iterations the estimation of unknown states enters into an asymptotic regime
that follows very closely twice the Gill-Massar bound for the MSE.

We also apply our previous result to the estimation of unknown unitary
transformations. This is an important application since the successful
realization of quantum information protocols and quantum devices requires the
use of efficient characterization tools. Among these, the most widely employed
is Quantum process tomography (QPT). This is based on a selection of probe
states that undergo the process to be estimated, which is followed by the
estimation of the states generated by the process. QPT has been applied to
multi-qubit processors [23], quantum communications channels [45], coherent
transport in biological mechanisms [25], ion traps [26], nuclear magnetic
resonance [27], superconducting circuits [28], nitrogen-vacancy color centers
[29], and few-photon linear-optical systems [30, 31, 32].

A quantum process is described as a completely-positive, trace preserving
(CPTP) map, which requires $d^{4}-d^{2}$ real numbers to be completely
characterized. If we know, however, that the process is unitary, then the
number of independent parameters can be further reduced. For instance,
$d^{2}+d$ measurement outcomes are necessary to distinguish among unitary
transformations [53]. Unitarity can be certified, for example, through
randomized benchmarking.

We can try to estimate a unitary transformation $U$ with the help of the
Uhlmann-Josza fidelity by applying $U$ to the elements of an orthonormal base
an estimating the generated states. This procedure, however, does not allow
obtaining a set of $d$ phases of the unitary transformation because the
Uhlmann-Josza fidelity is insensitive to global phases. Thereby, it is
necessary to increase the number of states onto which the unitary
transformation acts in order to obtain the missing complex phases. Thus, the
use of the Uhlmann-Josza fidelity for estimating unitary transformations
becomes equivalent to the problem of estimating $d$ pure states and $d$
unknown phases. Alternatively, it is possible to define an infidelity-guided
figure of merit to compare two unitary transformations $U_{1}$ and $U_{2}$.
This is given by $I(U_{1},U_{2})=1-|Tr(U_{1}U_{2}^{\dagger})|^{2}/d^{2}$ [34].
This can be shown to be equal to
$I(U_{1},U_{2})=1-|\langle\psi|(U^{\dagger}_{1}\otimes{\cal
I})(U_{2}\otimes{\cal I})|\psi\rangle|^{2}$, where $|\psi\rangle$ is the two-
qudit maximally entangled state
$\sum_{k}(1/\sqrt{d})|k\rangle\otimes|k\rangle$ [35]. Therefore, the
measurement of $I(U_{1},U_{2})$ requires the capability of preparing maximally
entangled two-qudit states and to project onto two-qudit states, which
separates the estimation of states from that of unitary transformations.
Instead, we apply our results on the estimation of pure states to the
estimation of unitary transformations via the optimization of SE. This allows
us to handle the estimation of states and unitary transformations within the
same theoretical framework, avoid the use of maximally entangled two-qudit
states, avoid increasing the number of measurements, and improve the
estimation accuracy given a fixed number of particles interacting with the
unknown unitary process. This is important, for instance, for measuring
biological samples [36] and materials [37] in scenarios where the number of
samples (photons) must be low to avoid sample damage. Our method estimates the
columns of an unknown unitary transformation separately, which after post-
processing leads to an estimation accuracy for the unitary transformations
close to $2d$ times the Gill-Massar bound for the MSE of a single unknown pure
state. Let us note that our estimation method for states and unitary
transformations, unlike the recent proposals [38], measures all photons in the
ensamble size independently, that is, no entanglement is used between photons.
Also, most methods employed to estimate unitary transformations employ
tomographic methods that are designed to estimate mixed states, which has an
estimation accuracy limited by the Gill-Massar lower bound for mixed states.
Thus, our estimation method for pure states leads to a clear improvement on
the estimation accuracy of unitary transformations.

##  II Method

The optimization of a real-valued function $f({\bm{x}})$ of a vector
${\bm{x}}$ of $d$ real variables can be implemented by means of the gradient
descent method. This is based on the iterative rule

$\bm{x}_{k+1}={\bm{x}}_{k}-a_{k}{\bm{g}}(x_{k}),$ (1)

where the function ${\bm{g}}_{k}$ is the gradient $\nabla f$ and $a_{k}$ is a
real gain coefficient. We are interested in the optimization of a real-valued
function $f({\bm{z}},{\bm{z}}^{*})$ with ${\bm{z}}$ a vector of $n$ complex
variables, where the function $f({\bm{z}},{\bm{z}}^{*})$ also depends on a set
of unknown fix parameters. In order to optimize $f({\bm{z}},{\bm{z}}^{*})$ we
might attempt to employ the iterative rule

$\bm{z}_{k+1}={\bm{z}}_{k}-a_{k}{\bm{g}}(\bm{z}_{k}).$ (2)

However, a few changes are needed. Real-valued functions of complex variables
violate the Cauchy-Riemann conditions, and consequently, they lack standard
complex derivative. This can be solved with the notion of Wirtinger
derivatives, which are defined by [39]

$\partial_{z_{i}}=\frac{1}{2}(\partial_{x_{i}}-i\partial_{y_{i}})~{}{\rm
and}~{}\partial_{z_{i}^{*}}=\frac{1}{2}(\partial_{x_{i}}+i\partial_{y_{i}}),$
(3)

where $x_{i}$ and $y_{i}$ are the real and imaginary parts of $z_{i}$,
respectively. These derivatives exist even if the function $f$ violates the
Cauchy-Riemann conditions. Minima and maxima of a real-valued function of
complex variables are completely characterized by the conditions
$\partial_{z_{i}^{*}}f=0~{}\forall~{}i=1,\dots,d$ or, equivalently,
$\partial_{z_{i}}f=0~{}\forall~{}i=1,\dots,d$ [40, 41, 42]. Thereby, we can
make the identification ${\bm{g}}=\partial_{{\bm{z}}^{*}}f$ with
$\partial_{{\bm{z}}^{*}}=(\partial_{z_{1}^{*}},\dots,\partial_{z_{d}^{*}})$.
Unfortunately, the function we seek to optimize depends on a set of unknown
fix parameters. Thus, even when we might be able to obtain an analytical
expression for the complex-valued gradient, we cannot evaluate it. In this
scenario the components of the complex-valued gradient are approximated as

${g}_{k,i}=\frac{f({\bm{z}}_{k+},{\bm{z}}_{k+}^{*})+\epsilon_{k,+}-f({\bm{z}}_{k-},{\bm{z}}_{k-}^{*})-\epsilon_{k,-}}{2c_{k}{\Delta}_{k,i}^{*}}.$
(4)

This expression can be evaluated as long as we have access to the values of
the function $f$ at the points ${\bm{z}}_{k\pm}={\bm{z}}_{k}\pm
c_{k}{\bm{\Delta}}_{k}$, where $c_{k}$ is a real gain coefficient and the
complex components ${\Delta}_{k,i}$ of vector ${\bm{\Delta}}_{k}$ provide a
direction for the approximation of the complex-valued gradient at each
iteration. The terms $\epsilon_{k,\pm}$ entering in Eq. (4) represent noise
affecting the evaluations, or measurements, of the function $f$.

Equations (2) and (4) are the core of the CSPSA method [17]. This allows us to
optimize real-valued functions of complex variables when the complex-valued
gradient cannot be evaluated. It is possible to solve this optimization
problem by resorting to a real parametrization of the complex variables.
However, it has been shown that CSPSA, which works on the field of the complex
numbers, provides a higher rate of convergence in the estimation of pure
states in comparison to similar methods working on the field of the real
numbers. It has been proven that, under suitable conditions, the sequence
${\bm{z}}_{k}$ of complex vectors generated by CSPSA converges in mean to the
minimizer of the function $f$. Furthermore, the approximation $\bm{g}$ of the
complex-valued gradient is asymptotically unbiased.

The gain coefficients $a_{k}$ and $c_{k}$ entering in Eqs. (2) and (4) and are
given by

$a_{k}=\frac{a}{(k+1+A)^{s}},~{}~{}c_{k}=\frac{b}{(k+1)^{r}}.$ (5)

These control the rate of convergence of CSPSA. The values of $a,A,s,b$ and
$r$ are tuned to increase the rate of convergence for each function $f$. The
complex coefficients $\Delta_{k,i}$ are independently and identically
generated as elements in the set $\\{\pm 1,\pm i\\}$ with uniform probability.
This choice leads to a boost in the rate of convergence of the algorithm.
These particular choices of gains and vector $\Delta_{k}$ obey the convergence
conditions of CSPSA.

Let us note that the CSPSA method requires twice the evaluation of the
function $f$ at each iteration. These values are employed to produce a new
estimate and are thereafter discarded. This procedure thus generate a large
amount of information about the function $f$ which can be also employed to
accelerate the convergence of CSPSA [18]. For instance, this data can be
processed via Maximum likelihood estimation (MLE). Given measurement settings
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}S}={m_{1},\dots,m_{M}}$
and corresponding data $D={d_{1},\dots,d_{M}}$, maximum likelihood estimation
seeks for a physical state $\rho$ that maximizes the likelihood functional

$P(D|\rho,S)=\Pi^{M}_{i=1}P(d_{i}|\rho,m_{i}),$ (6)

with $P(d_{i}|\rho,m_{i})$ the probability to obtain the measurement outcome
$d_{i}$ given the state $\rho$ and measurement setting $m_{i}$. The
optimization of the likelihood functional $P(D|\rho,S)$ requires an initial
guess. This is chosen as the estimate provided by CSPSA at each iteration. The
solution provided by MLE is then employed as the input estimate for the next
iteration with CSPSA. Let us note that set $MS$ contains all measurements
carried out until the $k$-th iteration. We employ the logarithmic multinomial
likelihood functional, which is optimized employing the Nelder-Mead (or
Simplex) algorithm.

To benchmark the estimation accuracy achieved by our method we will employ the
Gill-Massar lower bound for the mean-squared error. This is a fundamental
limit for the accuracy achieved in an estimation procedure. Let us consider a
real-valued metric $f$ for the accuracy that is a function of the covariance
matrix $\cal{C}$, the coefficients of which are given by

$\mathcal{C}_{ij}(\theta)=\mathbb{E}[(\tilde{\theta}_{i}-\theta_{i})(\tilde{\theta}_{j}-\theta_{j})|\tilde{\theta}],$
(7)

where $\theta$ is an unknown vector and $\tilde{\theta}$ its estimate. The
bound is the solution of the optimization problem

$\min f(\mathcal{C})$ (8)

with the constraint
$\text{tr}(\mathcal{J}^{-1}\mathcal{I}^{N_{T}})=N_{T}(d-1)$, where $N_{T}$ is
the total number of identical copies of the unknown state, $d$ is the
dimension of the state, and $\cal{J}$ and $\cal{I}$ are the quantum and
classical Fisher matrices of the whole ensemble, respectively. Then components
of the quantum Fisher matrix are given by

${\cal J}_{j,k}=\frac{1}{2}Tr(\rho[L_{j},L_{k}]),$ (9)

where where $L_{j}$ is the symmetric logarithm derivative with respect to the
parameters $S_{k}$ that define the state $\rho$. This derivative is implicitly
defined by $\partial\rho/\partial S_{j}=\\{\rho,L_{j}\\}/2$. The components of
the classical Fisher information matrix are given by

${\cal I}_{j,k}=\sum_{m}(1/p_{m})(\partial p_{m}/\partial S_{j})(\partial
p_{m}/\partial S_{k}).$ (10)

The covariance and classical and quantum Fisher matrices are through the
classical and quantum Cramer-Rao inequalities [43] ${\cal C}\geq{\cal
I}^{-1}\geq{\cal J}^{-1}$, which establish a lower limit for the Covariance
matrix. In general, it is still unknown whether or not this inequality can be
saturated. In the case of separable measurements, that is, measurements
carried out on each individual copy of the ensemble, the quantum Cramer-Rao
bound can not be attained and the fundamental limit is given by the Gill-
Massar lower bound [11]

If the function $f$ has the form $f(\cal{C})=\text{tr}(\cal{W}\cal{C})$ with
$\cal{W}$ being a weight matrix, then the minimum is given by [3]

$\text{tr}(\mathcal{WC})=\frac{\left[\text{tr}\left(\sqrt{\mathcal{J}^{-1/2}\mathcal{W}\mathcal{J}^{-1/2}}\right)\right]^{2}}{N_{T}(d-1)}.$
(11)

In the case of the mean-squared error we choose the weight matrix as the
identity, in which case we can identify the $MSE$ with $\text{tr}(\cal{C})$.
For a pure state the matrix $\cal{J}$ is written as [44]

$\mathcal{J}=4\mathbb{I}_{2(d-1)\times 2(d-1)}.$ (12)

Thereby, the bound yields

$MSE\geq MSE_{GM}(d,N_{T})=\frac{d-1}{N_{T}}.$ (13)

##  III Estimation of pure states via the optimization of the squared error

Our main aim is to obtain an estimate
$|\tilde{\psi}\rangle=\sum_{i}\tilde{z}_{i}|i\rangle$ of an unknown quantum
state $|\psi\rangle=\sum_{i}z_{i}|i\rangle$. In order to do this we consider
the squared error

$SE(\bm{z},\bm{\tilde{z}})=\sum_{i=1}^{d}|z_{i}-\tilde{z}_{i}|^{2},$ (14)

which is a function of the probability amplitudes $\tilde{z}_{i}$ of
$|\tilde{\psi}\rangle$. The probability amplitudes $z_{i}$ of $|\psi\rangle$
play the role of unknown fixed parameters. The unknown state $|\psi\rangle$
can be characterized as

$\bm{z}={\rm Arg}\\{\min_{\bm{\tilde{z}}}SE(\bm{z},\bm{\tilde{z}})\\},$ (15)

that is, SE achieves a global minimum when $\bm{\tilde{z}}=\bm{z}$.

To solve the minimization of $SE(\bm{z},\bm{\tilde{z}})$ we employ the CSPSA
method. This requires the capability to obtain the value of
$SE(\bm{z},\bm{\tilde{z}})$ for all $\bm{\tilde{z}}$. In the simplest case of
estimating the polarization state of a single photon, a Mach-Zehnder
interferometer whose arms are supplemented with unitary transformations acting
on the polarization degree of freedom allows us to infer the value of
$SE(\bm{z},\bm{\tilde{z}})$. In this setup, the initial state
$|\psi\rangle_{in}$ of a single photon before entering the interferometer is
given by $|\psi\rangle_{in}=|h\rangle_{1}$, where $|h\rangle_{1}$ describes a
horizontally polarized single photon. After the interaction of the photon with
the first beam splitter, the quantum state of the photon becomes
$(|h\rangle_{a}+|h\rangle_{b})/\sqrt{2}$, which corresponds to an equally
weighted coherent superposition of the two possible propagation paths $a$ and
$b$ for the photon inside the interferometer. The polarization state of the
photon changes conditional on the path, that is,

$U_{a}|h\rangle_{a}=z_{h}|h\rangle_{a}+z_{v}|v\rangle_{a}=|\psi\rangle,$ (16)

which is the unknown state to be estimated, and

$U_{b}|h\rangle_{b}=\tilde{z}_{h}|h\rangle_{b}+\tilde{z}_{v}|v\rangle_{b}=|\tilde{\psi}\rangle,$
(17)

which is the estimate of $|\psi\rangle$. Thereby, the state of the photon
before the second beam splitter becomes
$[(z_{h}|h\rangle_{a}+z_{v}|v\rangle_{a})+(\tilde{z}_{h}|h\rangle_{b}+\tilde{z}_{v}|v\rangle_{b})]/\sqrt{2}$.
After the second beam splitter, the state is given by

$\displaystyle|\psi\rangle_{out}$ $\displaystyle=$
$\displaystyle\frac{1}{2}[(z_{h}+\tilde{z}_{h})|h\rangle_{1}+(z_{v}+\tilde{z}_{v})|v\rangle_{1}]$
(18) $\displaystyle+$
$\displaystyle\frac{1}{2}[(z_{h}-\tilde{z}_{h})|h\rangle_{2}+(z_{v}-\tilde{z}_{v})|v\rangle_{2}],$

where the subindexes 1 and 2 indicates the output ports of the interferometer.
The probability $P_{2}$ of detecting a photon at output port 2 is

$P_{2}=\frac{1}{4}(|z_{h}-\tilde{z}_{h}|^{2}+|z_{v}-\tilde{z}_{v}|^{2}),$ (19)

which can be identified with the squared error of the complex probability
amplitudes as

$SE(\bm{z},\bm{\tilde{z}})=4P_{2}.$ (20)

Thus, in the setup above described the unitary transformation $U_{a}$ is
employed to create the unknown polarization state $|\psi\rangle$ defined by
the pair of complex probability amplitudes $(z_{h},z_{v})$. The unitary
transformation $U_{b}$ is employed to generate an estimate
$|\tilde{\psi}\rangle$, which is defined by the pair of complex probability
amplitudes $(\tilde{z}_{h},\tilde{z}_{v})$. Equation (20) indicates that the
transformation $U_{b}$ has to be changed in such a way that no photon is
detected at output port 2, in which case $|\tilde{\psi}\rangle=|\psi\rangle$.
We employ CSPSA and MLE to drive the sequence of choices of $U_{b}$ toward the
unknown state. Let us note that if the output ports of the interferometer are
supplemented with polarizing beam splitters and single-photon detectors it is
possible to measure independently the four combinations of coefficients
$|z_{h}-\tilde{z}_{h}|^{2}$, $|z_{v}-\tilde{z}_{v}|^{2}$,
$|z_{h}+\tilde{z}_{h}|^{2}$, and $|z_{v}+\tilde{z}_{v}|^{2}$.

The case of higher dimensions can be realized by considering a spatial qudit,
that is, a qudit encoded in the propagation paths of a single photon. The
initial state of the qudit is given by $|k\rangle$, where the state
$|k\rangle$ describes a single photon propagating along one of several
distinguishable paths $k=1,\dots,d$. On path $k$ a beam splitter transforms
the state $|k\rangle$ into the superposition
$(1/\sqrt{2})(|k_{1}\rangle+|k_{2}\rangle)$, where the subindexes $i=1,2$
distinguish the propagations paths at the exit ports of the beam splitter.
Thereafter, on paths $k_{1}$ and $k_{2}$ the unitary transformations $U$ and
$\tilde{U}$ are applied, respectively. These transformations create a
superposition of path states, that is,
$U|k_{1}\rangle=\sum_{k}z_{k}|k_{1}\rangle$ and
$\tilde{U}|k_{2}\rangle=\sum_{k}\tilde{z}_{k}|k_{2}\rangle$. This leads to the
state $(1/\sqrt{2})(U|k_{1}\rangle+\tilde{U}|k_{2}\rangle)$. Finally, paths
$k_{1}$ and $k_{2}$ for each $k=1,\dots,d$ are merged together by beam
splitters, which leads to the state
$(1/2)[\sum_{k}(z_{k}+\tilde{z}_{k})|k_{1^{\prime}}\rangle+\sum_{k}(z_{k}-\tilde{z}_{k})|k_{2^{\prime}}\rangle)]$
with $i=1^{\prime},2^{\prime}$ the output ports of each beam splitter. The
probability of detecting a photon on any path $k_{2^{\prime}}$ is given by

$P_{2}=\frac{1}{4}\sum_{k=1}^{d}|z_{k}-\tilde{z}_{k}|^{2}.$ (21)

Thereby, we have that

$SE(\bm{z},\bm{\tilde{z}})=4P_{2},$ (22)

which generalizes Eq. (20) to the case $d>2$.

We have formulated our proposal to measure the SE in terms of bulk optics
based setup. However, this proposal can easily be translated to other
experimental platforms, such as, for instance, integrated quantum photonics
[45] and space-division multiplexing optical fibres [46, 47]. On this
platforms unitary transformations can be implemented by means of sequences
composed of beam splitters and controlled phase transformations [48, 49].

The main steps of the MSE-based method for estimating pure quantum states are
summarized as pseudocode in Algorithm 1 above, where we have considered the
proposals to estimate the SE with the help of a multi-arm interferometer. An
implementation of the pseudocode in the Python programming language can be
found in the GitHub repository [50].

The minimization of SE via CSPSA requires an initial guess (or estimate) of
the unknown state. Since no a priori information about the unknown state is
available, the initial guess is also generated according to a uniform
distribution. At each iteration, CSPSA generates the $\bm{\Delta}$ vector the
components of which are randomly chosen. Also, at each iteration CSPSA uses
the value of SE on the states $\bm{z}_{k\pm}$. According to Eq. (20), the
value of the SE can be inferred from a probability, which requires an ensemble
of $N$ independently and identically prepared copies. Thereby, the total
number of copies employed after $k$ iterations of CSPSA is given by
$N_{T}=2Nk$. Since this ensemble is finite, the value of SE will be affected
by finite statistics effects. Thus, the estimation process for a fixed unknown
state has three sources of randomness: the choice of the initial guess, the
choice of the $\bm{\Delta}$ vector, and the measurement process of SE.
Thereby, each time that CSPSA is employed to obtain an estimate of a fixed
unknown state $\bm{z}$, a different estimate $\bm{\tilde{z}}$ is generated. In
this scenario the estimation accuracy for a fixed state $\bm{z}$ is given by

$MSE(\bm{z})=\mathbb{E}[SE(\bm{z},\bm{\tilde{z}})|\bm{\tilde{z}}],$ (23)

where the expectation is calculated over the set of all possible estimates
$\bm{\tilde{z}}$ of $\bm{z}$. The mean-squared error can also be cast as

$MSE(\bm{z})=\int\tilde{p}(\bm{\tilde{z}})SE(\bm{z},\bm{\tilde{z}})\bm{d\tilde{z}},$
(24)

where $\tilde{p}(\bm{\tilde{z}})$ probability density function of obtaining
the estimate $\bm{\tilde{z}}$ that characterizes the estimation procedure.

Algorithm 1 MSE-based estimation of pure states

1:Consider a known pure state $|\psi\rangle$, which is prepared on the upper
arm of the interferometer by means of the transformation $U$.

2:Choose an initial guess $|\tilde{\psi}_{0}\rangle$ and define
$\tilde{z}_{0,i}=\langle i|\tilde{\psi}_{0}\rangle$.

3:Set gain coefficients $a$, $A$, $s$, $b$ and $r$.

4:for $k=1,\ldots,k_{max}$ do

5: Set  $a_{k}=\frac{a}{(k+1+A)^{s}},\quad c_{k}=\frac{b}{(k+1)^{r}}.$

6: Choose $\Delta_{k,i}$ randomly in the set $\\{\pm 1,\pm i\\}$.

7: Calculate
$|\psi_{k\pm}\rangle=\sum_{i}\tilde{z}_{k\pm,i}|i\rangle/|\tilde{\bm{z}}_{k\pm}|$,
with $\tilde{\bm{z}}_{k\pm}=\tilde{\bm{z}}_{k}\pm c_{k}{\bm{\Delta}_{k}}$.

8: Prepare the states $|\psi_{k\pm}\rangle$ on the lower arm of the
interferometer by means of the transformation $U^{\prime}$

9: Estimate experimentally the square errors
$SE({\bm{z}},\tilde{{\bm{z}}}_{k\pm})$ with a sample of size $N$.

10: Estimate the gradient as
$\tilde{g}_{k,i}=\frac{SE({\bm{z}},\tilde{\bm{z}}_{k+})-SE({\bm{z}},\tilde{\bm{z}}_{k-})}{2c_{k}\Delta_{k,i}^{*}}.$

11: Actualize the guess
$\tilde{\bm{z}}_{k+1}=\tilde{\bm{z}}_{k}-a_{k}\tilde{\bm{g}}_{k}$.

12: Maximize the cumulative logarithmic likelihood function using
$|\phi\rangle=\sum_{i}\tilde{z}_{k+1,i}|i\rangle/|\tilde{\bm{z}}_{k+1}|$ as
starting point,  $|\tilde{\psi}_{k+1}\rangle=\arg\max_{|\phi\rangle}\log
P(D_{k};|\psi\rangle,S),\quad\text{s. t.}\quad\langle\phi|\phi\rangle=1,$ and
update the estimate as $\tilde{z}_{k+1,i}=\langle
i|\tilde{\psi}_{k+1}\rangle$.

13:end for

Figure 1: Left and right columns show mean and median of
$SE(\bm{z}_{i},\bm{\tilde{z}}_{i,j})$, correspondingly, with respect to
$\bm{\tilde{z}}_{i,j}$ as a function of the number of iterations for four
randomly chosen pure quantum states $\bm{z}_{i}$ in $d=2$ and for ensemble
size $N=10^{3}$ (light blue down threes), $10^{4}$ (solid red triangles), and
$10^{5}$ (solid yellow circles) per iteration. Shaded areas represent the
corresponding interquartile range. Figure 2: Mean (upper row) and median
(lower row) of $MSE(\bm{z}_{i})$ on $\Omega_{d}$ as a function of the
iteration number $k$, for dimension $d$=2, 4, 8, and 16 (from left to right)
and ensemble size $N=10^{3}$ (light blue down threes), $10^{4}$ (solid red
triangles), and $10^{5}$ (solid yellow circles) per iteration. Straight lines
depict the Gill-Massar lower bound $MSE_{GM}(2d,N_{T})$ for the respective
total ensemble size $N_{T}=2Nk$. Shaded areas represent interquartile range.

In order to study the properties of the estimation procedure we create the set
$\Omega_{d}=\\{\bm{z}_{i}\\}$ (with $i=1,\dots,m$) containing $m$ unknown pure
quantum states in dimension $d$. The states in $\Omega_{d}$ are independently
generated according to a uniform distribution. Each state in $\Omega_{d}$ is
estimated by minimizing SE by means of CSPSA concatenated to MLE. This creates
the set $\tilde{\Omega}_{i}=\\{\bm{\tilde{z}}_{i,j}\\}$ (with $j=1,\dots,n$)
for each $\bm{z}_{i}$, which is formed by $n$ estimates $\bm{\tilde{z}}_{i,j}$
of $\bm{z}_{i}$. The estimation accuracy of $\bm{z}_{i}$ is given by the
expectation value of $SE(\bm{z}_{i},\bm{\tilde{z}})$ over the set of all
estimates $\bm{\tilde{z}}$ of $\bm{z}_{i}$, which is approximated by the
expression

$MSE(\bm{z}_{i})=\frac{1}{n}\sum_{j=1}^{n}SE(\bm{z}_{i},\bm{\tilde{z}}_{i,j}).$
(25)

The left column of Fig. 1 shows $MSE(\bm{z}_{i})$, calculated according to Eq.
(25), as a function of the number of iterations $k$ with ensemble size
$N=10^{3},10^{4}$, and $10^{5}$, for four randomly chosen unknown states
$\bm{z}_{i}$ in $d=2$. The four states display very similar behaviors: a rapid
accuracy gain followed by an asymptotic linear regime. The latter emerges
approximately after the iteration $k=6$, where $MSE(\bm{z}_{i})$ reaches
values close to $5\times 10^{-6}$, $5\times 10^{-5}$ (from top to bottom on
each plot), and $5\times 10^{-4}$ for $N=10^{5},10^{4}$, and $10^{3}$,
respectively. The right column of Fig. 1 shows the median of
$SE(\bm{z}_{i},\bm{\tilde{z}}_{i,j})$, another central tendency indicator, in
the set $\tilde{\Omega}_{i}$ and the interquartile range for each one of the
four states. The median also exhibits a sharp accuracy gain followed by a
linear regime. However, the linear regime emerges approximately at iteration
$k=3$. Before this, the median exhibits a large interquartile range, which
indicates a large variation in the values of SE in the set
$\tilde{\Omega}_{i}$. Once the median of SE enters into the linear regime, the
interquartile range becomes extremely narrow. Thus, in the linear regime the
estimation procedure leads to a very sharp distribution of values of SE.
Furthermore, in each one the four inspected cases mean and median of the SE in
the linear regime are almost indistinguishable. In addition, in the linear
regime, the mean and median of SE are independent of the particular unknown
state.

Our estimation method is based on the optimization of the SE with the help of
CSPSA concatenated to MLE. Alternatively, we can reinterpret the method
considering SE optimization through CSPSA to accelerate the MLE convergence
rate. That is, MLE is calculated with the data obtained through SE
measurements. These are chosen with the help of CSPSA in such a way that they
increase the convergence of MLE towards the estimate. In this scenario, the
existence of a linear regime is consistent with the efficiency and asymptotic
normality of MLE estimators [51].

Figure 2 displays the expectation of $MSE(\bm{z})$ over the Hilbert space of
unknown states, that is,

$\overline{MSE}=\mathbb{E}[MSE(\bm{z})|\bm{z}],$ (26)

or equivalently,

$\overline{MSE}=\int p(\bm{z})MSE(\bm{z})\bm{dz},$ (27)

where $p(\bm{z})$ is the Haar-uniform probability density function for the
states in the Hilbert space. This quantity is approximated as the average of
$MSE(\bm{z}_{i})$ over $\Omega_{d}$, which is given by the expression

$\overline{MSE}=\frac{1}{m}\sum_{i=1}^{m}MSE(\bm{z}_{i}).$ (28)

Insets 2(a), 2(b), 2(c), and 2(d) show the behavior of $\overline{MSE}$ as a
function of the number $k$ of iterations for dimension $d=$2, 4, 8, and 16,
respectively, and for several ensemble sizes. As is apparent from these
figures, the expectation of $MSE(\bm{z})$ over the Hilbert space of unknown
states exhibits a rapid estimation accuracy gain followed by a linear regime.
The latter arises after a number of iterations that depends on the particular
dimension. In particular, the higher the dimension the more iterations are
needed for the emergence of the linear regimen. The insets also depict the
lower bound $MSE_{GM}(2d,N_{T})$ as a function of the iteration number $k$ and
the ensemble size $N$, for various dimensions. As the insets show, the
estimation accuracy characteristic of our method becomes very close to
$MSE_{GM}(2d,N_{T})$ as $N$ increases. In fact, the accuracy of the estimates
seems to be asymptotically close to $MSE_{GM}(2d,N_{T})$. Insets 2(e), 2(f),
2(g), and 2(h) illustrate the median of $MSE(\bm{z_{i}})$ over $\Omega_{d}$.
This exhibits a behavior similar to that of $\overline{MSE}$, but the linear
regime emerges earlier. Once the optimization method enters into the linear
regime, the mean and the median of $MSE(\bm{z}_{i})$ reach values that cannot
be distinguished. Furthermore, the interquartile range becomes extremely
narrow. This indicates that in the linear regime the minimization of SE via
the concatenation of CSPSA and MLE leads to an estimation procedure
characterized by a state-independent MSE. Insets 2(e), 2(f), 2(g), and 2(h)
also show twice the Gill-Massar lower bound for the MSE. As is apparent in all
insets, in the linear regime the median of the estimation accuracy is also
very close to $MSE_{GM}(2d,N_{T})$.

To study the linear regime we fitted the numerical data obtained from the
Monte Carlo simulations to the function $p/N_{T}^{a}$, as suggested by the
observation that behavior of $\overline{MSE}$ is close to
$MSE_{GM}(2d,N_{T})=(2d-1)/N_{T}$, where $N_{T}=2kN$, $p=2d-1$ and $a=1$. The
best fits for the values of $p$ and $a$ are shown in Table 1, where two sets
of $p$ and $a$ values are indicated for each dimension $d$ and ensemble size
$N$. The first set of values is obtained fitting data from iteration $k=10$
until $k=45$. The second set of values is obtained fitting data from iteration
$k=46$ until $k=100$. With the exception of the first dataset for $d=16$,
Table 1 indicates that the value of $a$ are in the range $[0.99,1.1]$ with an
average value $\bar{a}$ to 1.01. The anomalous behavior of the first dataset
for $d=16$ can be attributed to the fact that for $d=16$ more than 10
iterations are required for the onset of the linear asymptotic regime. The
first dataset exhibits values of $p$ that are larger than the values of $p$ of
the second dataset. In both datasets, the values of $p$ decrease with an
increase in the ensemble size $N$. This indicates a decrease in the rate at
which the algorithm approaches the minimizer. Finally, the values of $p$ for
$N=10^{4}$ and $10^{5}$ in the second dataset are close to the value of $2d$,
specially the latter. Thus, for larger values of $N$ and $k$, we can
approximate $\overline{MSE}$ as

$\overline{MSE}=\frac{2d+\alpha}{2kN},$ (29)

where $\alpha$ is a small quantity in comparison to $2d$.

| $N=10^{3}$ |  | $N=10^{4}$ |  | $N=10^{5}$ |
---|---|---|---|---|---|---
d=2 | p=6.34 | a=1.03 | p=9.97 | a=1.03 | p=6.52 | a=1.02
| p=7.41 | a=1.04 | p=4.21 | a=0.99 | p=3.57 | a=0.99
d=4 | p=10.56 | a=1.00 | p=9.12 | a=0.99 | p=8.77 | a=0.99
| p=10.18 | a=1.00 | p=8.12 | a=0.99 | p=8.50 | a=0.99
d=8 | p=56.69 | a=1.09 | p=25.76 | a=1.02 | p=19.66 | a=0.99
| p=21.30 | a=1.00 | p=20.16 | a=1.00 | p=17.41 | a=0.99
d=16 | p=4.80 | a=2.21 | p=1.18 | a=2.09 | p=8.84 | a=2.93
| p=147.62 | a=1.10 | p=39.13 | a=0.99 | p=33.87 | a=0.99
Table 1: Fit $\overline{MSE}=p/(2kN)^{a}$ in asymptotic regime.
$MSE_{GM}(2d,2kN)$ is obtained with $p=2d-1$ and $a=1$. For $d$ and $N$ fixed,
the first pair $(p,a)$ is obtained fitting iterations from $k=10$ until
$k=45$. The second pair is obtained fitting iterations from $k=46$ until
$k=100$. Figure 3: Right and left columns show the mean and median of
$SE(\bm{\tilde{z}})$, respectively, as a function of the number of iterations
for four randomly chosen unitary transformations in $d=2$ with ensemble size
$N=10^{3}$ (light blue down threes), $10^{4}$ (solid red triangles), and
$10^{5}$ (solid yellow circles) per iteration. Shaded areas represent the
corresponding interquartile range.

##  IV Estimation of unitary transformations via the optimization of the
squared error

The estimation of processes acting on quantum states is a much more demanding
problem than the estimation of quantum states. For instance, the estimation of
an unknown process acting onto a single qudit requires the characterization of
$d^{4}-d^{2}$ real parameters [52]. In the case of a unitary transformation,
only $d^{2}$ parameters must be determined.

In general, the estimation of a quantum process is carried out by carefully
choosing a set of states, letting the process act on them, and reconstructing
the output states by means of a quantum tomographic method [53]. We will
employ this strategy to estimate an unknown unitary transformation $U$. This
is suggested by Eq. (21), which can be cast in the form

$P_{2}=\frac{1}{4}\sum_{j=1}^{d}|U_{j,k}-\tilde{U}_{j,k}|^{2}.$ (30)

Thus, the probability $P_{2}$ is proportional to the squared error between the
k-th columns of the matrices $U$ and $\tilde{U}$, where $k$ is controlled by
the initial path state $|k\rangle$ followed by the single photon. Clearly, we
can reconstruct each column of $U$ by minimizing the squared error by CSPSA
and MLE. After estimating all columns of $U$ independently, we obtain an
estimate $\tilde{U}$ of $U$.

Figure 4: Mean (left column) and median (right column) of $MSE$ for randomly
generated unitary transformations as a function of the iteration number $k$
for dimension $d$=2 and ensemble size $N=10^{3}$ (light blue down threes),
$10^{4}$ (solid red triangles), and $10^{5}$ (solid yellow circles) per
iteration. Shaded areas represent interquartile range. From the top row to the
bottom row: estimates $\tilde{U}$ provided by CSPSA, estimates $\tilde{U}_{c}$
provided by CSPSA projected to the closest unitary transformation at each
iteration, and estimates $\tilde{U}_{gs}$ provided by CSPSA post-processed
with the Gramm-Schmidt orthogonalization procedure at each iteration.

However, the present estimation method cannot guarantee that the estimate
$\tilde{U}$ is really unitary. We consider two methods to obtain a unitary
estimate. We can obtain a unitary estimate from $\tilde{U}$ by means of the
expression [54]

$\tilde{U}_{c}=\tilde{U}(\tilde{U}\tilde{U}^{\dagger})^{-1/2},$ (31)

which is the closest unitary operator to the transformation $\tilde{U}$. In
order to quantify how close $\tilde{U}_{c}$ is from $U$ we employ the Hilbert-
Schmidt distance
$D(U,\tilde{U}_{c})=Tr[(U-\tilde{U}_{c})(U-\tilde{U}_{c})^{\dagger}]$. Another
method to generate a unitary estimate $\tilde{U}_{gs}$ consists in the
application of the Gramm-Schmidt orthogonalization procedure to the columns
$\tilde{U}_{j,k}$.

Figure 3 displays the mean and median of $SE(\bm{\tilde{z}})$ for the estimate
$\tilde{U}$ for four randomly chosen unitary transformations as a function of
the number of iterations $k$ for different ensemble sizes and $d=2$. The
overall behavior of the mean-squared error for the unitary transformations
resembles very closely the case of the estimation of pure states, that is, a
fast decrease of the MSE followed by a linear regime. In the case of the mean
MSE the linear regime arises at a higher number of iterations when compared to
the case of estimating pure states. In the case of the median, the linear
regime emerges at a similar number of iterations as in the case of the median
MSE for pure states. After 10 iterations, Fig. 3 shows that the mean and
median MSE achieve similar values, which are half order of magnitude higher
than the case of estimating pure states. This entail a loss of accuracy when
estimating unitary transformations with a method designed to estimate pure
states. This is, however, not dramatic. After 10 iterations the MSE reaches
values of the order of $0.5\times 10^{-3}$, $0.5\times 10^{-4}$ and $0.3\times
10^{-5}$ for ensamble sizes of $10^{3}$, $10^{4}$ and $10^{5}$,
correspondingly.

Figure 5: Mean (left column) and median (right column) of $MSE$ for randomly
generated unitary transformations as a function of the iteration number $k$
for dimension $d$=4 and ensemble size $N=10^{3}$ (light blue down threes),
$10^{4}$ (solid red triangles), and $10^{5}$ (solid yellow circles) per
iteration. Shaded areas represent interquartile range. From the top row to the
bottom row: estimates $\tilde{U}$ provided by CSPSA, estimates $\tilde{U}_{c}$
provided by CSPSA projected to the closest unitary transformation at each
iteration, and estimates $\tilde{U}_{gs}$ provided by CSPSA post-processed
with the Gramm-Schmidt orthogonalization procedure at each iteration. Figure
6: Mean (left column) and median (right column) of $MSE$ for randomly
generated unitary transformations as a function of the iteration number $k$
for dimension $d$=2 and ensemble size $N=10^{3}$ (light blue down threes),
$10^{4}$ (solid red triangles), and $10^{5}$ (solid yellow circles) per
iteration. Shaded areas represent interquartile range. From the top row to the
bottom row: estimates provided by CSPSA, estimates provided by CSPSA updated
at each iteration by projection to the closest unitary transformation, and
estimates provided by CSPSA updated at each iteration with the Gramm-Schmidt
orthogonalization procedure.

Figure 4 shows the mean and median mean-squared error achieved in process of
estimating unknown unitary transformations acting onto a 2-dimensional quantum
system. The left (right) column exhibits the mean (median) achieved with the
estimates $\tilde{U}$, $\tilde{U}_{c}$, and $\tilde{U}_{gs}$ from top to
bottom, respectively. The typical behavior of a rapidly increasing estimation
accuracy followed by a linear regimen is clearly present. This exhibits in the
linear regime a mean and a median that cannot be distinguished from each other
and an extremely narrow interquartile range, which indicates that after 10
iterations all unitary transformations are estimated with the same accuracy.
This is almost twice the accuracy obtained in estimating a 2-dimensional pure
state, as expected. The three estimates lead to very similar accuracies, but
the estimate $\tilde{U}_{c}$ generates a marginally better performance. Figure
5 exhibits similar results in the case $d=4$.

In the previous simulations we have considered that after each iteration the
estimate of $U$ is post-processed with the help of the Gram-Schmidt
orthogonalization procedure or the projection onto the closest unitary
transformation. This information has not been employed to modify the estimate
along the sequence of iterations. However, we can estimate the columns of the
unitary transformation, obtain a non-unitary estimate of the unitary, and
generate from it a unitary estimate. Thereafter, the columns of this unitary
estimate are employed as guesses for the next iteration round. Figure 6 shows
the effect on the estimation quality of performing such an update in the case
of $d=2$. The left (right) column exhibits the mean (median) achieved with the
estimate $\tilde{U}$ (for comparison purpose), the estimate $\tilde{U}_{c}$
updated by means of the projection onto the set of unitary transformations,
and the estimate $\tilde{U}_{gs}$ updated by means of the Gram-Schmidt
orthogonalization procedure from top to bottom, respectively. As it is
apparent from this figure, the update of the estimates leads to a modest
increase in estimation accuracy. However, in dimension $d=4$ there is a
significative improvement in the estimation accuracy. This is illustrated in
Fig. 7, where an improvement in half order of magnitude is achieved in
comparison to Fig. 5, where no update of the estimates is performed.
Nevertheless, the convergence rate to the linear regime is reduced by a few
iterations.

In Algorithm 2 we present a basic pseudocode for implementing the MSE-based
estimation of unitary transformations. An implementation of the pseudocode in
the Python programming language can be found in the GitHub repository [50]. We
consider the different choices for the post-processing of the estimates.
According to option 1 we project the possibly non-unitary estimate onto the
set of unitary transformations and evaluate the infidelity. These projections
are used in option 2 to provide a better update of the estimate.

Figure 7: Mean (left column) and median (right column) of $MSE$ for randomly
generated unitary transformations as a function of the iteration number $k$
for dimension $d$=4 and ensemble size $N=10^{3}$ (light blue down threes),
$10^{4}$ (solid red triangles), and $10^{5}$ (solid yellow circles) per
iteration. Shaded areas represent interquartile range. Estimates are updated
after each iteration. From the top row to the bottom row: estimates provided
by CSPSA, estimates provided by CSPSA updated at each iteration by projection
to the closest unitary transformation, and estimates provided by CSPSA updated
at each iteration with the Gramm-Schmidt orthogonalization procedure.

Since each of the $d$ columns of $U$ is estimated with an accuracy close to
$(2d+\alpha)/N_{T}$, where $N_{T}$ is the total ensamble size used in the
estimation of each column, we have that the unitary transformations are
estimated with an accuracy $MSE(U)$ given approximately by

$MSE(U)\approx\frac{d(2d+\alpha)}{N_{T}}.$ (32)

Thereby, the estimation accuracy of our procedure becomes

$MSE(U)\approx\frac{d^{2}(2d+\alpha)}{N^{*}_{T}},$ (33)

where $N^{*}_{T}=dN_{T}$ is the total number of copies used in the estimation
of all $d$ columns of $U$.

It is possible to follow a different approach to the estimation of a unitary
transformation. The multi-arm interferometer can be feed with a maximally
mixed state $I/d$. In this case the detection of photons at the output ports
leads to the SE of all the elements of the unknown unitary transformation.
This, however, requires the solution of the MLE problem considering all
elements of a density matrix. This procedure increases the computational cost
of the estimation process and entails a reduction of the estimation accuracy,
which becomes more severe in higher dimensions.

Algorithm 2 MSE-based estimation of unitary transformations

1:Consider a known unitary transformation $U$ on the upper arm of the
interferometer.

2:Choose initial estimate $\tilde{U}_{0}$ and define
$\tilde{z}_{0,i}^{j}=\tilde{U}_{0,ij}$.

3:Set gain coefficients $a$, $A$, $s$, $b$ and $r$.

4:for $k=1,\ldots,k_{max}$ do

5: Set  $a_{k}=\frac{a}{(k+1+A)^{s}},\quad c_{k}=\frac{b}{(k+1)^{r}}.$

6: for $j=1,\dots,d$ do

7: Choose $\Delta_{k,i}^{j}$ randomly in the set $\\{\pm 1,\pm i\\}$.

8: Calculate
$|\psi_{k\pm}^{j}\rangle=\sum_{i}\tilde{z}_{k\pm,i}^{j}|i\rangle/|\tilde{\bm{z}}_{k\pm}^{j}|$,
with $\tilde{\bm{z}}_{k\pm}^{j}=\tilde{\bm{z}}_{k}\pm c_{k}{\bm{\Delta}_{k}}$.

9: Feed the interferometer in the mode $|j\rangle$.

10: Prepare the states $|\psi_{k\pm}^{j}\rangle$ on the lower arm of the
interferometer with the transformation $U^{\prime}$.

11: Estimate experimentally the square errors
$SE({\bm{z}^{j}},\tilde{{\bm{z}}}^{j}_{k\pm})$ of the $j$-column of $U$ with a
sample of size $N$.

12: Estimate the gradient as
$\tilde{g}_{k,i}^{j}=\frac{SE({\bm{z}}^{j},\tilde{\bm{z}}^{j}_{k+})-SE({\bm{z}}^{j},\tilde{\bm{z}}^{j}_{k-})}{2c_{k}\Delta_{k,i}^{*}}.$

13: Actualize the guess
$\tilde{\bm{z}}^{j}_{k+1}=\tilde{\bm{z}}^{j}_{k}-a_{k}\tilde{\bm{g}}^{j}_{k}$.

14: Maximize the cumulative Likelihood function using
$|\phi^{j}\rangle=\sum_{i}\tilde{z}^{j}_{k+1,i}|i\rangle/|\tilde{\bm{z}}^{j}_{k+1}|$
as starting point,
$|\tilde{\psi}_{k+1}^{j}\rangle=\arg\max_{|\phi\rangle}\log
P(D_{k}^{j};|U|j\rangle,S),\quad\text{s. t.}\quad\langle\phi|\phi\rangle=1,$
and update the estimate as $\tilde{z}_{k+1,i}^{j}=\langle
i|\tilde{\psi}_{k+1}^{j}\rangle$.

15: end for

16: The estimated unitary matrix is
$\tilde{U}_{k+1,ij}=\tilde{z}_{k+1,i}^{j}$.

17: Option 1: In order to guarantee the unitarity of the estimator, we
consider two postprocessing methods:

  * •

$\tilde{U}_{c}$: Project $\tilde{U}$ into its closest unitary matrix.

  * •

$\tilde{U}_{gs}$: Apply the Gram-Schmidt procedure to the columns of
$\tilde{U}_{k+1}$.

18: Option 2: Re-update the estimates $\\{\tilde{\bm{z}}_{k+1}^{j}\\}$ with
the postprocessed unitary matrix,
$\tilde{z}_{k+1,i}^{j}=(\tilde{U}_{c})_{k+1,ij}\qquad\text{or}\qquad\tilde{z}_{k+1,i}^{j}=(\tilde{U}_{gs})_{k+1,ij}.$

19:end for

##  V Conclusions

We have studied the estimation of pure quantum states employing the mean-
squared error as accuracy metric. We have presented two setups within reach of
actual experimental techniques, one for a polarization qubit and other for a
path-encoded qudit, that allow one to measure the squared error. The mean-
squared error arises as a sampling of the squared error for a fixed unknown
state. The estimates of an unknown state are obtained by a combination of
stochastic optimization on the field of the complex numbers and maximum
likelihood estimation. The estimation of an unknown state is formulated as the
minimization of the infidelity with respecto to a known state. This
optimization problem is solved by means of CSPSA, which iteratively drives a
sequence of measurements in such a way that the infidelity approaches zero.
The rate of convergence of CPSA is increased by refining the estimates with
the help of maximum likelihood estimation applied onto the total sequence of
measurement results.

Monte Carlo numerical experiments show that the accuracy achieved in the
estimation of a fixed unknown state by the combination of CSPSA and MLE
exhibits, as a function of the number of iterations, two clearly defined
regimes: a fast decrease followed by a lineal behavior. The mean and median
squared error exhibit very close values, which indicates the absence of
outliers. The median squared error enters the lineal regime a few iterations
before than the mean-squared error. In addition, numerical experiments in a
wide range of dimension and ensemble size indicate that the achieved
estimation accuracy is nearly state independent. Furthermore, the estimation
accuracy reaches very close values to twice the Gill-Massar lower bound for
the mean-squared error.

We have also extended the estimation of pure states to the estimation of
unitary transformations. In first place we hace considered the estimation by
simply estimating each column of a unitary transformation, which leads to an
accuracy similar to obtained in the estimation of pure states. However, in
this case it is not possible to guaranty that the estimates fulfill the
condition of unitarity. In order to accomplish this, the estimates are
projected onto the closest unitary transformation or the estimated columns
undergo the Gram-Schmidt orthogonalization procedure. In both cases there is
an improvement in the estimation accuracy. A much more significative
improvement is obtained when incorporating the projection onto the closest
unitary transformation or the Gram-Schmidt orthogonalization procedure to the
iteration process that delivers the estimates. The present method also
inherits several of the characteristics exhibited in the estimation of states,
that is, it is independent of the unitary to be estimated, exhibits very close
values of mean and median accuracy, and a very narrow interquartile range.

###### Acknowledgements.

This work was supported by ANID – Millennium Science Initiative Program –
ICN17-012\. AD was supported by FONDECYT Grant 1180558. LP was supported by
ANID-PFCHA/DOCTORADO-BECAS-CHILE/2019-72200275.

## References

  * [1] D. H. Mahler, L. A. Rozema, A. Darabi, C. Ferrie, R. Blume-Kohout, and A. M. Steinberg, Phys. Rev. Lett. 111, 183601 (2013).
  * [2] S. Straupe, JETP Lett. 104, 510 (2016).
  * [3] Z. Hou, H. Zhu, G. Xiang, C.-F. Li, and G.-C. Guo, npj Quantum Inf. 2, 16001 (2016).
  * [4] L. Pereira, L. Zambrano, J. Cortés-Vega, S. Niklitschek, and A. Delgado, Phys. Rev. A 98, 012339 (2018).
  * [5] G. I. Struchalin, E. V. Kovlakov, S. S. Straupe, and S. P. Kulik, Phys. Rev. A 98, 032330 (2018).
  * [6] D. Goyeneche, G. Cañas, S. Etcheverry, E. S. Gomez, G. B. Xavier, G. Lima, and A. Delgado, Phys. Rev. Lett. 115, 090401 (2015).
  * [7] L. Zambrano, L. Pereira, and A. Delgado, Phys. Rev. A 100, 022340 (2019).
  * [8] R. D. Gill and S. Massar, Phys. Rev. A 61, 042312 (2000).
  * [9] R. Jozsa, J. Mod. Opt. 41, 2315 (1994).
  * [10] A. Uhlmann, Rep. Math. Phys. 9, 273 (1976).
  * [11] R. D. Gill and S. Massar, Phys. Rev. A 61, 042312 (2000).
  * [12] D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, Phys. Rev. A 64, 052312 (2001).
  * [13] R. T. Thew, K. Nemoto, A. G. White, and W. J. Munro, Phys. Rev. A 66, 012303 (2002).
  * [14] Z. Hradil, Phys. Rev. A 55, R1561 (1997).
  * [15] S. M. Tan, J. Mod. Opt. 44, 2233 (1997).
  * [16] C. Ferrie, Phys. Rev. Lett. 113, 190404 (2014).
  * [17] A. Utreras-Alarcón, M. Rivera-Tapia, S. Niklitschek, and A. Delgado, Sci. Rep. 9, 16143 (2019).
  * [18] L. Zambrano, L. Pereira, S. Niklitschek, and A. Delgado, Sci. Rep. 10, 12781 (2020).
  * [19] C. Benedetti, A. P. Shurupov, M. G. A. Paris, G. Brida, and M. Genovese, Phys. Rev. A 87, 052136 (2013).
  * [20] M. Bina, A. Mandarino, S. Olivares, and M. G. A. Paris, Phys. Rev. A 89, 012305 (2014).
  * [21] A. Mandarino, M. Bina, C. Porto, S. Cialdi, S. Olivares, and M. G. A. Paris, Phys. Rev. A 93, 062118 (2016).
  * [22] G. Casella and E. L. Lehmann, Theory of Point Estimation (Springer Texts in Statistics) (Springer-Verlag, New York, 1998).
  * [23] J. L. O’Brien, G. J. Pryde, A. Gilchrist, D. F. V. James, N. K. Langford, T. C. Ralph, and A. G. White, Phys. Rev. Lett. 93, 080502 (2004).
  * [24] J.-Y. Wang, B. Yang, S.-K. Liao, L. Zhang, Q. Shen, X.-F. Hu, J.-C. Wu, S.-J. Yang, H. Jiang, Y.-L. Tang, B. Zhong, H. Liang, W.-Y. Liu, Y.-H. Hu, Y.-M. Huang, B. Qi, J.-G. Ren, G.-S. Pan, J. Yin, J.-J. Jia, Y.-A. Chen, K. Chen, C.-Z. Peng, and J.-W. Pan, Nature Photon. 7, 387 (2013).
  * [25] J. Yuen-Zhou, J. J. Krich, M. Mohseni, and A. Aspuru-Guzik, Proc. Natl. Acad. Sci. USA 108, 17615 (2011).
  * [26] M. Riebe, K. Kim, P. Schindler, T. Monz, P. O. Schmidt, T. K. Körber, W. Hänsel, H. Häffner, C. F. Roos, and R. Blatt, Phys. Rev. Lett. 97, 220407 (2006).
  * [27] A. M. Childs, I. L. Chuang, and D. W. Leung, Phys. Rev. A 64, 012314 (2001).
  * [28] R. C. Bialczak, M. Ansmann, M. Hofheinz, E. Lucero, M. Neeley, A. D. O’Connell, D. Sank, H. Wang, J. Wenner, M. Steffen, A. N. Cleland, and J. M. Martinis, Nature Phys. 6, 409 (2010).
  * [29] M. Howard, J. Twamley, C. Wittmann, T. Gaebel, F. Jelezko, and J. Wrachtrup, New J. Phys. 8, 33 (2006).
  * [30] Y. Nambu, K. Usami, A. Tomita, S. Ishizaka, T. Hiroshima, Y. Tsuda, K. Matsumoto, and K. Nakamura, Proc. SPIE 4917, 13 (2002).
  * [31] F. De Martini, A. Mazzei, M. Ricci, and G. M. D’Ariano, Phys. Rev. A 67, 062307 (2003).
  * [32] J. B. Altepeter, D. Branning, E. Jeffrey, T. C. Wei, P. G. Kwiat, R. T. Thew, J. L. O’Brien, M. A. Nielsen, and A. G. White, Phys. Rev. Lett. 90, 193601 (2003).
  * [33] C. H. Baldwin, A. Kalev, and I. H. Deutsch, Phys. Rev. A 90, 012110 (2014).
  * [34] A. Acín, E. Jané, and G. Vidal, Phys. Rev. A 64, 050302(R) (2001).
  * [35] A. Acín, Phys. Rev. Lett. 87, 177901 (2001).
  * [36] M. A. Taylor, J. Janousek, V. Daria, J. Knittel, B. Hage, H.-A. Bachor, and W. P. Bowen, Nature Photon. 7, 229 (2013).
  * [37] F. Wolfgramm, C. Vitelli, F. A. Beduini, N. Godbout, and M. W. Mitchell, Nature Photon. 7, 28 (2012).
  * [38] Xiao-Qi Zhou, Hugo Cable, Rebeca Whittaker, Peter Shadbolt, Jeremy L. O’Brien, and Jonathan C. F. Matthews, Optica 2, 510 (2015).
  * [39] W. Wirtinger, Math. Ann. 97, 357 (1927).
  * [40] D. H. Brandwood, Proc. IEE-H 130, 11 (1983).
  * [41] Z. Nehari, Introduction to Complex Analysis (Allyn & Bacon, Boston, 1961).
  * [42] R. Remmert, Theory of Complex Functions (Springer-Verlag, New York, 1991).
  * [43] M. G. A. Paris, Int. J. Quantum. Inform. 7, 125 (2009).
  * [44] N. Li, C. Ferrie, J. A. Gross, A. Kalev, and C. M. Caves, Phys. Rev. Lett. 116, 180402 (2016).
  * [45] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, Nature Photon. 14, 273 (2020).
  * [46] G. B. Xavier and G. Lima, Commun. Phys. 3, 9 (2020).
  * [47] J. Cariñe, G. Cañas, P. Skrzypczyk, I. Šupić, N. Guerrero, T. Garcia, L. Pereira, M. A. S. Prosser, G. B. Xavier, A. Delgado, S. P. Walborn, D. Cavalcanti, and G. Lima, Optica 7, 542 (2020).
  * [48] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Phys. Rev. Lett. 73, 58 (1994).
  * [49] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optica 3, 1460 (2016).
  * [50] A Python programming language version of the pseudocode in Algorithm 1 can be found in the GitHub repository.
  * [51] L. E. Lehmann, Elements of Large-Sample Theory (Springer-Verlag, New York, 1999).
  * [52] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cambridge, UK, 2000).
  * [53] C. H. Baldwin, A. Kalev, and I. H. Deutsch, Phys. Rev. A 90, 012110 (2014).
  * [54] J. Keller, Mathematics Magazine 48, 192 (1975).