# ITSO: A novel Inverse Transform Sampling-based Optimization algorithm for
stochastic search

Nikolaos P. Bakas, Vagelis Plevris, Andreas Langousis, Savvas A.
Chatzichristofis Computation-based Science and Technology Research Center, The
Cyprus Institute, 20 Konstantinou Kavafi Street, 2121, Aglantzia Nicosia,
Cyprus. e-mail<EMAIL_ADDRESS>of Civil Engineering and Energy
Technology, OsloMet—Oslo Metropolitan University, Pilestredet 35, Oslo 0166,
Norway. e-mail<EMAIL_ADDRESS>of Civil Engineering, University
of Patras, 265 04 Patras, Greece. e-mail<EMAIL_ADDRESS>Systems Lab & Department of Computer Science, Neapolis University Pafos, 2
Danais Avenue, 8042 Pafos, Cyprus. e-mail<EMAIL_ADDRESS>
###### Abstract

Optimization algorithms appear in the core calculations of numerous Artificial
Intelligence (AI) and Machine Learning methods, as well as Engineering and
Business applications. Following recent works on the theoretical deficiencies
of AI, a rigor context for the optimization problem of a black-box objective
function is developed. The algorithm stems directly from the theory of
probability, instead of a presumed inspiration, thus the convergence
properties of the proposed methodology are inherently stable. In particular,
the proposed optimizer utilizes an algorithmic implementation of the
$n$-dimensional inverse transform sampling as a search strategy. No control
parameters are required to be tuned, and the trade-off among exploration and
exploitation is by definition satisfied. A theoretical proof is provided,
concluding that only falling into the proposed framework, either directly or
incidentally, any optimization algorithm converges in the fastest possible
time. The numerical experiments, verify the theoretical results on the
efficacy of the algorithm apropos reaching the optimum, as fast as possible.

Keywords: Stochastic Optimization, Inverse Transform Sampling, Black-box
Function, Global Convergence.

##  1 Introduction

Despite the numerous research works and industrial applications of Artificial
Intelligence algorithms, they have been criticized about lacking a solid
theoretical background [1]. The empirical results demonstrate impressive
performance, however, their theoretical foundation and analysis are often
vague [2]. Machine Learning (ML) models, frequently aim at identifying an
optimal solution [3, 4], which is computationally hard and attempts to explain
the procedure often based on the evaluation of the function’s Gradients [5,
6]. Accordingly, the identification of whether a stochastic algorithm will
work or not remains an open question for many research and real-world
applications [7, 8, 9, 10]. A vast number of optimization algorithms have been
developed and applied for the solution of the corresponding problems [11, 12,
13, 14], as well as mathematical proofs regarding their algorithmic
convergence [15, 16, 17], however, their basic formulation often stems from
nature-inspired procedures [18, 19] and not solid mathematical frameworks.
Accordingly, a vast number of research works have been published in order to
investigate the performance of black-box algorithms [20, 21, 22].

In its elementary form, the purpose of efficient optimization algorithms is to
find the argument yielding the minimum value of a black-box function $f(x)$,
defined on a set $A$, $f:A\to\mathbb{R}^{n}$. Accordingly, the inverse problem
of maximization is the minimization of the negated function $-f(x)$, while
problems with multiple objective functions often utilize single function
optimization algorithms to attain the best possible solution. $A$ is assumed a
compact subset of the Euclidean space ${{\mathbb{R}}^{n}}$, where $n$ is the
number of dimensions of the set $A$, however, the proposed method applies
similarly to discrete and continuous topological spaces. The unknown, black-
box function $f$, returns values for the given input
$x_{ij}={{\mathbf{x}}_{i}}=({{x}_{i1}},{{x}_{i2}},\cdots,{{x}_{in}})$ at each
computational discrete time step $i=\left\\{1,2,\ldots,{{f}_{e}}\right\\}$,
where ${{f}_{e}}$ is the number of maximum function evaluations. The sought
solution is a vector ${{\mathbf{x}}_{min}}\in A$, such that
$f({{\mathbf{x}}_{min}})\leq f(\mathbf{x}),\forall\mathbf{x}\in A$, which may
be written by

$\begin{split}{{\mathbf{x}}_{min}}={\operatorname{arg\,min}}\,f(\mathbf{x}):=\\{\mathbf{x}\in
A\subseteq\mathbb{R}^{n}\\\ \mid\forall\mathbf{y}\in A:f(\mathbf{y})\geq
f(\mathbf{x})\\}\end{split}$ (1)

The purpose of this work is to provide a rigor context for the optimization
problem of a black-box function, by adhering to the Probability Theory, aiming
at identifying the best possible solution $\mathbf{x}_{min}$, within the given
iterations $f_{e}$, during the execution of the algorithm.

The rest of the paper is organized as follows: In Section 2, the proposed
Inverse Transform Sampling Optimizer (ITSO) is presented. Additionally, the
same Section provides in details the theoretical formulation of the algorithm
as well as some programming and implementation techniques. Illustrative
examples of the optimization history are also comprised. In Section 3 the
theoretical proof of convergence is provided, as well as Lemma 1, deriving
that the suggested optimization framework is the fastest possible. The
numerical experiments are divided into three groups. Subsection 5 is about the
comparison with 13 nonlinear loss functions, 17 optimization methods, for 10
and 20 dimensions of search space, and 5000 and 10000 iterations per
dimension. Section 4, briefly presents the programming techniques that were
investigated, in order to implement the proposed method into a computer code.
Finally, the conclusions are drawn in Section 6.

##  2 Optimization by Inverse Transform Sampling

Let the probability distribution of the optimal values of $f$, be considered
as the product of some monotonically decreasing kernel $k$ over $f$ at some
time-step $i$ and hence given input ${{\mathbf{x}}_{i}}$. This assumption
stands in the foundation of the method and can be considered as rational,
instead of an arbitrary selection strategy, inspired by natural or other
phenomena. It is a straightforward application of the Probability theory to
the problem. The kernel may be considered as parametric, concerning time $i$.
The selection of the kernel $k$ should satisfy the condition that its limit is
the Dirac delta function centered at ${\operatorname{arg\,min}}\,f$,

${\mathop{\lim_{i\to\infty}}}\,k_{i}(f)={{\delta}_{m}}(x),$ (2)

where $\delta_{m}$ denotes the Dirac function centered at
${\operatorname{arg\,min}}\,f$. Although the selection of the kernel $k$ is
ambiguous, it can be chosen among a variety of functions satisfying Equation
2, such as the Gaussian:

$k_{i}(f)=exp(-{{f(\mathbf{x}_{i})}^{2}}g(i)),$ (3)

where $g(i)$ is a time increasing pattern. The function $g(i)$ controls the
shape of the kernel $k$, approximating numerically Equation 2. Additionally,
we may use:

$k_{i}({f})={\operatorname{max}}\,f-f(\mathbf{x}_{i}),$ (4)
$k_{i}({f})=1-\frac{f(\mathbf{x}_{i})-{\operatorname{min}}\,f+e_{i}}{{\operatorname{max}}\,f-{\operatorname{min}}\,f+e_{i}}\wedge{e_{i}\to
0},$ (5)

or any other non-negative Lebesgue-integrable function, which reverses the
order of the given set of all ${{f}_{{\hat{i}}}}$, where $\hat{i}$ is the
permutation of the indices $1,2,...,i$, such that the sequence of
${{f}_{{\hat{i}}}}$ being monotonically strictly decreasing. The duplicate
values of ${{f}_{{\hat{i}}}}$ should be extracted to avoid numerical
instabilities. These duplicates often appeared in the empirical calculations,
especially when the algorithm was close to a local or global stationary point.
A variety of kernel functions $k$ were investigated, and the results weren’t
affected significantly, even when distorting $k(f)$ with some random noise
${X}^{\prime}\sim\mathcal{U}(a,b)\forall k({{f}_{{\hat{i}}}})\in(a,b)$.

Accordingly, for each dimension $j$ of the vector space $A$, the marginal
probability density function ${P}_{{{X}_{j}}}$ is obtained numerically from
the kernel function:

${{P}_{{{X}_{j}}}}({x}_{ij})=k(f({{x}_{ij}}))\forall j\in\\{1,2,...n\\}.$ (6)

${{P}_{{{X}_{j}}}}({{x}_{ij}})dx$ is the probability that
${\operatorname{arg\,min}}\,f$ falls within the infinitesimal interval
$[{{x}_{ij}}-dx/2,{{x}_{ij}}+dx/2]$. The corresponding cumulative distribution
function (CDF) ${{F}_{{{X}_{j}}}}$, can be calculated by:

${{F}_{{{X}_{j}}}}({{x}_{ij}})=\int_{l{{b}_{j}}}^{{{x}_{ij}}}{{{P}_{{{X}_{j}}}}}({{\xi}})d\xi,$
(7)

where $l{{b}_{j}}$ is the lower bound of the $j^{th}$ dimension of the vector
space $A$ and can be numerically evaluated by some numerical integration rule,
such as the Riemann sum
$S_{j}=\sum\limits_{i=1}^{n}{{{P}_{{{X}_{j}}}}}(x_{ij}^{*})\Delta{{x}_{ij}}$,
with ${{{P}_{{{X}_{j}}}}}(x_{ij}^{*})$ computed by the application of the
kernel $k$ on some ${x}_{ij}$, such that
$f(x_{ij}^{*})=(f({{x}_{ij}})+f({{x}_{i-1,j}}))/2$, or another approximation
scheme. In the following pseudo-code 1, the algorithmic implementation of the
Inverse Transform Sampling method is demonstrated. The symbols are also noted
in the Nomenclature section.

A graphical demonstration of the evolution of the Probability Density, as well
as the corresponding Cumulative Distribution Functions, is presented in Figure
1, for $f(x)=(x-5)^{2}$, which is function with one extreme value and in
Figure 2, for $f(x)=sin(x+0.7)+0.01*(x+0.7)^{2}$, which has many extrema.
Interestingly, as the function evaluations increase, the CDF, is characterised
by sharp slopes, which are positioned in regions where the PDF exhibits high
values, and hence function $f(x)$ attends its lows. Figures 1 and 2, offer an
intuitive representation of the procedure for finding the minimum of the
function, within the suggested framework.

Data: $A,{{f}_{e}},n$

Result: ${\mathbf{x}^{*}={\operatorname{arg\,min}}\,f}$,
${{f}^{*}}=f(\mathbf{x}^{*})$

1 while  _$i\leq{{f}_{e}}$_ do

2 $j\leftarrow\mathcal{U}\left\\{1,n\right\\}$;

3 $r_{i}\leftarrow\mathcal{U}(0,1)$;

4 SORT $x_{ij}\hskip 5.69054pt\forall i$;

5 ${{x}_{ij}}\leftarrow{F_{j}}^{-1}(r_{i})$;

6 ${{f}_{i}}=f(\mathbf{x}_{i})$;

7 if  _${{f}_{i}}\leq{{f}^{*}}$_ then

8 ${{f}^{*}}\leftarrow{{f}_{i}}$;

9 ${{x}^{*}_{j}}\leftarrow{{x}_{ij}}$;

10

11 else

12 ${{x}_{ij}}\leftarrow x^{*}_{j}$;

13 end if

14

15 end while

return $\mathbf{x}^{*}={\operatorname{arg\,min}}\,f$

Algorithm 1 ITSO-Mathematical Framework

(a) PDF, $i=3$ (b) CDF, $i=3$ (c) PDF, $i=100$ (d) CDF, $i=100$

(e) PDF, $i=350$ (f) CDF, $i=350$ (g) PDF, $i=500$ (h) CDF, $i=500$

Figure 1: Probability Density and Cumulative Distribution Functions, utilizing
ITSO search strategy, for function $f(x)=(x-5)^{2}$. The shape sequentially
tends to the Heaviside function, centered at $x_{m}=5$

