# Efficient Depth Selection for the Implementation of Noisy Quantum
Approximate Optimization Algorithm

Yu Pan1,2<EMAIL_ADDRESS>Yifan Tong2 Shibei Xue3,4 Guofeng Zhang5 1State Key
Laboratory of Industrial Control Technology, Zhejiang University, Hangzhou,
310027, P. R. China  2Institute of Cyber-Systems and Control, College of
Control Science and Engineering, Zhejiang University, Hangzhou, 310027, P. R.
China  3Department of Automation, Shanghai Jiao Tong University, Shanghai
200240, P. R. China  4Key Laboratory of System Control and Information
Processing, Ministry of Education of China, Shanghai 200240, P. R. China
5Department of Applied Mathematics, The Hong Kong Polytechnic University, Hung
Hom, Kowloon, Hong Kong Special Administrative Region of China, P. R. China

###### Abstract

Noise on near-term quantum devices will inevitably limit the performance of
Quantum Approximate Optimization Algorithm (QAOA). One significant consequence
is that the performance of QAOA may fail to monotonically improve with depth.
In particular, optimal depth can be found at a certain point where the noise
effects just outweigh the benefits brought by increasing the depth. In this
work, we propose to use the model selection algorithm to identify the optimal
depth with a few iterations of regularization parameters. Numerical
experiments show that the algorithm can efficiently locate the optimal depth
under relaxation and dephasing noises.

##  I Introduction

In the absence of effective and scalable error correction techniques, the
current research on near-term Noisy Intermediate Scale Quantum (NISQ) Preskill
(2018) computing has been focused on Quantum Approximate Optimization
Algorithm (QAOA) which is expected to achieve quantum advantage for practical
application Farhi et al. (2014, 2015); Farhi and Harrow (2019). QAOA is a
hybrid quantum-classical algorithm for solving a special class of problems
such as combinatorial optimization. An approximate solution for such problems
would be acceptable as they have been proven to be NP-hard Zhou et al. (2020);
Borle et al. (2021).

QAOA consists of a quantum control protocol whose variational control
parameters are updated according to the measurement result of each iteration.
During each iteration, two noncommuting Hamiltonian controls are alternatively
applied on the quantum system, for which the control parameters can be
understood as angles of rotations or control durations. Each Hamiltonian
control is then implemented as a circuit made up of quantum gates. Due to the
high noise level and gate imperfection in the NISQ devices, the performance of
QAOA critically relies on the number of variational control parameters, which
is defined as the control depth in this paper. Nevertheless, the current
literature often assumes a noise-free model Farhi et al. (2014); Jiang et al.
(2017); Hastings (2018); Wang et al. (2018); Zhou et al. (2020); Larkin et al.
(2020); Bravyi et al. (2020); Phillip C. Lotshaw and Siopsis (2021). In
particular, the performance of QAOA will monotonically increase with the
control depth Farhi et al. (2014) if noise is not taken into consideration.
However, a large depth is not realistic for NISQ devices due to the prominent
noise effect and limited hardware capacity. Therefore, a very small depth,
which is usually determined by empirical rules, has to be chosen in these
works to validate the noise-free assumption.

QAOA will inevitably confront the need to increase the control depth as the
problem to solve becomes more and more complex. In that case, noise effect has
to be accounted for and mitigated to improve the overall performance. Although
circuit design, gate optimization and increasing the freedom of classical
control Alam et al. (2020, 2020); Lacroix et al. (2020); Herrman et al. (2021)
have been proposed to reduce the control depth and counteract the noise
effect, it is still not clear that how to select an optimal control depth
under the influence of a specific noise. Moreover, since the algorithm is
actually implemented on quantum hardware, the noise is often unknown and the
control parameters may have to be optimized via a data-driven way. It should
be noted that one does not need to know the underlying model of the quantum
system for the execution of QAOA. Hence, QAOA can be taken as a type of
iterative learning control Dong et al. (2020) which is naturally adaptive to
noise. More explicitly, during each iteration, the control parameters are
updated based on the measurement results (e.g., Farhi et al. (2014);
Guerreschi and Smelyanskiy (2017); Verdon et al. (2019); Sweke et al. (2020);
Liang et al. (2020)) to maximize the optimization accuracy no matter what
noise the quantum system has been exposed to.

Although extensive research efforts have been devoted to developing algorithms
for optimizing the control parameters, the automatic optimization of control
depth under quantum noise has not been studied so far. It is commonly known
that the performance of a learning or optimization model mainly lies on the
selection of hyperparameters, while control depth is the most important
hyperparameter for QAOA Bengtsson et al. (2020); Campos et al. (2021); Rebekah
Herrman and Siopsis (2021). In this paper, we propose an automatic algorithm
for optimizing the control depth of noisy QAOA. To be more precise, we adopt
the $l_{1}$ regularization technique to reduce the preselected large control
depth to an optimal value during the iteration. This approach has been widely
used to enforce sparse solutions James et al. (2013) in statistical learning.
The sparsity, or complexity of the $l_{1}$-regularized model, has to be
controlled by a proper criterion. For example, Akaike’s Information Criteria
(AIC) and Bayesian Information Criteria (BIC) Burnham and Anderson (2002) can
be used to find the balance between a good optimization accuracy and
complexity of the control depth. For noisy QAOA, instead of designing an
empirical criterion for model selection, the optimal depth is naturally
constrained. Due to the inevitable noise in NISQ devices, the computation
errors will accumulate over time and optimization accuracy will decrease with
the depth after reaching a peak. Particularly, it has been found by numerical
simulations and experiments in Marshall et al. (2020); Alam et al. (2020);
Harrigan et al. (2021) that the performance of QAOA may not monotonically
increase with depth under noise. As a result, the optimal depth can be simply
defined as the one achieving the highest accuracy. The optimal depth with the
highest accuracy is then found by solving an $l_{1}$-regularized optimization
problem with a fast algorithm, which is much more efficient than exhaustive
search.

This paper is organized as follows. A brief introduction to noisy QAOA and its
$l_{1}$ regularized variant is given in Section II. The algorithm for depth
optimization is given in Section III. In Section IV, the Max-Cut example,
metrics and numerical results are presented. Section V concludes this paper.

##  II Standard and Noisy QAOA

The noise-free QAOA is executed as a sequence of unitary control, which is
generated by applying a piece-wise Hamiltonian as

$H:=H_{o}\stackrel{{\scriptstyle\gamma_{1}}}{{\longrightarrow}}H_{c}\stackrel{{\scriptstyle\beta_{1}}}{{\longrightarrow}}\
\cdot\cdot\cdot\
H_{o}\stackrel{{\scriptstyle\gamma_{p}}}{{\longrightarrow}}H_{c}\stackrel{{\scriptstyle\beta_{p}}}{{\longrightarrow}}.$
(1)

Here $H_{o}$ is the problem-based Hamiltonian whose ground states encode the
solution to the optimization problem, and
$H_{c}=\sum_{n=1}^{N}\sigma_{x}^{(n)}$ is the noncommuting control
Hamiltonian. The variational parameters
$\beta=(\beta_{1},\cdot\cdot\cdot,\beta_{p}),\gamma=(\gamma_{1},\cdot\cdot\cdot,\gamma_{p})$
can be understood as control durations or angles of rotations. $2p$ is the
initial depth of control sequence. We have used the following notations

