# Magic state parity-checker
with pre-distilled components

Earl T. Campbell  Department of Physics & Astronomy, University of Sheffield,
Sheffield, S3 7RH, United Kingdom<EMAIL_ADDRESS>https://earltcampbell.com/ Mark Howard  Department of Physics & Astronomy,
University of Sheffield, Sheffield, S3 7RH, United Kingdom.
<EMAIL_ADDRESS>https://markhoward.info

(September 3, 2024)

###### Abstract

Magic states are eigenstates of non-Pauli operators. One way of suppressing
errors present in magic states is to perform parity measurements in their non-
Pauli eigenbasis and postselect on even parity. Here we develop new protocols
based on non-Pauli parity checking, where the measurements are implemented
with the aid of pre-distilled multiqubit resource states. This leads to a two
step process: pre-distillation of multiqubit resource states, followed by
implementation of the parity check. These protocols can prepare single-qubit
magic states that enable direct injection of single-qubit axial rotations
without subsequent gate-synthesis and its associated overhead. We show our
protocols are more efficient than all previous comparable protocols with
quadratic error reduction, including the protocols of Bravyi and Haah.

Error corrected quantum computers require additional gadgets and tricks to
enable fully universal, fault-tolerant quantum computation (see Ref. [Campbell
et al.(2017)Campbell, Terhal, and Vuillot] for a review). Of the various
competing approaches, magic state distillation is especially efficient. At
first, magic state distillation was conceived of as a way to implement the $T$
gate (also called the $\pi/8$ phase gate), with successive generations of
protocols making ever greater improvements [Bravyi and Kitaev(2005), Meier et
al.(2013)Meier, Eastin, and Knill, Bravyi and Haah(2012), Jones(2013a), Haah
et al.(2017)Haah, Hastings, Poulin, and Wecker]. Any desired unitaries could
be approximated by efficient synthesis of $T$ gates and Clifford gates, with
recent years also bringing new, optimised synthesis methods [Kliuchnikov et
al.(2013)Kliuchnikov, Maslov, and Mosca, Gosset et al.(2014)Gosset,
Kliuchnikov, Mosca, and Russo, Ross and Selinger(2016), Paetznick and
Svore(2014), Bocharov et al.(2015)Bocharov, Roetteler, and Svore]. However, we
can circumvent the need for synthesis if we instead prepare magic states
tailored for injecting specific gates. Preparing tailored single qubit magic
states can directly provide single qubit rotations other than the $T$-gate and
distillation routines for these states have been proposed [Landahl and
Cesare(2013), Jones(2013b), Duclos-Cianci and Poulin(2015), Campbell and
O’Gorman(2016), Haah(2017)]. Complex multi-qubit circuits can also be directly
injected once multi-qubit magic states have been distilled [Jones(2013c),
Eastin(2013), Campbell and Howard(2017a), Campbell and Howard(2017b), Haah et
al.(2017)Haah, Hastings, Poulin, and Wecker].

Here we propose a family of two-step protocols for distilling magic states
that are tailored for injecting a desired single-qubit $Z$-axis rotation. A
distinctive feature is that the first step creates a multi-qubit magic state
using the synthillation protocol [Campbell and Howard(2017a), Campbell and
Howard(2017b)]. In the second step, this multi-qubit resource is then used to
fault-tolerantly perform a parity check in a non-Pauli basis. Combined, the
protocol takes in single qubit states and outputs single qubit states, with
multi-qubit magic states appearing only fleetingly.

A single round of our protocol will quadratically reduce errors using fewer
inputs per output than any previous protocol with quadratic error reduction.
Higher order error reductions can be achieved by concatenation of our
protocols. Higher order reductions in noise are possible without concatenation
by using sophisticated protocols that give a large error reduction in a single
round [Jones(2013a), Haah et al.(2017)Haah, Hastings, Poulin, and Wecker, Haah
et al.(2018)Haah, Hastings, Poulin, and Wecker, Hastings and Haah(2018)], but
there are important practical considerations for why one might favour a
concatenated approach (see Sec. 6.2 for a discussion).

We begin by covering some basic notation. Sec. 2 gives an overview of our
approach. Note that the first step is our previously proposed synthillation
method [Campbell and Howard(2017a), Campbell and Howard(2017b)], so the
details will not be repeated here. Rather we focus on how non-Pauli parity
checking is possible given these pre-distilled resources, giving a detailed
explanation in Sec. 3. The protocol’s performance is analysed in Sec. 4. In
Sec. 5, we present a small bonus result that proves equivalent performance is
unlikely to be possible using codes with conventional transversal gate
constructions.

##  1 Notation

We denote axial rotations about the Pauli-$Z$ axis as

$R(\theta)=\exp(i\theta Z)=\cos(\theta){\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+i\sin(\theta)Z.$ (1)

If the angle is $\theta=\pi/2^{\ell}$ for integer $\ell$ then $R(\theta)$
belongs in the $\ell^{\mathrm{th}}$ level of the Clifford hierarchy [Gottesman
and Chuang(1999)]. Therefore, $R(\pi/8)$ is the $\pi/8$-phase gate, also known
at the $T$ gate. Unitaries inside the Clifford hierarchy are special because
they can be realised using state-injection and a bounded number of appropriate
magic states. However, all the analysis in this paper holds for any $\theta$,
even values corresponding to unitaries and magic states not connected to the
Clifford hierarchy.

We use $W(\theta)$ for the Hermitian operator
$W(\theta):=R(\theta)XR(\theta)^{\dagger}=R(2\theta)X$. Note that $W(\pi/8)$
is a Clifford and plays a similar role to the Hadamard; it interchanges $X$
and $Y$ whereas the Hadamard interchanges $X$ and $Z$. The relevant magic
states are eigenstates of $W(\theta)$ and sit on the equator of the Bloch
sphere

$|R(\theta)\rangle=W(\theta)|R(\theta)\rangle=R(\theta)|+\rangle.$ (2)

More generally, when $U$ is a diagonal gate (acting on $n$ qubits) we use
$|U\rangle:=U(|+\rangle^{\otimes n})$. In this notation, the familiar $T$
state is $|R(\pi/8)\rangle$. We use CZ for control-Z and CCZ for control-
control-Z.

##  2 Overview of new protocols

###  2.1 Protocols for $\pi/8$ phase gates

The protocols presented here can be used both for distillation of $T$-states
that can implement $\pi/8$ phase gates, or more generally for distillation of
$|R(\theta)\rangle$ magic states for smaller angle rotations. We begin by
sketching the simple case of $T$-state distillation. We say a protocol is an
$n\rightarrow k$ protocol if it takes $n$ inputs and outputs $k$ magic states
with some success probability (typically this probability approaches unity in
the low noise limit). Our two-step protocols for $T$-state distillation are
$3k+4\rightarrow k$ protocols for even $k$ with quadratic reduction of noise.
The resource overhead of the protocol is roughly captured by
$n/k=3+\frac{4}{k}$, which approaches 3 for large $k$. For $k=2$ we have a
$10\rightarrow 2$ protocol and so the protocol is very similar to the MEK
proposal proposed by Meier, Eastin and Knill [Meier et al.(2013)Meier, Eastin,
and Knill]. Therefore, compared to MEK, we reduce the $n/k$ overhead from 5 to
3 by going to larger blocks sizes (larger $k$). Another class of protocols was
proposed by Bravyi and Haah [Bravyi and Haah(2012)], which are
$3k+8\rightarrow k$ protocols for even $k$ with quadratic reduction of noise.
The Bravyi-Haah protocols also have $n/k\rightarrow 3$ in the large protocol
limit. Both our protocols and the Bravyi-Haah protocols have the same
asymptotic limit, but our protocols approach this limit faster. For example,
to achieve $n/k=4$ we can use modest size $16\rightarrow 4$ protocols, whereas
the comparable Bravyi-Haah protocol is $32\rightarrow 8$ with double the block
size. This effect becomes more pronounced as the protocol is concatenated.

Figure 1: An illustration of how the two steps of our protocol are chained
together for the special case $N=1$. The combined process is described in Eq.
(9) for general $N$.

Furthermore, our work presents a different approach to distillation as we
break the process up into two steps, making use of a multi-qubit magic state
resource. Let us describe how this two-step process works, first considering
the $10\rightarrow 2$ case. In the first step we prepare a “pre-distilled”
magic state that can inject a CCZ (or Toffoli) gate. It has been known for
several years that a single CCZ magic state can be prepared from 8 $T$-states
with quadratic error reduction [Jones(2013c), Eastin(2013)]. It is no longer
appropriate to describe this process in simple $n\rightarrow k$ notation and
so we introduce the more detailed description of these protocols as being

$\left\\{\begin{array}[]{c}8\\\ |R(\pi/8)\rangle\\\
\epsilon_{\pi/8}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}1\\\
|CCZ\rangle\\\ O(\epsilon_{\pi/8}^{2})\end{array}\right\\}.$ (3)

In this notation, the left hand side gives the input resources and the right
hand size gives the output resources. The top line gives the quantity of
resources, the middle describes the species of magic state and the bottom line
gives the infidelity. Our second step, which was not previously known, is to
notice that a single $|CCZ\rangle$ resource can be used to check the parity in
a pair of $|R(\pi/8)\rangle$ states, which implements the transform