(a) PDF, $i=3$ (b) CDF, $i=3$ (c) PDF, $i=150$ (d) CDF, $i=150$

(e) PDF, $i=250$ (f) CDF, $i=250$ (g) PDF, $i=500$ (h) CDF, $i=500$

Figure 2: Probability Density and Cumulative Distribution Functions, utilizing
ITSO search strategy, for function $f(x)=sin(x+0.7)+0.01*(x+0.7)^{2}$. The
shape sequentially tends to the Heaviside function, centered at $x_{m}=-2.24$

##  3 Convergence Properties

During the optimization process, the algorithm generates some input variables
$x_{ij}$ as arguments for the black-box function $f$. By utilizing the values
of $f$, we may compute the values of the distribution
${{P}_{{{X}_{j}}}}(x_{ij})$, by Equation 6, and
${{F}_{{{X}_{j}}}}({{x}_{ij}})$ by Equation 7, for all $x_{ij}$.

######  Definition 1.

We define the argument of $f$ within $f_{e}$ iterations, corresponding to the
minimum of $f$ as $x_{m}={\operatorname{arg\,min}}\,f$.

In iteration $i$, the algorithm will have searched the space $\hat{A}\subset
A$, within the limits $lb_{j},ub_{j}$. By definition, $P$ is the probability
density (likelihood) that the optimum occurs in a region
$\mathbf{x}\pm\mathbf{dx}$. Hence,

$E[\mathbf{X}]=\begin{Bmatrix}E[X_{1}]\\\ E[X_{2}]\\\ \ldots\\\
E[X_{n}]\end{Bmatrix}=\begin{Bmatrix}\int_{lb_{1}}^{ub_{1}}{\xi}P_{X_{1}}(\xi)d\xi\\\
\int_{lb_{2}}^{ub_{2}}{\xi}P_{X_{2}}(\xi)d\xi\\\ \ldots\\\
\int_{lb_{n}}^{ub_{n}}{\xi}P_{X_{n}}(\xi)d\xi\end{Bmatrix},$ (8)

where $E[\cdot]$, denotes the expectation of a random variable or vector.

######  Theorem 1.

If $P_{X_{j}}$ tends to Dirac $\delta_{m}$, then ITSO will converge to $x_{m}$

###### Proof.

For each dimension $j$, we may write

$E[X_{j}]=\int_{lb_{j}}^{ub_{j}}{\xi}P_{X_{j}}(\xi)d\xi=\int_{lb_{j}}^{ub_{j}}{\xi}{F^{\prime}_{X_{j}}}(\xi)d\xi,$
(9)

and integrating by parts, we obtain

$\begin{split}E[X_{j}]=[\xi
F_{X_{j}}(\xi)]_{lb_{j}}^{ub_{j}}-\int_{lb_{j}}^{ub_{j}}{1}F_{X_{j}}(\xi)d\xi=ub_{j}*1-lb_{j}*0\\\
-\int_{lb_{j}}^{ub_{j}}{1}F_{X_{j}}(\xi)d\xi.\end{split}$ (10)

If we apply the theorem of the antiderivative of inverse functions [23], we
obtain

$\int_{0}^{1}{{{F_{X_{j}}}^{-1}}}(y)dy+\int_{lb_{j}}^{ub_{j}}{F_{X_{j}}}(x)dx=ub_{j}*1-lb_{j}*0.$
(11)

Hence, for Equations 10 and 11 we deduce that

$E[X_{j}]=\int_{0}^{1}{{{F_{X_{j}}}^{-1}}}(y)dy.$ (12)

With subscript $m$ denoting that the Dirac function is centered at the
argument that minimizes $f$, i.e.