$\sigma_{z}=\left(\begin{array}[]{cc}1&0\\\ 0&-1\end{array}\right),\
\sigma_{x}=\left(\begin{array}[]{cc}0&1\\\ 1&0\end{array}\right),\
\sigma_{-}=\left(\begin{array}[]{cc}0&0\\\ 1&0\end{array}\right).$ (2)

The computational basis states of a single qubit are written as
$|0\rangle=(0\quad 1)^{T}$ and $|1\rangle=(1\quad 0)^{T}$. During each
iteration, the final state is generated by the Hamiltonian control as

$|\psi(\beta,\gamma)\rangle=U(H_{c},\beta_{p})U(H_{o},\gamma_{p})\cdot\cdot\cdot
U(H_{c},\beta_{1})U(H_{o},\gamma_{1})|s\rangle,\\\ $

where we have defined the unitary evolutions as
$U(H_{c},\beta_{j})=e^{-\mbox{i}H_{c}\beta_{j}}$ and
$U(H_{o},\gamma_{j})=e^{-\mbox{i}H_{o}\gamma_{j}}$. The initial state is
$|s\rangle=|+\rangle_{1}\cdot\cdot\cdot|+\rangle_{N}$, with
$|+\rangle=1/\sqrt{2}(|0\rangle+|1\rangle)$. The standard QAOA attempts to
find the approximate solution to the following problem

$\min_{x=(\beta,\gamma)}f(x)=\min_{x}\langle\psi(x)|H_{o}|\psi(x)\rangle,$ (3)

where $f(x)$ is the expectation of $H_{o}$ taken on the final state
$|\psi(x)\rangle$. Since $H_{o}$ can be decomposed as a combination of quantum
projective measurements, $f(x)$ can be estimated by repeating the quantum
control procedures for many times and taking the average of measurement
results. Clearly, the state preparation, evolution and measurement can be
implemented on quantum hardware, while the control parameters are optimized
using a classical algorithm.

As QAOA is expected to run on noisy quantum hardware, the practical
optimization accuracy may not monotonically increase with the depth under the
influence of quantum noise. Specifically, the state of an open quantum system
is described by a density matrix instead of a pure state vector. In order to
study the influence of the noise, we simulate the dynamical evolution of the
open quantum system using quantum master equation. The general formulation of
quantum master equation in the Lindblad form reads Breuer and Petruccione
(2007)

$\dot{\rho}=-\mbox{i}[H,\rho]+\sum_{n=1}^{N}\gamma_{n}(L_{n}\rho
L_{n}^{\dagger}-\frac{1}{2}L_{n}^{\dagger}L_{n}\rho-\frac{1}{2}\rho
L_{n}^{\dagger}L_{n}),$ (4)

with $\gamma_{n}$ being the coupling strength between the $n$-th qubit and its
environment. We assume that each qubit is independently coupled to the
environment via the coupling operator $L_{n}$. $H$ is the piece-wise control
Hamiltonian defined in (1) that alternates between $H_{o}$ and $H_{c}$. It is
clear from (4) that the decay of a quantum state will accumulate as the total
control time increases.

The initial state is then written as $\rho_{0}=|s\rangle\langle s|$. The state
generated at the end of each iteration is denoted as $\rho(x)$. It should be
noted that in this case $H_{o}$ is measured against $\rho(x)$ instead of a
pure state, and the measurement result is calculated by
$\mbox{tr}(H_{o}\rho(x))$. The regularized model for noisy QAOA is thus given
by

$\displaystyle\min_{x}\mbox{tr}(H_{o}\rho(x))+\lambda||x||_{1},$ (5)
$\displaystyle||x||_{1}=\sum_{j=1}^{p}(|\beta_{j}|+|\gamma_{j}|),$

where $||\cdot||_{1}$ refers to the $l_{1}$ norm, and $\lambda>0$ is a
regularization parameter.

##  III Algorithm

###  III.1 Proximal Gradient Descent

The unregularized version of (5) can be solved simply by calculating its
numerical gradient as Guerreschi and Smelyanskiy (2017); Wang et al. (2018);
Sweke et al. (2020); Liang et al. (2020)

$[\nabla\mbox{tr}(H_{o}\rho(x))]_{i}\approx\frac{\mbox{tr}(H_{o}\rho(x_{i}-\epsilon))-\mbox{tr}(H_{o}\rho(x_{i}+\epsilon))}{2\epsilon}$
(6)

with a small $\epsilon>0$. Here $[\nabla\mbox{tr}(H_{o}\rho(x))]_{i}$ denotes
the $i$-th element of the gradient vector. The expectation
$\mbox{tr}(H_{o}\rho(x_{i}\pm\epsilon))$ is obtained by perturbing the control
parameter $x_{i}$ by a small amount of $\epsilon$ and measuring $H_{o}$ with
respect to the generated states. The standard gradient update is then given by

$x_{k+1}=x_{k}-\eta\nabla\mbox{tr}(H_{o}\rho(x_{k})),$ (7)

where $\eta$ is the step size, and $k$ is the current iteration step. However,
since the regularization term in (5) is not differentiable, the standard
gradient descent is not directly applicable to this problem.

Algorithm 1 Proximal Gradient Descent

1: for $k=1,\cdots,K$ do

2: $\quad z_{k}=x_{k}-\eta\nabla\mbox{tr}(H_{o}\rho(x_{k}));$

3: $\quad x_{k+1}=S_{\lambda\eta}(z_{k});\quad\backslash\backslash$ Soft
thresholding

4: end for

5: return $x_{K+1}$;

Proximal Gradient (PG) descent Beck and Teboulle (2009) provides a solution to
the regularized problem (5) by minimizing the quadratic approximation to
$\mbox{tr}(H_{o}\rho(x_{k}))$ around $x_{k}$ as

$\displaystyle x_{k+1}=$
$\displaystyle\mathop{\arg\min}_{y}\mbox{tr}(H_{o}\rho(x_{k}))+\nabla\mbox{tr}(H_{o}\rho(x_{k}))^{T}(y-x_{k})$
$\displaystyle+\frac{1}{2\eta}||y-x_{k}||_{2}^{2}+||y||_{1},$ (8)

in which the quadratic gradient $\nabla^{2}f(\cdot)$ is approximated by
$I/\eta$ while leaving the regularization term $||y||_{1}$ alone. The solution
to (8) can be obtained using PG descent detailed in Algorithm 1, which
introduces an additional soft-thresholding operation after the standard
gradient descent.

The soft-thresholding operation is defined by