$\left\\{\begin{array}[]{cc}2&1\\\ |R(\pi/8)\rangle&|CCZ\rangle\\\
\epsilon_{\pi/8}&\epsilon_{\\#}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2\\\
|R(\pi/8)\rangle\\\
O(\epsilon_{\pi/8}^{2}+\epsilon_{\\#})\end{array}\right\\}.$ (4)

Chaining these steps together so that $\epsilon_{\\#}=O(\epsilon_{\pi/8}^{2})$
yields

$\left\\{\begin{array}[]{c}10\\\ |R(\pi/8)\rangle\\\
\epsilon_{\pi/8}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{cc}2&1\\\
|R(\pi/8)\rangle&|CCZ\rangle\\\
\epsilon_{\pi/8}&O(\epsilon_{\pi/8}^{2})\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2\\\
|R(\pi/8)\rangle\\\ O(\epsilon_{\pi/8}^{2})\end{array}\right\\},$ (5)

or simply $10\rightarrow 2$ for short.

Next we outline the two-step process for larger block protocols. Instead of
using $|CCZ\rangle$ resources, larger block protocols will use a resource that
we denote as $|CCZ_{\\#N}\rangle$. This resource can inject a gate
$CCZ_{\\#N}$, which is $N$ copies of the CCZ gate all sharing one control
qubit in common, and the relevant magic state is simply

$|CCZ_{\\#N}\rangle:=CCZ_{\\#N}(|+\rangle^{\otimes 2N+1}).$ (6)

In the first step, we borrow results on synthillation [Campbell and
Howard(2017a), Campbell and Howard(2017b)] that provide protocols implementing

$\left\\{\begin{array}[]{c}4N+4\\\ |R(\pi/8)\rangle\\\
\epsilon_{\pi/8}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}1\\\
|CCZ_{\\#N}\rangle\\\ O(\epsilon_{\pi/8}^{2})\end{array}\right\\}.$ (7)

This protocol is described in the section "$U_{N\\#}$ family" of Ref.
[Campbell and Howard(2017b)]. To avoid repetition, here we simply treat this
synthillation routine as a black box with known properties. Rather, here we
focus on the second step — the main technical contribution of this work — by
showing that a pre-distilled $|CCZ_{\\#N}\rangle$ can be used to parity check
on $2N$ magic states

$\left\\{\begin{array}[]{cc}2N&1\\\ |R(\pi/8)\rangle&|CCZ_{\\#N}\rangle\\\
\epsilon_{\pi/8}&\epsilon_{\\#}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2N\\\
|R(\pi/8)\rangle\\\
O(\epsilon_{\pi/8}^{2}+\epsilon_{\\#})\end{array}\right\\}.$ (8)

Chaining steps one and two, so that $\epsilon_{\\#}=O(\epsilon_{\pi/8}^{2})$,
we obtain a family of protocols for all integer $N$

$\left\\{\begin{array}[]{c}6N+4\\\ |R(\pi/8)\rangle\\\
\epsilon_{\pi/8}\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{cc}2N&1\\\
|R(\pi/8)\rangle&|CCZ_{\\#N}\rangle\\\
\epsilon_{\pi/8}&O(\epsilon_{\pi/8}^{2})\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2N\\\
|R(\pi/8)\rangle\\\ O(\epsilon_{\pi/8}^{2})\end{array}\right\\},$ (9)

or simply $6N+4\rightarrow 2N$ for short. Alternatively, using $k=2N$ we have
a family of $3k+4\rightarrow k$ protocols with even $k$. So we have returned
to the choice of symbols used by Bravyi and Haah who found $3k+8\rightarrow k$
protocols with even $k$. Fig. 1 illustrates this for the $N=1$ ($k=2$) case.
This completes the outline of our protocols for $T$-state distillation.

###  2.2 Protocols for general phase gates

Next, we describe a family of protocols for other equatorial magic states
$|R(\theta)\rangle$. Step one will remain the same, again making use of
synthillation of $|CCZ_{\\#N}\rangle$ resource states. Step two generalises to

$\displaystyle\left\\{\begin{array}[]{ccc}2N&1&N\\\
|R(\theta)\rangle&|CCZ_{\\#N}\rangle&R(2\theta)\\\
\epsilon_{\theta}&\epsilon_{\\#}&\eta\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2N\\\
|R(\theta)\rangle\\\
O(\epsilon_{\theta}^{2}+\epsilon_{\\#}+\eta)\end{array}\right\\}.$ (16)

When $R(2\theta)$ appears uncluttered by a ket, as it does above, it refers to
inputting a phase gate $R(2\theta)$ rather than a magic state. Most
interesting is the case when $R(2\theta)$ sits in the Clifford hierarchy and
so can be injected using appropriate magic states. Proving the validity of the
above mapping is at the core of this paper. Chaining this with synthillation
yields an overall protocol

$\displaystyle\left\\{\begin{array}[]{@{\hspace{0pt}}ccc@{\hspace{0pt}}}2N&4N+4&N\\\
|R(\theta)\rangle&|R(\pi/8)\rangle&R(2\theta)\\\
\epsilon_{\theta}&\epsilon_{\pi/8}&\eta\end{array}\right\\}$
$\displaystyle\rightarrow\left\\{\begin{array}[]{@{\hspace{0pt}}ccc@{\hspace{0pt}}}2N&1&N\\\
|R(\theta)\rangle&|CCZ_{\\#N}\rangle&R(2\theta)\\\
\epsilon_{\theta}&O(\epsilon_{\pi/8}^{2})&\eta\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{@{\hspace{0pt}}c@{\hspace{0pt}}}2N\\\
|R(\theta)\rangle\\\
O(\epsilon_{\theta}^{2}+\epsilon_{\pi/8}^{2}+\eta)\end{array}\right\\}.$ (26)

It is important to notice that there is no noise reduction in $\eta$.
Therefore, it is crucially important they are predistilled (e.g.
$\eta\sim\epsilon^{2}$) and so we refer to the rotations $R(2\theta)$ as
pivotal. Despite the need for high fidelity pivotal rotations, other protocols
have used pivotal rotations and found significant reductions of resource costs
compared against using gate-synthesis. For instance, the $N=1$ two-step
protocol has an identical resource cost to the protocols introduced by
Campbell and O’Gorman [Campbell and O’Gorman(2016)]. Pivotal rotations (though
not under this name) played a similar role in the protocols proposed by
Duclos-Cianci and Poulin [Duclos-Cianci and Poulin(2015)], which in our
notation can be described as

$\displaystyle\left\\{\begin{array}[]{ccc}2N&8N+4&N\\\
|R(\theta)\rangle&|R(\pi/8)\rangle&R(2\theta)\\\
\epsilon_{\theta}&\epsilon_{\pi/8}&\eta\end{array}\right\\}\rightarrow\left\\{\begin{array}[]{c}2N\\\
|R(\theta)\rangle\\\
O(\epsilon_{\theta}^{2}+\epsilon_{\pi/8}^{2}+\eta)\end{array}\right\\}.$ (33)

Note that for the majority of their paper, Duclos-Cianci and Poulin only
discuss the $N=1$ case, but they do sketch the higher $N$ case later in the
paper.

One could say our protocols are compressed as they essentially give a slight
compression in $T$-cost of the Duclos-Cianci and Poulin protocols [Duclos-
Cianci and Poulin(2015)]. Our protocols also have very different inner
workings. This provides a new perspective on magic state distillation, but
also the two-step feature has a potentially significant practical advantage.
The exotic resources are higher value since the $R(2\theta)$ pivotal rotation
and $|R(\theta)\rangle$ resources are more difficult to prepare than standard
$|R(\pi/8)\rangle$ magic states. However, in the two-step protocols, one does
not risk using the exotic resources until the first step has succeeded. This
contrasts, with both the Duclos-Cianci-Poulin protocols and Campbell-O’Gorman
protocols for which all the resources are committed at the same time, with a
single error anywhere leading to loss of all resources.

Figure 2: Gadgets for measuring parity in the $W(\theta)$ basis. In (a) we
illustrate how control-$W(\theta)$ gates are used to measure parity. In (b) we
show how the control-$W(\theta)$ gates can be decomposed into CNOT and
control-phase gates. In the text we describe this phase-gate circuit by $V$
(see Eq. (35)), which can be broken up into a product of $U_{j}$ gates (see
Eq. (38)).

##  3 Implementing step two

Here we show how to use $CCZ$ circuits, implemented using synthillation, to
perform a parity check. We will construct circuits that measure the parity of
$2N$ qubits in the basis of $W(\theta)$. Using an ancilla and a
control-$W(\theta)^{\otimes 2N}$ gate would achieve this, but a smaller
resource overhead is needed if we instead use the parity check gadgets in Fig.
2. We use a combination of control-$W(\theta)$ that triggers when the control
bit is in the $|1\rangle$ state with an unconventional control-$W(\theta)$
(shown with open circles) that triggers when the control bit is in the
$|0\rangle$ state. This parity check gadget measures in the $W(\theta)$ basis,
but with an additional $W(\theta)$ rotation on half the qubits. This resulting
Kraus operator for the desired even parity outcome is

$\displaystyle K$ $\displaystyle=\frac{1}{2}({\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}^{\otimes N}\otimes W(\theta)^{\otimes N}+W(\theta)^{\otimes
N}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}^{\otimes N})$
$\displaystyle=\frac{1}{2}({\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}^{\otimes N}\otimes
W(\theta)^{\otimes N})({\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}^{\otimes
2N}+W(\theta)^{\otimes 2N}).$ (34)

We assume the noisy $|R(\theta)\rangle$ magic states are diagonal in the
$W(\theta)$ basis, and so the additional $W(\theta)$ rotations have no effect.
The assumption of diagonal noise is mild and can be removed at the expense of
hideously complicating the noise analysis (see e.g. App D of Ref. [Campbell
and O’Gorman(2016)]).

Using $W(\theta)=R(2\theta)X$ it follows that this parity measurement circuit
can be split into a sequence of control-$X$ gates followed by a phase gate
circuit (see Fig. 2b). Algebraically, this phase gate circuit is

$\displaystyle V$ $\displaystyle=\prod_{j=1,\ldots,N}\left[|0\rangle\langle
0|_{0}R_{2j}(2\theta)+|1\rangle\langle 1|_{0}R_{2j-1}(2\theta)\right],$ (35)

where the subscripts denote which qubits the rotations act on, with qubit
labels running from $0$ to $2N$. Using the shorthand

$\displaystyle U_{j}$ $\displaystyle:=|0\rangle\langle
0|_{0}R_{2j}(2\theta)+|1\rangle\langle 1|_{0}R_{2j-1}(2\theta),$ (36)

we have

$\displaystyle V$ $\displaystyle=\prod_{j=1,\ldots,N}U_{j}.$ (37)

To recap, given $|R(\theta)\rangle$ of error rate $\epsilon$ and the ability
to implement $V$ we can parity check in the $W(\theta)$ basis, outputting
$|R(\theta)\rangle$ of error rate $O(\epsilon^{2})$ with some probability
$p=1-O(\epsilon)$.

Next, we show how to implement $V$ with some pre-distilled resources. First,
we use the decomposition $R(2\theta)=\cos(2\theta){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}+i\sin(2\theta)Z$ to expand out $U_{j}$ as

$\displaystyle U_{j}=$ $\displaystyle|0\rangle\langle
0|_{0}(\cos(2\theta){\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+i\sin(2\theta)Z_{2j})+|1\rangle\langle
1|_{0}(\cos(2\theta){\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+i\sin(2\theta)Z_{2j-1}).$ (38)

Collecting the $\cos$ and $\sin$ terms, we have

$\displaystyle U_{j}=$ $\displaystyle\cos(2\theta){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}+i\sin(2\theta)[|0\rangle\langle 0|_{0}Z_{2j}+|1\rangle\langle
1|_{0}Z_{2j-1}]$ (39) $\displaystyle=$
$\displaystyle\cos(2\theta){\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+i\sin(2\theta)M_{j},$

where we have introduced further shorthand

$\displaystyle M_{j}$ $\displaystyle:=|0\rangle\langle
0|_{0}Z_{2j}+|1\rangle\langle 1|_{0}Z_{2j-1}$
$\displaystyle=Z_{2j}(|0\rangle\langle 0|_{0}+|1\rangle\langle
1|_{0}Z_{2j}Z_{2j-1})$ $\displaystyle=Z_{2j}CZ_{0,2j}CZ_{0,2j-1},$ (40)

which is unitary, Hermitian and Clifford.

Figure 3: In (a) we show how $U_{j}$ (a pair of control phase gates) can be
implemented using an ancilla, a control-$M_{j}$ gate and a pivotal rotation.
An algebraic proof of this equivalence starts at Eq. (38) and ends after Eq.
(44). In (b) we show how a pair of CCZ gates, which share a pair of controls,
is Clifford equivalent to only a single CCZ gate. In (c) we uses this identity
twice to simplify a more complex circuit down to a $CCZ_{\\#2}$. In general,
one finds that a sequence of $N$ control-$M_{j}$ can be simplified to a
$CCZ_{\\#N}$ circuit (see Eq. (45) and Eq. (47)).

Next we show that each $U_{j}$ can be implemented with access to a single
$|+\rangle$ ancilla, a control-$M_{j}$ unitary and a $R(2\theta)$ rotation
(see Fig. 3a). We prepare the ancilla in the $|+\rangle$ state and use it as
the control qubit for the control-$M_{j}$ unitary, which gives the state
$|0\rangle|\psi\rangle+|1\rangle(M_{j}|\psi\rangle)$. Next we rotate the
ancilla by $HR(2\theta)H$ and measure in the computational basis. This is
equivalent to measuring with projections

$\displaystyle\langle 0|HR(2\theta)H$ $\displaystyle=\cos(2\theta)\langle
0|+i\sin(2\theta)\langle 1|$ (41) $\displaystyle\langle 1|HR(2\theta)H$
$\displaystyle=\cos(2\theta)\langle 1|+i\sin(2\theta)\langle 0|.$ (42)

In the eventuality of a “+1” outcome, we find

$\displaystyle|0\rangle|\psi\rangle+|1\rangle(M_{j}|\psi\rangle)\rightarrow(\cos(2\theta){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}+i\sin(2\theta)M_{j})|\psi\rangle=U_{j}|\psi\rangle,$ (43)

as desired. However, when a “1” outcome is measured we have

$\displaystyle|0\rangle|\psi\rangle+|1\rangle(M_{j}|\psi\rangle)\rightarrow(i\sin(2\theta){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}+\cos(2\theta)M_{j})|\psi\rangle=M_{j}U_{j}|\psi\rangle.$ (44)

We see that an $M_{j}$ gate will correct for the different measurements
outcomes, and since $M_{j}$ is Clifford this does not contribute to the
resource cost. This completes the proof of the identity in Fig. 3a.

Similarly we chain together $N$ such circuits, assuming we have access to $N$
copies of $R(2\theta)$ and the circuit

$\displaystyle\tilde{V}$ $\displaystyle=\prod_{j=1,\ldots,N}(|0\rangle\langle
0|_{-j}+|1\rangle\langle 1|_{-j}M_{j}),$ (45)
$\displaystyle=\prod_{j=1,\ldots,N}CZ_{-j,2j}CCZ_{-j,0,2j}CCZ_{-j,0,2j-1}.$
(46)

This is composed of a sequence of control-$M_{j}$ gates, where each gate is
controlled from a different ancillary qubit. We have labelled these new
ancilla with negative integers from $-1$ to $-N$. Each control-$M_{j}$
consists of a CZ gate and two CCZ gates. On first inspection this seems to
imply $2N$ CCZ gates are needed. However, using the identity in Fig. 3b we see
a pair of CZZ gates can sometimes be realised using a single CCZ. Note this
identity only works because the pair of CCZ gates share two control bits in
common. Applying this identity repeatedly (see e.g. Fig. 3c) can reduce the
circuit to one using only $N$ CCZ gates. Algebraically, the identity is

$\displaystyle\tilde{V}$
$\displaystyle=\left(\prod_{j=1}^{N}CX_{2j,2j-1}\right)\left(\prod_{j=1}^{N}CCZ_{-j,0,2j-1}\right)\left(\prod_{j=1}^{N}CX_{2j,2j-1}\right)\left(\prod_{j=1}^{N}CZ_{-j,2j}\right),$
(47)

where $CX_{c,t}$ is a control-$X$ gate with control qubit $c$ and target qubit
$t$. The resource intensive part is the non-Clifford component of $N$ CCZ
gates all sharing one single control qubit in common (qubit $0$). Here we
denote such a circuit as CCZ#N, which has been elsewhere called Tof#N. For
these gates, the problem of optimal synthesis into CNOT + $T$ gates has been
solved and the circuit requires $4N+3$ $T$ gates (see Example IV.2. of Ref.
[Campbell and Howard(2017b)]). Recall that we are requiring that CCZ#N is
predistilled to a higher fidelity. The most efficient known method to achieve
this is to use the synthillation protocol that can perpare CCZ#N using only
$(4N+4)$ $T$-states of $\epsilon_{\pi/8}$ error rate, and so this is the
first-step of our two-step protocol.

We have demonstrated how step two works using a series of circuit identities
(for any integer $N$). For completeness, we show in Fig. 4 how these circuit
identities plug together for $N=1$ and $N=2$. The circuits could be further
expanded by replacing the non-Clifford gate $CCZ_{\\#N}$ with the magic state
$|CCZ_{\\#N}\rangle$ and the appropriate Clifford injection circuit.

Figure 4: An overview of the second-step of our protocol for $N=1$ and $N=2$,
which can be extended to any integer $N$. The first-step is to use
synthillation to inject the non-Clifford gates shown in the pink section of
the circuit. The only other non-Clifford elements of the circuits are the
pivotal rotations $R(\theta)$. Input to the circuit are $2N$ mixed states
$\rho(\theta)$, which are $|R(\theta)\rangle$ magic states with $\epsilon$
phase noise. In the absence of noise, the circuit measures the parity in the
$W(\theta)$ bases. Therefore, when we see a SUCCESS event, the $\rho(\theta)$
states are output with quadratically reduced noise. In the event of a FAILURE,
we discard the qubits and attempt again.

##  4 Noise analysis

This section presents a performance analysis for one round of our protocols.
Subsection 4.1 reports some results of numerical simulations for smaller size
protocols with $N=1$ and $N=2$. Subsection 4.2 focuses on providing a simple,
yet rigorous derivation, of analytic upper bounds on output noise. The
analytic results hold for all $N$, but are loose and actual performance will
be much better than analytically bounded.

###  4.1 Numerical analysis

We performed full state vector numerical simulations using IBM’s QISKit (code
available as ancillary file). We simulated the effect of leading order errors
for circuits with $N=1$ and $N=2$. The output error probabilities are for a
single qubit with other output qubits traced out. Numerical results were
independent of which output qubit is chosen and independent of $\theta$ (up-to
numerical accuracy of 2 significant figures).

For $N=1$, we found that

$\displaystyle\epsilon^{\prime}$
$\displaystyle=8\epsilon_{\pi/8}^{2}+\epsilon_{\theta}^{2}+\frac{1}{4}\eta+O(\epsilon_{\pi/8}^{3},\epsilon_{\theta}^{2},\eta^{2},\ldots).$
(48)

The leading order coefficients for output error are identical to those for the
MEKk protocols proposed by Campbell and O’Gorman (themselves a modified form
of MEK) and so it seems that our new protocols (with $N=1$) perform
identically in this regard. We have a slight performance advantage in terms of
success probability due to the two step nature. We found

$\displaystyle p_{\mathrm{synth}}$
$\displaystyle=1-8\epsilon_{\pi/8}+O(\epsilon_{\pi/8}^{2}),$ (49)
$\displaystyle p_{\mathrm{parity}}$
$\displaystyle=1-2\epsilon_{\theta}-\frac{1}{2}\eta+O(\epsilon_{\pi/8}^{2},\epsilon_{\theta}^{2},\eta^{2},\ldots),$

where $p_{\mathrm{synth}}$ is the probability of step one succeeding and
$p_{\mathrm{parity}}$ is the probability of step two succeeding. To leading
order, the previous MEKk protocols had a success probability equal to
$p_{\mathrm{mek}}=p_{\mathrm{synth}}p_{\mathrm{parity}}$. Here, we don’t
commit to the second step until the first step is successful, which will lead
to superior rates of generating magic states. For the setting $\theta=\pi/8$,
the protocol simplifies to a $10\rightarrow 2$ protocol with
$\epsilon^{\prime}=9\epsilon_{\pi/8}^{2}+O(\epsilon_{\pi/8}^{3})$ and overhead
$n/k=10/2=5$, very similar to the original MEK protocol.

For $N=2$, we found that

$\displaystyle\epsilon^{\prime}$
$\displaystyle=16\epsilon_{\pi/8}^{2}+3\epsilon_{\theta}^{2}+\frac{1}{4}\eta+O(\epsilon_{\pi/8}^{3},\epsilon_{\theta}^{2},\eta^{2},\ldots),$
(50) $\displaystyle p_{\mathrm{synth}}$
$\displaystyle=1-12\epsilon_{\pi/8}+O(\epsilon_{\pi/8}^{2}),$ (51)
$\displaystyle p_{\mathrm{parity}}$
$\displaystyle=1-6\epsilon_{\theta}-\eta+O(\epsilon_{\pi/8}^{2},\epsilon_{\theta}^{2},\eta^{2},\ldots),$
(52)

By going to $N=2$, we incur only mildly worse constant prefactors, but gain a
significant efficiency improvement in terms of magic states output per input.
For the setting $\theta=\pi/8$, the protocol simplifies to a $16\rightarrow 4$
protocol with
$\epsilon^{\prime}=19\epsilon_{\pi/8}^{2}+O(\epsilon_{\pi/8}^{3})$ and
overhead $n/k=16/4=4$. To obtain the same $n/k$ overhead using Bravyi-Haah
protocols (which are limited to $\theta=\pi/8$), we need to go to a larger
size $32\rightarrow 8$ protocol with
$\epsilon^{\prime}=25\epsilon_{\pi/8}^{2}+O(\epsilon_{\pi/8}^{3})$. This
confirms that our protocols can obtain similar resource overheads with a
smaller scale quantum computer, and without any sacrifice in terms of error
suppression or success probability.

###  4.2 Algebraic analysis

Here we take an analytic approach. We do not know of an analytic method of
determining the exact expressions for $\epsilon^{\prime}$, but can prove a
rigourous upper bound using standard norm inequalities. Actual performance
will be much better than proven here. We begin by considering the effect of
noise on the input states

$\rho(\theta):=(1-\epsilon_{\theta})|R(\theta)\rangle\langle
R(\theta)|+\epsilon_{\theta}Z|R(\theta)\rangle\langle R(\theta)|Z.$ (53)

We extend later to account for CCZ#N noise and pivotal rotation noise, but
when these components work perfectly the circuit implements

$\rho(\theta)^{\otimes 2N}\rightarrow\mathcal{E}(\rho(\theta)^{\otimes
2N})=K\rho(\theta)^{\otimes 2N}K^{\dagger}$ (54)

where $K$ is the parity projecting Kraus operator introduced in Eq. (3).
Rather than the whole multi-qubit output, we are interested in the fidelity of
a single output qubit, and so introduce the channel

$\mathcal{E}_{i}(\rho(\theta)^{\otimes
2N})=\mathrm{tr}_{i}[K\rho(\theta)^{\otimes 2N}K^{\dagger}],$ (55)

where $\mathrm{tr}_{i}[\cdots]$ is the partial trace over all but the
$i^{\mathrm{th}}$ qubit. The output of this channel is the unnormalised state

$\mathcal{E}_{i}(\rho(\theta)^{\otimes
2N})=p_{\mathrm{good}}|R(\theta)\rangle\langle
R(\theta)|+p_{\mathrm{bad}}Z|R(\theta)\rangle\langle R(\theta)|Z.$ (56)

The renormalisation constant $p_{\mathrm{good}}+p_{\mathrm{bad}}$ is the
probability of the parity check yielding a “+1" outcome. When the parity check
process is error-free, this occurs whenever the input states contains no
errors or an even number of errors, and so

$\displaystyle p_{\mathrm{good}}+p_{\mathrm{bad}}$
$\displaystyle=\frac{1}{2}(1+(1-2\epsilon_{\theta})^{2N})$ (57)

The term $p_{\mathrm{bad}}$ is the probability of an error on
$i^{\mathrm{th}}$ qubit and an odd number of errors on the remaining $2N-1$
qubits

$\displaystyle p_{\mathrm{bad}}$
$\displaystyle=\epsilon_{\theta}\sum_{j=1}^{N}\genfrac{(}{)}{0.0pt}{1}{2N-1}{2j-1}\epsilon_{\theta}^{2j-1}(1-\epsilon_{\theta})^{2N-2j},$
(58)
$\displaystyle=\frac{\epsilon_{\theta}}{2}(1-(1-2\epsilon_{\theta})^{2N-1}),$
$\displaystyle<(2N-1)\epsilon_{\theta}^{2},$

where in the first line we have a binomial coefficient. The inequality follows
from Bernoulli’s inequality. This shows quadratic noise suppression in
$\epsilon_{\theta}$.

Next, we account for $\epsilon_{\\#}$ noise in the CCZ#N gate. We can write
the corresponding magic state as

$\rho_{\\#N}=(1-\epsilon_{\\#})\Psi_{\\#N}+\epsilon_{\\#}\sigma_{\\#N},$ (59)

where

$\Psi_{\\#N}=|CCZ_{\\#N}\rangle\langle CCZ_{\\#N}|,$ (60)

and $\sigma_{\\#N}$ carries some $Z$ noise. We define $\mathcal{F}$ as the
channel describing the action of the whole circuit (including implicit
injection gadget for $CCZ_{\\#N}$), assuming ideal pivotal rotations, acting
on $\rho_{\\#N}$ and $\rho(\theta)^{\otimes 2N}$. We will use that
$\Psi_{\\#N}$ leads to a parity

$\mathcal{F}(\Psi_{\\#N}\otimes\rho(\theta)^{\otimes
2N})=\mathcal{E}(\rho(\theta)^{\otimes 2N})$ (61)

and by linearity of $\mathcal{F}$ we deduce

$\displaystyle\mathcal{F}(\rho_{\\#N}\otimes\rho(\theta)^{\otimes
2N})=(1-\epsilon_{\\#})\mathcal{E}(\rho(\theta)^{\otimes
2N})+\epsilon_{\\#}\mathcal{F}(\sigma_{\\#N}\otimes\rho(\theta)^{\otimes
2N}).$ (62)

Again, we are interested in only the single output qubit, and so introduce
$\mathcal{F}_{i}=\mathrm{tr}_{i}\mathcal{E}_{i}$, which straightforwardly
yields

$\mathcal{F}_{i}(\rho_{\\#N}\otimes\rho(\theta)^{\otimes
2N})=(1-\epsilon_{\\#})\mathcal{E}_{i}(\rho(\theta)^{\otimes
2N})+\epsilon_{\\#}\mathcal{F}_{i}(\sigma_{\\#N}\otimes\rho(\theta)^{\otimes
2N}).$ (63)

This yields a single qubit state of the form in Eq. (56) with new parameters
$p^{\prime}_{\mathrm{good}}$ and $p^{\prime}_{\mathrm{bad}}$, which are tricky
to exactly calculate but can again be bounded. The joint probability
$p^{\prime}_{\mathrm{good}}+p^{\prime}_{\mathrm{bad}}$ can be lower bounded by
assuming $\mathcal{F}_{i}(\sigma_{\\#N}\otimes\rho(\theta)^{\otimes 2N})=0$
and so

$\displaystyle p^{\prime}_{\mathrm{good}}+p^{\prime}_{\mathrm{bad}}$
$\displaystyle\geq(1-\epsilon_{\\#})(p_{\mathrm{good}}+p_{\mathrm{bad}})$ (64)

The $p^{\prime}_{\mathrm{bad}}$ term can be upper bounded by considering the
worst-case scenario that
$\mathcal{F}_{i}(\sigma_{\\#N}\otimes\rho(\theta)^{\otimes 2N})$ leads to a
logical error with unit probability, and so

$\displaystyle p^{\prime}_{\mathrm{bad}}$
$\displaystyle\leq(1-\epsilon_{\\#})p_{\mathrm{bad}}+\epsilon_{\\#}<(2N-1)(1-\epsilon_{\\#})\epsilon_{\theta}^{2}+\epsilon_{\\#},$
(65)

where the second inequality follows from Eq. (58). These bounds are very loose
and overestimate $p^{\prime}_{\mathrm{bad}}$ by quite a lot. Nevertheless,
they are simple to obtain and rigorous.

Next, we further consider phase noise on the pivotal rotation, each failing
with probability $\eta$. In other words, all pivotal rotations act perfectly
with probability $(1-\eta)^{N}$. Therefore, the channel implemented is not
$\mathcal{F}_{i}$ but something of the form

$\mathcal{G}_{i}=(1-\eta)^{N}\mathcal{F}_{i}+(1-(1-\eta)^{N})\mathcal{F}^{\prime}_{i},$
(66)

where $\mathcal{F}^{\prime}_{i}$ is the noisy part of the channel with diamond
norm not exceeding unity. Therefore,

$\displaystyle\mathcal{G}_{i}(\rho_{\\#N}\otimes\rho(\theta)^{\otimes 2N})=$
$\displaystyle(1-\eta)^{N}(p^{\prime}_{\mathrm{good}}|R(\theta)\rangle\langle
R(\theta)|+p^{\prime}_{\mathrm{bad}}Z|R(\theta)\rangle\langle R(\theta)|Z)$
(67)
$\displaystyle+(1-(1-\eta)^{N})\mathcal{F}^{\prime}_{i}(\rho_{\\#N}\otimes\rho(\theta)^{\otimes
2N}).$

The worst case scenario is that $\mathcal{F}^{\prime}_{i}$ always generates an
error, adding a $(1-(1-\eta)^{N})$ contribution to the error term. Therefore,
after renormalising the error probability is bounded by

$\displaystyle\epsilon^{\prime}$
$\displaystyle\leq\frac{(1-\eta)^{N}p^{\prime}_{\mathrm{bad}}+(1-(1-\eta)^{N})}{(1-\eta)^{N}(p_{\mathrm{good}}^{\prime}+p_{\mathrm{bad}}^{\prime})+(1-(1-\eta)^{N})}$
$\displaystyle\leq\frac{(1-\eta)^{N}((2N-1)(1-\epsilon_{\\#})\epsilon_{\theta}^{2}+\epsilon_{\\#})+(1-(1-\eta)^{N})}{(1-\eta)^{N}(1-\epsilon_{\\#})(p_{\mathrm{good}}+p_{\mathrm{bad}})+(1-(1-\eta)^{N})}$

The result scales as $O(\epsilon_{\\#})$, but this error rate is itself the
output of performing the synthillation protocol using noisy
$|R(\pi/8)\rangle$-states of error rate $\epsilon_{\pi/8}$. In particular, Eq.
(128) of Ref. [Campbell and Howard(2017b)] shows that

$\epsilon_{\\#}\leq
1-\frac{2(1-\epsilon_{\pi/8})^{4N+4}}{1+(1-2\epsilon_{\pi/8})^{4N+4}}~=(6+14N+8N^{2})\epsilon_{\pi/8}^{2}+O(\epsilon_{\pi/8}^{3})$
(68)

This suffices to conclude that

$\epsilon^{\prime}\leq(2N-1)\epsilon_{\theta}^{2}+(6+14N+8N^{2})\epsilon_{\pi/8}^{2}+N\eta+O(\epsilon_{\theta}^{3},\epsilon_{\pi/8}^{3},\eta^{2},\ldots)$
(69)

where $O(\epsilon^{3})$ collects all higher order terms. For instance, for
$N=1$ and $N=2$ this yields

$\displaystyle\epsilon^{\prime}$
$\displaystyle\leq\begin{cases}\epsilon_{\theta}^{2}+28\epsilon_{\pi/8}^{2}+\eta+O(\epsilon_{\theta}^{3},\epsilon_{\pi/8}^{3},\eta^{2},\ldots)&\mbox{
for }N=1\\\
3\epsilon_{\theta}^{2}+66\epsilon_{\pi/8}^{2}+2\eta+O(\epsilon_{\theta}^{3},\epsilon_{\pi/8}^{3},\eta^{2},\ldots)&\mbox{
for }N=2\\\ \end{cases}$ (70)

Comparing this with the numerical expressions Eq. (48) and Eq. (50), we see
the analytic upper bound is very loose and grossly overestimates the
prefactors.

##  5 No small triorthogonal codes

Many other distillation protocols are based upon projections into codespaces
with a transversal non-Clifford, with transversality proofs typically using
some notion of triorthogonal matrices [Bravyi and Haah(2012)]. While the
protocols proposed here do not manifestly have this form, it is natural to ask
whether there is some codespace projection with equivalent performance.
Indeed, Jones’ first-level distiller protocol [Jones(2013a)] is effectively
equivalent to projecting onto the codespace of Bravyi-Haah triorthogonal codes
[Bravyi and Haah(2012)], and Haah has recently introduced level-lifting as a
general methodology for finding such equivalences [Haah(2017)]. Furthermore,
it has long been known that for any distillation protocol there exists a
codespace projection that achieves the same error suppression [Campbell and
Browne(2009)], though it may not achieve the same success probability or admit
a transversal non-Clifford gate.

In this section, we show that there exist no triorthogonal codes with fewer
than 14 qubits. This bound is tight since the smallest Bravyi-Haah code is a
14 qubit triorthogonal construction. It follows that the $10\rightarrow 2$ MEK
protocol is not equivalent to a projection onto a triorthogonal code. What
distinguishes MEK is that it is a highly compressed circuit that is obtained
from taking a larger circuit and cancelling some $T$-gates. This suggests that
something happens during the compression process of eliminating extraneous
$T$-gates that breaks the equivalence to triorthogonal codes. Since our
protocols can be understood as a generalisation of MEK protocols, it seems
unlikely similar performance parameters will be achievable using projections
onto codes with exotic transversality properties.

We present the definition of triorthogonality

###### Definition 1

(Def 1. of Ref. [Bravyi and Haah(2012)])A binary matrix $G$ of size $m\times
n$ is called triorthogonal iff the supports of any pair and any triple of its
rows have even overlap, that is,

$\sum_{j=1}^{n}G_{a,j}G_{b,j}=0\pmod{2}$ (71)

for all pairs of rows $1\leq a<b\leq m$ and

$\sum_{j=1}^{n}G_{a,j}G_{b,j}G_{c,j}=0\pmod{2}$ (72)

for all triples of rows $1\leq a<b<c\leq m$.

The definition of triorthogonality allows a matrix to have either odd or even
rows, and it is standard to use a horizontal line to demarcate the split

$G=\left(\frac{G_{1}}{G_{0}}\right),$ (73)

so $G_{1}$ contains odd weight rows and $G_{0}$ contains even weight rows.
Assuming $G$ is row-wise linearly independent, it describes an $[[n,k,d]]$
quantum code where: $n$ is the number of columns in $G$; $k$ is the number of
rows in $G_{1}$; and $d\geq 2$ if and only if $G_{0}$ is non-trivially
supported on every column. We also use the notion of a biorthogonal matrix,
which obeys the constraint for pairs of rows but not for triples of rows.

Let $G$ be a triorthogonal matrix with block matrix form

$G=\left(\frac{G_{1}}{G_{0}}\right)=\left(\begin{array}[]{cc}A&B\\\ \hline\cr
C&D\\\ \mathbf{1}&\mathbf{0}\end{array}\right),$ (74)

where $\mathbf{1}$ and $\mathbf{0}$ are the all-1 and all-0 row vectors of
appropriate width. Using column permutation the matrix can always be brought
into this form. Without loss of generality, we assume that the last row has
weight $w$ and is the highest weight row in the span of $G_{0}$. Let $B$ and
$D$ be width $u$, so that the total matrix width is $n=w+u$.

From triorthogonality of $G$, it follows that the submatrix

$L=\left(\begin{array}[]{c}A\\\ C\end{array}\right),$ (75)

is biorthogonal with all rows being even weight and that

$R=\left(\begin{array}[]{c}B\\\ \hline\cr D\end{array}\right),$ (76)

is biorthogonal with $B$ containing odd weight rows. Since $B$ contains odd
weight rows, the matrix $D$ cannot contain the all-1 vector as this would
violate biorthogonality. However, since the code is distance 2, the matrix $D$
must be supported on every column. Therefore, there must exist at least 2 non-
trivial rows in $D$. The smallest possible width for $D$ is then achieved by

$D=\left(\begin{array}[]{cccccc}1&1&1&1&0&0\\\ 1&1&0&0&1&1\end{array}\right),$
(77)

which has width $u=6$ and contains weight 4 vectors, and so we can infer that
$w\geq 4$. Other $D$ are possible, but one cannot obtain smaller parameters:
There are at least two rows and for any pair of rows they must overlap on at
least 2 columns and also have support on at least 2 other non-overlapping
columns. From this we see that $n=w+u\geq 4+6=10$. So this is already enough
to prove there are no $[[n,k,2]]$ triorthogonal codes with $k\leq 1$ and $n$
less than 10.

However, the bound $w\geq 4$ was obtained based on the rows of $D$ only, but
$w$ is the max weight across all rows in the span of $G_{0}$. If $C$ contains
a row of weight $w_{c}$, then we know that $w\geq w_{c}+4$. For any row of $C$
we can add the last row of $G_{0}$, which generates a row of weight
$w_{c}^{\prime}=w-w_{c}$, entailing that $w\geq w_{c}^{\prime}+4=w-w_{c}+4$
and so $w_{c}\geq 4$. Putting this together yields $w\geq 8$ and $u\geq 6$ so
that $n\geq 14$. There are no triorthogonal codes with fewer than 14 qubits.

##  6 Discussion

###  6.1 Variation of the two-step protocol

This subsection discusses one possible variant of the two-step protocol.
Consider a quantum algorithm that needs many magic states of the form
$|R(\theta)\rangle$, but with different values of $\theta$. As presented, our
main protocol cannot be used to full effect as the very large $N$ limit
assumes that we need many magic states with the same $\theta$. However, it is
straightforward to check that one can distill pairs of states with the same
$\theta_{j}$. That is, we may input states of the form
$\bigotimes_{j=1}^{N}\rho(\theta_{j})^{\otimes 2}$ and use $N$ pivotal
rotations with corresponding angles $2\theta_{j}$.

###  6.2 Quadratic vs higher order error suppression

In our introduction, we remarked on the existence of distillation routines
that offer much larger reductions in errors without using concatenation
[Jones(2013a), Haah et al.(2017)Haah, Hastings, Poulin, and Wecker, Haah et
al.(2018)Haah, Hastings, Poulin, and Wecker, Hastings and Haah(2018)]. The
appeal of these protocols is better asymptotic performance in the limit of
large quantum computers and large error reduction. The analysis underpinning
these results assumes that it is appropriate to quantify resource costs by the
ratio of input to output states. But a more realistic picture is given by an
involved full space time analysis; also accounting for the cost of Clifford
gates and quantum error correction. In such an analysis, it is possible to
scale the size of error correction codes between rounds of magic state
distillation [Raussendorf et al.(2006)Raussendorf, Harrington, and Goyal,
Fowler et al.(2012)Fowler, Mariantoni, Martinis, and Cleland, O’Gorman and
Campbell(2017)]. This scaling trick is extremely effective, and is arguably
the most important tool in the arsenal of magic state distillation techniques.
Although the idea has been known for some time, it has gone without a name. In
an effort to popularise this trick, O’Gorman and Campbell recently proposed
the phrase “balanced investment” [O’Gorman and Campbell(2017)].

Balanced investment relies on distinct rounds of magic state distillation with
successive error reduction. Therefore, balanced investment is more compatible
with protocols giving quadratic error reduction, such as presented here, than
with the protocols of Refs. [Jones(2013a), Haah et al.(2017)Haah, Hastings,
Poulin, and Wecker]. This argument is qualitative, and we need detailed
resource analyses to make concrete quantitative statements. However, such full
resource investigations are difficult and time-consuming and have only been
undertaken for the Reed-Muller and Bravyi-Haah protocols. Naturally, such a
numerical investigation also falls outside the scope of this paper.

##  7 Conclusions

We presented a new two-step method of magic state distillation that is very
competitive at preparing single-qubit magic states, offering a way to
circumvent the need for costly gate-synthesis of single-qubit rotations. An
important aspect of these new protocols are the preparation of multi-qubit
magic states using synthillation. Pressing open questions include how these
competing approach fare when all resource costs are considered, though such an
analysis will depend heavily on the architecture considered.

We would also like to explore whether the synthillation driven techniques
proposed here could be extended to protocols with larger than quadratic error
reduction [Jones(2013a), Haah et al.(2017)Haah, Hastings, Poulin, and Wecker].
After completing this work, Hastings and Haah proposed some new approaches to
synthillation that may provide a starting point for attacking this problem
[Haah and Hastings(2017)].

##  8 Acknowledgements

This research was supported by the EPSRC (grant ref EP/M024261/1). We thank
IBM and developers of the QISKit, which was used for numerical simulations.

## References

  * [Campbell et al.(2017)Campbell, Terhal, and Vuillot] Earl T Campbell, Barbara M Terhal, and Christophe Vuillot.  Roads towards fault-tolerant universal quantum computation.  _Nature_ , 549(7671):172, 2017.  doi: 10.1038/nature23460.
  * [Bravyi and Kitaev(2005)] Sergey Bravyi and Alexei Kitaev.  Universal quantum computation with ideal Clifford gates and noisy ancillas.  _Phys. Rev. A_ , 71:022316, 2005.  doi: 10.1103/PhysRevA.71.022316.
  * [Meier et al.(2013)Meier, Eastin, and Knill] Adam M. Meier, Bryan Eastin, and Emanuel Knill.  Magic-state distillation with the four-qubit code.  _Quant. Inf. and Comp._ , 13:195, 2013.
  * [Bravyi and Haah(2012)] Sergey Bravyi and Jeongwan Haah.  Magic-state distillation with low overhead.  _Phys. Rev. A_ , 86:052329, Nov 2012.  doi: 10.1103/PhysRevA.86.052329.
  * [Jones(2013a)] Cody Jones.  Multilevel distillation of magic states for quantum computing.  _Phys. Rev. A_ , 87:042305, Apr 2013a.  doi: 10.1103/PhysRevA.87.042305.
  * [Haah et al.(2017)Haah, Hastings, Poulin, and Wecker] Jeongwan Haah, Matthew B. Hastings, D. Poulin, and D. Wecker.  Magic state distillation with low space overhead and optimal asymptotic input count.  _Quantum_ , 1:31, October 2017.  ISSN 2521-327X.  doi: 10.22331/q-2017-10-03-31.
  * [Kliuchnikov et al.(2013)Kliuchnikov, Maslov, and Mosca] Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca.  Asymptotically optimal approximation of single qubit unitaries by clifford and $t$ circuits using a constant number of ancillary qubits.  _Phys. Rev. Lett._ , 110:190502, May 2013.  doi: 10.1103/PhysRevLett.110.190502.
  * [Gosset et al.(2014)Gosset, Kliuchnikov, Mosca, and Russo] David Gosset, Vadym Kliuchnikov, Michele Mosca, and Vincent Russo.  An algorithm for the t-count.  _Quant. Inf. & Comp._, 14(15-16):1261–1276, 2014.
  * [Ross and Selinger(2016)] Neil J Ross and Peter Selinger.  Optimal ancilla-free clifford+ t approximation of z-rotations.  _Quant. Inf. and Comp._ , 16:901, 2016.
  * [Paetznick and Svore(2014)] Adam Paetznick and Krysta M Svore.  Repeat-until-success: Non-deterministic decomposition of single-qubit unitaries.  _Quant. Inf. & Comp._, 14(15-16):1277–1301, 2014.
  * [Bocharov et al.(2015)Bocharov, Roetteler, and Svore] Alex Bocharov, Martin Roetteler, and Krysta M. Svore.  Efficient synthesis of probabilistic quantum circuits with fallback.  _Phys. Rev. A_ , 91:052317, May 2015.  doi: 10.1103/PhysRevA.91.052317.
  * [Landahl and Cesare(2013)] Andrew J Landahl and Chris Cesare.  Complex instruction set computing architecture for performing accurate quantum $z$ rotations with less magic.  _arXiv preprint arXiv:1302.3240_ , 2013.  URL https://arxiv.org/pdf/1302.3240.pdf.
  * [Jones(2013b)] Cody Jones.  Distillation protocols for fourier states in quantum computing.  _arXiv preprint arXiv:1303.3066_ , 2013b.  URL https://arxiv.org/pdf/1303.3066.pdf.
  * [Duclos-Cianci and Poulin(2015)] Guillaume Duclos-Cianci and David Poulin.  Reducing the quantum-computing overhead with complex gate distillation.  _Phys. Rev. A_ , 91:042315, Apr 2015.  doi: 10.1103/PhysRevA.91.042315.
  * [Campbell and O’Gorman(2016)] Earl T Campbell and Joe O’Gorman.  An efficient magic state approach to small angle rotations.  _Quantum Science and Technology_ , 1(1):015007, 2016.  doi: doi:10.1088/2058-9565/1/1/015007.
  * [Haah(2017)] Jeongwan Haah.  Towers of generalized divisible quantum codes.  _arXiv preprint arXiv:1709.08658_ , 2017.  URL https://arxiv.org/pdf/1709.08658.pdf.
  * [Jones(2013c)] Cody Jones.  Low-overhead constructions for the fault-tolerant toffoli gate.  _Phys. Rev. A_ , 87:022328, Feb 2013c.  doi: 10.1103/PhysRevA.87.022328.
  * [Eastin(2013)] Bryan Eastin.  Distilling one-qubit magic states into toffoli states.  _Phys. Rev. A_ , 87:032321, Mar 2013.  doi: 10.1103/PhysRevA.87.032321.
  * [Campbell and Howard(2017a)] Earl T. Campbell and Mark Howard.  Unifying gate synthesis and magic state distillation.  _Phys. Rev. Lett._ , 118:060501, Feb 2017a.  doi: 10.1103/PhysRevLett.118.060501.
  * [Campbell and Howard(2017b)] Earl T. Campbell and Mark Howard.  Unified framework for magic state distillation and multiqubit gate synthesis with reduced resource cost.  _Phys. Rev. A_ , 95:022316, Feb 2017b.  doi: 10.1103/PhysRevA.95.022316.
  * [Haah et al.(2018)Haah, Hastings, Poulin, and Wecker] Jeongwan Haah, Matthew B Hastings, D Poulin, and D Wecker.  Magic state distillation at intermediate size.  _Quant. Inf. and Comp._ , 18:0114, 2018.
  * [Hastings and Haah(2018)] Matthew B. Hastings and Jeongwan Haah.  Distillation with sublogarithmic overhead.  _Phys. Rev. Lett._ , 120:050504, Jan 2018.  doi: 10.1103/PhysRevLett.120.050504.
  * [Gottesman and Chuang(1999)] Daniel Gottesman and Isaac L. Chuang.  Demonstrating the viability of universal quantum computation using teleportation and single-qubit operations.  _Nature_ , 402:390, 1999.  doi: 10.1038/46503.
  * [Campbell and Browne(2009)] Earl T. Campbell and Dan E. Browne.  On the structure of protocols for magic state distillation.  _Lecture Notes in Computer Science_ , 5906:20, 2009.  doi: 10.1007/978-3-642-10698-9_3.  arXiv:0908.0838.
  * [Raussendorf et al.(2006)Raussendorf, Harrington, and Goyal] R. Raussendorf, J. Harrington, and K. Goyal.  A fault-tolerant one-way quantum computer.  _Annals of Physics_ , 321(9):2242 – 2270, 2006\.  ISSN 0003-4916.  doi: 10.1016/j.aop.2006.01.012.
  * [Fowler et al.(2012)Fowler, Mariantoni, Martinis, and Cleland] Austin G. Fowler, Matteo Mariantoni, John M. Martinis, and Andrew N. Cleland.  Surface codes: Towards practical large-scale quantum computation.  _Phys. Rev. A_ , 86:032324, Sep 2012.  doi: 10.1103/PhysRevA.86.032324.
  * [O’Gorman and Campbell(2017)] Joe O’Gorman and Earl T. Campbell.  Quantum computation with realistic magic-state factories.  _Phys. Rev. A_ , 95:032338, Mar 2017.  doi: 10.1103/PhysRevA.95.032338.
  * [Haah and Hastings(2017)] Jeongwan Haah and Matthew B Hastings.  Codes and protocols for distilling $t$, controlled-$s$, and toffoli gates.  _arXiv preprint arXiv:1709.02832_ , 2017.  URL https://arxiv.org/pdf/1709.02832.pdf.


# Deep Neural Nets with Fixed
Bias Configuration

Harbir Antil1, Thomas S. Brown1,2, Rainald Löhner2, Fumiya Togashi3, Deepanshu
Verma4 1Center for Mathematics and Artificial Intelligence (CMAI), College of
Science, George Mason University, Fairfax, VA 22030-4444, USA.
hantil<EMAIL_ADDRESS>2Center for Computational Fluid Dynamics,
College of Science, George Mason University, Fairfax, VA 22030-4444, USA.
<EMAIL_ADDRESS>3Applied Simulations, Inc., 1211 Pine Hill Road, McLean, VA
22101, USA<EMAIL_ADDRESS>4Department of Mathematics, Emory
University, Atlanta, GA 30322, USA<EMAIL_ADDRESS>
###### Abstract.

For any given neural network architecture a permutation of weights and biases
results in the same functional network. This implies that optimization
algorithms used to ‘train’ or ‘learn’ the network are faced with a very large
number (in the millions even for small networks) of equivalent optimal
solutions in the parameter space. To the best of our knowledge, this
observation is absent in the literature. In order to narrow down the parameter
search space, a novel technique is introduced in order to fix the bias vector
configurations to be monotonically increasing. This is achieved by augmenting
a typical learning problem with inequality constraints on the bias vectors in
each layer. A Moreau-Yosida regularization based algorithm is proposed to
handle these inequality constraints and a theoretical convergence of this
algorithm is established. Applications of the proposed approach to standard
trigonometric functions and more challenging stiff ordinary differential
equations arising in chemically reacting flows clearly illustrate the benefits
of the proposed approach. Further application of the approach on the MNIST
dataset within TensorFlow, illustrate that the presented approach can be
incorporated in any of the existing machine learning libraries.

This work is partially supported by the Defense Threat Reduction Agency (DTRA)
under contract HDTRA1-15-1-0068 where Jacqueline Bell served as the technical
monitor, and by NSF grants DMS-2110263, DMS-1913004, and the Air Force Office
of Scientific Research under Award NO: FA9550-19-1-0036.

##  1\. Introduction

Background. A typical neural network can be represented as a function
$\mathcal{F}:\mathbb{R}^{n_{0}}\to\mathbb{R}^{n_{L}}$ that consists of the
composition of layer functions $\\{f_{\ell}\\}_{\ell=0}^{L-1}$ and can written
as

$\mathcal{F}=f_{L-1}\circ f_{L-2}\circ\cdots\circ f_{0}\,.$ (1)

Each layer function is parameterized by a weight matrix
$W_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell+1}}$, a bias vector
$b_{\ell}\in\mathbb{R}^{n_{\ell+1}}$, and incorporates a nonlinear activation
function $\sigma$, for instance, ReLU [5].

The weights $W_{\ell}$ and biases $b_{\ell}$ are determined in a process known
as training the network. Training a neural network can be written in the
framework of constrained optimization as follows: for training data
$\\{u^{i},S(u^{i})\\}_{i=1}^{N}$ (input/output pairs), solve

$\displaystyle\min_{\\{W_{\ell}\\}_{\ell=0}^{L-1},\\{b_{\ell}\\}_{\ell=0}^{L-2}}J(\\{(y_{L}^{i},S(u^{i}))\\}_{i},\\{W_{\ell}\\}_{\ell},\\{b_{\ell}\\}_{\ell})$
(2) $\displaystyle\mbox{subject to }\quad
y_{L}^{i}=\mathcal{F}(u^{i};(\\{W_{\ell}\\},\\{b_{\ell}\\}))\qquad
i=1,\dots,N,$

where the function $J$, known as the loss function, measures the error of the
approximation of $S(u^{i})$ by the network output $y_{L}^{i}$ in some way.

Figure 1. A neural network with 2 hidden layers of width 10. Switching the
bias values and edges connected to node A (in red) with the those of node B
(in magenta) results in a different ordering of the parameters, but the same
neural network function.

It is easy to see that there are many permutations of the parameters that
result in the same network. This is illustrated in Figure 1 where a network is
represented graphically with the nodes of the graph representing the biases
and the edges of the graph representing the weights. Switching the bias values
A (red) and B (magenta) results in the same network as long as the
corresponding weights (also in red and magenta) are also switched. This
implies that it is possible to change the order of the parameters without
changing the action of the function $\mathcal{F}$. In fact, for networks where
each layer consists of $n_{\ell}$ neurons, there are
$\sum_{\ell=1}^{L-2}n_{\ell}!$ ways to rearrange the parameters and obtain the
same network. For the relatively small network displayed in Figure 1 there are
7,257,600 different ways to permute the parameters and still obtain the same
network. In the context of training, this means that for this particular
example there are 7,257,600 different solutions to the optimization problem
(2) that will result in the same neural network function $\mathcal{F}$. This
level of non-uniqueness is troubling and highly unsatisfactory. Apparently,
this simple looking observation has not received any attention so far.

Problem formulation. The parameter search space may be narrowed down by fixing
the bias vector configurations to be monotonically increasing. This is
achieved by augmenting the learning problem (2) with inequality constraints on
the bias vectors in each layer. The resulting optimization problem that
describes the training of the network is

$\displaystyle\min_{\\{W_{\ell}\\}_{\ell=0}^{L-1},\\{b_{\ell}\\}_{\ell=0}^{L-2}}$
$\displaystyle
J(\\{(y_{L}^{i},S(u^{i}))\\}_{i},\\{W_{\ell}\\}_{\ell},\\{b_{\ell}\\}_{\ell})$
(3a) subject to $\displaystyle
y_{L}^{i}=\mathcal{F}(u^{i};(\\{W_{\ell}\\},\\{b_{\ell}\\}))\,,\qquad
i=1,\dots,N,$ (3b) $\displaystyle b_{\ell}^{j}\leq b_{\ell}^{j+1}\,,\qquad
j=1,\dots,n_{\ell+1}-1,\quad\ell=0,\dots,L-2,$ (3c)

where for a fixed $\ell$, the quantities $b_{\ell}^{j}$ are the entries of the
bias vector $b_{\ell}$. The added inequality constraints fix the configuration
of the network, so that any permutation of the parameters (excluding the case
when two adjacent bias values are equal) either violates the constraints or
results in a different neural network.

Outline of the paper. Section 2 introduces a Moreau-Yosida regularization
based algorithm to handle inequality constraints given in (3c). This is
followed by a convergence result of the Moreau-Yosida regularized problem to
the original problem in Section 3. Section 4 describes the ResNets used for
the examples shown in Section 5. In the first example the proposed algorithm
is applied to a standard trigonometric function. This is followed by an
application to a realistic application in chemically reacting flows, which are
governed by stiff ordinary differential equations. The final example on the
MNIST dataset within TensorFlow, illustrates that the presented approach can
be incorporated in any of the existing machine learning libraries.

##  2\. Regularized problem

The inequality constraints in (3c) are difficult to implement directly.
However, it is possible to implement them implicitly by augmenting the loss
function $J$ with a penalty term similar to a Moreau-Yosida regularization.
The well-known Moreau-Yosida regularization is frequently used to implement
inequality constraints in the context of optimization problems with partial
differential equations as constraints, see [9, 12, 8, 1]. Before introducing
the loss function with bias order regularization, a more precise definition of
the loss function is given. For training data
$\\{u^{i},S(u^{i})\\}_{i=1}^{N}$, consider

$J:=\frac{1}{2N}\sum_{i=1}^{N}\|y_{L}^{i}-S(u^{i})\|_{2}^{2}+\frac{\lambda}{2}\sum_{\ell=0}^{L-1}\big{(}\|W_{\ell}\|_{1}+\|b_{\ell}\|_{1}+\|W_{\ell}\|_{2}^{2}+\|b_{\ell}\|_{2}^{2}\big{)},$
(4)

where $\lambda>0$ is a regularization parameter. In order to fit $J$ into the
framework introduced above, $b_{L-1}$ is taken to be the zero vector in
$\mathbb{R}^{n_{L}}$. Even though a mean squared error term is used above to
measure the approximation error of the neural network, this is easily
generalizable to other terms such as cross-entropy, likelihood, etc. [5].

Using the notation above to represent the entries of each bias vector, namely
$b_{\ell}=(b_{\ell}^{j})_{j=1}^{n_{\ell+1}}\in\mathbb{R}^{n_{\ell+1}}$, the
new loss function is defined as

$J_{\gamma}:=J+\frac{\gamma}{2}\sum_{\ell=0}^{L-2}\sum_{j=1}^{n_{\ell+1}-1}\|\min\\{b_{\ell}^{j+1}-b_{\ell}^{j},0\\}\|_{2}^{2},$
(5)

where $\gamma$ is the so-called penalization parameter. The regularized
optimization problem can now be written as

$\displaystyle\min_{\\{W_{\ell}\\}_{\ell=0}^{L-1},\\{b_{\ell}\\}_{\ell=0}^{L-2}}J_{\gamma}(\\{(y_{L}^{i},S(u^{i}))\\}_{i},\\{W_{\ell}\\}_{\ell},\\{b_{\ell}\\}_{\ell})$
(6a) $\displaystyle\mbox{subject to }\quad
y_{L}^{i}=\mathcal{F}(u^{i};(\\{W_{\ell}\\},\\{b_{\ell}\\}))\qquad
i=1,\dots,N.$ (6b)

Note that, even though it is not explicitly written in the formulation above,
all of the variables $W_{\ell}$, $b_{\ell}$, and $y^{i}_{L}$ depend on the
parameter $\gamma$.

In Appendix A the first order optimality conditions for this problem are
derived where the DNN used is a Deep Residual Neural Net (ResNet).

##  3\. Convergence of $J_{\gamma}$

In order to show that as $\gamma\to\infty$, the minimum value of $J_{\gamma}$
converges to the minimum value of $J$ and the constraints (3c) are also
satisfied, let $\theta$ represent the concatenation of all of the parameters
which are being minimized, i.e. $\theta$ contains all of the entries of
$\\{W_{\ell}\\}$ and $\\{b_{\ell}\\}$. Furthermore, assume that a fixed set of
training data is being used and so the loss function $J$ given in (4) can be
rewritten as

$J(\theta):=\frac{1}{2N}\sum_{i=1}^{N}\|\mathcal{F}(u^{i};\theta)-S(u^{i})\|_{2}^{2}+\frac{\lambda}{2}\sum_{\ell=0}^{L-1}\big{(}\|\theta\|_{1}+\|\theta\|_{2}^{2}\big{)}.$

Introducing the notation

$g(\theta):=\sum_{\ell=0}^{L-2}\sum_{j=1}^{n_{\ell+1}-1}\|\min(b_{\ell}^{j+1}-b_{\ell}^{j},0)\|_{2}^{2}\,,$

the regularized loss function in (5) may be rewritten as

$J_{\gamma}(\theta^{\gamma}):=J(\theta^{\gamma})+\frac{\gamma}{2}g(\theta^{\gamma}).$

Now, the framework of [11, Section 10.11] can be used to show the following
convergence results. The proof follows exactly as in [11] after a
transformation of notation.

######  Proposition 3.1.

Let $J_{0}$ be the minimum value attained from solving (3). For each $\gamma$,
let $\theta^{\gamma}$ be a minimizer of $J_{\gamma}$. The following hold

  1. (a)

$J_{\gamma}(\theta^{\gamma})\geq J_{\tilde{\gamma}}(\theta^{\tilde{\gamma}})$,
for $\gamma\geq\tilde{\gamma}$;

  2. (b)

$J_{0}\geq J_{\gamma}(\theta^{\gamma})$ for each $\gamma$;

  3. (c)

$\lim_{\gamma\to\infty}\frac{\gamma}{2}g(\theta^{\gamma})=0$.

In particular, part $(c)$ shows that as $\gamma\to\infty$ the inequalities in
(3c) are satisfied.

######  Remark 3.2.

Typically, when a Moreau-Yosida regularization is implemented, a path-
following technique is used to increase the size of $\gamma$ gradually. This
means that a sequence of optimization problems is solved for subsequently
larger values of $\gamma$. The initial $\gamma$ value is taken to be small,
and the solution to the problem, $\theta^{\gamma}$, is used as the initial
guess for the next optimization problem with a larger value of $\gamma$. This
path-following process continues until $\gamma$ is sufficiently large. For the
numerical examples presented below, a path-following technique was not used,
and it was sufficient (as the results show) to solve each problem for a single
value of $\gamma$.

##  4\. ResNets

As some of the examples shown below are obtained for Deep Residual Neural
Networks (ResNets), a small description follows. Recall the definition of
$\mathcal{F}$ in (1). In the sequel, the layer functions will be denoted as
$f_{\ell}=f_{\ell}(y_{\ell};(W_{\ell},b_{\ell})):\mathbb{R}^{n_{\ell}}\to\mathbb{R}^{n_{\ell+1}}$,
where the dependence of $f_{\ell}$ on $\sigma$ is not explicitly written. With
this representation the neural network can be viewed as an iterative
progression of updating the output of each layer as

$y_{\ell+1}=f_{\ell}(y_{\ell};(W_{\ell},b_{\ell}))\qquad\ell=0,\dots,L-1,$

with initial input $y_{0}$ and final output $y_{L}$. As before, in order to
preserve consistency, $b_{L-1}$ is taken to be the zero vector. One way to
define the layer functions is

$\displaystyle f_{0}(y_{0})$ $\displaystyle:=\sigma(W_{0}y_{0}+b_{0}),$
$\displaystyle f_{\ell}(y_{\ell})$
$\displaystyle:=P_{\ell}y_{\ell}+\tau\sigma(W_{\ell}y_{\ell}+b_{\ell})\qquad\ell=1,\dots,L-2,$
$\displaystyle f_{L-1}$ $\displaystyle:=W_{L-1}y_{L-1},$

for matrices $P_{\ell}\in\mathbb{R}^{n_{\ell}\times n_{\ell+1}}$ and a scalar
$\tau>0$. With this definition and for $L>2$, $\mathcal{F}$ is termed as a
Deep Residual Neural Network (ResNet). If each $P_{\ell}$ is taken to be the
identity matrix, which requires that the hidden layers have a uniform width,
then this network can be viewed as a forward Euler discretization of an ODE.
For more on these kinds of networks see [7, 14, 3, 2, 6], among others. Notice
that if $\tau=1$ and $\\{P_{\ell}\\}_{\ell=0}^{L-1}$ contains only zero
entries, then $\mathcal{F}$ is a standard feedforward deep neural network [5].

##  5\. Numerical Results

In this section several examples are given that not only show the efficacy of
this method, but also the advantages of using the method. In Section 5.1, a
first example compares a single ResNet to learn the function $\sin(x)$ with
and without bias ordering. This simple example shows that the method performs
as desired, and in fact outperforms the same network trained without bias
ordering. Following this, Section 5.2, a more complicated experiment is
reported that use parallel ResNets to learn a model related to chemically
reacting flows [4]. This example shows that the proposed method is a useful
technique for practical problems in machine learning. Finally, in in Section
5.3, an example is described where the bias ordering regularization is applied
to a classification problem using MNIST data and implemented in Keras. This
example shows that the proposed method is also suitable for Convolutional
Neural Networks.

In Sections 5.1 and 5.2 a BFGS optimization routine with Armijo line search is
used during training to solve the optimization problem with or without bias
ordering. In order to avoid overfitting, validation data is used with a
patience of 400 iterations. This means that the training data is separated
into two sets: training and validation data. The validation data is not used
to update the weights and biases, rather it is used to measure the error of
the network on unseen data (the validation data) during training. If the
validation error increases, the patience iterations provide a buffer during
which training continues. If during these iterations, the validation error
reaches a new minimum, then training continues as before, otherwise the
training routine is terminated.

###  5.1. ResNet to learn $\sin(x)$

For the first example a simple ResNet with 2 hidden layers of width 50 is used
to learn the function $\sin(x)$. The skip parameter $\tau$ is taken to be 1.
For data, 1000 evenly spaced points from the interval $[0,2\pi]$ are generated
and then randomly split into training data (400 points), testing data (400
points), and validation data (200 points).

The same experiment was performed twice, once using the loss function $J$, and
once with the loss function $J_{\gamma}$ with $\gamma=100$. The resulting
trained networks are shown in Figure 2, where the bias values are represented
by the neurons of the network, and the weights are represented by the
connections between neurons. The color for the input and output are set to
zero. In the case where the network was trained with loss function
$J_{\gamma}$, the resulting biases were perfectly ordered and so the
inequality constraints were satisfied. This can also be seen in the bottom
panel in Figure 2. Note that for the network with unconstrained bias values,
there are on the order of $10^{64}$ permutations of the parameters that will
give the same ResNet approximation. For the ResNet with ordered bias values,
however, any permutation of the parameters will either result in a different
ResNet approximation, or violate the bias ordering.

Figure 2.  Network visualizations for the ResNets trained with loss function
$J$ (top) and $J_{\gamma}$ (bottom). The network consists of 2 hidden layers
with 50 neurons in each layer to learn $\sin(x)$. It is clear from the panels
that the ordering of the bias values is only enforced in the bottom panel.

In Figure 3 the output of the two networks is compared on the 400 test points,
with exact values in blue and ResNet output in red. The left plot of Figure 3
shows the results for a network trained with a standard loss function $J$,
while the plot on the right shows results from a network trained with the
augmented loss function $J_{\gamma}$ which orders the biases. It is evident
from the plots that the ResNet that implements bias order produced more
accurate results. This is quantified above the plots with the relative error
measured in the 2-norm.

Figure 3. Results from training ResNets to learn the function $\sin(x)$. Exact
values in blue, ResNet output in red. The left panel shows results for the
standard ResNet with $J$ and the right panel shows results where bias ordering
is enforced using $J_{\gamma}$. Both visual inspection and quantitative
inspection of the error confirms that the proposed approach works better in
this example.

###  5.2. Applications to chemically reacting flows

In this section experiments are presented involving a ResNet approximation of
a system of stiff ODEs that model a reduced H2-O2 reaction. This reduced model
(see [13]) tracks 8 species and temperature as they interact over time and is
completely separated from any advection and diffusion in space. For more
information on this problem and more experiments using a parallel ResNet
approximation see [4].

For the experiments included in this work, nine parallel ResNets with input
dimension 10 and output dimension 1 (9 total output quantities) and 8 hidden
layers of width 30 were trained on H2-O2 reaction data created by solving the
stiff ODE system using CHEMKIN [10]. Given an input vector representing the
data at time $t_{k}$, each ResNet is trained to learn a single quantity
(temperature, for example) at time $t_{k+1}$. More details can be found in
[4]. The generated data used to train and test these networks consists of
thirteen subsets corresponding to initial conditions with a fixed equivalence
ratio of 1 and different temperatures varying from 1200K to 2400K in
increments of 100K. The parallel ResNets were trained on the data sets
corresponding to initial temperatures 1200K, 1500K, 1800K, 2100K, and 2400K.

The experiment above was performed twice, once with loss functions $J$ (one
for each parallel ResNet) as described in Section 2, and once with regularized
loss functions $J_{\gamma}$ with $\gamma=1000$. All other network
hyperparameters including the initial values of the weights and biases prior
to training are kept the same for the two experiments. In Table 1 the number
of BFGS iterations used in training the parallel networks for both experiments
described above is compared. The network with Moreau-Yosida regularization to
order the bias values used fewer iterations to train in five of the nine
networks.

Summary of BFGS iterations used during training
---
| with loss functions $J$ | with loss functions $J_{\gamma}$
Network 1 (temperature) | 2638 | 786
Network 2 (O) | 2506 | 5125
Network 3 (H) | 2772 | 1123
Network 4 (OH) | 5957 | 4444
Network 5 (HO2) | 230 | 540
Network 6 (H2O2) | 886 | 1738
Network 7 (H2O) | 7267 | 583
Network 8 (O2) | 8872 | 2547
Network 9 (H2) | 3719 | 4570
Table 1. A comparision of the number of BFGS iterations used during training
for the two sets of parallel ResNets. The quantity that the ResNet is learning
(output) is written in parentheses.

Recall from Proposition 3.1, that the Moreau-Yosida regularization approach
will be enforcing the ordering (3c) approximately.

The loss functions $J_{\gamma}$ penalize the violation of the ordering, but do
not strictly implement the order itself. Even so, the biases in eight of the
nine parallel networks were ordered perfectly for the networks that trained
with the loss functions $J_{\gamma}$. The only violation of the bias ordering
occured in Network 9, the network used to learn H2. In the first hidden layer
of this network neurons six and seven were ordered incorrectly. In the trained
network, the value of neuron 6 was approximately -0.0175652 and the value of
neuron 7 was approximately -0.0175660, and therefore the size of the order
violation was -8e-7, which is negligible. To reiterate, in each of the 9
networks, there were 240 bias values (30 per hidden layer) for a total of
2,160 different biases. In all of these biases only a single pair (negligibly)
violated the monotonic ordering using the proposed method.

Figure 4 compares the results of the two sets of parallel ResNets tested with
initial conditions with initial temperature 1400K (left set of plots) and
2000K (right set of plots). Note that the networks were not trained on this
data. To test the networks only the initial condition comes from the CHEMKIN
data. The output of the ResNets from the initial condition are then combined
and used as the input to the parallel ResNets for the next timestep. This
process is repeated for the duration of the reaction. To compare the results,
the known CHEMKIN data are represented in blue, the results from the networks
trained with $J$ are represented with dashed red lines and the results for the
networks trained with $J_{\gamma}$ are represented with dash/dotted black
lines. It can clearly be seen that the results of the parallel ResNets trained
with $J_{\gamma}$ match the data more closely. Furthermore, note that the
$x-$axis (time axis) is scaled logarithmically in order to display the details
of H2O2 which happen quickly and early in the reaction. Therefore, while the
results in black anticipate the reaction, they only do so slightly, on the
order of $10^{-6}$ seconds.

Figure 4. Results from training parallel ResNets on a reduced H2-O2 reaction
model. Values from known data in blue, results from ResNets trained with
traditional loss functions in red, and results from ResNets trained with loss
functions with a Moreau-Yosida penalization term to order the biases in black.
Clearly, the proposed approach outperforms the existing one.

###  5.3. A Convolutional Neural Network classification problem

For this example, a convolutional neural network is constructed using Keras to
solve a classification problem using the MNIST data set. The purpose of this
example is to show that the flexibility of the bias ordering technique and
that it can be incorporated into existing neural network software. Unlike the
other examples presented, the neural network is not a ResNet. Instead the
network consists of 4 convolutional layers followed by two dense layers. The
diagram in Figure 1 can also be used to describe a convolutional layer, where
the nodes of the graph still represent the bias value, but now the edges
represent a convolutional filter or convolutional kernel. In this way, the
bias ordering technique easily applies to these networks as well.

The network used for this example consists of four convolutional layers
followed by two dense layers. The first two convolutional layers consist of 32
convolutional filters of width 3. Convolutional layers three and four both
consist of 64 convolutional filters of width 3. The first dense layer consists
of 512 neurons, and the second dense layer (the output layer) consists of 10
neurons since there are 10 digits (or classes) in the MNIST data. Unlike the
other examples, bias values are included on the output layer of this network.
A softmax activation function is used for the output layer, while a ReLU
activation function is used for the other five layers. Max pooling is used
after layers 2 and 4. In order to avoid overfitting, dropout is used in the
network.

Two versions of the network described above are constructed using Keras with a
Tensorflow backend. Both of the networks use a categorical cross entropy loss
function, and the difference between the two is that one of the networks also
implements the Moreau-Yosida regularization term through a custom bias
regularizer. The regularizer is defined using the following lines in Python

⬇

def MY_regularizer(bias, gamma=100):

bias_length = tf.size(bias)

bias_diff = bias[1:bias_length]-bias[0:bias_length-1]

bias_min = tf.math.minimum(bias_diff,0)

bias_norm = tf.math.reduce_sum(tf.math.square(bias_min))

return 0.5*gamma*bias_norm

In the above code , the value of the regularization parameter is taken to be
100, but this is easily customized. Once the regularizer is defined, it can be
implemented into any layer with the bias_regularizer option.

Both of the networks were trained on 60,000 samples from the MNIST data set.
The loss functions were minimized using stochastic gradient descent with a
learning rate of 0.01. Batch normalization was used during training with a
batch size of 32. Each network trained for 5 epochs before the networks were
tested on 10,000 MNIST samples. In Figure 5 the bias values for each layer of
the two networks are plotted. It can clearly be seen that the regularization
resulted in successfully ordering the bias values in each layer. The accuracy
of the two trained networks were comparable. For the networks corresponding to
the plots in Figure 5, the network with the bias regularization had an
accuracy of 98.7% and the network without regularization had an accuracy of
98.66%.

Figure 5. The bias values by layer for the two neural networks described in
Section 5.3. These plots show the successful implementation of the Moreau-
Yosida bias order regularization in Keras with a Tensorflow backend.

##  6\. Conclusions

A method to reduce the very large search space of equivalently optimal neural
nets has been introduced. The key idea is to try to enforce the biases to be
monotonically increasing in each layer of neurons. This is accomplished by a
Moreau-Yosida regularization-based approach to solve the resulting
optimization problem. The convergence of the regularized problem has been
proven. The benefit of this new approach has been demonstrated not only on
simple approximation cases but also on a realistic problem arising in
chemically reacting flows. The numerical experiments presented also show that
the regularization technique is effective for relatively small values of the
penalization parameter $\gamma$.

##  Appendix A Derivation of the first order optimality conditions

In the context of constrained optimization, the problem (6) is typically
solved by using a gradient based method. Indeed this is the approach taken in
the numerical experiments. Each evaluation of the gradient requires solving
the state equation or forward problem and solving the adjoint equation
problem. In what follows, the state and adjoint equations will be derived for
the problem in (6) where the neural network being used is a Deep Residual
Neural Network (ResNet), as introduced in Section 4.

The state and adjoint equations as well as the gradient will be derived by
using the Lagrangian approach. For brevity, the following will be written for
a single input $u$, rather than for a set of training data. Similar to Section
3, $\theta$ is used to represent the concatenation of the weights and bias,
i.e., the parameters being optimized. For appropriate adjoint variables
$\psi=(\psi_{j})_{j=1}^{L}$ the Lagrangian functional corresponding to (6) is

$\mathcal{L}(u,\theta,\psi):=J_{\gamma}(\theta)+\langle
y_{1}-\tau\sigma(W_{0}u+b_{0}),\psi_{1}\rangle\\\
+\sum_{j=2}^{L-1}\Bigg{\langle}y_{j}-P_{j-1}y_{j-1}-\tau\sigma(W_{j-1}y_{j-1}+b_{j-1}),\psi_{j}\Bigg{\rangle}+\langle
y_{L}-W_{L-1}y_{L-1},\psi_{L}\rangle,$

From here, the state and adjoint equations result from evaluating the
derivatives of the Lagrangian with respect to $y_{j}$ and $\psi_{j}$ at a
stationary point. Furthermore, the gradient is derived by taking the
derivatives of the Lagrangian with respect to $\theta$.

  1. (i)

State Equation.

$\displaystyle\begin{aligned} y_{1}&=\tau\sigma(W_{0}u+b_{0}),\\\
y_{j}&=P_{j-1}y_{j-1}+\tau\sigma(W_{j-1}y_{j-1}+b_{j-1}),\qquad 2\leq j\leq
L-1,\\\ y_{L}&=W_{L-1}y_{L-1}.\end{aligned}$ (7a)

  2. (ii)

Adjoint Equation.

$\displaystyle\psi_{j}$
$\displaystyle=P_{j}^{T}\psi_{j+1}-\tau\left[-W_{j}^{T}\left(\psi_{j+1}\odot\sigma^{\prime}\left(W_{j}y_{j+1}+b_{j}\right)\right)\right]\qquad
j=L-2,\ldots,1$ (7b) $\displaystyle\psi_{L-1}$
$\displaystyle=-W_{L-1}^{T}\psi_{L},$ $\displaystyle\psi_{L}$
$\displaystyle=-\partial_{y_{L}}J_{\gamma}(\theta).$

  3. (iii)

Derivative with respect to $\theta$.

$\displaystyle\partial_{W_{L-1}}\mathcal{L}=$
$\displaystyle-\psi_{L}\;y_{L-1}^{T}+\partial_{W_{L-1}}J_{\gamma}(\theta)$
(7c) $\displaystyle=$
$\displaystyle\partial_{y_{L}}J_{\gamma}(\theta)\;y_{L-1}^{T}+\partial_{W_{L-1}}J_{\gamma}(\theta),$
$\displaystyle\partial_{W_{j}}\mathcal{L}=$ $\displaystyle-
y_{j}\>\left(\psi_{j+1}\odot\sigma^{\prime}(W_{j}y_{j}+b_{j})\right)^{T}+\partial_{W_{j}}J_{\gamma}(\theta)\qquad$
$\displaystyle j=0,...,L-2,$ $\displaystyle\partial_{b_{j}}\mathcal{L}=$
$\displaystyle-\;\psi_{j+1}^{T}\;\sigma^{\prime}(W_{j}y_{j}+b_{j})+\partial_{b_{j}}J_{\gamma}(\theta)\qquad$
$\displaystyle j=0,...,L-2\,.$

In the jargon of machine learning the equations in (7b) are called back
propagation. The gradient is represented by the right hand side of (7c), where
the contributions from the Moreau-Yosida regularization terms enter in the
term $\partial_{b_{j}}J_{\gamma}(\theta)$.

## References

  * [1] Harbir Antil, Thomas S. Brown, Deepanshu Verma, and Mahamadi Warma.  Optimal control of fractional PDEs with state and control constraints.  To appear in Pure and Applied Functional Analysis, 2021.
  * [2] Harbir Antil, Howard C Elman, Akwum Onwunta, and Deepanshu Verma.  Novel deep neural networks for solving bayesian statistical inverse.  arXiv preprint arXiv:2102.03974, 2021.
  * [3] Harbir Antil, Ratna Khatri, Rainald L Lohner, and Deepanshu Verma.  Fractional deep neural network via constrained optimization.  Machine Learning: Science and Technology, 2020.
  * [4] Thomas S Brown, Harbir Antil, Rainald Löhner, Fumiya Togashi, and Deepanshu Verma.  Novel DNNs for stiff ODEs with applications to chemically reacting flows.  arXiv preprint arXiv:2104.01914, 2021.
  * [5] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2016\.
  * [6] Stefanie Günther, Lars Ruthotto, Jacob B. Schroder, Eric C. Cyr, and Nicolas R. Gauger.  Layer-parallel training of deep residual neural networks.  SIAM J. Math. Data Sci., 2(1):1–23, 2020.
  * [7] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 770–778, 2016.
  * [8] Michael Hintermüller and Michael Hinze.  Moreau-Yosida regularization in state constrained elliptic control problems: error estimates and parameter adjustment.  SIAM J. Numer. Anal., 47(3):1666–1683, 2009.
  * [9] Kazufumi Ito and Karl Kunisch.  Lagrange multiplier approach to variational problems and applications, volume 15 of Advances in Design and Control.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008.
  * [10] R. J. Kee, F. M. Rupley, J.A. Miller, M.E. Coltrin, J.F. Grcar, E. Meeks, H.K. Moffat, A.E. Lutz, G. Dixon-Lewis, M.D. Smooke, J. Warnatz, G.H. Evans, R. S. Larson, R.E. Mitchell, L.R. Petzold, W.C. Reynolds, M. Caracotsios, W.E. Stewart, P. Glarborg, C. Wang, and O. Adigun.  Chemkin collection, release 3.6, 2000.
  * [11] David G. Luenberger.  Optimization by vector space methods.  John Wiley & Sons, Inc., New York-London-Sydney, 1969.
  * [12] Ira Neitzel and Fredi Tröltzsch.  On regularization methods for the numerical solution of parabolic control problems with pointwise state constraints.  ESAIM Control Optim. Calc. Var., 15(2):426–453, 2009.
  * [13] Eric L. Petersen and Ronald K. Hanson.  Reduced kinetics mechanisms for ram accelerator combustion.  Journal of Propulsion and Power, 15(4):591–600, 1999.
  * [14] Lars Ruthotto and Eldad Haber.  Deep neural networks motivated by partial differential equations.  J. Math. Imaging Vision, 62(3):352–364, 2020.