${{\delta}_{m}}(x)=\left\\{\begin{array}[]{*{35}{l}}+\infty,&x={\operatorname{arg\,min}}\,f\\\
0,&x\neq{\operatorname{arg\,min}}\,f\\\ \end{array}\right.,$ (13)

and

$H_{m}(x):=\int_{-\infty}^{x}{\delta_{m}(s)}ds.$ (14)

As $P_{X_{j}}$ tends to Dirac function, $F_{X_{j}}$ tends to the Heaviside
step function centered at $x_{m}$, hence by Equation 12 we deduce that

$E[X_{j}]\xrightarrow[]{i\to f_{e}}\int_{0}^{1}{{{H_{m}}^{-1}}}(y)dy,$ (15)

and by Equation 11

$E[X_{j}]\xrightarrow{}ub_{j}-\int_{0}^{1}{{H_{m}}}(y)dy=ub_{j}-(ub_{j}-x_{m})*1=x_{m}$
(16)

∎

######  Lemma 1.

_(ITSO Convergence Speed)_ ITSO is the fastest possible optimization framework

###### Proof.

Any distribution that doesn’t tend to Dirac, could be considered as a Dirac
plus a positive function of $x$. In this case, the algorithmic framework would
search through a strategy that produces a ${P}^{\prime}$ over $\mathbf{x}$, as

${P}^{\prime}={{P}^{*}}\pm{{\delta}_{m}}.$ (17)

Hence, with $F^{*}$ indicating the CDF corresponding to $P^{*}$, Equations 15,
and 12, result in

$E[X]\xrightarrow[]{i\to
f_{e}}\int_{0}^{1}{{{H_{m}}^{-1}}}(y)dy+\int_{lb}^{ub}{F^{*}(x)dx},$ (18)

.

and thus, by Equations 17, and 16, we obtain

$E[X]\xrightarrow{}x_{m}\pm\epsilon.$ (19)

Hence the algorithm would converge to a point different than $x_{m}$

∎

##  4 Programming techniques

A variety of programming techniques were investigated, in order to implement
the proposed method into a computer code. To keep the algorithm simple and
reduce the computational time, we applied inverse transform sampling by
keeping in each iteration $i$ the best function evaluations, and randomly
sampling among them. This is equivalent to a kernel function that vanishes
over the worst function evaluations, and distributing the probability mass to
the best performing ones. Accordingly, similar programming techniques may be
investigated in future works, within the suggested framework.

1 Initialize:
$\mathbf{x}=\text{rand}(n),\quad\mathbf{opti\\_x}=\mathbf{x},\quad{opti\\_f}=f(\mathbf{\mathbf{opti\\_x}})$;

2 for  _$i=1:f_{e}$_ do

3 for  _$j=1:n$_ do

4 $\mathbf{indsBEST}=\text{sortperm}({opti\\_evals})[1:\alpha]$;

5 $\mathbf{inp\\_x=x\\_evals}[\mathbf{indsBEST},j]$;

6 $\mathbf{inp\\_y=opti\\_evals}[\mathbf{indsBEST}]$;

7 ${rr}=\text{min}\\{\mathbf{inp\\_x\\}}$

8
$\quad+\text{rand}(0,1)(\text{max}\\{\mathbf{inp\\_x\\}}-\text{min}\\{\mathbf{inp\\_x\\}})$;

9 $\mathbf{x}[j]={rr}$;

10 $f_{i}=f(\mathbf{x})$;

11 if  _${{f}_{i}}\leq$ opti_f_ then

12 ${opti\\_f}={{f}_{i}}$;

13 $\mathbf{opti\\_x}=\mathbf{x}$;

14

15 else

16 $\mathbf{x}=\mathbf{opti\\_x}$;

17 end if

18

19 end for

20

21 end for

return $\quad{opti\\_f},\quad\mathbf{opti\\_x}$

Algorithm 2 ITSO-Short Pseudocode

The Algorithm 2 represents a simple version of the supplementary code in the
appendix, which may easily be programmed. The variable $\mathbf{opti\\_evals}$
is a dynamic vector, containing all values of the objective function returned
during the optimization history, until step $i$, and $\mathbf{x\\_evals}$ is
an $i\times n$ matrix, containing all the design vectors $\mathbf{x}_{1:i}$,
corresponding to $\mathbf{opti\\_evals}$. The integer $\alpha$ is a parameter
regarding how many instances of the optimization history are kept in order to
randomly sample among them; for example if $f_{e}=10^{4}$, we may select
$\alpha=10^{2}$. In this variation of the code, the Inverse Transform Sampling
is approximately implemented in line 7, by calculating the random variable
$rr$, among the extrema of the vector $\mathbf{inp\\_x}$, corresponding to the
range where: for the $j^{th}$ dimension of all $\mathbf{x}_{1:i}$, the
$\alpha$ best function values were returned by the black-box function $f$. The
sought solution is the vector $\mathbf{opti\\_x}$, and the mimimum attained
value of the objective function $opti\\_f$.

##  5 Numerical Experiments

In this section, we present the results obtained by running the ITSO
algorithm, as well as Adaptive Differential Evolution (rand 1 bin),
Differential Evolution (rand 1 bin), and Differential Evolution (rand 2 bin)
with and without radius limited, Compass Coordinate Direct Search,
Probabilistic Descent Direct search, Random Search, Resampling Inheritance
Memetic Search, Resampling Memetic Search, Separable Natural Evolution
Strategies, Simultaneous Perturbation Stochastic Approximation, and
Exponential Natural Evolution Strategies from the Julia Package
BlackBoxOptim.jl [24], and Nelder-Mead, Particle Swarm, and Simulated
Annealing from Optim.jl [25]. To demonstrate the performance of each optimizer
in attaining the minimum, we firstly run the algorithm $r=10$ times, obtain
$f_{k}(\mathbf{x}_{i})$ for $k=\\{1,2,\dots,r\\}$ and all iterations
$i=\left\\{1,2,\ldots,{{f}_{e}}\right\\}$, and average the results

$\hat{f}(\mathbf{x}_{i})=\frac{\sum_{k=1}^{r}f_{k}(\mathbf{x}_{i})}{r}.$ (20)

Then, we normalize the vector of obtained function evaluations
$\mathbf{v}=\\{\hat{f}(\mathbf{x}_{1}),\hat{f}(\mathbf{x}_{1}),\dots,\hat{f}(\mathbf{x}_{f_{e}})\\}$
in the domain $\left[0,1\right]$ through

$\hat{f}_{n}(\mathbf{x}_{i})=\frac{\hat{f}(\mathbf{x}_{i})-\operatorname{min}\mathbf{v}}{\operatorname{max}\mathbf{v}-\operatorname{min}\mathbf{v}},$
(21)

and finally use the optimization history $h$

$h(\mathbf{x}_{i})=\left(\frac{\sum\limits_{l=1}^{m}{\hat{f}_{n}^{l}(\mathbf{x}_{i})}}{m}\right)^{\frac{1}{10}},$
(22)

where $m$ indicates the number of black-box functions used for the evaluation.
Equation 22 was selected as a performance metric, in order to obtain a clear
representation of the various optimizers utilized, as powers smaller than one
(in our case $\frac{1}{10}$), have the property to magnify the attained values
at the final steps of the optimization history. In Figure 3 the numerical
experiments for $m=13$ functions are presented. Each line corresponds to the
normalized average optimization history $h$ (Equation 22, for all functions,
which were repeated 10 times). We may see a clear prevalence of the proposed
framework, in terms of convergence performance, as expected by the theoretical
investigation. The numerical experiments for the comparison with other
optimizers, may be reproduced by running the file __run.jl.

(a) 10 variables and 5000 evaluations

(b) 20 variables and 5000 evaluations

(c) 20 variables and 10000 evaluations

(d) Utilized Optimizers

Figure 3: Average Optimization history $h$ (Equation 22, for all 13 Functions, repeated 10 times. Table 1: Average Minimum Values of $h$ Attained by each Optimizer | $n=10$ | $n=20$ | $n=20$
---|---|---|---
| $f_{e}=5000$ | $f_{e}=5000$ | $f_{e}=10000$
| Adapt. Diff. Evol. (rand 1 bin)
---
0.00991 | 0.02210 | 0.00852
| Adapt. Diff. Evol. (rand 1 bin, radious limited)
---
0.00734 | 0.01762 | 0.00482
| Diff. Evol. (rand 1 bin)
---
0.01348 | 0.02196 | 0.01161
| Diff. Evol. (rand 1 bin, radious limited)
---
0.00760 | 0.01377 | 0.00564
| Diff. Evol. (rand 2 bin)
---
0.02140 | 0.03712 | 0.02687
| Diff. Evol. (rand 2 bin, radious limited)
---
0.02072 | 0.03110 | 0.02008
| Compass Coordinate Direct Search
---
0.00084 | 0.00118 | 0.00045
| Probabilistic Descent Direct search
---
0.00984 | 0.01424 | 0.01239
Random Search | 0.05678 | 0.08307 | 0.07936
| Resampling Inheritance Memetic Search
---
0.00238 | 0.00513 | 0.00243
| Resampling Memetic Search
---
0.00129 | 0.00236 | 0.00223
| Separable Natural Evolution Strategies
---
0.02038 | 0.01165 | 0.01127
| Simultaneous Perturbation Stochastic Approximation
---
0.13967 | 0.13998 | 0.13227
| Exponential Natural Evolution Strategies
---
0.01972 | 0.01415 | 0.01286
Nelder-Mead | 0.01858 | 0.03307 | 0.03452
Particle Swarm | 0.00254 | 0.00166 | 0.00136
| Simulated Annealing
---
0.03825 | 0.03812 | 0.03721
Proposed-ITSO | 0.00001 | 0.00020 | 0.00017

##  6 Discussion and Conclusions

In this work a novel approach was presented for the well known problem of
finding the argument that minimizes a black-box, function or system. A vast
volume of approximation algorithms have been proposed, mainly heuristic, such
as genetic, evolutionary, particle swarm, as well as their variations.
However, they stem from nature-inspired procedures, and hence their converge
is investigated a-posteriori. Despite their efficiency, they are often
deprecated by researchers, due to the lack of rigorous mathematical
formulation, as well as complexity of implementation. To the contrary, the
proposed algorithm, initiates its formulation from well established
probabilistic definitions and theorems, and its implementation demands a few
lines of computer code. Furthermore, the convergence properties were found
stable, as a proof that the suggested framework attains the best possible
solution in the fewest possible iterations. The numerical examples validate
the theoretical results and may be reproduced by the provided computer code.
We consider the suggested method as a powerful framework which may easily be
adopted to the sought solution of any problem involving the minimization of a
black-box function.

##  Appendix A Programming Code

The corresponding computer code is available on GitHub
https://github.com/nbakas/ITSO.jl. The examples of Figure 3 may be reproduced
by running __run.jl. The sort version of the Algorithm 1 is available in Julia
[26] (file ITSO-short.jl), Octave [27] (ITSOshort.m), and Python [28] (ITSO-
short.py). The implementation of the framework is integrated in a few lines of
computer code, which can be easily adapted for case specific applications with
high efficiency.

##  Appendix B Black-Box Functions

The following functions were used for the numerical experiments. Equations 23,
24 (Elliptic, Cigar), were utilized from [29], Cigtab (Eq. 25), Griewank 26
from [30], Quartic (Eq. 27 from [31], Schwefel (Eq. 28), Rastrigin (Eq. 29),
Sphere (Eq. 30), and Ellipsoid (Eq. 31) from [32, 24], and Alpine (Eq. 32)
from [33]. Equations 33, 34, 35, were developed by the authors. The code
implementation for the selected equations appears in file functions_opti.jl in
the supplementary computer code.

The exact variation used in this work is as follows, where we have adopted the
notation presented in the Nomenclature section, where $i$ denotes the step of
the optimization history, and $j$ the dimension of the design variable
$x_{ij}$.

$\displaystyle
f_{elliptic}(\mathbf{x}_{i})=\sum_{j=1}^{n}{c_{j}(x_{ij}+\frac{3}{2})^{2}},\text{where}$
(23) $\displaystyle\mathbf{c}=10^{3}\\{0,\frac{1}{n-1},\dots,1\\}.$
$f_{cigar}(\mathbf{x}_{i})=x_{1}^{2}+\sum_{j=2}^{n}{|x_{ij}|}.$ (24)
$f_{cigtab}(\mathbf{x}_{i})=x_{1}^{2}+\sum_{j=2}^{n-1}{|x_{ij}|}+x_{n}^{2}.$
(25)
$f_{griewank}(\mathbf{x}_{i})=1+{\frac{1}{4000}}\sum_{{j=1}}^{n}x_{ij}^{2}-\prod_{{j=1}}^{n}\cos\left({\frac{x_{ij}}{{\sqrt{j}}}}\right).$
(26) $f_{quartic}(\mathbf{x}_{i})=\sum_{j=1}^{n}{j(x_{ij}-2)^{4}}.$ (27)
$\displaystyle
f_{schwefel}(\mathbf{x}_{i})=\sum_{j=1}^{n}{c_{j}^{2}},\text{where}$ (28)
$\displaystyle{c_{j}}=\sum_{k=1}^{j}{(x_{ik}-9)}.$ $\displaystyle
f_{rastrigin}(\mathbf{x}_{i})=10n+\sum_{j=1}^{n}{(x_{ij}+\frac{7}{10})^{2}}$
(29) $\displaystyle-10\sum_{j=1}^{n}{\cos(2\pi(x_{ij}+\frac{7}{10})^{2})}.$
$f_{sphere}(\mathbf{x}_{i})=\sum_{j=1}^{n}{(x_{ij}-\frac{13}{10})^{2}}.$ (30)
$\displaystyle
f_{elipsoid}(\mathbf{x}_{i})=\sum_{j=1}^{n}{(x_{ij}-\sqrt{2})^{2}}.$ (31)
$f_{alpine}(\mathbf{x}_{i})=\sum_{j=1}^{n}{|x_{ij}\sin{x_{ij}}+\frac{1}{10}x_{ij}|}.$
(32) $f_{x\\_j}(\mathbf{x}_{i})=\sum_{j=1}^{n}{(x_{ij}-j-\frac{21}{10})^{2}}.$
(33) $f_{x\\_5}(\mathbf{x}_{i})=\sum_{j=1}^{n}{(x_{ij}-5)^{2}}-5.$ (34)
$f_{sin\\_x}(\mathbf{x}_{i})=\sum_{j=1}^{n}{(\sin(x_{ij}+\frac{7}{10})+\frac{(x_{ij}+\frac{7}{10})^{2}}{100}}).$
(35)

## References

  * [1] M. Hutson, “AI researchers allege that machine learning is alchemy,” _Science_ , may 2018. [Online]. Available: http://www.sciencemag.org/news/2018/05/ai-researchers-allege-machine-learning-alchemy
  * [2] D. Sculley, J. Snoek, A. Rahimi, and A. Wiltschko, “Winner’s Curse? On Pace, Progress, and Empirical Rigor,” _ICLR Workshop track_ , 2018.
  * [3] S. Sra, S. Nowozin, and S. J. Wright, _Optimization for machine learning_. Mit Press, 2012.
  * [4] L. Bottou, F. E. Curtis, and J. Nocedal, “Optimization methods for large-scale machine learning,” _Siam Review_ , vol. 60, no. 2, pp. 223–311, 2018.
  * [5] C. Rudin, “Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead,” _Nature Machine Intelligence_ , vol. 1, no. 5, pp. 206–215, 2019. [Online]. Available: https://doi.org/10.1038/s42256-019-0048-x
  * [6] J. Wu, M. Poloczek, A. G. Wilson, and P. Frazier, “Bayesian optimization with gradients,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 5267–5278.
  * [7] B. Doerr, “Probabilistic tools for the analysis of randomized optimization heuristics,” in _Theory of Evolutionary Computation_. Springer, 2020, pp. 1–87.
  * [8] C. Doerr, “Complexity theory for discrete black-box optimization heuristics,” in _Theory of Evolutionary Computation_. Springer, 2020, pp. 133–212.
  * [9] M. Papadrakakis, N. D. Lagaros, and V. Plevris, “Design optimization of steel structures considering uncertainties,” _Engineering Structures_ , vol. 27, no. 9, pp. 1408–1418, 2005.
  * [10] ——, “Optimum design of space frames under seismic loading,” _International Journal of Structural Stability and Dynamics_ , vol. 1, no. 01, pp. 105–123, 2001.
  * [11] N. Moayyeri, S. Gharehbaghi, and V. Plevris, “Cost-based optimum design of reinforced concrete retaining walls considering different methods of bearing capacity computation,” _Mathematics_ , vol. 7, no. 12, p. 1232, 2019.
  * [12] V. Plevris and M. Papadrakakis, “A hybrid particle swarm—gradient algorithm for global structural optimization,” _Computer-Aided Civil and Infrastructure Engineering_ , vol. 26, no. 1, pp. 48–68, 2011.
  * [13] N. D. Lagaros, N. Bakas, and M. Papadrakakis, “Optimum Design Approaches for Improving the Seismic Performance of 3D RC Buildings,” _Journal of Earthquake Engineering_ , vol. 13, no. 3, pp. 345–363, mar 2009. [Online]. Available: http://www.tandfonline.com/doi/full/10.1080/13632460802598594
  * [14] N. D. Lagaros, M. Papadrakakis, and N. P. Bakas, “Automatic minimization of the rigidity eccentricity of 3D reinforced concrete buildings,” _Journal of Earthquake Engineering_ , vol. 10, no. 4, pp. 533–564, jul 2006\. [Online]. Available: http://www.tandfonline.com/doi/full/10.1080/13632460609350609
  * [15] A. D. Bull, “Convergence rates of efficient global optimization algorithms,” _Journal of Machine Learning Research_ , vol. 12, no. Oct, pp. 2879–2904, 2011.
  * [16] G. Rudolph, “Convergence analysis of canonical genetic algorithms,” _IEEE transactions on neural networks_ , vol. 5, no. 1, pp. 96–101, 1994\.
  * [17] M. Clerc and J. Kennedy, “The particle swarm-explosion, stability, and convergence in a multidimensional complex space,” _IEEE transactions on Evolutionary Computation_ , vol. 6, no. 1, pp. 58–73, 2002.
  * [18] K. R. Opara and J. Arabas, “Differential evolution: A survey of theoretical analyses,” _Swarm and evolutionary computation_ , vol. 44, pp. 546–558, 2019\.
  * [19] N. Razmjooy, V. V. Estrela, H. J. Loschi, and W. Fanfan, “A comprehensive survey of new meta-heuristic algorithms,” _Recent Advances in Hybrid Metaheuristics for Data Clustering, Wiley Publishing_ , 2019.
  * [20] M. A. Muñoz and K. A. Smith-Miles, “Performance analysis of continuous black-box optimization algorithms via footprints in instance space,” _Evolutionary computation_ , vol. 25, no. 4, pp. 529–554, 2017.
  * [21] C. Audet and W. Hare, _Derivative-free and blackbox optimization_. Springer, 2017.
  * [22] C. Audet and M. Kokkolaras, “Blackbox and derivative-free optimization: theory, algorithms and applications,” 2016.
  * [23] F. D. Parker, “Integrals of Inverse Functions,” _The American Mathematical Monthly_ , vol. 62, no. 6, p. 439, jun 1955. [Online]. Available: http://www.jstor.org/stable/2307006?origin=crossref
  * [24] R. Feldt, “Blackboxoptim.jl,” 2013-2018. [Online]. Available: https://github.com/robertfeldt/BlackBoxOptim.jl
  * [25] J. Myles White, T. Holy, O. Contributors (2012), P. Kofod Mogensen, J. Myles White, T. Holy, P. Contributors (2016), Other. Kofod Mogensen, A. Nilsen Riseth, J. Myles White, T. Holy, and O. Contributors (2017), “Optim.jl,” 2012,2016,2017. [Online]. Available: https://github.com/JuliaNLSolvers/Optim.jl
  * [26] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” _SIAM review_ , vol. 59, no. 1, pp. 65–98, 2017.
  * [27] Contributors, “Gnu octave,” http://hg.savannah.gnu.org/hgweb/octave/file/tip/doc/interpreter/contributors.in, 28 Feb 2020.
  * [28] ——, “Python 3.8.2,” https://www.python.org/, 2020.
  * [29] J. Liang, B. Qu, P. Suganthan, and A. G. Hernández-Díaz, “Problem definitions and evaluation criteria for the cec 2013 special session on real-parameter optimization,” _Computational Intelligence Laboratory, Zhengzhou University, Zhengzhou, China and Nanyang Technological University, Singapore, Technical Report_ , vol. 201212, no. 34, pp. 281–295, 2013.
  * [30] C.-K. Au and H.-F. Leung, “Eigenspace sampling in the mirrored variant of (1, $\lambda$)-cma-es,” in _2012 IEEE Congress on Evolutionary Computation_. IEEE, 2012, pp. 1–8.
  * [31] M. Jamil and X.-S. Yang, “A literature survey of benchmark functions for global optimization problems,” _arXiv preprint arXiv:1308.4008_ , 2013.
  * [32] S. Finck, N. Hansen, R. Ros, and A. Auger, “Real-parameter black-box optimization benchmarking 2009: Presentation of the noiseless functions,” Citeseer, Tech. Rep., 2010.
  * [33] K. Hussain, M. N. M. Salleh, S. Cheng, and R. Naseem, “Common benchmark functions for metaheuristic evaluation: A review,” _JOIV: International Journal on Informatics Visualization_ , vol. 1, no. 4-2, pp. 218–223, 2017.


# Simulating the dynamics of braiding of Majorana zero modes using an IBM
quantum computer

John P. T. Stenger  Department of Physics and Astronomy, University of
Pittsburgh, Pittsburgh, PA 15260, USA  Nicholas T. Bronn  IBM Quantum, IBM
T.J. Watson Research Center, Yorktown Heights, New York 10598, USA  Daniel J.
Egger  IBM Quantum, IBM Research – Zurich, Säumerstrasse 4, 8803 Rüschlikon,
Switzerland  David Pekker  Department of Physics and Astronomy, University of
Pittsburgh, Pittsburgh, PA 15260, USA

(August 27, 2024)

###### Abstract

We simulate the dynamics of braiding Majorana zero modes on an IBM Quantum
computer. We find the native quantum gates introduce too much noise to observe
braiding. Instead, we use Qiskit Pulse to develop scaled two-qubit quantum
gates that better match the unitary time evolution operator and enable us to
observe braiding. This work demonstrates that quantum computers can be used
for simulation, and highlights the use of pulse-level control for programming
quantum computers and constitutes the first experimental evidence of braiding
via dynamical Hamiltonian evolution.

##  I Introduction

Quantum computers may significantly outperform classical ones in the area of
simulation of quantum systems Feynman (1982); Lloyd (1996) and other
specialized algorithms Boixo  _et al._ (2018); Shor (1997). However, we are
currently in the era of noisy quantum computing Preskill (2018) where only a
small number of qubits, with relatively short coherences times, can be
entangled. Fortunately, recent results demonstrating hardware specific
optimization Shi  _et al._ (2019); Glaser  _et al._ (2015); Leung  _et al._
(2017) and quantum advantage with short-depth noisy circuits Bravyi  _et al._
(2020a) offer hope for the usefulness of near-term systems. Here, we simulate
a quantum topological condensed matter system on an IBM Quantum processor
within the qubits’ coherence times using pulse-level instructions provided by
Qiskit Pulse McKay  _et al._ (2018); Alexander  _et al._ (2020). Ideally, the
simulation would occur by a continuous time-evolution of the qubits under the
appropriate spin Hamiltonian obtained from a transformation of the fermion
Hamiltonian. Practically, this “analog” simulation must be decomposed and
mapped onto the calibrated native basis gates of a quantum computer, making it
“digital”. This digital implementation on noisy quantum hardware limits the
flexibility needed to avoid the accumulation of unnecessary errors. Here, we
demonstrate a “semi-analog” approach to noise reduction: a pulse-scaling
technique that, without additional calibration Alexander  _et al._ (2020);
Garion  _et al._ (2020), gets us closer to the ideal analog simulation.

Topologically protected quantum computation works by moving nonabelian anyons,
such as Majorana zero modes (MZMs), around each other in two dimensions to
form three dimensional braids in space-time Kitaev (2001); Nayak  _et al._
(2008). This approach is advantageous as it offers protection from local
perturbations. The quest for topological quantum computation has focused on
hybrid superconducting-semiconducting Alicea (2012); Beenakker (2013); Chen
_et al._ (2017); Deng  _et al._ (2016); Kitaev (2001); Lutchyn  _et al._
(2010); Mourik  _et al._ (2012); Oreg  _et al._ (2010); Sau  _et al._ (2010,
2012); Stenger  _et al._ (2019); Ben-Shach  _et al._ (2015) and fractional
quantum hall devices Camino  _et al._ (2007); Goldman and Su (1995); de C.
Chamon  _et al._ (1997). While trivial-topological phase transition Xiao  _et
al._ (2020) and measurement-based braiding Wootton (2017) have been observed
on a (non-topological) quantum computer, thus far there has been no definitive
experimental evidence of braiding due to dynamical state evolution Bartolomei
_et al._ (2020); Camino  _et al._ (2007); Ofek  _et al._ (2010); McClure  _et
al._ (2012); Willett  _et al._ (2013); Nakamura  _et al._ (2019); Willett  _et
al._ (2019).

In this work, we simulate the key part of a topological quantum computer: the
dynamics of braiding of a pair of MZMs on a tri-junction. Specifically, we use
the results of Ref. Backens  _et al._ (2017) to map a minimum model of a
topological tri-junction Alicea  _et al._ (2011); Halperin  _et al._ (2012);
Hyart  _et al._ (2013); van Heck  _et al._ (2012); Hassler  _et al._ (2011);
Stenger  _et al._ (2019) onto a three-qubit Hamiltonian. Braiding is
implemented by parametrically adjusting the Hamiltonian parameters; the time
evolution is implemented using the Suzuki-Trotter decomposition, with each
time step implemented by one- and two-qubit gates. We significantly boost the
fidelity of our quantum time-evolution code by using pulse-level control to
scale cross resonance (CR) gates Chow  _et al._ (2011) derived from those pre-
calibrated on the backend, thereby enabling coupling of qubit pairs with
shorter CR gate times. Specifically, we observe that using native CNOT gates,
we can move a MZM from one arm of the tri-junction device to another arm, thus
performing 1/6 of a full braid. However, using the scaled gates, we are able
observe a complete braid. We remark that this result can be interpreted as an
experimental observation of braiding in a quantum system.

Figure 1: Topological tri-junctions. (a) Schematic depiction of a tri-junction
device as proposed in Refs. Alicea  _et al._ (2011); Stenger  _et al._ (2019);
Halperin  _et al._ (2012); Hyart  _et al._ (2013); van Heck  _et al._ (2012);
Hassler  _et al._ (2011). The arms of the device are represented by blue bars,
Majorana modes by red circles and the tunable couplings between pairs of inner
Majorana modes by dashed lines. (b) Minimum model composed of three qubits,
each qubit contributes a pair of Pauli operators to the description of the
Majoranas. The representation is written in the blue ovals. Figure 2: Gate
sets that implement a single Suzuki-Trotter step $t_{i}$ of the time
evolution. (a) Implementation with 2-qubit basis gates; each set of gates is
labeled by the corresponding Hamiltonian terms above. (b) Implementation with
analog scaled 2-qubit ECR gates $A(\theta)$, where
$ZX(\theta)=R_{x}(\pi/2)A(\theta)R_{x}(\pi/2)$, as shown in Fig. 4c.

A minimal model of a topological tri-junction – Our model is constructed using
the six Majorana operators depicted in Fig. 1a and is described by the
Hamiltonian 111See Section A of the supplement for a model that describes both
high energy modes and MZMs and the corresponding quantum code.

$H=i\alpha\sum_{a=0}^{2}\gamma_{a}^{x}\gamma_{a}^{y}+i\frac{1}{2}\sum_{a\neq
b}J_{ab}\gamma_{a}^{x}\gamma_{b}^{x}.$ (1)

Here, $\gamma_{a}^{x}$, $\gamma_{a}^{y}$ are Majorana operators on the inner
and outer ends of arm $a$ respectively, $\alpha$ describes the coupling
between Majorana modes on the same arm, and $J_{ab}=-J_{ba}$ couples Majorana
modes at the tri-junction. Initially, we set $J_{01}=J_{\text{max}}$ and
$J_{12}=J_{20}=0$. This setting results in $\gamma_{0}^{x}$ and
$\gamma_{1}^{x}$ fusing into one complex fermion and $\gamma_{2}^{x}$ and
$\gamma_{2}^{y}$ fusing into another complex fermion. The low-energy sector is
spanned by the operators $\gamma_{0}^{y}$ and $\gamma_{1}^{y}$, which are the
two MZMs that we would like to braid. Braiding is preformed by moving the two
MZMs around the arms of the tri-junction, which is accomplished by modulating
each $J_{ab}(t)$ as a function of time from $0$ to some maximum value which we
call $J_{\text{max}}$ and then back to $0$ following the protocol enumerated
in Table 1. This is the same braiding procedure that would be used on a real
tri-junction Stenger  _et al._ (2019). This braiding procedure relies on a
separation of timescales: braiding should be slow compared to the timescale of
the fused MZMs, but fast compared to the splitting of the two “free” MZMs
222See Section B of the supplement for a discussion of topological protection.

step | time | couplings
---|---|---
1 | $t:0\rightarrow\tau$ |  $J_{01}:J_{\text{max}}\rightarrow 0$, $J_{12}:0\rightarrow J_{\text{max}}$
2 | $t:\tau\rightarrow 2\tau$ |  $J_{12}:J_{\text{max}}\rightarrow 0$, $J_{20}:0\rightarrow J_{\text{max}}$
3 | $t:2\tau\rightarrow 3\tau$ |  $J_{20}:J_{\text{max}}\rightarrow 0$, $J_{01}:0\rightarrow J_{\text{max}}$
4-6 | $t:3\tau\rightarrow 6\tau$ | repeat steps 1-3
Table 1: Procedure for braiding Majorana zero modes depicted in Fig. 1 by
modulating the couplings $J_{01}$, $J_{12}$, and $J_{20}$. Each step moves a
Majorana from one arm to another. In our calculations, we set
$J_{\text{max}}=1$, $\alpha=3\leavevmode\nobreak\ J_{\text{max}}$, and
$\tau=3.3\leavevmode\nobreak\ 1/J_{\text{max}}$.

To simulate braiding on a quantum computer we map the fermionic Hamiltonian
Eq. (1) onto a Hamiltonian

$H(t)=\alpha\sum_{a=0}^{2}\sigma_{a}^{z}+J_{01}(t)\sigma_{0}^{y}\sigma_{1}^{x}+J_{12}(t)\sigma_{1}^{y}\sigma_{2}^{x}+J_{20}(t)\sigma_{0}^{y}\sigma_{1}^{z}\sigma_{2}^{x}$
(2)

that acts on qubits. Here, $\sigma_{a}^{\nu}$ is the Pauli $\nu$-matrix acting
on qubit $a$. The Hilbert space of the minimal model of Eq. (1) is equivalent
to the Hilbert space of three qubits, depicted in Fig. 1b. The Majorana
operators are related to the qubit operators via the Jordan-Wigner
transformation, see Fig. 1(b).

To simulate the dynamics of Eq. (2) we must find the quantum gates that
approximate the time evolution operator
$U=\mathcal{T}\prod_{i}\exp{\left(-iH(t_{i})dt\right)}$, where $\mathcal{T}$
is the time ordering operator and $dt$ is a small time step. Throughout, we
rely on the second order Suzuki-Trotter approximation to decompose the time
evolution into manageable pieces. We choose the time step $dt$ to maximize the
fidelity of the braiding procedure which is determined by a competition
between (1) the Suzuki-Trotter error that is minimized by making a large
number of small time-steps $dt$ and (2) the error from the imperfect quantum
gates which is minimized by reducing the gate count by making $dt$ large.

Quantum simulation with basis gates – We split the time evolution of a single
Suzuki-Trotter step into four parts. Each part is decomposed into single-qubit
rotations labeled $R_{\nu,i}(\theta)$ for a $\theta$-rotation around the
$\nu$-axis of qubit $i$, and CNOT gates labeled $C_{x,ij}$, as depicted in
Fig. 2a. We initialize the qubits to

$\psi_{\pm}(0)=R_{y,2}(\pi)R_{y,0}(\mp\pi/2)C_{x,01}R_{x,0}(-\pi/2)\ket{000},$
(3)

where the states $\psi_{\pm}(0)=\ket{e}\pm\ket{g}$ with $\ket{g}$ and
$\ket{e}$ being the ground and first excited state of the tri-junction
corresponding to the even $\ket{e}$ and odd $\ket{g}$ parity states of the two
MZMs. Because there is a finite coupling between the MZMs, given by $\alpha$
in Eq. (2), $\ket{g}$ and $\ket{e}$ are at slightly different energies. Any
state outside of the low energy subspace defined by $\ket{g}$ and $\ket{e}$
can be understood as a linear combination of higher energy states. After a
single braid we expect:
$\psi_{\pm}\underset{\text{braid}}{\rightarrow}\psi_{\mp}$. To check if
braiding was successful, we invert the initialization gates so that
$\ket{000}$ corresponds to successful braiding.

Figure 3: Panel (a): Comparing the probability distributions obtained using
noiseless classical simulations and quantum computer calculations (see legend)
obtained after moving a MZM from one arm of a tri-junction to another arm (as
depicted in panel (b)). The $\psi_{+}(\tau)$ and $\psi_{-}(\tau)$ states form
the low-energy subspace, with $\psi_{+}(\tau)$ corresponding to the target
state in the absence of time-step and gate errors; in constructing the plot we
normalized probabilities in this subspace so that
$P_{\psi_{+}(\tau)}+P_{\psi_{-}(\tau)}=1$. The probability to end up in all
other states is grouped together and labeled ‘High Energy States’.

A full braid with only three time steps per MZM swap (i.e. $\tau$ is broken
into three time steps) requires 96 CNOT gates and cannot be simulated on
ibmq_athens as the CNOT gates (the dominate source of error) have an error of
$1.3\%$ and $0.9\%$ between qubits $(0,1)$ and $(1,2)$, respectively. The CNOT
gate for qubits $(0,2)$ is built out of the other two. However, moving a
single MZM from one arm of the tri-junction to another arm (thus performing
1/6th of a complete braid), requires only $3\times 4=12$ two qubit gates. This
corresponds to step 1 in table 1 with $J_{20}=0$. The set of gates depicted in
Fig. 2 was translated into the basis gates using Qiskit’s transpile function
at optimization level 3 and 1024 shots were measured. Figure 3 shows the
resulting probability distribution 333See section C of the supplement for
additional data and details of the measurement procedure.. The probabilities
labeled ‘noiseless’ were calculated using the QASM simulator in Qiskit to
perform fully coherent quantum evolution and represent the size of time-step
errors. We observe a clear preference for the $\psi_{+}(\tau)$ state when
running on ibmq_athens, indicating that we have successfully moved a MZM from
one arm of the tri-junction to another. However, when performing the full
braid, there is no longer a clear preference as we measure
$P_{-}(\tau)/P_{+}(\tau)=1.04\pm 0.12$ where $P_{\pm}$ is the probability of
system ending in the $\ket{\pm}state$

Designing scaled quantum gates –  We decrease two-qubit gate errors by
designing scaled gates which implement smaller rotations in the two-qubit
Hilbert subspace than the basis CNOT gates, which apply a full $\pi/2$
rotation in the Hilbert space. We focus on the operator
$U_{yx}(\theta)=e^{-i\theta\sigma^{y}_{1}\sigma^{x}_{2}/2}$, which arises from
the Jordan-Wigner transformation, see Fig. 2a. $U_{yx}(\theta)$ can be
implemented with 2 CNOT gates as in Fig. 4(a). The basis CNOT gates are
created from $U_{zx}(\pi/2)$ implemented by CR. These CNOTs are echoed CR
pulses Sheldon  _et al._ (2016) on the control qubit with rotary echoes
Sundaresan  _et al._ (2020) on the target qubit, combined with single-qubit
gates before and after the echoed CR gate that generate the correct direction
of CNOT. Since errors mainly arise during CR pulses, we minimize their
duration in the braiding algorithm. This is achieved with pulse-level control
enabled by Qiskit Pulse McKay  _et al._ (2018). We implement $U_{yx}(\theta)$
using $U_{zx}(\theta)$ operations (see Fig. 4b) derived from the highly-
calibrated CNOT pulse schedules Gokhale  _et al._ (2020); Alexander  _et al._
(2020). Modifying the CNOT pulse schedules avoids additional calibrations
which is paramount when running jobs through a queue on cloud-based quantum
computers 444see section D of the Supplement for details on the implementation
of the scaled gates. The rotation angle $\theta$ depends on the area under the
pulses, and is often considerably less than $2\times\pi/2$ for the Majorana
braiding simulation schedule (Fig. 2). This allows us to build
$U_{zx}(\theta)$ gates with considerably shorter CR pulses, and hence
introduce less error per Suzuki-Trotter step.

Figure 4: Quantum process tomography of $U_{yx}(\theta)$ as a function of
$\theta$ on ibmq_athens. Both the double-CNOT circuit (a) and the parametric
$ZX(\theta)$-based circuit (b) create $U_{yx}(\theta)$. The circuit in (b) at
$\theta=\pi/2$ is implemented by the pulse schedule shown in (c). Here, each
sample lasts $0.222\leavevmode\nobreak\ {\rm ns}$. (d) Fidelity $F_{(a)}$
(blue up-triangles) and $F_{(b)}$ (orange down triangles) of quantum circuits
(a) and (b), respectively, obtained with QPT. (e) Relative error reduction of
the pulse-efficient circuit in (b) over the circuit in (a), i.e.
$1-E_{(b)}/E_{(a)}$ where $E_{(x)}=1-F_{(x)}$.

To benchmark the $U_{yx}$ analog quantum gate, we compare it to the double-
CNOT basis gate implementation (see Fig. 4a and b). We measure the gate
fidelity with quantum process tomgraphy (QPT) for 15 values of $\theta$
linearly spaced between 0 and $\pi$. Each measurement is done with 2048 shots
and repeated four times to gather statistics. We mitigate readout errors by
preparing each of the four basis states and measuring the outcome which we use
to correct the QPT data Bravyi  _et al._ (2020b); Barron and Wood (2020).
Fidelity measurements of the benchmark circuit are interleaved with those of
the scaled CR circuit to mitigate biases in our comparison due to drifts. We
observe that the scaled CR pulses systematically have a higher fidelity than
the double-CNOT implementation at all measured values of $\theta$, see Fig.
4d. The double-CNOT benchmark should have a constant fidelity as the rotation
angle $\theta$ is implemented by a virtual $Z$ gate McKay  _et al._ (2017)
once the schedule is transpiled to ibmq_athens. We therefore attribute the
fidelity fluctuations in Fig. 4d to drifts as the data were acquired over a
three-day period. Such fluctuations are also observed in the fidelity of the
scaled CR pulses. We observe that the analog circuit is strongly advantageous
to the basis gate circuit at small $\theta$, see Fig. 4e, which is the
relevant case for quantum simulation. We attribute the decrease in fidelity of
the analog circuit at larger $\theta$ to an increase in errors caused by the
longer duration CR pulses.

Quantum simulation with scaled gates –  The Hamiltonian (2) contains a three-
qubit coupling term which induces a unitary time evolution that cannot be
efficiently encoded using our analog two-qubit quantum gates. We perform one
more basis transformation, which results in the three-qubit gate being
replaced by a pair of two-qubit gates 555see section E of the supplement for
details of the basis transformation. The resulting quantum circuit, after hand
optimization to combine single-qubit rotations, is depicted in Fig. 2b.

Figure 5: Bias towards braiding as a function of delay time added to the ZX
rotation gate. We plot $(P_{-}-P_{+})/(P_{-}+P_{+})$ so that 1 (-1)
corresponds to perfect braiding for the plus (minus) state and 0 corresponds
to an equal distribution of shots in the braiding target state and in the
trivial (non-braiding) state.

We now return to Fig. 3 to compare the fidelity of the simulations implemented
with basis and scaled gates. We note that quantum simulations performed with
basis and scaled gates use a different wave function basis and therefore have
slightly different time-step errors as seen from the high energy states in Fig
3. We observe a larger leakage out of the low energy subspace for the scaled
computation than for the basis gates, but the fidelity in the low energy basis
is better. Crucially, the scaled-gate computation on the quantum computer is
much closer to its noiseless ideal than is the basis-gate calculation on the
quantum computer to its own noiseless ideal, highlighting the improvement in
the fidelity achieved with analog quantum gates.

In Fig. 5 we plot the results of running the braiding pulse schedule on
ibmq_athens with a variable amount of incoherent error added in the form of an
additional delay time appended to all 2-qubit gates 666See Section F of the
supplement for additional data. The probabilities $P_{+}$ and $P_{-}$ are
found by counting the number of shots that end in the $\psi_{+}$ and
$\psi_{-}$ states, respectively, after the braiding procedure. Each data point
is the average of four trials with 8192 shots in each trial after applying the
same readout error mitigation as with the QPT Bravyi  _et al._ (2020b); Barron
and Wood (2020). We use the ratio $(P_{-}-P_{+})/(P_{+}+P_{-})$ to measure
success in braiding, this ratio should be positive if the final state is
closer to $\psi_{-}(0)$ and negative if it closer to $\psi_{+}(0)$. In the
absence of additional noise, we observe a strong preference for successful
braiding. As we introduce additional noise this preference slowly diminishes,
disappearing completely when the delay exceeds 150 ns or about $25\%$ of the
duration of a CNOT gate. The observation that additional noise washes out the
braiding signal (a) supports our interpretation that we are indeed observing a
quantum coherent process of braiding and (b) explains why we were unable to
observe braiding with CNOT gates which extend the pulse schedule significantly
past the 150 ns per two-qubit gate at which the braiding signal disappears
777See Section G of the supplement for a noise model.

In summary – We have demonstrated that pulse-level control of quantum
computers enables us to simulate the braiding of Majorana zero modes, thus
expanding the “digital” capabilities offered by the native basis gates on IBM
Quantum backends. These two-qubit operations were derived from the highly-
calibrated basis CNOT of the backend and required no further calibration. Our
demonstration shows that we have reached the point at which quantum computers
can perform interesting quantum simulations, but to achieve a sufficient
quality it is crucial to understand the performance of the hardware and
compose software that respects its limitations. Looking towards the future,
the ease of programming obtained from digital abstractions is outweighed by
the increase in performance obtained by programming the “analog” pulse
schedules in a noise-aware method. We therefore argue that the path forward
for quantum simulation in the noisy quantum era is the use of abstraction-free
programming of quantum computers. This allows for continuous time evolution on
part of of the Hamiltonian and could be a path towards fully continuous
Hamiltonian evolution.

##  II Acknowledgements

The authors acknowledge use of the IBM Quantum Experience devices for this
work. The authors also thank S.M. Frolov, N. Earnest-Noble, D.T. McClure, N.
Kanazawa and E. Pritchett for insightful discussions. J.S. and D.P.
acknowledge support from NSF PIRE-1743717.

## References

  * Feynman (1982) R. P. Feynman, Int. J. Theor Phys. 21, 467 (1982).
  * Lloyd (1996) S. Lloyd, Science (New York, N.Y.) 273, 1073 (1996).
  * Boixo  _et al._ (2018) S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595–600 (2018).
  * Shor (1997) P. W. Shor, SIAM Journal on Computing 26, 1484–1509 (1997).
  * Preskill (2018) J. Preskill, Quantum 2, 79 (2018).
  * Shi  _et al._ (2019) Y. Shi, N. Leung, P. Gokhale, Z. Rossi, D. I. Schuster, H. Hoffmann, and F. T. Chong, Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems (2019), 10.1145/3297858.3304018.
  * Glaser  _et al._ (2015) S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, and et al., The European Physical Journal D 69 (2015).
  * Leung  _et al._ (2017) N. Leung, M. Abdelhafez, J. Koch, and D. Schuster, Physical Review A 95 (2017).
  * Bravyi  _et al._ (2020a) S. Bravyi, D. Gosset, R. König, and M. Tomamichel, Nature Physics 16, 1040 (2020a).
  * McKay  _et al._ (2018) D. C. McKay, T. Alexander, L. Bello, M. J. Biercuk, L. Bishop, J. Chen, J. M. Chow, A. D. Córcoles, D. Egger, S. Filipp, J. Gomez, M. Hush, A. Javadi-Abhari, D. Moreda, P. Nation, B. Paulovicks, E. Winston, C. J. Wood, J. Wootton, and J. M. Gambetta, “Qiskit backend specifications for openqasm and openpulse experiments,”  (2018), arXiv:1809.03452 .
  * Alexander  _et al._ (2020) T. Alexander, N. Kanazawa, D. J. Egger, L. Capelluto, C. J. Wood, A. Javadi-Abhari, and D. C McKay, Quantum Science and Technology 5, 044006 (2020).
  * Garion  _et al._ (2020) S. Garion, N. Kanazawa, H. Landa, D. C. McKay, S. Sheldon, A. W. Cross, and C. J. Wood, “Experimental implementation of non-Clifford interleaved randomized benchmarking with a controlled-S gate,”  (2020), arXiv:2007.08532 .
  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Alicea (2012) J. Alicea, Rep. Prog. Phys. 75, 076501 (2012).
  * Beenakker (2013) C. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * Chen  _et al._ (2017) J. Chen, P. Yu, J. Stenger, M. Hocevar, D. Car, S. R. Plissard, E. P. A. M. Bakkers, T. D. Stanescu, and S. M. Frolov, Sci. Adv. 3, e1701476 (2017).
  * Deng  _et al._ (2016) M. T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus, Science 354, 1557 (2016).
  * Lutchyn  _et al._ (2010) R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Mourik  _et al._ (2012) V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
  * Oreg  _et al._ (2010) Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Sau  _et al._ (2010) J. D. Sau, S. Tewari, R. M. Lutchyn, T. D. Stanescu, and S. Das Sarma, Phys. Rev. B 82, 214509 (2010).
  * Sau  _et al._ (2012) J. D. Sau, S. Tewari, and S. Das Sarma, Phys. Rev. B 85, 064512 (2012).
  * Stenger  _et al._ (2019) J. P. T. Stenger, M. Hatridge, S. M. Frolov, and D. Pekker, Phys. Rev. B 99, 035307 (2019).
  * Ben-Shach  _et al._ (2015) G. Ben-Shach, A. Haim, I. Appelbaum, Y. Oreg, A. Yacoby, and B. I. Halperin, Physical Review B 91 (2015), 10.1103/physrevb.91.045403.
  * Camino  _et al._ (2007) F. E. Camino, W. Zhou, and V. J. Goldman, Phys. Rev. Lett. 98, 076805 (2007).
  * Goldman and Su (1995) V. J. Goldman and B. Su, Science 267, 1010 (1995).
  * de C. Chamon  _et al._ (1997) C. de C. Chamon, D. E. Freed, S. A. Kivelson, S. L. Sondhi, and X. G. Wen, Phys. Rev. B 55, 2331 (1997).
  * Xiao  _et al._ (2020) X. Xiao, J. K. Freericks, and A. F. Kemper, “Topological quantum computing on a conventional quantum computer,” (2020), arXiv:2006.05524 [quant-ph] .
  * Wootton (2017) J. R. Wootton, Quantum Science and Technology 2, 015006 (2017).
  * Bartolomei  _et al._ (2020) H. Bartolomei, M. Kumar, R. Bisognin, A. Marguerite, J.-M. Berroir, E. Bocquillon, B. Plaçais, A. Cavanna, Q. Dong, U. Gennser, and et al., Science 368, 173–177 (2020).
  * Ofek  _et al._ (2010) N. Ofek, A. Bid, M. Heiblum, A. Stern, V. Umansky, and D. Mahalu, Proceedings of the National Academy of Sciences 107, 5276–5281 (2010).
  * McClure  _et al._ (2012) D. T. McClure, W. Chang, C. M. Marcus, L. N. Pfeiffer, and K. W. West, Physical Review Letters 108 (2012).
  * Willett  _et al._ (2013) R. L. Willett, C. Nayak, K. Shtengel, L. N. Pfeiffer, and K. W. West, Physical Review Letters 111 (2013).
  * Nakamura  _et al._ (2019) J. Nakamura, S. Fallahi, H. Sahasrabudhe, R. Rahman, S. Liang, G. C. Gardner, and M. J. Manfra, Nature Physics 15, 563–569 (2019).
  * Willett  _et al._ (2019) R. L. Willett, K. Shtengel, C. Nayak, L. N. Pfeiffer, Y. J. Chung, M. L. Peabody, K. W. Baldwin, and K. W. West, “Interference measurements of non-abelian e/4 & abelian e/2 quasiparticle braiding,”  (2019), arXiv:1905.10248 .
  * Backens  _et al._ (2017) S. Backens, A. Shnirman, Y. Makhlin, Y. Gefen, J. E. Mooij, and G. Schön, Phys. Rev. B 96, 195402 (2017).
  * Alicea  _et al._ (2011) J. Alicea, Y. Oreg, G. Refael, F. von Oppen, and M. P. A. Fisher, Nat. Phys. 7, 412 (2011).
  * Halperin  _et al._ (2012) B. I. Halperin, Y. Oreg, A. Stern, G. Refael, J. Alicea, and F. von Oppen, Phys. Rev. B 85, 144501 (2012).
  * Hyart  _et al._ (2013) T. Hyart, B. van Heck, I. C. Fulga, M. Burrello, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev. B 88, 035121 (2013).
  * van Heck  _et al._ (2012) B. van Heck, A. R. Akhmerov, F. Hassler, and C. W. J. Beenakker, New J. Phys. 14, 035019 (2012).
  * Hassler  _et al._ (2011) F. Hassler, A. R. Akhmerov, and C. W. J. Beenakker, New J. Phys. 13, 095004 (2011).
  * Chow  _et al._ (2011) J. M. Chow, A. D. Córcoles, J. M. Gambetta, C. Rigetti, B. R. Johnson, J. A. Smolin, J. R. Rozen, G. A. Keefe, M. B. Rothwell, M. B. Ketchen, and M. Steffen, Phys. Rev. Lett. 107, 080502 (2011).
  * Note (1) See Section A of the supplement for a model that describes both high energy modes and MZMs and the corresponding quantum code.
  * Note (2) See Section B of the supplement for a discussion of topological protection.
  * Note (3) See section C of the supplement for additional data and details of the measurement procedure.
  * Sheldon  _et al._ (2016) S. Sheldon, E. Magesan, J. M. Chow, and J. M. Gambetta, Physical Review A 93 (2016).
  * Sundaresan  _et al._ (2020) N. Sundaresan, I. Lauer, E. Pritchett, E. Magesan, P. Jurcevic, and J. M. Gambetta, “Reducing unitary and spectator errors in cross resonance with optimized rotary echoes,”  (2020), arXiv:2007.02925 .
  * Gokhale  _et al._ (2020) P. Gokhale, A. Javadi-Abhari, N. Earnest, Y. Shi, and F. T. Chong, “Optimized quantum compilation for near-term algorithms with openpulse,”  (2020), arXiv:2004.11205 .
  * Note (4) See section D of the Supplement for details on the implementation of the scaled gates.
  * Bravyi  _et al._ (2020b) S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, “Mitigating measurement errors in multi-qubit experiments,”  (2020b), arXiv:2006.14044 .
  * Barron and Wood (2020) G. S. Barron and C. J. Wood, “Measurement error mitigation for variational quantum algorithms,”  (2020), arXiv:2010.08520 .
  * McKay  _et al._ (2017) D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Physical Review A 96 (2017).
  * Note (5) See section E of the supplement for details of the basis transformation.
  * Note (6) See Section F of the supplement for additional data.
  * Note (7) See Section G of the supplement for a noise model.
  * Jurcevic  _et al._ (2020) P. Jurcevic, A. Javadi-Abhari, L. S. Bishop, I. Lauer, D. F. Bogorin, M. Brink, L. Capelluto, O. Günlük, T. Itoko, N. Kanazawa, A. Kandala, G. A. Keefe, K. Kruslich, W. Landers, E. P. Lewandowski, D. T. McClure, G. Nannicini, A. Narasgond, H. M. Nayfeh, E. Pritchett, M. B. Rothwell, S. Srinivasan, N. Sundaresan, C. Wang, K. X. Wei, C. J. Wood, J.-B. Yau, E. J. Zhang, O. E. Dial, J. M. Chow, and J. M. Gambetta, “Demonstration of quantum volume 64 on a superconducting quantum computing system,”  (2020), arXiv:2008.08571 .
  * Magesan and Gambetta (2020) E. Magesan and J. M. Gambetta, Phys. Rev. A 101, 052308 (2020).

##  Appendix A Quantum computer code for simulating a tri-junction with longer
topological superconducting wires

Figure 6: A ten qubit setup. On the left we depict the fermion model and on
the right the corresponding qubit “device”. Note that there is an extra qubit
with operators $\sigma_{0}^{\nu}$ which does not correspond to a site in the
fermion model. This qubit is used to enforce fermion commutation relations
while preserving locality Backens  _et al._ (2017).

The three qubit model is convenient due to its simplicity, however, the
overlap of the Majorana modes is controlled by hand through the parameter
$\alpha$ in Eq. (1). Instead, we would like to control the overlap by the
number of qubits which separate the end modes – thus obtain topological
protection. In this model, every arm of the T-junction is composed of multiple
qubits, see Fig. 6. The Hamiltonian becomes,

$\begin{split}H(t)=&i\sum_{a=1}^{3}\left(\mu\sum_{i=0}^{L-1}\gamma_{ia}^{y}\gamma_{ia}^{x}+\Delta\sum_{i=0}^{L-2}\gamma_{ia}^{y}\gamma_{(i+1)a}^{x}\right)\\\
+&i\frac{1}{2}\sum_{a,b,c}J_{ab}(t)\epsilon_{abc}\gamma_{0a}^{x}\gamma_{0b}^{x},\end{split}$
(4)

where $\epsilon_{abc}$ is the Levi-Civita symbol, $\mu$ is the on-site
potential, $\Delta$ is the p-wave coupling term between sites and
$J_{ab}=J_{ba}$ couples the arms of the tri-junction. $\gamma_{ia}^{x}$ and
$\gamma_{ia}^{y}$ are the two Majorana operators that act on site $i$ of arm
$a$. In other words, we have three Kitaev chains coupled at the first site of
each chain. We can model each fermion operator as a string of spin operators
acting on qubits:

$\begin{split}&\gamma_{ia}^{x}=\sigma_{0}^{a}\prod_{j=0}^{i-1}\sigma_{ja}^{z}\sigma_{ia}^{x},\\\
&\gamma_{ia}^{y}=\sigma_{0}^{a}\prod_{i=0}^{i-1}\sigma_{ja}^{z}\sigma_{ia}^{y},\end{split}$
(5)

where $\sigma_{ia}^{\nu}$ is the Pauli $\nu$-matrix acting on the qubit
labeled $ia$. The auxiliary qubit with Pauli matrix $\sigma_{0}^{a}$ is added
to keep string operators local Backens  _et al._ (2017). In terms of these
spin operators, the Hamiltonian is

$\begin{split}H(t)=&\sum_{a=1}^{3}\left(\mu\sum_{i=0}^{L-1}\sigma_{i,a}^{z}+\Delta\sum_{i=0}^{L-2}\sigma_{i,a}^{x}\sigma_{i+1,a}^{x}\right)\\\
+&\frac{1}{2}\sum_{a\neq b\neq
c}J_{ab}(t)\sigma_{0,a}^{x}\sigma_{0,b}^{x}\sigma_{0}^{c}.\end{split}$ (6)

Notice that the auxiliary qubit $\sigma^{c}_{0}$ appears only in the junction
coupling term.

Once again we write the sums in the exponential for the time evolution
operator $U=\mathcal{T}\prod_{i}\exp{\left(-iH(t_{i})dt\right)}$ as products
of exponentials. In this way, time evolution is defined by the gate sets

$e^{-i\mu\sigma_{ia}^{z}dt/2}=R_{z,ia}(\mu dt)$ (7)
$e^{-i\Delta\sigma_{ia}^{x}\sigma_{(i+1)a}^{x}dt/2}=C_{x,ia,(i+1)a}R_{x,ia}(\Delta
dt)C_{x,ia,(i+1)a}$ (8)
$\begin{split}&e^{-iJ_{ab}(t)\sigma_{0,a}^{x}\sigma_{0,b}^{x}\sigma_{0}^{c}dt/2}\\\
&=C_{c,0,0a}C_{x,0b,0a}R_{x,0a}(J_{ab}(t)dt)C_{x,0b,0a}C_{c,0,0a}\end{split}$
(9)

where $C_{\nu,\alpha,\beta}$ is a controlled $\nu\in\\{x,y,z\\}$ gate acting
on the qubits labeled by $\alpha$ and $\beta$.

The qubits are initialized to the two degenerate ground states of the
Hamiltonian written in the Jordan-Wigner basis (Eq. (6)) for $\mu=0$,
$J_{12}=J_{20}=0$, and $J_{01}=\Delta$. These ground states can be reached
from the empty state $\psi_{0}=\ket{000\ldots}$ by applying a particular gate
set. First we initialize the auxiliary third arm of the tri-junction:

$\begin{split}&\psi_{a=3}=\\\
&C_{x,23,23}R_{y,23}\left(\frac{\pi}{2}\right)R_{x,13}(\pi)C_{x,03,13}R_{y,03}\left(\frac{\pi}{2}\right)R_{x,0}(\pi)\psi_{0}.\end{split}$
(10)

This is the same for both ground states. They can be reached by:

$\begin{split}\psi_{\pm}(0)=&\ldots
R_{y,12}\left(\mp\frac{\pi}{2}\right)R_{y,02}\left(\pm\frac{\pi}{2}\right)\\\
&\ldots
R_{y,11}\left(\mp\frac{\pi}{2}\right)R_{y,01}\left(\pm\frac{\pi}{2}\right)\psi_{a=3}\end{split}$
(11)

If we start in the $\psi_{+}(0)$ state and apply the conjugate transpose of
the gate set that defines $\psi_{-}(0)$ at the end of the braiding procedure
then $\psi_{0}$ corresponds to successful braiding.

Although there are currently quantum devices with enough qubits to simulate
the nine fermion device, the extra two qubit gates (6 per time step) make this
simulation impractical. However, as gate fidelities improve checking the
topological protection of Majoranas on the nine fermion device will become
possible.

##  Appendix B Hints of topological protection from the three qubit model

Topological protection from braiding manifests itself as a tolerance to local
perturbations. However, braiding MZMs at a finite rate results in the
topological protection being imperfect. In our case the braiding time is
dictated by the separation of time scales associated with the strength of the
MZM couplings in the tri-junction. To study tolerance to local perturbations
we rewrite the Hamiltonian in Eq. (1) of the main text so that each arm of the
tri-junction has a unique coupling $\alpha_{a}$.

$H=i\sum_{a=0}^{2}\alpha_{a}\gamma_{a}^{x}\gamma_{a}^{y}+i\frac{1}{2}\sum_{a\neq
b}J_{ab}\gamma_{a}^{x}\gamma_{b}^{x}$ (12)

In Fig. 7 we plot the braiding fidelity as a function of a small perturbation
to $\alpha_{0}$. As the procedure becomes more adiabatic, the system becomes
more protected. Each curve in the figure corresponds to a different protocol
time $\tau$. As $\tau$ increases the curve flattens out indicating a higher
tolerance to local perturbations. This is consistent with the expectation that
there is additional protection to local perturbations as the MZM couplings in
the tri-junction arms become weaker, e.g. as the number of qubits that compose
each arm of the tri-junction increase.

To truly get the benefit of topological protection, one must use a tri-
junction model with arms composed of multiple qubits each coupled by local
parameters. The parameter $\alpha_{0}$ can be thought of as a coarse graining
of these local parameters. The system becomes truly topologically protected
only as $\alpha_{0}$ approaches zero. In this case, each step of the braiding
procedure is achieved by tuning the local parameters so that the edge between
topological and trivial superconductors (qubits) can slide towards the tri-
junction, the tri-junction couplings are adjusted, and then the topological-
trivial edge is brought back out to the the end of the tri-junction arm.

Figure 7: Braiding fidelity as a function of a local perturbation. The
probability $P_{-}$ to be in the target quantum state is shown against a shift
$d\alpha_{0}$ of parameter $\alpha_{0}=\alpha+d\alpha_{0}$ on qubit 0 where
$\alpha_{1}=\alpha_{2}=\alpha$. The fidelity is plotted as a function of the
change of $\alpha_{0}$ on the first qubit to the optimized $\alpha_{0}$. The
protection of the braiding procedure depends on the protocol time $\tau$. The
parameter $\alpha$ is optimized for each protocol time.

##  Appendix C Tracking the wave function trajectory during braiding

Figure 8: Probability distributions in the computational subspace during the
first three steps of the braiding procedure. The yellow (Noiseless) and orange
(ibmq_athens) bars are for the state which originates in the $\ket{+}$ state
and the purple (Noiseless) and blue (ibmq_athens) bars are for the state which
originates in the $\ket{-}$ state. The labels on the horizontal axis indicate
where each state would end up during each step of the process if there was no
trotterization error.

The first three steps of the braiding procedure swap the two MZMs. While the
fidelity of the swap is characterized by the target state on the third step,
it is instructive to do an experiment which ends at each of the three steps to
see that the states are following the correct path. To do this we unwind the
basis in which braiding occurs so that the target state after each step is
mapped onto a specific occupation state of the device. This is accomplished by
applying

$\begin{split}R_{y,2}\left(\frac{\pi}{2}\right)R_{x,1}\left(\frac{\pi}{2}\right)R_{y,0}\left(\frac{\pi}{2}\right)&\text{
after step 1,}\\\
R_{y,2}\left(\frac{\pi}{2}\right)R_{y,1}\left(\frac{\pi}{2}\right)R_{x,0}\left(\frac{\pi}{2}\right)&\text{
after step 2, and}\\\ R_{y,0}\left(\frac{\pi}{2}\right)&\text{ after step
3.}\end{split}$ (13)

Steps 4 and 5 require several CNOT gates to unwind the basis, therefore, we
end this analysis at step 3. The $\ket{+}$ and $\ket{-}$ states stay separated
throughout the process, as seen by the separation of probability distributions
in Fig. 8. This observation supports our interpretation that the quantum
device is truly following the expected Hamiltonian evolution and is not ending
in the correct state by other means.

##  Appendix D Implementation with Qiskit Pulse

IBM quantum computers provide highly calibrated single-qubit and CNOT gates
Jurcevic  _et al._ (2020). Virtual single-qubit $Z$-gates are created by phase
shifting subsequent microwave drives McKay  _et al._ (2017). In a basis gate
quantum computation where the only two-qubit operation is the CNOT gate, the
operator $U_{yx}(\theta)=e^{-i\theta\sigma^{y}_{1}\sigma^{x}_{2}/2}$ is built
with an $R_{y,1}(\theta)$ rotation sandwiched between two CNOT gates, see Fig.
4a of the main text. However, the same operation can be implemented in a
pulse-efficient manner by moving the parameter $\theta$ into a single two-
qubit CR gate which ideally implements
$U_{zx}(\theta)=e^{-i\theta\sigma^{z}_{1}\sigma^{x}_{2}/2}$ Chow  _et al._
(2011). The required circuit is thus $R_{x,1}(\pi/2)\cdot U_{zx}(\theta)\cdot
R_{x,1}(-\pi/2)$, see Fig. 4b of the main text.

Pulse-level control of IBM backends is enabled by Qiskit Pulse McKay  _et al._
(2018); Alexander  _et al._ (2020). The CNOT basis gates are built from echoed
CR pulses which consist of the $CR(\pm\pi/4)$ pulses sandwiching an
$X$-rotation echo pulse applied to the control qubit to cancel undesired $ZI$
and $IX$ terms of the CR Hamiltonian Sheldon  _et al._ (2016). The
compensation rotary pulses $C(\pm\pi/4)$ applied to the target qubit suppress
the remaining non-commuting $ZZ$ and $IY$ terms Sundaresan  _et al._ (2020).

The $CR(\pm\pi/4)$ pulses are calibrated to have the shortest duration while
minimizing leakage outside of the computational basis, as determined by
randomized benchmarking Sheldon  _et al._ (2016). The effect of decoherence is
thus minimized and the pulse amplitude that retains the qubits in the
computational subspace is maximized. Since the largest errors arise when
performing CR pulses we wish to minimize their duration in the braiding
algorithm. This is achieved by (a) using the $U_{zx}(\theta)$ gate to
implement $U_{yx}(\theta)$ instead of two CNOT gates and (b) by creating the
pulse schedules for $U_{zx}(\theta)$ by modifying the highly-calibrated CNOT
pulse schedules Gokhale  _et al._ (2020); Alexander  _et al._ (2020). The
rotation angle $\theta$ depends on the area under the pulses, and is often
considerably less than $2\times\pi/2$ for the Trotterization that braids MZM.
This allows us to build $U_{zx}(\theta)$ gates with considerably shorter
duration CR pulses than two CNOT gates, hence introducing less error per
Trotter step. Since $CR(\pm\theta)$ and $C_{\pm}$ are all implemented with
flat-top Gaussian pulses and that $\theta$ depends non-linearly on the pulse
amplitude $A$ Magesan and Gambetta (2020) we modify the pulse area by first
stretching and compressing the width $w$ of the flat-tops. We only scale the
amplitude when $w=0$. This avoids any additional calibration as the relation
$\theta(w)$ is linear. Avoiding additional calibration on cloud-based quantum
computers is paramount when running jobs through a queue. The $CR(\pm\pi/4)$
pulses are given in the parametric form GaussianSquare$(d,A,\sigma,w)$ with an
area given by

$\displaystyle\alpha^{*}=|A^{*}|w^{*}+|A^{*}|\sigma\sqrt{2\pi}\,\text{erf}(n_{\sigma}).$
(14)

Here, $n_{\sigma}$ is the number of standard deviations $\sigma$ contained in
the pulse with total duration $d$ and flat-top width $w$ and amplitude $A$.
The quantities $d$, $w$, and $\sigma$ are all specified in units of Arbitrary
Waveform Generator (AWG) samples which last $0.222\leavevmode\nobreak\ {\rm
ns}$ on ibmq_athens. The star superscript indicates that we are referring to
the parameters of the CNOT schedule. To scale the CR gates we first calculate
the target area of each pulse

$\displaystyle\alpha(\theta)=\frac{\theta}{\pi/2}\alpha^{*}.$ (15)

As long as $\alpha(\theta)>|A^{*}|\sigma\sqrt{2\pi}\,\text{erf}(n_{\sigma})$
we change the width of the pulse following

$\displaystyle
w(\theta)=\frac{\alpha(\theta)}{|A^{*}|}-\sigma\sqrt{2\pi}\,\text{erf}(n_{\sigma}).$
(16)

When $\alpha(\theta)<|A^{*}|\sigma\sqrt{2\pi}\,\text{erf}(n_{\sigma})$, i.e.
when the flat-top vanishes, we instead scale the amplitude of the remaining
Gaussian pulse according to

$\displaystyle|A(\theta)|=\frac{\alpha(\theta)}{\sigma\sqrt{2\pi}\,\text{erf}(n_{\sigma})}.$
(17)

The phase of the pulse, i.e. $\arg(A)$, is left unchanged to implement the
$\pm$ rotations in the echo. Since the AWGs can only load pulses if their
duration is a multiple of $m=16$ samples we set the duration of our pulses to

$\displaystyle d=\left\lceil\frac{w(\theta)+n_{\sigma}\sigma}{m}\right\rceil
m\quad\text{samples}.$ (18)

The pulse schedule implementing $U_{yx}(\theta)$, shown in Fig. 4c of the main
text, has three single-qubit pulses. The first pulse is the $R_{x,1}(\pi/2)$
seen in Fig. 4b of the main text. The second pulse is the $R_{x,1}(\pi)$
needed in the echoed CR gate. The third pulse is $R_{x,1}(\pi/2)$ which
corresponds to the second $R_{x,1}(\pi)$ pulse in the echoed CR gate together
with the $R_{x,1}(-\pi/2)$ seen in Fig. 4b of the main text.

Figure 9: Pulse schedules corresponding to the (a) scaled circuit and (b)
basis circuit, corresponding to gates given in Fig. 2, with the end of the
scaled schedule marked by the vertical gray dashed line. Here, the ‘d’ and ‘u’
correspond to the (on-resonant) DriveChannels and (off-resonant)
ControlChannels of Qiskit Pulse, respectively. The pulses are shaded according
to phase, with light pulses in-phase and dark pulses in the quadrature phase.
The circular arrows correspond to $Z$-rotations executed in software by
changing the phase of subsequent pulses McKay  _et al._ (2017). The duration
is given in units of sample time, $dt=0.222$ ns.

A side-by-side comparison of the pulse schedules generated by a highly-
transpiled (optimization level 3) circuit using the basis CNOTs to generate
$U_{yx}(\theta)$ interactions as in Fig. 4a and the scaled pulses of Fig. 4b-c
is presented in Fig. 9 to highlight their differences. These pulse schedules
correspond to the circuits of Fig. 2 for a single Trotter step, in this case
the second ($n_{t}=1$) step of the first Braiding step, as $J_{01}\rightarrow
0$ and $J_{12}\to J_{\rm max}$. The duration of the scaled circuit is
$5952/9200=62\%$ of the basis circuit, leading to a reduction in the amount of
decoherence the would occur regardless of the applied pulses. Furthermore, the
bulk of errors on current noisy quantum hardware occur during the two-qubit
operations, as depicted by the yellow pulses on the ControlChannels in both
schedules. The area of the two-qubit pulses is substantially less when scaling
the CR pulses than when using CNOT gates which indicates why the scaled pulses
yield a successful MZM braiding that is impossible to observe with CNOT-based
circuits. For the Trotter step in Fig. 9 the ratio of areas shows that the
scaled CR is $463.51/1544.87=30\%$ that of the basis CNOTs.

##  Appendix E Rotated basis

The Hamiltonian of the three dot system (Eq. (2) in the main text) has a thee-
qubit coupling term $J_{20}(t)\sigma_{0}^{y}\sigma_{1}^{z}\sigma_{2}^{x}$
whose evolution operator cannot be continuously generated from scaled CR
gates. We therefore rotate the Hamiltonian into a basis where there are no
three-qubit coupling terms.

Take $U=(a+b)/\sqrt{2}$ as a general unitary operator. We want $U$ to
transform all terms in the Hamiltonian into terms with less than three qubits.
That is $U^{\dagger}\sigma_{0}^{y}\sigma_{1}^{z}\sigma_{2}^{x}U$,
$U^{\dagger}\sigma_{1}^{y}\sigma_{2}^{x}U$,
$U^{\dagger}\sigma_{0}^{y}\sigma_{1}^{x}U$, and $U^{\dagger}\sigma_{a}^{z}U$
must contain fewer than three Pauli matrices for all $a$. Take
$U=(a+b)/\sqrt{2}$ as a general unitary operator and $\lambda$ to be a general
Hermitian operator. Then we have that

$U^{\dagger}\lambda U=\begin{cases}\lambda&\text{if}\quad a\lambda(\lambda
a)^{-1}=+b\lambda(\lambda b)^{-1}\\\ ab\lambda&\text{if}\quad a\lambda(\lambda
a)^{-1}=-b\lambda(\lambda b)^{-1}\end{cases}$ (19)

In other words, we need to find an $a$ and $b$ such that $a$ and $b$ have
different commutation relations with
$\sigma_{0}^{y}\sigma_{1}^{z}\sigma_{2}^{x}$ and $ab\lambda$ is a one or two
qubit operator. Additionally, for all other operators in the Hamiltonian
$\lambda^{\prime}$, we need that either $ab\lambda^{\prime}$ is a one or two
qubit operator or for $a$ and $b$ to have the same commutation relation with
$\lambda^{\prime}$.

Notice that $a=\sigma^{z}_{1}$ and $b=\sigma^{y}_{0}\sigma^{x}_{1}$ have all
of the above properties. Applied to the qubit Hamiltonian (Eq. (2) of the main
text), we get:

$\begin{split}\bar{H}(t)&=(\sigma^{z}_{1}+\sigma^{y}_{0}\sigma^{x}_{1})H(t)(\sigma^{z}_{1}+\sigma^{y}_{0}\sigma^{x}_{1})/2\\\
&=\alpha\left(\sigma_{0}^{x}\sigma_{1}^{y}+\sigma_{0}^{y}\sigma_{1}^{x}+\sigma_{2}^{z}\right)\\\
&+J_{01}(t)\sigma_{1}^{z}+J_{12}(t)\sigma_{1}^{y}\sigma_{2}^{x}+J_{20}(t)\sigma_{1}^{x}\sigma_{2}^{x}.\end{split}$
(20)

We have exchanged the three-qubit term for two two-qubit terms. Since there
are more multi-qubit terms overall, there is no reason to expect that this
basis will have less error if the braiding procedure is implemented using
basis gates. However, since there are no three-qubit gates, we can simulate
the entire braiding process using the scaled CR gates.

To initialize the state in this basis we use the initialization in Eq. (3) of
the main text and apply the basis rotation gate set:

$R_{z,1}(\pi)C_{y,01}R_{y,0}(-\pi/2)C_{y,01}$ (21)

The last step to implementing the Hamiltonian evolution in this basis is to
translate the evolution operators into quantum gates. To do this, we
Trotterize the evolution operator and use the generalized unitary-to-gate
transcription.

$\begin{split}&e^{-i\phi\sigma^{a}_{i}/2}=R_{a,i}(\phi)\\\
&e^{-i\phi\sigma^{a}_{i}\sigma^{b}_{j}/2}=C_{a,ij}R_{b,j}(\phi)C_{a,ij}\end{split}$
(22)

which holds for $b\neq z$ which is always the cases in this basis.
Alternatively, we can use the scalled $ZX(\theta)$ gate described in the main
text.

$\begin{split}&e^{-i\phi\sigma^{a}_{i}\sigma^{b}_{j}/2}=\\\
&R_{bj}(\pi/2)R_{ai}(\pi/2)ZX(\theta)R_{ai}(-\pi/2)R_{bj}(-\pi/2)\end{split}$
(23)

##  Appendix F Noise induced drift in the braiding experiment

In Fig. 10 we plot all of the data accumulated to construct Fig. 5 of the main
text, including data at large delay times which was cut off in the figure in
the main text as it shows no signs of braiding. The data is acquired in
sequential order with all of the low time-delay data taken before the long
time-delay data.

Figure 10: Bias towards braiding as a function of delay time added after each
$ZX$ rotation. Similar to Fig. 5 in the main text but the data has been
extended to longer time delays.

When the delay is greater than $\sim 200\leavevmode\nobreak\ {\rm ns}$ the
probability to end up in $\ket{+}$ or $\ket{-}$ is essentially independent of
the initial state, indicating that the added noise overwhelms the quantum
braiding dynamics. We attribute the delay dependency of
$(P_{-}-P_{+})/(P_{+}+P_{-})$ to drifts in the quantum device over the course
of the experiment. Similar drifts are also seen in Fig. 4d of the main text,
where the fidelity of the double-CNOT circuit oscillates with $\theta$.

##  Appendix G Error Model

Figure 11: (a) Ratio of the probability of successfully braiding to the total
probability of being in the computational space. (b) Probability of being in
the computational space. The blue and orange curves show the scaled CR gates
and the basis CNOT gates, respectively. The gray shaded region is the typical
range of CNOT errors in real devices. While the basis gate procedure has
little or no bias towards successfully braiding, the scaled CR procedure is
clearly biased towards successfully braiding. The difference in fidelity ratio
at zero CNOT error is due to an increased number of single-qubit gates in the
scaled CR procedure. For these plots, we have $\alpha=0.2\leavevmode\nobreak\
J_{\rm max}$ and $\tau=3.3\leavevmode\nobreak\ J_{\rm max}$. We find that the
best balance between Suzuki-Trotter error and gate error is to perform each
protocol in thee time steps.

To illustrate the advantage of scaled gates over basis ones for quantum
simulation we construct a simple error model that describes errors accumulated
in the CR procedure. The error model consists of single bit flip errors which
are proportional to the duration of the CR gate, described by the
superoperator

$\begin{split}&\text{Error}_{1}=(1-\epsilon)\text{SuperOp}(\hat{I})+\epsilon\text{SuperOp}(\hat{X})\\\
&\text{Error}_{2}=\text{Error}_{1}\otimes\text{Error}_{1}\end{split}$ (24)

where $\text{Error}_{1}$ refers to a single-qubit error and $\text{Error}_{2}$
refers to a two-qubit error, $\text{SuperOp}(\hat{O})$ denotes the
superoperator of operator $\hat{O}$ and $\epsilon$ is determined from the
device.

For the ibmq_athens device, the single-qubit gate errors range from $2.2\times
10^{-4}$ to $2.8\times 10^{-4}$. For the CR error, we take the CNOT error and
scale it to the phase shift we want to apply. In other words,
$\epsilon=(\phi/\pi)\epsilon_{CNOT}$, were $\epsilon_{CNOT}$ is the CNOT error
and $\phi$ is the desired rotation angle. The CNOT errors in the ibmq_athens
device range from $6.9\times 10^{-3}$ to $9.4\times 10^{-3}$. In Fig. 11, we
plot probabilities for both the basis gate procedure using controlled gates
and the scaled procedure using CR gates for a range of CNOT errors. The top
panel plots the probability of braiding if we project onto the computational
basis $(P_{+}-P_{-})/(P_{+}+P_{-})$ while the bottom panel shows the
probability of being in the the computational basis $P_{+}+P_{-}$. The blue
curve is generated using the scaled CR procedure while the orange curve is for
the basis gate procedure. Typical values of the CNOT error are shaded in gray.
There is no bias for the basis gate procedure in the shaded region, but there
is a bias for the scaled CR gates. For the real pulse experiment, the braiding
probability (Fig. 5 of the main text) is in the expected range, however, there
is only a small bias for the computational subspace ($P_{+}+P_{-}=0.28\pm
0.01$) which suggest that there are sources of error that we do not include in
our model.