$S_{\lambda\eta}({[z_{k}]}_{i})=\left\\{\begin{array}[]{rcl}&{[z_{k}]}_{i}-\lambda\eta,&{[z_{k}]}_{i}>\lambda\eta\\\
&0,&-\lambda\eta\leq{[z_{k}]}_{i}\leq\lambda\eta\\\
&{[z_{k}]}_{i}+\lambda\eta,&{[z_{k}]}_{i}<-\lambda\eta\end{array}\right.$ (9)

The values of the control parameters are shifted towards zero by an amount of
$\lambda\eta$ after soft thresholding. Moreover, once the condition
$|{[z_{k}]}_{i}|\leq\lambda\eta$ is satisfied, the $i$-th control parameter
will be penalized to exactly zero, which in turn reduces the control depth by
removing the corresponding control action from the control sequence.

###  III.2 Choice of $\lambda$

The choice of regularization parameter $\lambda$ is very important for
determining the depth. Generally speaking, a large $\lambda$ leads to a large
threshold value $\lambda\eta$, and thus the control parameters are more likely
to be penalized to zero during the iteration. On the contrary, a small
$\lambda$ will keep most of the parameters in the model. Therefore, a set of
values have to be tested in order to determine an appropriate regularization
parameter that yields the optimal depth. To do this, for each $\lambda$, an
optimized accuracy is obtained using the regularized QAOA with reduced depth.
The optimized accuracies for different $\lambda$ are compared, and the optimal
depth is the one that achieves the best accuracy.

The common practice to solve the $l_{1}$-regularized optimization problem is
to start with a relatively large $\lambda$, which will shrink by a constant
factor after each round of experiment. As $\lambda$ decreases, the selected
depth in general will increase to improve the optimization accuracy. However,
due to the noise effect, the optimization accuracy will reach a peak with a
certain $\lambda$. Continuing decreasing $\lambda$ from that point is not
necessary since it will introduce overwhelming noise, which makes this
$\lambda$ optimal.

It should be noted that only a small number of values have to be tested for
$\lambda$, which makes this algorithm very efficient. To be more specific, the
regularization term and objective function $\mbox{tr}(H_{o}\rho(\cdot))$
should be comparable for the optimization to proceed. As a result, any
effective $\lambda$ will be confined within a small range. After a few rounds
of shrinkage, the regularization strength will soon become too small for depth
selection.

##  IV Numerical Results

###  IV.1 Max-Cut and approximation ratio

The performance of depth optimization is demonstrated on Max-Cut problem Farhi
et al. (2015). Consider an $N$-node non-directed and weighted graph $G=(V,E)$.
Max-Cut is the partition of $V$ into two subsets $V_{1}$ and $V_{2}$, for
which the sum of weights of edges between the nodes of two disjoint subsets is
maximized. By assigning $1$ to the nodes of one subset and $-1$ to the nodes
of the other subset, Max-Cut can be formulated as a binary optimization
problem which is equivalent to the minimization of expectation of the
corresponding Hamiltonian

$H_{o}=\sum_{(i,j)\in E}\omega_{ij}\sigma_{z_{i}}\sigma_{z_{j}}.$ (10)

We define the approximation ratio as follows

$r=1-\frac{\min_{x}\mbox{tr}(H_{o}\rho(x))-C_{\min}}{C_{\max}-C_{\min}}\in[0,1],$
(11)

with $C_{\min}$ being the theoretical minimum value and $C_{\max}$ being the
theoretical maximum value. The approximation ratio $r$ is a measure of how
close the final state is to the optimal solution, and a larger $r$ indicates a
better solution.

Figure 1: The randomly generated Max-Cut example used in the numerical
experiment.

The Max-Cut problem to test in this paper is randomly generated as described
in Fig. 1. The weighted graph is made up of 5 vertices and 8 edges. For
simplicity, we set $\gamma=\gamma_{1}=...=\gamma_{N}$ in the simulation
definition.

Figure 2: The performance comparison between (a) noisy QAOA and (b) PG
algorithm under the influence of relaxation with $\gamma=0.2$. Each algorithm
has been iterated for 300 times to obtain one approximation ratio. The
approximation ratios achieved by noisy QAOA for different number of control
parameters are depicted in (a), with an unconspicuous peak at $p=4$ followed
by a slight downtrend till $p=7$. The PG algorithms all start with the same
initial depth but different $\lambda$. (b) shows the finally selected number
of parameters after optimization. The ideal number of control parameters lies
in the region $p=[4,7]$, which is consistent with (a).

###  IV.2 Depth Selection with Relaxation Process

The relaxation process is simulated by letting

$L_{1}=...=L_{n}=\sigma_{-}.$ (12)

That is, the effect of noise for each qubit is to relax the system from the
excited state to the ground state, which constitutes a major source of quantum
decoherence.

The initial value of each control parameter is set as 0.1. The learning rate
is 0.008. In addition, $H_{o}$ and $H_{c}$ are multiplied by a factor of 6 to
accelerate the simulation of dynamics. As shown in Fig. 2, increasing the
control depth for noisy QAOA does not always lead to performance improvement.
On one hand, the noise effect accumulates with the growth of $p$, as the the
total control duration increases with the depth when the initial value of each
control parameter is fixed. On the other hand, it is increasingly hard to
optimize the model with more parameters, which in general will slow down the
pace of improvement once the approximation ratio has reached a certain level.
At some point, the negative effects will outweigh the benefits of increasing
the depth. In Fig. 2(a), the peak value of the approximation ratio is 0.8702,
which has been achieved with 8 control parameters, corresponding to $p=4$. In
addition, the curve shows a plateau after reaching the peak. Due to the
relatively low noise level, it is possible that the increased noise brought by
a few additional control steps can be counteracted with more freedoms in
optimization.

The results of PG algorithm applied on the regularized model (5) are shown in
Fig. 2(b). The initial depth is $p=8$, which corresponds to 16 control
parameters. The finally selected numbers of parameters are obtained with
different values of $\lambda$. It can be seen that the algorithm can
accurately identify the same plateau as Fig. 2(a) that indicates the ideal
number of control parameters. In particular, the leftmost point on the plateau
can be chosen as the solution, since it corresponds to the smallest model
complexity. In practice, when such a plateau exists, we can always stop the
experiment early if the difference between the last approximation ratio and
the current value is less than a given threshold.

Figure 3: The initial control sequence, the finally selected control sequences
with $\lambda=1$ and $\lambda=0.4$ are drawn. The numbers of control
parameters for each sequence are 16, 8 and 14, respectively. The achieved
approximation ratio is marked on the right. Different types of control
operations are represented by different colors. Adjacent control operations of
the same type are combined into a single operation. Figure 4: The performance
comparison between (a) noisy QAOA and (b) PG algorithm under the influence of
relaxation with $\gamma=0.5$. The approximation ratios achieved by 300
iterations are depicted in (a), with respect to different number of control
parameters. (b) shows the finally selected number of parameters using the PG
algorithms and 300 iterations. The approximation ratios are obtained with
different $\lambda$, while the initial depths are the same. The selected
number of control parameters is consistent with the best control depth of
plain QAOA.

The reduced depths and durations for $\lambda=1$ and $\lambda=0.4$ are
visualized in Fig. 3. It can be seen that although the approximation ratios
are similar for the two cases, the control sequences are significantly
different. In particular, the control sequence with 8 selected parameters is
approximately $2/3$ the length of the control sequence with 14 selected
parameters. It should be noted that if one control parameter is penalized to
zero, then the corresponding control operation can be removed from the
sequence. The control operations before and after the removed operation are of
the same type, which means that they can be combined into one operation.
Therefore, 8 selected parameters yield only 6 control operations, which is a
significant reduction in depth when compared to the initial control sequence
that involves 16 control operations.

Figure 5: The orange curve of the finally selected number of parameters is
obtained under the high noise level $\gamma=0.5$ by varying the value of
$\lambda$ from 0.5 to 2, while the blue line is obtained under the low noise
level $\gamma=0.2$.

As the noise level is increased, it becomes increasingly difficult to
counteract the noise effect brought by the additional control steps.
Therefore, once the peak is reached, the approximation ratio is expected to
descend faster than the low-noise case. As shown in Fig. 4(a), an obvious peak
can be found when executing the QAOA with different number of control
parameters. In particular, the optimal number of control parameters is either
5 or 6, which can be precisely identified using the PG algorithm.

As an intelligent algorithm, the PG algorithm is capable of efficiently
determining the best depth with a small set of regularization parameters. That
is, the selection result is robust to small variations in $\lambda$, such that
a sparsely distributed set of regularization parameters are sufficient for
finding the near-optimal solution. As shown in Fig. 5, the selection result is
the same for any $\lambda$ that is confined in the region $\lambda\in(1,2)$.
As a result, as long as we have tried one $\lambda$ in $(1,2)$, the near-
optimal number of control parameters can always be found. For example, if the
initial value of $\lambda$ is chosen as 6 which will be shrinked by a constant
factor of 0.6 after each experiment, i.e., $\lambda=6\rightarrow
3.6\rightarrow 2.16\rightarrow 1.296\rightarrow\cdot\cdot\cdot$, the optimal
depth can be found with around 5 experiments.

Figure 6: The noise level is $\gamma=0.2$. The orange curve shows the selected
numbers of control parameters using 200 iterations of PG descent followed by
100 iterations of vanilla gradient descent. The blue curve shows the selected
numbers of control parameters using 300 iterations of PG descent.

A more precise way to implement the depth selection method is by using the PG
algorithm for the initial iterations, and then fixing the control depth for
further optimization with vanilla gradient descent. By this way, the ultimate
performance of QAOA with the selected depth can be precisely estimated, since
the final accuracy is obtained without the regularization term. Here we test
this implementation by using PG descent for the first 200 iterations and then
employing vanilla gradient descent for the next 100 iterations. As shown in
Fig. 6, an obvious peak appears on the curve of selected parameters which is
in sharp contrast with the results obtained by 300 iterations of PG descent,
and the peak exactly matches the best depth in Fig. 2(a). That is, the updated
selection method could accurately identify the best depth instead of a
plateau.

Figure 7: The performance comparison between (a) noisy QAOA and (b) PG
algorithm under the influence of pure dephasing noise with the coupling
strength being 0.4. Numbers on the right of the enlarged image in (b) are the
finally selected parameters for different $\lambda$.

###  IV.3 Depth Selection with Dephasing Process

The dephasing process is simulated by letting

$L_{1}=...=L_{n}=\sigma_{z}.$ (13)

The experimental setup is the same as the relaxation process except that the
coupling strength is set as 0.4. As shown in Fig. 7(a), the dephasing process
has a strong influence on the approximation ratio, since the approximation
ratio decreases rapidly after reaching the peak. In this case, the optimal
depth is 6. As seen from Fig. 7(b), PG algorithm has successfully identified
the optimal depth with $\lambda=1$. As the regularization strength decreases,
the selected number of parameters increases in general, but the achieved
approximation ratio will begin to decline if the reserved parameters are too
many. These results demonstrate that depth selection algorithm maintains
stable performance under various types of quantum noises.

##  V Conclusion

In this work, we have demonstrated that model selection algorithms such as PG
can be applied to find the optimal depth of QAOA that yields the highest
accuracy under realistic noises. These algorithms are efficient in the sense
that only a few regularization parameters have to be tested before the optimal
depth is found, which significantly reduces the number of experiments as
compared to random search. In particular, it has been shown that the
$l_{1}$-regularized algorithm is capable of determining the best depth with a
sparsely distributed set of regularization parameters which will shrink by a
constant factor after each round of experiment.

Although PG algorithm is specifically considered in this work, the proposed
framework is compatible with any other model selection algorithms and
regularization terms von Luxburg and Schölkopf (2011). Moreover, the
optimization process does not rely on the information of the underlying
dynamics of quantum systems, which makes it easy to implement on near-term
quantum devices.

## Acknowledgements

This research was supported by the National Natural Science Foundation of
China under Grant Nos. 62173296 and 61873162. G. Zhang also acknowledges
support from Hong Kong Research Grant Council (Grants Nos. 1520841, 15203619,
and 15506619), Shenzhen Fundamental Research Fund, China, under Grant No.
JCYJ20190813165207290, and the CAS AMSS-polyU Joint Laboratory of Applied
Mathematics.

## References

  * Preskill [2018] John Preskill.  Quantum Computing in the NISQ era and beyond.  _Quantum_ , 2:79, 2018.  ISSN 2521-327X.
  * Farhi et al. [2014] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  _arXiv:1411.4028_ , 2014.
  * Farhi et al. [2015] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm applied to a bounded occurrence constraint problem.  _arXiv:1412.6062_ , 2015.
  * Farhi and Harrow [2019] Edward Farhi and Aram W Harrow.  Quantum supremacy through the quantum approximate optimization algorithm.  _arXiv:1602.07674_ , 2019.
  * Zhou et al. [2020] Leo Zhou, Sheng-Tao Wang, Soonwon Choi, Hannes Pichler, and Mikhail D. Lukin.  Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices.  _Phys. Rev. X_ , 10:021067, 2020.
  * Borle et al. [2021] Ajinkya Borle, Vincent E. Elfving, and Samuel J. Lomonaco.  Quantum approximate optimization for hard problems in linear algebra.  _arXiv:2006.15438_ , 2021.
  * Jiang et al. [2017] Zhang Jiang, Eleanor G. Rieffel, and Zhihui Wang.  Near-optimal quantum circuit for grover’s unstructured search using a transverse field.  _Phys. Rev. A_ , 95:062317, 2017.
  * Hastings [2018] Matthew Hastings.  Classical and quantum bounded depth approximation algorithms.  _2019 Quantum Information & Computation_, 19:1116–1140, 2018.
  * Wang et al. [2018] Zhihui Wang, Stuart Hadfield, Zhang Jiang, and Eleanor G. Rieffel.  Quantum approximate optimization algorithm for maxcut: A fermionic view.  _Phys. Rev. A_ , 97:022304, 2018.
  * Larkin et al. [2020] Jason Larkin, Matías Jonsson, Daniel Justice, and Gian Giacomo Guerreschi.  Evaluation of QAOA based on the approximation ratio of individual samples.  _arXiv:2006.04831_ , 2020.
  * Bravyi et al. [2020] Sergey Bravyi, Alexander Kliesch, Robert Koenig, and Eugene Tang.  Obstacles to variational quantum optimization from symmetry protection.  _Phys. Rev. Lett._ , 125:260505, 2020.
  * Phillip C. Lotshaw and Siopsis [2021] Rebekah Herrman James Ostrowski Phillip C. Lotshaw, Travis S. Humble and George Siopsis.  Empirical performance bounds for quantum approximate optimization.  _Quantum Information Processing_ , 20(403), 2021.
  * Alam et al. [2020] Mahabubul Alam, Abdullah Ash-Saki, and Swaroop Ghosh.  Design-space exploration of quantum approximate optimization algorithm under noise.  In _2020 IEEE Custom Integrated Circuits Conference (CICC)_ , pages 1–4, 2020.
  * Alam et al. [2020] Mahabubul Alam, Abdullah Ash-Saki, and Swaroop Ghosh.  Circuit compilation methodologies for quantum approximate optimization algorithm.  In _2020 53rd Annual IEEE/ACM International Symposium on Microarchitecture (MICRO)_ , pages 215–228, 2020.
  * Lacroix et al. [2020] Nathan Lacroix, Christoph Hellings, Christian Kraglund Andersen, Agustin Di Paolo, Ants Remm, Stefania Lazar, Sebastian Krinner, Graham J. Norris, Mihai Gabureac, Johannes Heinsoo, Alexandre Blais, Christopher Eichler, and Andreas Wallraff.  Improving the performance of deep quantum optimization algorithms with continuous gate sets.  _PRX Quantum_ , 1:110304, 2020.
  * Herrman et al. [2021] Rebekah Herrman, Phillip C. Lotshaw, James Ostrowski, Travis S. Humble, and George Siopsis.  Multi-angle quantum approximate optimization algorithm.  _arXiv:2109.11455_ , 2021.
  * Dong et al. [2020] Daoyi Dong, Xi Xing, Hailan Ma, Chunlin Chen, Zhixin Liu, and Herschel Rabitz.  Learning-based quantum robust control: Algorithm, applications, and experiments.  _IEEE Transactions on Cybernetics_ , 50(8):3581–3593, 2020.
  * Guerreschi and Smelyanskiy [2017] Gian Giacomo Guerreschi and Mikhail Smelyanskiy.  Practical optimization for hybrid quantum-classical algorithms.  _arXiv:1701.01450_ , 2017.
  * Verdon et al. [2019] Guillaume Verdon, Michael Broughton, and Jacob Biamonte.  A quantum algorithm to train neural networks using low-depth circuits.  _arXiv:1712.05304_ , 2019.
  * Sweke et al. [2020] Ryan Sweke, Frederik Wilde, Johannes Meyer, Maria Schuld, Paul K. Faehrmann, Barthélémy Meynard-Piganeau, and Jens Eisert.  Stochastic gradient descent for hybrid quantum-classical optimization.  _Quantum_ , 4:314, 2020.  ISSN 2521-327X.
  * Liang et al. [2020] Daniel Liang, Li Li, and Stefan Leichenauer.  Investigating quantum approximate optimization algorithms under bang-bang protocols.  _Phys. Rev. Research_ , 2:033402, 2020.
  * Bengtsson et al. [2020] Andreas Bengtsson et al.  Improved success probability with greater circuit depth for the quantum approximate optimization algorithm.  _Phys. Rev. Applied_ , 14:034010, 2020.
  * Campos et al. [2021] E. Campos, D. Rabinovich, V. Akshay, and J. Biamonte.  Training saturation in layerwise quantum approximate optimization.  _Phys. Rev. A_ , 104:L030401, 2021.
  * Rebekah Herrman and Siopsis [2021] Travis S. Humble Rebekah Herrman, James Ostrowski and George Siopsis.  Lower bounds on circuit depth of the quantum approximate optimization algorithm.  _Quantum Information Processing_ , 20(59), 2021.
  * James et al. [2013] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.  _An Introduction to Statistical Learning_.  Springer-Verlag New York, 2013.
  * Burnham and Anderson [2002] Kenneth P. Burnham and David R. Anderson.  _Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach_.  Springer-Verlag New York, 2002.
  * Marshall et al. [2020] Jeffrey Marshall, Filip Wudarski, Stuart Hadfield, and Tad Hogg.  Characterizing local noise in QAOA circuits.  _IOP SciNotes_ , 1(2):025208, 2020.
  * Harrigan et al. [2021] M. P. Harrigan, K. J. Sung, M. Neeley, et al.  Quantum approximate optimization of non-planar graph problems on a planar superconducting processor.  _Nat. Phys._ , 17:332–336, 2021.
  * Breuer and Petruccione [2007] Heinz Peter Breuer and Francesco Petruccione.  _The Theory of Open Quantum Systems_.  Oxford University Press, 2007.
  * Beck and Teboulle [2009] Amir Beck and Marc Teboulle.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM Journal on Imaging Sciences_ , 2(1):183–202, 2009.
  * von Luxburg and Schölkopf [2011] Ulrike von Luxburg and Bernhard Schölkopf.  Statistical learning theory: Models, concepts, and results.  In Dov M. Gabbay, Stephan Hartmann, and John Woods, editors, _Inductive Logic_ , volume 10 of _Handbook of the History of Logic_ , pages 651–706. North-Holland, 2011.


In order to transform this equation into a more manageable form, we look for a
function that, multiplied by the left-hand side, transforms it into the
derivative of a product. That is, we look for a function $\mu(t)$ such that

$\mu(t)\left[\frac{\partial P}{\partial t}+\lambda
P\right]=\frac{\partial}{\partial t}\bigl{[}\mu(t)P\bigr{]}$ (A.59)

It is easy to verify that $\mu(t)=\exp(\lambda{t})$ fits the bill. Equation
(A.58) therefore becomes

$\frac{\partial}{\partial
t}\Bigl{[}e^{\lambda{t}}P(n;t)\Bigr{]}=e^{\lambda{t}}\lambda P(n-1;t)$ (A.60)

For $n=1$ we have

$\frac{\partial}{\partial
t}\Bigl{[}e^{\lambda{t}}P(1;t)\Bigr{]}=e^{\lambda{t}}\lambda
e^{-\lambda{t}}=\lambda$ (A.61)

That is, integrating both sides and multiplying by $e^{-\lambda{t}}$

$P(1;t)=\lambda te^{-\lambda{t}}$ (A.62)

For arbitrary $n$, I’ll show by induction that

$P(n;t)=\frac{(\lambda{t})^{n}}{n!}e^{-\lambda{t}}$ (A.63)

We have already derived the result for $n=0$ and for $n=1$. For arbitrary $n$,
we have

$\begin{array}[]{lcll}\displaystyle\frac{\partial}{\partial
t}\Bigl{[}e^{\lambda{t}}P(n+1;t)\Bigl{]}&=&e^{\lambda{t}}\lambda P(n;t)&\\\
\displaystyle&=&e^{\lambda{t}}\lambda\frac{(\lambda{t})^{n}}{n!}e^{-\lambda{t}}&\mbox{(induction
hypothesis)}\\\ \displaystyle&=&\lambda\frac{(\lambda{t})^{n}}{n!}\end{array}$
(A.64)

So, integrating

$e^{\lambda{t}}P(n+1;t)=\frac{\lambda}{n!}\int(\lambda{t}^{n})dt=\frac{(\lambda{t})^{n+1}}{(n+1)!}+C$
(A.65)

where $C=0$ because of the initial conditions, so

$P(n+1;t)=e^{-\lambda{t}}\frac{(\lambda{t})^{n+1}}{(n+1)!}$ (A.66)

(end of example)

The distribution that results from this example:

$P_{X}(x)=e^{-x}\frac{x^{n}}{n!}$ (A.67)

is the Poisson distribution that, in the example, gives us the probability
that $n$ events take place in a time $x$. Figure A.4 shows the shape of this
distribution as a function of $x$ for various values of $n$.

$0$ $0.05$ $0.1$ $0.15$ $0.2$ $0.25$ $0.3$ $0$ $5$ $10$ $15$ $20$ $25$ $30$
$n=1$$n=5$$n=10$Poisson PDF Figure A.4: The Poisson PDF for various values of
$n$.

Example VII:
The Poisson distribution can also be seen as a limiting case of the binomial
distribution. If $p$ is the probability of success, then $\nu=Np$ is the
expected number of successful trials, as per (A.50). This approximation is
valid for large $N$. In this case, we have

$P(n;N)=\frac{N!}{n!(N-n)!}\left(\frac{\nu}{N}\right)^{n}\left(1-\frac{\nu}{N}\right)^{N-n}$
(A.68)

Taking $N\rightarrow\infty$, we have

$\displaystyle P_{\nu}(n)$ $\displaystyle=\lim_{N\rightarrow\infty}P(n;N)$
(A.69)
$\displaystyle=\lim_{N\rightarrow\infty}\frac{N\cdot(N-1)\cdots(N-n+1)}{n}\frac{\nu^{n}}{N^{n}}\left(1-\frac{\nu}{N}\right)^{N}\left(1-\frac{\nu}{N}\right)^{-n}$
$\displaystyle=\lim_{N\rightarrow\infty}\frac{N\cdot(N-1)\cdots(N-n+1)}{N^{n}}\frac{\nu^{n}}{n!}\left(1-\frac{\nu}{N}\right)^{N}\left(1-\frac{\nu}{N}\right)^{-n}$
$\displaystyle=1\cdot\frac{\nu^{n}}{n!}e^{-\nu}\cdot 1$
$\displaystyle=\frac{\nu^{n}}{n!}e^{-\nu}$

So, once again, we find that the number of successes has a Poisson
distribution.

(end of example)

The characteristic function of the distribution (A.67) is

$\tilde{P}(\omega)=e^{\lambda(e^{i\omega}-1)}$ (A.70)

from which we obtain

$\langle X\rangle=\lambda$ (A.71)

####  A.1.2 Functions of Random Variables

If $X$ is a random variable on $\Omega$, and
$f:\Omega\rightarrow\Omega^{\prime}$, then $Y=f(X)$ is a random variable on
$\Omega^{\prime}$. Here I’ll consider, for the sake of simplicity, the case
$\Omega=\Omega^{\prime}={\mathbb{R}}$ (all our considerations can be
generalized to arbitrary continua $\Omega$ under fairly general conditions,
essentially that $\Omega$ be a metric space). In order to determine the
distribution of $y$, I begin with a preliminary observation. For a random
variable $X$, let ${\mathbb{P}}_{X}[x,x+\Delta{x}]$ the probability that the
value of $X$ falls in $[x,x+\Delta{x}]$. Then, for small $\Delta{x}$,

$\displaystyle{\mathbb{P}}_{X}[x,x+\Delta{x}]$ $\displaystyle=P(X\leq
x+\Delta{x})-P(X\leq x)$ (A.72) $\displaystyle=\frac{\partial}{\partial
x}P(X\leq x)\Delta{x}+O(\Delta{x}^{2})$
$\displaystyle=P_{X}(x)\Delta{x}+O(\Delta{x}^{2})$

Let now $f$ be invertible, and $g=f^{-1}$. Then

$\displaystyle P_{Y}\Delta{y}$ $\displaystyle={\mathbb{P}}_{Y}[y,y+\Delta{y}]$
(A.73) $\displaystyle={\mathbb{P}}_{X}[g(y),g(y+\Delta{y})]$
$\displaystyle\approx{\mathbb{P}}_{X}\Bigl{[}g(y),g(y)+\left|\frac{dg}{dy}\right|\Delta{y})\Bigr{]}$
$\displaystyle=P_{X}(g(y))\left|\frac{dg}{dy}\right|\Delta{y}$

from which we get

$P_{Y}(y)=P_{X}(g(y))\left|\frac{dg}{dy}\right|$ (A.74)

Note that equivalently one could have defined

$P_{Y}(y)=\int\delta(y-f(x))P_{X}\,dx=\langle\delta(y-f(x))\rangle_{X}$ (A.75)

where the subscript on the average reminds us that we are taking the average
with respect to the distribution of $X$. From this, we can determine the
characteristic function of $Y$:

$\displaystyle\tilde{P}_{Y}(\omega)$ $\displaystyle=\int
e^{i\omega{y}}P_{Y}(y)\,dy$ (A.76) $\displaystyle=\int P_{X}(x)\Bigl{[}\int
e^{i\omega{y}}\delta(y-f(x))\,dy\Bigr{]}\,dx$ $\displaystyle=\int
e^{i\omega{f(x)}}P_{X}(x)\,dx$ $\displaystyle=\langle\exp\bigl{[}i\omega
f(x)\bigr{]}\rangle_{X}$

If $Y=aX$, then

$\tilde{P}_{Y}(\omega)=\langle\exp\bigl{[}i\omega
aX\bigr{]}\rangle_{X}-\tilde{P}_{X}(a\omega)$ (A.77)

* * *

Consider now the sum of two random variables: $Z=X+Y$. Each value of $Z$ can
be obtained through an infinity of events: each time $X$ takes an arbitrary
value $x$, and $y$ takes a value $z-x$, $Z$ takes the same value, namely $z$.
Summing up all these possible events we obtain

$P_{Z}(z)=\int_{-\infty}^{\infty}P_{X}(x)P_{Y}(z-x)\,dx$ (A.78)

This is known as the convolution of $P_{X}$ and $P_{Y}$, often indicated as
$P_{Z}=P_{X}*P_{Y}$. The properties of the Fourier transform entail that the
corresponding relation between characteristic functions is

$\tilde{P}_{Z}(\omega)=\tilde{P}_{X}(\omega)\tilde{P}_{Y}(\omega)$ (A.79)

* * *

Let $Y=\\{y_{1},\ldots,y_{n}\\}$ be a set of independent and identically
distributed (i.i.d.) variables with cumulative distribution
${\mathcal{P}}_{Y}$ and density $P_{Y}$. Consider the function $\min(Y)$: we
are interested in finding its density $P_{\min}$ and cumulative distribution
${\mathcal{P}}_{\min}$. We have:

${\mathcal{P}}_{Y}(x)={\mathbb{P}}\bigl{[}\min(Y)\leq{x}\bigr{]}=1-{\mathbb{P}}\bigl{[}\min(Y)\geq{x}\bigr{]}$
(A.80)

We have $\min(Y)\geq{x}$ iff we have $y_{i}\geq{x}$ for all $i$, that is

$\displaystyle{\mathcal{P}}_{\min}(x)$
$\displaystyle=1-{\mathbb{P}}\Bigl{[}\forall y\in{Y}.y\geq{x}\Bigr{]}$ (A.81)
$\displaystyle=1-{\mathbb{P}}\bigl{[}y\geq{x}\bigr{]}^{n}$
$\displaystyle=1-\Bigl{(}1-{\mathbb{P}}\bigl{[}y\leq{x}\bigr{]}\Bigr{)}^{n}$
$\displaystyle=1-\Bigl{(}1-{\mathcal{P}}_{Y}(x)\Bigr{)}^{n}$

The density is

$\displaystyle P_{\min}(x)$
$\displaystyle=\frac{d}{dx}{\mathcal{P}}_{\min}(x)$ (A.82)
$\displaystyle=n\Bigl{(}1-{\mathcal{P}}_{Y}(x)\Bigr{)}^{n-1}\frac{d}{dx}{\mathcal{P}}_{Y}(x)$
$\displaystyle=n\Bigl{(}1-{\mathcal{P}}_{Y}(x)\Bigr{)}^{n-1}P_{Y}(x)$

For the function $\max(Y)$, working in a similar way, we have

$\displaystyle{\mathcal{P}}_{\max}(x)$
$\displaystyle=({\mathcal{P}}_{Y}(x))^{n}P_{\max}(x)$
$\displaystyle=n({\mathcal{P}}_{Y}(x))^{n-1}P_{Y}(x)$ (A.83)

####  A.1.3 The Central Limit Theorem

The Central Limit Theorem (important enough to be granted its own acronym:
CLT) is one of the fundamental results in basic probability theory and the
main reason why the Gaussian distribution is so important and so common in
modeling natural events. In a nutshell, the theorem tells us the following: if
we take a lot of random variables, independent and identically distributed
(i.i.d.), and add them up, the result will be a random variable with Gaussian
distribution. So, for example, if we repeat an experiment many times and take
the average of the results that we obtain (the average is, normalization
apart, a sum), no matter what the characteristics of the experiment are, the
resulting average will have (more or less) a Gaussian distribution.

But, ay, there’s the rub! The theorem works only in the assumption that the
moments of the distributions involved be finite. We shall see shortly what
happens if this assumption is not satisfied.

Let $X_{1},\ldots,X_{n}$ be a set of i.i.d. random variables with distribution
$P_{X}$, zero mean, and (finite) variance $\sigma^{2}$. Note that
$Y=\sum_{i}X_{i}$ has zero mean and variance $n\sigma^{2}$, while
$Y=(\sum_{i}X_{i})/n$ has zero mean and variance $\sigma^{2}/n$. It is
therefore convenient to work with the variable

$Z_{n}=\frac{1}{\sqrt{n}}\sum_{i}X_{i}$ (A.84)

which has zero mean and variance $\sigma^{2}$ independently of $n$.

######  Theorem A.1.1.

For any distribution $P_{X}$ with finite mean and variance, and
$X_{1},\ldots,X_{n}$ i.i.d. with distribution $P_{X}$, for
$n\rightarrow\infty$, we have $Z_{n}\rightarrow{Z_{\infty}}$, where
$Z_{\infty}$ is a Gaussian random variable with zero mean and variance
$\sigma^{2}$ equal to the variance of $P_{X}$.

###### Proof.

Consider the first terms of the expansion of the characteristic function of
$P_{X}$:

$\tilde{P}_{X}(\omega)=\int
e^{i\omega{x}}P_{X}(x)\,dx=1-\frac{1}{2}\sigma^{2}\omega^{2}+O(\omega^{3})$
(A.85)

The characteristic function of $Y=\sum_{i}X_{i}$ is given by (A.79):

$\tilde{P}_{Y}(\omega)=\prod_{i}\tilde{P}_{X_{i}}(\omega)=\bigl{[}\tilde{P}_{X}(\omega)\bigr{]}^{n}$
(A.86)

(the second equality holds because the $X$s have the same distribution) while
(A.77) with $a=1/\sqrt{n}$ gives

$\tilde{P}_{Z}(\omega)=P_{Y}\left(\frac{\omega}{\sqrt{n}}\right)=\left[P_{Y}\left(\frac{\omega}{\sqrt{n}}\right)\right]^{n}\approx\left(1-\frac{\sigma^{2}\omega^{2}}{2n}\right)^{n}\stackrel{{\scriptstyle
n\rightarrow\infty}}{{\longrightarrow}}\exp(-\frac{1}{2}\sigma^{2}\omega^{2})$
(A.87)

Finally, from (A.16) we have the inverse transform

$P_{Z}(z)=\frac{1}{\sigma\sqrt{2\pi}}\exp\Bigl{(}-\frac{z^{2}}{2\sigma^{2}}\Bigr{)}$
(A.88)

∎

This theorem is true, in the form in which we have presented it, only for
distributions $X$ with finite mean and variance222I have assumed zero mean
since, if the mean of the $X$ is non-zero, the mean of $Z$ goes to infinity;
this doesn’t represent a major hurdle for the theorem, which can easily be
generalized by subtracting the mean from the variables $X$ and then adding it
back.. However, the key to the theorem is an invariance property of the
characteristic function of the Gaussian. Consider the equality (A.86); we can
split it up as:

$\tilde{P}_{Z}(\omega;n)=\bigl{[}\tilde{P}_{X}(\omega)\bigr{]}^{n}=\bigl{[}\tilde{P}_{X}(\omega)\bigr{]}^{n/2}\bigl{[}\tilde{P}_{X}(\omega)\bigr{]}^{n/2}=\tilde{P}_{Z}(\omega;n/2)\tilde{P}_{Z}(\omega;n/2)$
(A.89)

Taking the limit $n\rightarrow\infty$, this gives us
$P_{Z}(\omega)=P_{Z}(\omega)P_{Z}(\omega)$. That is: the condition for a
distribution to be a central limit is that the product of two characteristic
functions have the same functional form as the original distributions. As we
have seen in (A.20), the Gaussian distribution does have this property. Nay:
it is the only distribution with finite moments that has this property, hence
its appearance in the theorem in the finite moments case, and hence its great
importance in application as a model of many processes resulting from the sum
of identical sub-processes.

If we abandon the finite moment hypothesis, however, there is a more general
distribution to which (A.89) applies: the stable Levy distribution. So, a more
general form of the CLT can be enunciated as:

######  Theorem A.1.2.

For any distribution $P_{X}$, and $X_{1},\ldots,X_{n}$ i.i.d. with
distribution $P_{X}$, for $n\rightarrow\infty$, we have

$\lim_{n\rightarrow\infty}\frac{1}{\sqrt{n}}\sum_{i=1}^{n}X_{i}=Z_{\infty}$
(A.90)

where $Z_{\infty}$ is a random variable with Levy distribution. If the
variance of $P_{X}$ is finite and equal to $\sigma^{2}$, then $Z_{\infty}$ has
a Gaussian distribution with variance $\sigma^{2}$.

####  A.1.4 Stochastic Processes

A stochastic process is a set of random variables $X(t)$ indexed by a variable
$t$ (commonly identified with time) that takes value either in ${\mathbb{N}}$
or ${\mathbb{R}}^{+}$ (less frequently in ${\mathbb{R}}$). We indicate with
$P(x,t)$ the probability that the process take value $x$ at time $t$ (the
probability density if $t$ is continuous; I shall omit the subscript $X$ to
avoid complicating the notation), and with $P(x_{2},t_{2};x_{1},t_{1})$ the
joint probability density for the two variables $X(t_{1})$ and $X(t_{2})$. The
multiple joint probability density $P(x_{1},t_{1};\ldots,x_{n},t_{n})$ is
defined analogously. In the following, whenever possible, I shall use the
joint probability $P(x_{2},t_{2};x_{1},t_{1})$ to simplify the notation, but
all considerations hold for the more general multiple joint probability.

Just as a stochastic variable is instantiated to a specific value $x\in\Omega$
with a certain probability, so a stochastic process is instantiated as a
trajectory $X:{\mathbb{R}}\rightarrow\Omega$ (or with a discrete series
$X:{\mathbb{N}}\rightarrow\Omega$ if the process is discrete). Each $X(t)$,
for fixed $t$, is a stochastic variable with a probability distribution that,
in general, depends on $t$. A stochastic process is stationary if all these
distributions are the same, that is, $P(x,t)\equiv{P(x)}$, or, equivalently,
if

$P(x_{1},t_{1};x_{2},t_{2})=P(x_{1},t_{1}+\tau;x_{2},t_{2}+\tau)$ (A.91)

In a stochastic process, there are two ways of computing averages: one can
compute the ensemble average $\langle{X(t)}\rangle$, that is, the average of
the random variable $X(t)$, or the mean value along a trajectory

$\bar{X}=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{0}^{T}x(t)\,dt$ (A.92)

A process is ergodic if the two coincide

$\langle X\rangle=\bar{X}$ (A.93)

Ergodicity is an important property for random walks: many times we are
interested in the characteristics of the motion of one individual, but many of
the equations that we shall use involve ensemble probabilities based on a
whole population. Ergodicity allows us to switch from one to the other with
impunity.

Note that in a stationary process the correlation
$\langle{X(t_{1})}X(t_{2})\rangle$ does not depend on $t_{1}$ and $t_{2}$
individually, but only on their difference $\tau=t_{2}-t_{1}$. Joint
probabilities are positive, symmetric
($P(x_{1},t_{1};x_{2},t_{2})=P(x_{2},t_{2};x_{1},t_{1})$) and normalized:

$\int\\!\\!\\!\int_{\Omega^{2}}P(x_{1},t_{1};x_{2},t_{2})\,dx_{1}\,dx_{2}=1$
(A.94)

Joint probabilities can be reduced by integration

$P(x_{1},t_{1})=\int_{\Omega}P(x_{1},t_{1};x_{2},t_{2})\,dx_{2}$ (A.95)

and Bayes theorem can be extended to stochastic processes

$P(x_{1},t_{2}=\int_{\Omega}P(x_{2},t_{2}|x_{1},t_{1})P(x_{1},t_{1})\,dx_{1}$
(A.96)

A process if Markov if, for all $t_{1}<t_{2}<\cdots<t_{n}$,

$P(x_{n},t_{n}|x_{n-1},t_{n-1};\ldots;x_{1},t_{1})=P(x_{n},t_{n}|x_{n-1},t_{n-1})$
(A.97)

this entails that, at any time, the status of the process encodes all the
information necessary to make predictions about its future: it is not
necessary to know how the process reached that status. The Markov property can
be chained:

$\displaystyle P(x_{1},t_{1};x_{1},t_{2};x_{3},t_{3})$
$\displaystyle=P(x_{2},t_{2};x_{3},t_{3}|x_{1},t_{1})P(x_{1},t_{2})$ (A.98)
$\displaystyle=P(x_{3},t_{3}|x_{2},t_{2})P(x_{2},t_{2}|x_{1},t_{1})P(x_{1},t_{1})$

Finally, it can be shown [19] that Markov processes must satisfy the Chapman-
Kolmogorov equation:

$P(x_{3},t_{3}|x_{1},t_{1})=\int_{\Omega}P(x_{3},t_{3}|x_{2},t_{2})P(x_{2},t_{2}|x_{1},t_{1})\,dx_{2}$
(A.99)

The characteristics of the Markov process is evidenced by the fact that the
probability of the transition
$(x_{1},t_{1})\rightarrow(x_{2},t_{2})\rightarrow(x_{3},t_{3})$ is the product
of the probabilities of the transitions
$(x_{1},t_{1})\rightarrow(x_{2},t_{2})$ and
$(x_{2},t_{2})\rightarrow(x_{3},t_{3})$, that is, the two transitions are
statistically independent.

####  A.1.5 Gaussian and Wiener processes

I shall provide here some details on two types of processes of considerable
importance for random walks and diffusion.

A stochastic process $X(t)$ is Gaussian with zero mean if
$\langle{X(t)}\rangle=0$ and

$P(x_{i},t_{i})=\sqrt{\frac{A_{ii}}{2\pi}}\exp\Bigl{(}-\frac{1}{2}A_{ii}x^{2}\Bigr{)}$
(A.100)

($A_{ii}>0$). The joint probability $P(x_{1},t_{1};\ldots;x_{n},t_{n})$ then
follows a multivariate Gaussian distribution

$P(x_{1},t_{1};\ldots;x_{n},t_{n})=\frac{\mbox{det}(\mathbf{A})^{1/2}}{(2\pi)^{n/2}}\exp\left[-\frac{1}{2}\sum_{i,j=1}^{n}x_{i}A_{ij}x_{j}\right]$
(A.101)

Where $\mathbf{A}\in{\mathbb{R}}^{n\times{n}}$ is symmetric (strictly)
positive definite. The matrix $\mathbf{A}$ is a measure of the covariance
between two variables of the Gaussian process

$\langle X(t_{i})X(t_{j})\rangle=(\mathbf{A}^{-1})_{ij}$ (A.102)

(this is true since we assume zero mean). A process is uncorrelated if
$\langle{X(t_{i})}X(t_{j})\rangle=D\delta(t_{i}-t_{j})$, in which case
$A_{ij}=D^{-1}\delta_{ij}$.

A Wiener process $W$ is a process in which the variables $W(t)$ are real and
with independent increments $W(t_{2})-W(t_{1})$ that follow a Gaussian
distribution. That is, they define a conditional probability

$P(w_{2},t_{2}|w_{1},t_{1})=\frac{1}{\sigma\sqrt{2\pi(t_{2}-t_{1})}}\exp\left[-\frac{(w_{2}-w_{1})^{2}}{2\sigma^{2}(t_{2}-t_{1})}\right]$
(A.103)

from which the covariance can be computed

$\displaystyle\langle(W(t_{2})-\langle{W}\rangle)(W(t_{1})-\langle{W}\rangle)\rangle$
$\displaystyle=\langle(W(t_{2})-W(0))(W(t_{1})-W(0))\rangle$ (A.104)
$\displaystyle=\int_{-\infty}^{\infty}(w_{2}-w_{0})\,dw_{2}\int_{-\infty}^{\infty}\,dw_{1}(w_{1}-w_{0})P(w_{2},t_{2};w_{1},t_{1})$
$\displaystyle=\sigma^{2}\min(t_{1},t_{2})+w_{0}^{2}$

From this we get

$\langle W(t)^{2}\rangle=\sigma^{2}t+w_{0}^{2}$ (A.105)

Wiener processes are related to Gaussian processes, in particular to
uncorrelated (white) Gaussian processes. Let $X(t)$ be a Gaussian process with
$\langle{X(t_{1})}X(t_{2})\rangle=\sigma^{2}\delta(t_{2}-t_{1})$, and define a
new stochastic process as the integral of $X(t)$:

$Y(t)=\int_{0}^{t}X(u)\,du$ (A.106)

then

$\displaystyle\langle Y(t_{2})Y(t_{1})\rangle$
$\displaystyle=\int_{0}^{t_{2}}\\!\\!\\!\\!\\!du_{2}\int_{0}^{t_{1}}\\!\\!\\!\\!\\!du_{1}\langle{X(u_{1})}X(u_{2})\rangle$
(A.107)
$\displaystyle=\int_{0}^{t_{2}}\\!\\!\\!\\!\\!du_{2}\int_{0}^{t_{1}}\\!\\!\\!\\!\\!du_{1}\delta(u_{2}-u_{1})$

By the properties of the Dirac function

$\int_{0}^{t_{1}}\\!\\!\\!\\!\\!du_{1}\delta(u_{2}-u_{1})=\begin{cases}1&0<u_{2}<t_{1}\\\
0&\mbox{otherwise}\end{cases}$ (A.108)

Then

$\langle Y(t_{2})Y_{(}t_{1})\rangle=\sigma^{2}\min(t_{2},t_{1})$ (A.109)

which coincides with (A.105) for $w_{1}=0$. That is, the integral of a
Gaussian process is a Wiener process.