# A Different Perspective
On The Stochastic Convex Feasibility Problem

James Renegar and Song Zhou  School of Operations Research and Information
Engineering, Cornell University, U.S.

###### Abstract.

We investigate the stochastic convex feasibility problem, the problem of
finding $x$ satisfying a system of constraints $f_{\omega}(x)\leq 0$
($\omega\in\Omega$), where the functions $f_{\omega}$ are convex, and where
$\Omega$ is the sample space in a probability triple
$(\Omega,{\mathcal{F}},P)$. Our aim is to solve the problem approximately, in
the sense that given values $\epsilon>0$ and $0<\Gamma<1$, a point $x$ is
found for which the probability measure $P$ satisfies $P\\{\omega\in\Omega\mid
f_{\omega}(x)\leq\epsilon\\}\geq 1-\Gamma$. Our primary vehicle for obtaining
such a point $x$ is a subgradient method taking “Polyak steps”, for which the
only input is the size of minibatches to be employed in the effort. We prove
stochastic iteration bounds for the general setting, and improved bounds when
a notion of “Hölderian growth” is satisfied, showing, for example, linear
convergence when the growth is linear (as happens, say, for a finite system of
linear inequalities). We also derive deterministic iteration bounds for
knowing “with high confidence” that an iterate $x_{k}$ is in hand which
fulfills the desired inequality, $P\\{\omega\in\Omega\mid
f_{\omega}(x_{k})\leq\epsilon\\}\geq 1-\Gamma$.

Research supported NSF grant DMS-1812904.

##  1\. Introduction

The convex feasibility problem – the problem of finding a solution to a system
of convex constraints – is an integral component of convex optimization, and
from some perspectives the two are equivalent, in the deterministic setting.
In the stochastic setting the feasibility problem is less studied than the
problem that often is referred to as “stochastic optimization”, the problem of
minimizing a function of the form
$\mathbb{E}_{\omega\in\Omega}[f_{\omega}(x)]$, where $\Omega$ is the sample
space in a probability triple $(\Omega,{\mathcal{F}},P)$.

For us the stochastic convex feasibility problem is that of finding, or
approximating, a point $x$ satisfying a system of functional constraints
$f_{\omega}(x)\leq 0$ ($\omega\in\Omega$) where the functions $f_{\omega}$ are
convex.111A close cousin is the problem of finding, or approximating, a point
in the intersection of closed convex sets, $\cap_{\omega\in\Omega}X_{\omega}$,
for which (deterministic and randomized) projection methods have been
extensively studied (see Necoara, Richtárik and Patrascu [4] for discussion
that fits well with our paper). The intersection problem can be
(superficially) viewed as a special case of the functional constraint problem,
by defining $f_{\omega}(x)=\mathrm{dist}(x,X_{\omega})$ (Euclidean distance
from $x$ to $X_{\omega}$) – then the full gradient step at $x$ corresponds to
projecting $x$ onto $X_{\omega}$. However, this viewpoint does not begin to
capture the rich body of geometric ideas that have arisen and matured in the
literature on projection methods, and so the intersection problem is best
considered as a cousin rather than a special case of the functional constraint
problem. Beginning with Polyak [5] (and perhaps most recently, Necoara and
Nedić [3]), the prominent approach for designing iterative methods is to view
the feasibility problem, under mild assumptions, as being equivalent to the
stochastic optimization problem in which the goal is to minimize
$\mathbb{E}[f_{\omega}^{+}(x)]$ (where
$f_{\omega}^{+}(x):=\max\\{0,f_{\omega_{\ell}}(x)\\}$). Then if an algorithm
employs minibatches in making iterations, it is natural to take a
(sub)gradient step at the current iterate $x_{k}$ for the sample-expectation
function $x\mapsto\frac{1}{L}\sum_{\ell=1}^{L}f_{\omega_{\ell}}^{+}(x)$, where
$\omega_{1},\ldots,\omega_{L}$ are independently drawn from
$(\Omega,{\mathcal{F}},P)$.222The size of minibatches, $L$, can be chosen
“large but not too large”, so as to fully yet efficiently utilize available
parallelism in making the step. Ideally, the number of iterations is reduced
by a factor of $1/L$. A manner in which algorithmic complexity is measured is
by the expected number of iterations sufficient to reach an iterate $x_{k}$
for which $\mathbb{E}[f_{\omega}^{+}(x_{k})]\leq\epsilon$, where $\epsilon>0$
is a parameter (perhaps chosen by the user).

We take a different route. For us the minibatch problem is to take a
subgradient step at the current iterate $x_{k}$ for the sample-maximum
function $x\mapsto\max_{\omega_{1},\ldots,\omega_{L}}f_{\omega_{\ell}}(x)$.
Our manner of measuring complexity is by the expected number of steps
sufficient to reach an iterate $x_{k}$ for which “most” of the functions
satisfy $f_{\omega}(x)\leq\epsilon$, in the sense that $P\\{\omega\mid
f_{\omega}(x_{k})\leq\epsilon\\}\geq 1-\Gamma$, where $0<\Gamma<1$ is a
parameter signifying the proportion of functions that can be ignored. (The
only input to the main algorithm is $L$, the size of minibatches, not
$\Gamma$, $\epsilon$, or any other parameter.)

This alternative approach to the stochastic convex feasibility problem results
in surprisingly tractable analyses. For example, bounding the expected number
of iterations for the main algorithm is done by computing the hitting time for
a compound Bernoulli process, elementary probability. Besides bounding the
expected number of iterations quite generally, we easily establish
concentration bounds on the number of iterations, and more interesting, we
establish much-improved bounds when a notion akin to “Hölderian growth” is
satisfied, for all degrees of growth, not just the linear growth of piecewise-
linear convex functions or the quadratic growth of strongly convex functions.
Finally, we establish the analogous results under a slight modification to the
algorithm which results in the user knowing with “high confidence” an iterate
$x_{k}$ is in hand that satisfies the desired inequality, $P\\{\omega\mid
f_{\omega}(x_{k})\leq\epsilon\\}\geq 1-\Gamma$. Perhaps surprisingly, the
iteration bounds here are deterministic – all of the probability gets wrapped
into the confidence level (albeit at the expense of potentially large
minibatches).

In the following three sections, a progression of three theorems is
established. Immediately following the statement of a theorem, we provide
discussion, including when appropriate, relations to the existing literature.

##  2\. Basics

For a convex function $h$ which attains its minimum value $h^{*}$, the
subgradient method with “Polyak step size” is the iterative algorithm given by
$x_{+}=x-\frac{h(x)-h^{*}}{\|g\|^{2}}g$ where $g\in\partial h(x)$
(subdifferential) and $\|\,\,\|$ is the Euclidean norm. Under mild conditions,
the iterates converge to the set of minimizers, $\\{x\mid h(x)=h^{*}\\}$.
Likewise, for any constant $c\geq h^{*}$, the iterative algorithm
$x_{+}=x-\frac{\max\\{0,h(x)-c\\}}{\|g\|^{2}}g$ converges to $\\{x\mid
h(x)\leq c\\}$. We also refer to the steps made by this algorithm as being
Polyak steps.

Let $X_{\Omega}:=\\{x\mid f_{\omega}(x)\leq 0\textrm{ for all $\omega$}\\}$,
which we assume is nonempty. Input to the following algorithm is $x_{0}$, the
initial iterate, and $L$, the size of minibatches. Let $z_{0}$ denote the
point in $X_{\Omega}$ which is closest to $x_{0}$. We assume for all $\omega$
that $f_{\omega}$ is $M$-Lipschitz333$|f_{\omega}(x)-f_{\omega}(y)|\leq
M\|x-y\|$. on an open neighborhood of the closed ball
${\mathcal{B}}:=\\{x\mid\|x-z_{0}\|\leq\|x_{0}-z_{0}\|\\}$, a consequence
being that if $x\in{\mathcal{B}}$ and $g_{\omega}\in\partial f_{\omega}(x)$,
then $\|g_{\omega}\|\leq M$.

We are especially interested in the case where $\Omega$ is finite but contains
a vast number of elements, and the case that $\Omega$ is infinite. We assume
samples are drawn independently from $\Omega$ with replacement.

$\mathtt{PolyakFM}$ (Polyak Feasibility Method)
0) Input: $x_{0}$ (initial iterate), and integer $L\geq 1$ (size of minibatch)
Initialization: $k=1$
1) Independently select samples $\omega_{k,1},\ldots,\omega_{k,L}$ from
$(\Omega,{\mathcal{F}},P)$.
2) Determine $\ell_{k}:=\arg\max_{\ell=1,\ldots,L}f_{k,\ell}(x_{k-1})$, then
let $\epsilon_{k-1}:=f_{k,\ell_{k}}(x_{k-1})$ and compute

$x_{k}=\begin{cases}x_{k-1}&\textrm{if $\epsilon_{k-1}\leq 0$ }\\\
x_{k-1}-\frac{\epsilon_{k-1}}{\|g_{\omega_{k,\ell_{k}}}\|^{2}}\,g_{\omega_{k,\ell_{k}}}&\textrm{otherwise}\;,\end{cases}$

where $g_{\omega_{k,\ell_{k}}}\negthinspace\negthinspace\in\partial
f_{\omega_{k,\ell_{k}}}(x_{k-1})$.
3) Let $k+1\rightarrow k$ and return to Step 1.

To ease notation, let $\bar{\omega}_{k}:=(\omega_{k,1},\ldots,\omega_{k,L})$.
So as to view the algorithm as being entirely determined by $x_{0}$ and the
sequence of sample vectors $\bar{\omega}_{1},\bar{\omega}_{2},\ldots$, let us
think of $g_{\omega}=g_{\omega}(x)$ as being a fixed element of $\partial
f_{\omega}(x)$, that is, $x\mapsto g_{\omega}(x)$ is a subgradient selection.
This causes no loss of generality in that if
$g_{\omega_{k,\ell_{k}}}(x_{k-1})$ is needed by the algorithm (i.e., if
$f_{\omega_{k,\ell}}(x_{k-1})>0$ for some $\ell$), the subgradient can be
chosen on the fly, then used and immediately wiped from memory, because as
will become clear, the point $x_{k-1}$ will never again be visited by the
algorithm, and hence $g_{\omega_{k,\ell_{k}}}(x_{k-1})$ will not be needed a
second time. Likewise if index $\ell_{k}$ is not uniquely determined, then
choose any one of the candidate indices, use it, and discard from memory.
Henceforth, we consider the algorithm to be entirely determined by $x_{0}$ and
whatever sample vectors are drawn, $\bar{\omega}_{1},\bar{\omega}_{2},\ldots$.

For points $x$ and scalars $\epsilon>0$, let

$\Omega(x,\epsilon):=\\{\omega\in\Omega\mid f_{\omega}(x)\leq\epsilon\\}\;.$

Given $\epsilon$ and $0<\Gamma<1$, the goal is to compute a point $x$ for
which $P(\Omega(x,\epsilon))\geq 1-\Gamma$. (Clearly, the algorithm itself is
independent of $\epsilon$ and $\Gamma$.)

In the following theorem we assume
$0<\epsilon<M\mathrm{dist}(x_{0},X_{\Omega})$ (Euclidean distance from $x_{0}$
to $X_{\Omega}$), because if $\epsilon\geq M\mathrm{dist}(x_{0},X_{\Omega})$
then $\Omega(x_{0},\epsilon)=\Omega$, so $x_{0}$ achieves the goal of
obtaining an iterate $x_{k}$ satisfying $P(\Omega(x_{k},\epsilon))\geq
1-\Gamma$.

######  Theorem 2.1.

For any scalars $0<\Gamma<1$ and
$0<\epsilon<M\,\mathrm{dist}(x_{0},X_{\Omega})$, the algorithm
$\mathtt{PolyakFM}$ reaches an iterate $x_{k}$ satisfying
$P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$ within an expected number of
iterations not exceeding

$\mathbb{E}={\scriptstyle{\textstyle{\frac{1}{p}}}}\left({\scriptstyle{\textstyle{\frac{M\,\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)^{2}\quad\textrm{where
}p=1-(1-\Gamma)^{L}\;.$

Moreover, we have the concentration inequalities

$k\geq 2\mathbb{E}\Rightarrow\mathrm{Prob}\big{(}\,x_{k}\,\textrm{is first
iterate satisfying }P(\Omega(x_{k},\epsilon))\geq
1-\Gamma\,\big{)}\,\leq\,{\scriptstyle{\textstyle{\frac{1}{2}}}}\bigg{(}\frac{1}{1+\frac{1}{2}\frac{p}{1-p}}\bigg{)}^{k-\lceil
2\,\mathbb{E}\rceil}$ (2.1)

Remarks:

  1. (1)

For integers $L$ satisfying $1\leq L\leq 1/\Gamma$, it is not difficult to
show
$\frac{1}{L\Gamma}\leq\frac{1}{p}<{\scriptstyle{\textstyle{\frac{2}{L\Gamma}}}}$,
and thus the theorem indicates that the expected number of iterations scales
almost ideally with the size of minibatches.

  2. (2)

If $\Omega$ is finite, say $\Omega=\\{1,\ldots,m\\}$, with $P(i)=\frac{1}{m}$
for all $i$, then the proof of the theorem can readily be modified to allow
the minibatch samples to be drawn from $\Omega$ without replacement (assuming
$1\leq L\leq m$). In this case the expected number of iterations scales
perfectly, and choosing $L=m$ results in a deterministic algorithm and
iteration bound. However, if $m$ is astronomical in size, choosing $L\approx
m$ is computationally untenable.

  3. (3)

In general it is impossible to know with certainty whether
$P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$. Nonetheless, “high confidence”
embellishments to the algorithm can be developed, the subject of §4. There,
the values $\epsilon_{k}$ computed by $\mathtt{PolyakFM}$ will play the
fundamental role of being inferred values for which
$P(\Omega(x_{k},\epsilon_{k}))\geq 1-\Gamma$ with high probability (but this
depends on choosing $L$ appropriately large, unlike the present section where
$L$ is allowed to be any positive integer).

  4. (4)

As for a number of papers on the feasibility problem, we can slightly
generalize the setting by assuming the region of interest is $X_{\Omega}\cap
Y$, where $Y$ is a closed convex set for which, given $x$, it is easy to
compute the orthogonal projection $\Pi_{Y}(x)$, the closest point in $Y$ to
$x$. The appropriate modification for $\mathtt{PolyakFM}$ is simply to project
onto $Y$ the point computed in Step 2, that is, replace $x_{k}$ by
$\Pi_{Y}(x_{k})$. As anyone familiar with analyses of subgradient methods can
easily verify, all of our results remain valid for this slightly more general
setting. Since our notation already is cluttered, we choose to leave $Y$ out
of the picture.

The following simple proposition plays a significant role in the subgradient
literature. The proposition allows us to reduce the proof of Theorem 2.1 to
establishing the expected time for a compound Bernoulli process to reach a
specific value, that is, the proof reduces to establishing a hitting time.

######  Proposition 2.2.

Assume $h$ is a convex function, and assume $x,z$ are points satisfying
$h(x)>0$ and $h(z)\leq 0$. For $g\in\partial h(x)$, the Polyak step

$x_{+}=x-{\scriptstyle{\textstyle{\frac{h(x)}{\|g\|^{2}}}}}\,g\;$

satisfies

$\|x_{+}-z\|^{2}\leq\|x-z\|^{2}-({\scriptstyle{\textstyle{\frac{h(x)}{\|g\|}}}})^{2}\;.$
(2.2)

Proof: We have

$\displaystyle\|x_{+}-z\|^{2}$
$\displaystyle=\|x-{\scriptstyle{\textstyle{\frac{h(x)}{\|g\|^{2}}}}}g-z\|^{2}$
$\displaystyle=\|x-z\|^{2}+2{\scriptstyle{\textstyle{\frac{h(x)}{\|g\|^{2}}}}}\,\langle
g,z-x\rangle+({\scriptstyle{\textstyle{\frac{h(x)}{\|g\|}}}})^{2}$
$\displaystyle\leq\|x-z\|^{2}+2{\scriptstyle{\textstyle{\frac{h(x)}{\|g\|^{2}}}}}(h(z)-h(x))+({\scriptstyle{\textstyle{\frac{h(x)}{\|g\|}}}})^{2}$
(because $g$ is a subgradient at $x$)
$\displaystyle\leq\|x-z\|^{2}-({\scriptstyle{\textstyle{\frac{h(x)}{\|g\|}}}})^{2}$
$\displaystyle\qquad\textrm{(using $h(z)-h(x)\leq-h(x)$)}\;,$

completing the proof. $\Box$

Remark: The proof is trivially extended to allow extrapolated subgradient
steps, $x_{+}=x=\delta{\scriptstyle{\textstyle{\frac{h(x)}{\|g\|^{2}}}}}g$ for
$0<\delta<2$, the importance of which is that choosing $\delta>1$ often
results empirically in improved performance, and often results in obtaining a
point in the interior of $X_{\Omega}$ if the interior is not empty (we do not
assume the interior is nonempty). For extrapolated steps, the inequality (2.2)
becomes
$\|x_{+}-z\|^{2}\leq\|x-z\|^{2}-\delta(2-\delta)({\scriptstyle{\textstyle{\frac{h(x)}{\|g\|}}}})^{2}$.
This can be used exactly as we rely on (2.2) throughout the paper. We choose
not to include this layer of detail so as to avoid the additional notation.

An easy corollary to the proposition is that for $\mathtt{PolyakFM}$, we have
$\|x_{k}-z_{0}\|<\|x_{k-1}-z_{0}\|$ unless $x_{k}=x_{k-1}$. In particular, all
iterates lie in the ball ${\mathcal{B}}$.

Also an important consequence of the proposition is that the iterates of
$\mathtt{PolyakFM}$ satisfy

$\mathrm{dist}(x_{k},X_{\Omega})^{2}\leq\mathrm{dist}(x_{k-1},X_{\Omega})^{2}-(\epsilon_{k-1}/M)^{2}\quad\textrm{if
$\epsilon_{k-1}\geq 0$}.$ (2.3)

Proof of Theorem 2.1: We begin by inductively defining a variable which keeps
track of whether a subgradient iterate $x_{k}$ has been reached for which the
desired inequality $P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$ holds – in this
case we rely on the symbol “$\bullet$” – and which also keeps track, when such
an iterate has yet to be reached, of the number of iterations $k$ for which
$\epsilon_{k-1}\geq\epsilon$. In this it is useful to define
$j+\bullet=\bullet$ for integers $j$.

Let

$X_{0}=\begin{cases}\bullet&\textrm{if $P(\Omega(x_{0},\epsilon))\geq
1-\Gamma$}\\\ 0&\textrm{otherwise}\;,\end{cases}$

and for sequences of vectors $\bar{\omega}_{1},\bar{\omega}_{2},\ldots$ (where
$\bar{\omega}_{k}=(\omega_{k,1},\ldots,\omega_{k,L})$), inductively define

$X_{k}(\bar{\omega}_{1},\ldots,\bar{\omega}_{k})=\begin{cases}\bullet&\textrm{if
$P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$}\\\ X_{k-1}+1&\textrm{if
$P(\Omega(x_{k},\epsilon))<1-\Gamma$ and $\epsilon_{k-1}\geq\epsilon$}\\\
X_{k-1}&\textrm{if $P(\Omega(x_{k},\epsilon))<1-\Gamma$ and
$\epsilon_{k-1}<\epsilon$}\;,\end{cases}$

where $x_{k}=x_{k}(\bar{\omega}_{1},\ldots,\bar{\omega}_{k})$. (Note that if
state $\bullet$ is reached, it remains the state thereafter.)

The goal is to understand, probabilistically, when state $\bullet$ is reached,
assuming all samples $\omega_{k,\ell}$ are independently drawn from
$(\Omega,{\mathcal{F}},P)$.

According to Proposition 2.2, the distances $\|x_{k}-z_{0}\|$ ($k=0,1,\ldots$)
form a decreasing sequence (strictly decreasing except when $x_{k}=x_{k-1}$).
Moreover, the number of indices $k$ satisfying $\epsilon_{k-1}\geq\epsilon$
cannot exceed $N:=\lfloor(\frac{M\,\|x_{0}-z_{0}\|}{\epsilon})^{2}\rfloor$,
else by induction, (2.3) implies $\mathrm{dist}(x_{k},X_{\Omega})$ becomes
negative, a contradiction. Consequently, if indices $k_{1}<k_{2}<\cdots<k_{N}$
all satisfy $\epsilon_{k-1}\geq\epsilon$, then necessarily
$\epsilon_{k_{N}+1}<\epsilon$ – thus, necessarily
$\Omega(x_{k_{N}},\epsilon)=\Omega$ – and hence
$X_{k_{N}}(\bar{\omega}_{1},\ldots,\bar{\omega}_{k_{N}})=\bullet$. Hence, any
state occurring is either $\bullet$ or one of the values $0,1,\ldots,N-1$.

Until reaching state $\bullet$, the values $X_{k}$ change in a simple
stochastic manner expressed by the conditional probability

$\mathrm{Prob}_{\omega_{k,1},\ldots,\omega_{k,L}}(X_{k}=X_{k-1}+1\mid\bar{\omega}_{1},\ldots,\bar{\omega}_{k-1})=1-P(\Omega(x_{k-1},\epsilon))^{L}\geq
1-(1-\Gamma)^{L}=p\;.$ (2.4)

Consequently, starting at $x_{0}$, the expected number of iterations
sufficient to reach state $\bullet$ is no greater than the expected number of
steps for the compound Bernoulli process

$Y_{0}=0,\quad Y_{k}=\begin{cases}Y_{k-1}+1&\textrm{with probability $p$}\\\
Y_{k-1}&\textrm{with probability $1-p$}\end{cases}$

to reach value $N=\lfloor(\frac{M\,\|x_{0}-z_{0}\|}{\epsilon})^{2}\rfloor$, an
expectation which is easily shown to equal $\mathbb{E}:=N/p$.

For the concentration inequalities, consider any integer $k\geq 2\mathbb{E}$
and observe

$\displaystyle{\scriptstyle{\textstyle{\frac{\mathrm{Prob}\big{(}Y_{k+1}=N\textrm{
and }Y_{k}=N-1\big{)}}{\mathrm{Prob}\big{(}Y_{k}=N\textrm{ and
}Y_{k-1}=N-1\big{)}}}}}$
$\displaystyle={\scriptstyle{\textstyle{\frac{p\,{k\choose
N-1}p^{N-1}(1-p)^{k-(N-1)}}{p\,{k-1\choose N-1}p^{N-1}(1-p)^{(k-1)-(N-1)}}}}}$
(2.5)
$\displaystyle={\scriptstyle{\textstyle{\frac{k(1-p)}{k-(N-1)}}}}\leq{\scriptstyle{\textstyle{\frac{\frac{2N}{p}(1-p)}{\frac{2N}{p}-(N-1)}}}}={\scriptstyle{\textstyle{\frac{1-p}{1-p+\frac{1}{2}p+\frac{p}{N}}}}}<{\scriptstyle{\textstyle{\frac{1}{1+\frac{1}{2}\frac{p}{1-p}}}}}\;.$

The inequality (2.1) now follows by induction, the base case being $k=\lceil
2\mathbb{E}\rceil$ for which $\mathrm{Prob}(Y_{k}=N\textrm{ and
}Y_{k-1}=N-1)\leq{\scriptstyle{\textstyle{\frac{1}{2}}}}$ (an immediate
consequence of $\mathbb{E}$ being an upper bound on the expected number of
iterations). $\Box$

##  3\. Growth

In recent years a substantial literature has emerged that is focused on
establishing improved iteration bounds of first-order methods when applied to
to minimizing a convex function $h$ possessing Hölderian growth, this being
the property that there exist constants $\mu>0$ and $d\geq 1$ such that for
$x$ in the domain of $h$, we have the lower bound

$h(x)\geq\mu\,\mathrm{dist}(x,X^{*})^{d}$

where $X^{*}$ is the set of minimizers of $h$. When $d=1$, this is sometimes
referred to as “linear” growth or “sharpness”, whereas for $d=2$, the function
has “quadratic” growth. Strongly convex functions possess quadratic growth but
so do many other convex functions.

In this section we establish improved iteration bounds for $\mathtt{PolyakFM}$
when the collection of functions $f_{\omega}$ ($\omega\in\Omega$) possesses a
property exactly analogous to Hölderian growth. For motivation, consider a
collection of affine functions
$f_{\omega}(x)=\alpha_{\omega}^{T}x+\beta_{\omega}$ for which $X_{\Omega}$ is
nonempty. Assuming $\Omega$ is finite – say, $\Omega=\\{1,\ldots,m\\}$ – a
classic result of Hoffman [2] asserts that there exists a constant $\mu>0$
(now known as the Hoffman constant) such that for every $x\notin X_{\Omega}$,
$\max_{\omega}f_{\omega}(x)\geq\mu\,\mathrm{dist}(x,X_{\Omega})$, that is,
there is at least one function $f_{\omega}$ for which
$f_{\omega}(x)\geq\mu\,\mathrm{dist}(x,X_{\Omega})$.

However, when $m$ is astronomical in size (and especially when $\Omega$ is
infinite), there is little use in the property of having, for each $x$, at
least one of the functions satisfying
$f_{\omega}(x)\geq\mu\,\mathrm{dist}(x,X_{\Omega})$. Indeed, if there is only
one such function $f_{\omega}$ for an iterate $x_{k}$, what are the chances
that a randomized algorithm like $\mathtt{PolyakFM}$ will make use of that
function in computing $x_{k+1}$? Instead, what is needed is that a
“considerable portion” of samples $\omega\in\Omega$ satisfy
$f_{\omega}(x)\geq\mu\,\mathrm{dist}(x,\mathrm{Feas})$. This leads us to the
following assumption.

> Growth Assumption: There exist constants $\mu>0$, $d\geq 1$ and
> $0<\Delta\leq 1$ such that for every $x\in{\mathcal{B}}\setminus
> X_{\Omega}$,
>
> $P\\{\omega\in\Omega\mid
> f_{\omega}(x)\geq\mu\,\mathrm{dist}(x,X_{\Omega})^{d}\\}\geq\Delta\;.$

In the following theorem, as in Theorem 2.1, we assume
$\epsilon<M\mathrm{dist}(x_{0},X_{\Omega})$ since otherwise $x_{0}$ achieves
the goal of obtaining an iterate $x_{k}$ satisfying
$P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$.

######  Theorem 3.1.

Make the Growth Assumption. For $0<\epsilon<M\mathrm{dist}(x_{0},X_{\Omega})$
and $0<\Gamma<\Delta$, $\mathtt{PolyakFM}$ reaches an iterate $x_{k}$
satisfying $P(\Omega(x,\epsilon))\geq 1-\Gamma$ – and
$\mathrm{dist}(x_{k},X_{\Omega})\leq(\epsilon/\mu)^{1/d}$ – within an expected
number of iterations not exceeding

$\mathbb{E}^{\prime}={\scriptstyle{\textstyle{\frac{4}{p}}}}\left(1+\left({\scriptstyle{\textstyle{\frac{M}{\mu^{1/d}\epsilon^{1-1/d}}}}}\right)^{2}\min\left\\{{\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}-1}}}},\log_{2}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)\right\\}\right)\textrm{
where }p=1-(1-\Gamma)^{L}\;.$ (3.1)

Moreover, (2.1) provides concentration inequalities if $\mathbb{E}$ is
replaced by $\mathbb{E}^{\prime}$, and the requirement $k\geq 2\mathbb{E}$ is
replaced by $k\geq 2\mathbb{E}^{\prime}$.

Remarks:

  1. (1)

When $d=1$, the dependence on $1/\epsilon$ is only logarithmic. Necoara and
Nedić [3] established a similar iteration bound for an algorithm which takes a
subgradient step for the sample expectation,
$x\mapsto\frac{1}{L}\sum_{\ell=1}^{L}f_{\omega_{\ell}}(x)$. Their growth
assumption roughly translates as being that there exists $\tilde{\mu}>0$ such
that for all $x\in{\mathcal{B}}$, there holds the lower bound
$\mathbb{E}[f_{\omega}^{+}(x)]\geq\tilde{\mu}\,\mathrm{dist}(x,X_{\Omega})$.
While their minibatches require computing $L$ subgradients, their bounds have
no dependence on a quantity similar to $1/p$, and in fact, essentially replace
$1/p$ by a value that under non-trivial conditions decreases proportionally to
$1/L$, assuming that subgradient steps are extended to a particular length
that in some important cases can be computed in practice. They do not consider
growth of degree $d>1$, nor do we see how their proofs might be extended to
provide bounds of the quality of Theorem 3.1, especially with regards to
$\epsilon$. Nonetheless, we would be remiss not to acknowledge that the
interesting results of [3] were the primary catalyst for our paper.

  2. (2)

For all degrees of growth, the dependence on $\epsilon$ improves on the worst-
case complexity of subgradient methods in the deterministic setting,
$O(1/\epsilon^{2})$.

  3. (3)

With additional bookkeeping in the proof, the theorem can be improved if the
growth assumption is satisfied both for, say, degree $1$ with constant
$\mu_{1}$, and for degree $2$ with constant $\mu_{2}$. Then for $x$ farther
than distance $\mu_{1}/\mu_{2}$ away from $X_{\Omega}$, the quadratic lower
bound exceeds the linear lower bound, that is, the quadratic lower bound
dominates. Assuming the quadratic growth dominates at $x_{0}$, the manner of
proof of the theorem – and the fact that the algorithm does not require $d$ or
$\mu$ as input – allow (3.1) to be replaced by a sum of two expressions, one
obtained by substituting into (3.1) the values $d=2$, $\mu=\mu_{2}$ and the
choice $\epsilon=\mu_{1}^{2}/\mu_{2}$, and the other obtained by substituting
$d=1$, $\mu=\mu_{1}$ and replacing $\mathrm{dist}(x_{0},X_{\Omega})$ by
$\mu_{1}/\mu_{2}$. Similarly if the growth assumption is satisfied by any two
degrees $d<d^{\prime}$ (or any three (or more) degrees, but then the notation
becomes unwieldy).

Proof of Theorem 3.1: For non-negative integers $i$, define
$e_{i}=(\frac{1}{2})^{i}M\mathrm{dist}(x_{0},X_{\Omega})$, and let
$I:=\lceil\log_{2}(\frac{M\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon})\rceil$,
the smallest integer $i$ for which $e_{i}\leq\epsilon$.

Theorem 2.1 implies for each $i$ that with probability 1, the algorithm
reaches an iterate $x_{k}$ for which $P(\Omega(x_{k},e_{i}))\geq 1-\Gamma$.
Let $x_{k_{i}}$ ($=x_{k_{i}}(\bar{\omega}_{1},\ldots,\bar{\omega}_{k_{i}})$)
be the first such iterate. Note that $k_{0}=0$ because
$\Omega(x_{0},e_{0})=\Omega$ due to the choice
$e_{0}=M\,\mathrm{dist}(x_{0},X_{\Omega})$.

For $i\geq 1$, if in Theorem 2.1 we substitute $x_{k_{i-1}}$ for $x_{0}$ and
$e_{i}$ for $\epsilon$, we find that the expected number of iterations to
reach $x_{k_{i}}$ from $x_{k_{i-1}}$ does not exceed
$\frac{1}{p}(\frac{M\,\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{k_{i}}})^{2}$.
Consequently, beginning at $x_{0}$, the expected number of iterations to reach
an iterate $x_{k}$ satisfying $P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$ does
not exceed

${\scriptstyle{\textstyle{\frac{1}{p}}}}\sum_{i=1}^{I}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}}\right)^{2}\,\,=\,\,{\scriptstyle{\textstyle{\frac{1}{p}}}}\left(4+\sum_{i=2}^{I}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}}\right)^{2}\right).$
(3.2)

Since for $i\geq 1$, the point $x_{k_{i}}$ satisfies
$P(\Omega(x_{k_{i}},e_{i}))\geq 1-\Gamma$, our assumption $\Gamma<\Delta$
implies there exists $\omega_{i}\in\Omega(x_{k_{i}},e_{i})$ for which
$f_{\omega_{i}}(x_{k_{i}})\geq\mu\,\mathrm{dist}(x_{k_{i}},X_{\Omega})^{d}$.
Consequently, we have both

$e_{i}\geq f_{\omega_{i}}(x_{k_{i}})\textrm{ and
}f_{\omega_{i}}(x_{k_{i}})\geq\mu\,\mathrm{dist}(x_{k_{i}},X_{\Omega})^{d},\textrm{
and thus }\mathrm{dist}(x_{k_{i}},X_{\Omega})\leq(e_{i}/\mu)^{1/d}\;.$ (3.3)

Since $e_{i}={\scriptstyle{\textstyle{\frac{1}{2}}}}e_{i-1}$, this gives for
$i\geq 2$ that

${\scriptstyle{\textstyle{\frac{\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}}\leq{\scriptstyle{\textstyle{\frac{2}{e_{i-1}}}}}\left({\scriptstyle{\textstyle{\frac{e_{i-1}}{\mu}}}}\right)^{1/d}={\scriptstyle{\textstyle{\frac{2}{\mu^{1/d}e_{i-1}^{1-1/d}}}}}={\scriptstyle{\textstyle{\frac{2}{\mu^{1/d}(2^{I-i}e_{I-1})^{1-1/d}}}}}\leq{\scriptstyle{\textstyle{\frac{2}{\mu^{1/d}(2^{I-i}\epsilon)^{1-1/d}}}}}\;.$

Substituting into (3.2) shows the expected time to reach an iterate $x_{k}$
satisfying $P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$ does not exceed

${\scriptstyle{\textstyle{\frac{1}{p}}}}\left(4+\left({\scriptstyle{\textstyle{\frac{2M}{\mu^{1/d}\epsilon^{1-1/d}}}}}\right)^{2}\,\sum_{i=1}^{I-1}\left({\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}}}}}\right)^{I-i}\right)\;.$

Finally, observe

$\sum_{i=1}^{I-1}\left({\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}}}}}\right)^{I-i}\leq\min\\{I-1,\sum_{j=1}^{\infty}\left({\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}}}}}\right)^{j}\\}=\min\\{I-1,{\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}-1}}}}\\}\;,$

completing the proof that $\mathbb{E}^{\prime}$ is an upper bound on the
expected number of iterations to reach an iterate $x_{k}$ satisfying
$P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$. That the iterate also satisfies
$\mathrm{dist}(x_{k},X_{\Omega})\leq(\epsilon/\mu)^{1/d}$ is established
exactly like (3.3).

To verify the concentration inequalities, first observe we could have
established that $\mathbb{E}^{\prime}$ is an upper bound on the expected
number of iterations by relying on first principles as in the proof of Theorem
2.1, where the focus was on determining the expected number of iterations for
the compound Bernoulli compound $Y_{0},Y_{1},\ldots$ to first reach value
$N=\lfloor(\frac{M\|x_{0}-z_{0}\|}{\epsilon})^{2}\rfloor$. In the present
context the relevant target is instead

$\displaystyle N^{\prime}$
$\displaystyle=\sum_{i=1}^{I}\Bigl{\lfloor}({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}})^{2}\Bigr{\rfloor}$
$\displaystyle\leq
4\left(1+\left({\scriptstyle{\textstyle{\frac{M}{\mu^{1/d}\epsilon^{1-1/d}}}}}\right)^{2}\min\left\\{{\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}-1}}}},\log_{2}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)\right\\}\right)$
$\displaystyle=p\,\mathbb{E}^{\prime}\;.$

The desired concentration inequalities follow, by substituting $N^{\prime}$
for $N$ in (2.5) (and substituting $\mathbb{E}^{\prime}$ for ${\mathbb{E}}$).
$\Box$

For general convex functions possessing Hölderian growth, the trick of using
summations like (3.2) in establishing improved iteration bounds was first
developed extensively in [6], albeit in a non-stochastic setting, yet in a
manner that applies to a wide range of first order methods, not just
subgradient methods.

##  4\. Confidence

While $\mathtt{PolyakFM}$ will, with probability 1, compute an iterate $x_{k}$
satisfying $P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$, in general it is
impossible or impractical to be certain that such an iterate has been reached
(an exception being when $\Omega$ consists of only a modest number of
samples). In this section we present a modification to $\mathtt{PolyakFM}$
resulting in an algorithm which determines, with high confidence, that an
iterate satisfying $P(\Omega(x_{k},\epsilon))\geq 1-\Gamma$ is in hand.

The only change to $\mathtt{PolyakFM}$ is that now the user inputs $x_{0}$ and
values $0<\Gamma<1$, $0<\alpha<1$, the latter being an allowable probability
of error. (The algorithm computes the batch size based on $\Gamma$ and
$\alpha$.) Here the iterate pairs $(x_{k},\epsilon_{k})$ satisfy
$P(\Omega(x_{k},\epsilon_{k})\geq 1-\Gamma$ “with confidence $1-\alpha$”, by
which is meant that across all sequences of samples
$\bar{\omega}_{1},\bar{\omega}_{2},\ldots$, the probability that even one
error is made does not exceed $\alpha$ (an error occurs when a pair
$(x_{k},\epsilon_{k})$ fails to satisfy the inequality
$P(\Omega(x_{k},\epsilon_{k}))\geq 1-\Gamma$). For $\epsilon>0$, we establish
bounds on the number of iterations needed to reach a pair
$(x_{k},\epsilon_{k})$ for which $\epsilon_{k}\leq\epsilon$. Perhaps
surprisingly, these iteration bounds are deterministic – all of the
probability gets wrapped into the level of confidence, $1-\alpha$ (but at the
expense of large minibatches if $\Gamma\approx 0$).

Following is the modified algorithm, in which the batch size $L$ slowly
increases with the number of iterations.

$\mathtt{confidentPFM}$ (Confident Polyak Feasibility Method)
0) Input: $x_{0}$, $0<\Gamma<1$, $0<\alpha<1$ (allowable probability of
error).
Initialization: $k=1$
1) Compute $L=\lceil\frac{1}{\Gamma}\ln(\frac{2k^{2}}{\alpha})\rceil$ and
independently select samples $\omega_{k,1},\ldots,\omega_{k,L}$ from
$(\Omega,{\mathcal{F}},P)$.
2) Determine
$\ell_{k}:=\arg\max_{\ell=1,\ldots,L}f_{\omega_{k,\ell}}(x_{k-1})$, then let
$\epsilon_{k-1}:=f_{\omega_{k,\ell_{k}}}(x_{k-1})$ and compute

$x_{k}=\begin{cases}x_{k-1}&\textrm{if $\epsilon_{k-1}\leq 0$}\\\
x_{k-1}-\frac{\epsilon_{k-1}}{\|g_{\omega_{k,\ell_{k}}}\|^{2}}\,g_{\omega_{k,\ell_{k}}}&\textrm{otherwise}\;,\end{cases}$

where $g_{\omega_{k,\ell_{k}}}\negthinspace\negthinspace\in\partial
f_{\omega_{k,\ell_{k}}}(x_{k-1})$.
3) Let $k+1\rightarrow k$ and return to Step 1.

In the following theorem, as in Theorems 2.1 and 3.1, we assume the desired
accuracy $\epsilon$ satisfies $\epsilon<M\,\mathrm{dist}(x_{0},X_{\Omega})$
since otherwise all samples $\omega$ satisfy $f_{\omega}(x_{0})\leq\epsilon$
and thus $\epsilon_{0}$, computed in the first iteration, satisfies
$\epsilon_{0}\leq\epsilon$.

######  Theorem 4.1.

The probability that ${\tt confidentPFM}$ makes no errors is at least
$1-\alpha$.

If $0<\epsilon<M\mathrm{dist}(x_{0},X_{\Omega})$, then ${\tt confidentPFM}$
computes $\epsilon_{k}$ satisfying $\epsilon_{k}\leq\epsilon$

$\textrm{within
}\,1+\Bigl{\lfloor}\left({\scriptstyle{\textstyle{\frac{M\,\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)^{2}\Bigr{\rfloor}\textrm{
iterations.}$ (4.1)

Moreover, if (1) the Growth Assumption holds, (2) $\Gamma<\Delta$ and (3)
${\tt confidentPFM}$ makes no errors, then the number of iterations does not
exceed

$5+4\left({\scriptstyle{\textstyle{\frac{M}{\mu^{1/d}\epsilon^{1-1/d}}}}}\right)^{2}\min\left\\{{\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}-1}}}},\log_{2}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)\right\\}\;.$
(4.2)

Remarks:

  1. (1)

We know of no similar general result in the literature on subgradient methods,
although there certainly have been stronger high-confidence results in
optimization settings focused on narrower classes of functions (e.g., the
results for strongly convex stochastic composite optimization in the seminal
work of Ghadimi and Lan [1, §4.2]).

  2. (2)

The main shortcoming of the setup is that the number of samples per iteration
is proportional to $1/\Gamma$, a large number of samples if $\Gamma\approx 0$.
On the other hand, due to the slow increase in sample size per iteration, the
total number of samples drawn in iterations $1,\ldots,k$ is
$O(\frac{k}{\Gamma}\log(\frac{k}{\alpha})$, ideal in $\alpha$ and off by at
most a factor of $\log(k)$ from being ideal in $k$.

Proof of Theorem 4.1: First we verify that the probability of error is no
greater than $\alpha$. Let
$x_{k-1}=x_{k-1}(\bar{\omega}_{1},\ldots,\bar{\omega}_{k-1})$ be an iterate
(possibly $k-1=0$). Let
$L_{k}=\lceil\frac{1}{\Gamma}\ln(\frac{2k^{2}}{\alpha})\rceil$, the number of
independent samples $\omega_{k,\ell}$ drawn for iteration $k$. Observe that by
definition of $\epsilon_{k-1}$, all of the samples satisfy
$\omega_{k,\ell}\in\Omega(x_{k-1},\epsilon_{k-1})$.

Assume the algorithm makes an error on iterate $k$, that is, assume
$P(\Omega(x_{k-1},\epsilon_{k-1}))<1-\Gamma$. Then the probability of having
chosen all samples from $\Omega(x_{k-1},\epsilon_{k-1})$ is

$P(\Omega(x_{k-1},\epsilon))^{L_{k}}<(1-\Gamma)^{L_{k}}=(1-\Gamma)^{\lceil\frac{1}{\Gamma}\ln(\frac{2k^{2}}{\alpha})\rceil}\leq{\scriptstyle{\textstyle{\frac{\alpha}{2k^{2}}}}}$

– the final inequality holds because an integer $n$ satisfies
$(1-\Gamma)^{n}\leq\frac{\alpha}{2k^{2}}$ if and only if
$n\geq{\scriptstyle{\textstyle{\frac{\ln(\frac{\alpha}{2k^{2}})}{\ln(1-\Gamma)}}}}={\scriptstyle{\textstyle{\frac{\ln(2k^{2}/\alpha)}{-\ln(1-\Gamma)}}}}$,
and hence if $n\geq\frac{1}{\Gamma}\ln(\frac{2k^{2}}{\alpha})$. Consequently,
the probability of $\mathtt{confidentPFM}$ ever making an error does not
exceed $\sum_{k=1}^{\infty}\frac{\alpha}{2k^{2}}<\alpha$.

We next establish (4.1). This is easily accomplished. Indeed, whenever
$\epsilon_{k-1}>\epsilon$, we have from (2.3) that

$\mathrm{dist}(x_{k},X_{\Omega})^{2}<\mathrm{dist}(x_{k-1},X_{\Omega})^{2}-(\epsilon/M)^{2}\;,$

an amount of decrease that can occur for at most
$\lfloor(\frac{M\,\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon})^{2}\rfloor$
iterations. Since $\epsilon_{k}$ is computed one iteration later than $x_{k}$,
(4.1) is thus a valid upper bound.

Now assume the Growth Assumption holds and $\Gamma<\Delta$. As in the proof of
Theorem 3.1, for integers $i\geq 0$ let
$e_{i}=(\frac{1}{2})^{i}M\mathrm{dist}(x_{0},X_{\Omega})$. Define
$I=\lceil\frac{M\,\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}\rceil$, the
smallest integer $i$ for which $e_{i}\leq\epsilon$.

For all $i$ we know from (4.1) that for some iterate $x_{k}$, the value
$\epsilon_{k}$ computed in the following iteration satisfies $\epsilon_{k}\leq
e_{i}$. Let $x_{k_{i}}$ be the first such iterate. Of course $k_{0}=0$ since
$\Omega(x_{0},e_{0})=\Omega$ due to
$e_{0}=M\,\mathrm{dist}(x_{0},X_{\Omega})$.

For $i\geq 1$, if in (4.1) we substitute $x_{k_{i-1}}$ for $x_{0}$ and $e_{i}$
for $\epsilon$, we find that the number of iterations to reach $x_{k_{i}}$
from $x_{k_{i-1}}$ does not exceed
$\lfloor(\frac{M\,\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}})^{2}\rfloor$
(one iteration less than (4.1) because $x_{k_{i}}$ is computed one iteration
before $e_{k_{i}}$). Consequently, beginning at $x_{0}$, the number of
iterations to reach an iterate $x_{k}$ for which $\epsilon_{k}\leq\epsilon$
does not exceed

$\sum_{i=1}^{I}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}}\right)^{2}\,\,=\,\,4+\sum_{i=2}^{I}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{k_{i-1}},X_{\Omega})}{e_{i}}}}}\right)^{2}.$
(4.3)

Observe that this deterministic iteration bound is the exact analogue of the
expected iteration bound (3.2) – only the factor $1/p$ is now missing. From
this point onwards the proof follows that of Theorem 3.1, assuming
$\mathtt{confidentPFM}$ makes no errors. In particular, for the iterates
$x_{k_{i}}$ we need that $P(\Omega(x_{k_{i}},\epsilon_{k_{i}}))\geq 1-\Gamma$,
implying $P(\Omega(x_{k_{i}},e_{i}))\geq 1-\Gamma$.

Thus, assuming $\mathtt{confidentPFM}$ makes no errors, and proceeding exactly
as in the proof of Theorem 3.1 beginning at (3.2) (but disregarding the factor
$1/p$), we find that the number of iterations required to reach $x_{k_{I}}$
from $x_{0}$ does not exceed

$4\left(1+\left({\scriptstyle{\textstyle{\frac{M}{\mu^{1/d}\epsilon^{1-1/d}}}}}\right)^{2}\min\left\\{{\scriptstyle{\textstyle{\frac{1}{4^{1-1/d}-1}}}},\log_{2}\left({\scriptstyle{\textstyle{\frac{M\mathrm{dist}(x_{0},X_{\Omega})}{\epsilon}}}}\right)\right\\}\right)\;.$

Since $\epsilon_{k_{I}}$ ($\leq\epsilon$) is not computed until the following
iteration, (4.1) is thus a valid deterministic bound. $\Box$

## References

  * [1] Saeed Ghadimi and Guanghui Lan.  Optimal stochastic approximation algorithms for strongly convex stochastic composite optimization i: A generic algorithmic framework.  SIAM Journal on Optimization, 22(4):1469–1492, 2012.
  * [2] Alan J Hoffman.  On approximate solutions of systems of linear inequalities.  In Selected Papers Of Alan J Hoffman: With Commentary, pages 174–176. World Scientific, 2003.
  * [3] Ion Necoara and Angelia Nedić.  Minibatch stochastic subgradient-based projection algorithms for feasibility problems with convex inequalities.  Computational Optimization and Applications, pages 1–32, 2021.
  * [4] Ion Necoara, Peter Richtárik, and Andrei Patrascu.  Randomized projection methods for convex feasibility: Conditioning and convergence rates.  SIAM Journal on Optimization, 29(4):2814–2852, 2019.
  * [5] BT Polyak.  Random algorithms for solving convex inequalities.  In Studies in Computational Mathematics, volume 8, pages 409–422. Elsevier, 2001.
  * [6] James Renegar and Benjamin Grimmer.  A simple nearly optimal restart scheme for speeding up first-order methods.  Foundations of Computational Mathematics, pages 1–46, 2021.


# Improving qubit coherence using closed-loop feedback

Antti Vepsäläinen<EMAIL_ADDRESS>Research Laboratory of Electronics,
Massachusetts Institute of Technology  Roni Winik  Research Laboratory of
Electronics, Massachusetts Institute of Technology  Amir H. Karamlou
Department of Electrical Engineering and Computer Science, Massachusetts
Institute of Technology  Research Laboratory of Electronics, Massachusetts
Institute of Technology  Jochen Braumüller  Research Laboratory of
Electronics, Massachusetts Institute of Technology  Agustin Di Paolo  Research
Laboratory of Electronics, Massachusetts Institute of Technology  Youngkyu
Sung  Department of Electrical Engineering and Computer Science, Massachusetts
Institute of Technology  Bharath Kannan  Department of Electrical Engineering
and Computer Science, Massachusetts Institute of Technology  Morten Kjaergaard
Research Laboratory of Electronics, Massachusetts Institute of Technology
Center for Quantum Devices, University of Copenhagen  David K. Kim  MIT
Lincoln Laboratory  Alexander J. Melville  MIT Lincoln Laboratory  Bethany M.
Niedzielski  MIT Lincoln Laboratory  Jonilyn L. Yoder  MIT Lincoln Laboratory
Simon Gustavsson  Research Laboratory of Electronics, Massachusetts Institute
of Technology  William D. Oliver  Department of Electrical Engineering and
Computer Science, Massachusetts Institute of Technology  MIT Lincoln
Laboratory

###### Abstract

Superconducting qubits are a promising platform for building a larger-scale
quantum processor capable of solving otherwise intractable problems. In order
for the processor to reach practical viability, the gate errors need to be
further suppressed and remain stable for extended periods of time. With recent
advances in qubit control, both single- and two-qubit gate fidelities are now
in many cases limited by the coherence times of the qubits. Here we
experimentally employ closed-loop feedback to stabilize the frequency
fluctuations of a superconducting transmon qubit, thereby increasing its
coherence time by 26% and reducing the single-qubit error rate from
$($8.5$\pm$2.1$)\times 10^{-4}$ to $($5.9$\pm$0.7$)\times 10^{-4}$.
Importantly, the resulting high-fidelity operation remains effective even away
from the qubit flux-noise insensitive point, significantly increasing the
frequency bandwidth over which the qubit can be operated with high fidelity.
This approach is helpful in large qubit grids, where frequency crowding and
parasitic interactions between the qubits limit their performance.

##  I Introduction

High fidelity single- and two-qubit gates are a prerequisite for high-depth
circuits and quantum error correction. For many qubit modalities, including
superconducting qubits, the qubit frequencies and their controls are subject
to temporally correlated noise - most notably $1/f$-type noise - resulting in
correlated errors and slow drifts in frequency. Here, we use active feedback
control to stabilize the frequency drift of a flux-tunable superconducting
transmon qubit and thereby suppress its temporal frequency fluctuations. This
results in improved coherence times and gate fidelities, which remain stable
over extended periods of time.

Closed-loop feedback control is ubiquitously used in a wide variety of
engineering applications. There is an increasing number of experiments where
feedback is used to modify the evolution of quantum systems, such as
stabilizing the motion of the atoms in optical cavities [1], cooling quantum
mechanical resonators [2], or stabilizing Rabi oscillations in superconducting
qubits [3]. These experiments are based on continuous monitoring of the
system, which inevitably decoheres the quantum state. In Ref. [4], a method
based on interleaving the probing sequences with separate periods of time used
for computation was introduced and demonstrated to mitigate the effect of slow
magnetic field fluctuations. In Ref. [5], it was further shown that this
method can be used to decouple an electron spin from low frequency magnetic
field fluctuations, resulting in improved coherence times. Ref. [6] employs a
slightly different approach in a trapped-ion system, where spectator qubits
are used to probe spatially correlated errors in the control laser amplitude
and targeting. Here, we demonstrate that closed-loop feedback can be used to
mitigate decoherence and improve gate fidelities in a frequency-tunable
superconducting qubit system.

The dominant source of decoherence in superconducting qubits is typically
either charge noise or flux noise but depending on the qubit design also
photon shot noise [7] or Bogoliubov quasiparticles may contribute [8]. These
noise sources are primarily intrinsic and local to the device [9, 10] — as
opposed to noise in the control electronics — though recently there has been
some evidence of correlated noise between the qubits [11, 12, 13]. In our
experiment, we employ flux-tunable transmon qubits [14], which are widely used
in contemporary superconducting quantum processors [15, 16]. Due to their
design, transmons are mostly insensitive to charge-noise [14], but suffer from
noise in magnetic flux, which is used to tune their frequency. In order to
protect the qubits from flux noise, these qubits are typically operated at
bias points where their frequency is first-order insensitive to small changes
in flux, colloquially referred to as a sweet spot, see Fig. 1a. However, in
many architectures the qubits cannot be operated at the sweet spot
indefinitely while performing logic operations [17, 18]. Additionally, it is
necessary to operate some of the qubits away from their sweet spots, in part
due to parasitic couplings to two-level fluctuators with frequencies near the
sweet spot [19, 20, 21, 22] or couplings to other qubits or their higher
excited states in the quantum processor [16]. We demonstrate that using the
feedback protocol to suppress low frequency noise enables gate fidelities
exceeding 99.9% even far away from the flux sweet spot, greatly extending the
range of available operation frequencies for such qubits.

Figure 1: Qubit frequency power spectral density a) The spectrum of the qubit
as a function of flux bias. The orange dot shows the flux bias point at which
the qubit is operated. b) Schematic of the frequency estimation pulse sequence
(grey) interleaved with the sequence used for computation (light blue). The
qubit frequency is adjusted in between (green). c) Schematic of the feedback
loop. The measurement record $q_{i}$ of a repeated Ramsey experiment is used
to estimate the qubit frequency offset $\hat{\delta}_{\rm q}[n]$ subject to a
noisy environment $i_{\rm noise}$. The feedback signal $p[n]$ controls an AWG
that produces a current $i[n]$ to cancel the fluctuations in the qubit
frequency. d) Power spectral density of the qubit frequency fluctuations. The
blue dots (line) show the measured (simulated) spectral density of the qubit
frequency fluctuations estimated from $N=20$ Ramsey experiments
$S_{\hat{\delta}_{\rm q}\hat{\delta}_{\rm q}}(f)$, limited by the statistical
sampling noise (dashed black line). The spectrum with the sampling noise
suppressed using cross-correlation between consecutive samples is shown with
green dots along with a fit (green line). Orange dots (line) are the measured
(simulated) power spectral density of the error signal $\hat{\delta_{\rm
q}}[n]$ with the feedback activated. The simulated power spectral density of
the actual qubit frequency fluctuations $f_{\rm q}(t)$ is shown with solid red
line. The gray background describes the noise power removed by the feedback.
The frequency response of the feedback signal is overlaid in the plot with
solid black line.

##  II Results

The feedback protocol we employ consists of three phases that are continuously
repeated, see Fig 1b. In the probing phase, the qubit frequency is estimated
using a simple single-qubit frequency estimation algorithm. After the
frequency of the qubit is estimated, the magnetic flux through the qubit SQUID
loop is adjusted to set the qubit frequency to its target value. This is
followed by a computation phase, where an algorithm can be run with a freshly
stabilized qubit.

The probing phase consists of $N$ repeated Ramsey measurements. For each of
those measurements, the qubit is first prepared in a superposition state
$|\psi\rangle=\left(|0\rangle+|1\rangle\right)/\sqrt{2}$ using a $\pi/2$
rotation around the $y$-axis of the Bloch sphere. The state preparation is
followed by a period of free evolution for a duration $\tau$, during which the
qubit state acquires a phase $\phi=2\pi\int_{0}^{\tau}\delta_{\rm
q}(t)\mathrm{d}t$, where $\delta_{\rm q}(t)=f_{\rm d}-f_{\rm q}(t)$ is the
detuning between a microwave drive frequency $f_{\rm d}$ defining a rotating
reference frame, and $f_{\rm q}(t)$ is the fluctuating qubit frequency in the
presence of noise. A second $\pi/2$ pulse is then applied around the $x$-axis,
and the state of the qubit is measured using dispersive readout [23]. To
simplify the feedback protocol, we make a quasi-static approximation and
assume that the qubit frequency remains constant within one frequency
estimation experiment – $N$ Ramsey measurements – but may fluctuate between
the experiments. With this assumption, the probability of measuring the qubit
in the excited state is given by

$p_{1}=\frac{1}{2}+\frac{1}{2}\cos(2\pi\delta_{\rm q}\tau-\pi/2),$ (1)

which can be inverted to yield the frequency shift

$\delta_{\rm q}=\frac{\pm\arccos(2p_{1}-1)+2\pi k+\pi/2}{2\pi\tau},$ (2)

where $k$ is an integer. Equation (2) is a one-to-one mapping from $p_{1}$ to
the frequency detuning $\delta_{\rm q}$ over the domain $\delta_{\rm
q}\in[-\frac{1}{4\tau},\frac{1}{4\tau}]$. This implies the fluctuations in
$\delta_{\rm q}$ need to be within $\pm\frac{1}{4\tau}$ between the estimation
steps, approximately $70\text{\,}\mathrm{\SIUnitSymbolMicro s}$ in our
experiment.

Figure 2: Improvement in qubit coherence and stability a) The coherence time
$T_{2}$ of the qubit is measured by interleaving a Ramsey measurement with
feedback sequences (orange dots), compared to when the feedback is off (blue
dots). The solid lines are a fit to the data. Dashed lines show a simulation
of the decay envelope assuming the noise spectral densities shown in Fig. 1d.
b) The measured coherence time of the qubit as a function of the duration of
the Ramsey experiment is shown with the feedback (orange dots) and without
(blue dots). The solid lines are the expected coherence times based on the
measured noise spectral density. The inset shows the measured deviation of the
qubit frequency from the target frequency during the experiment. Each point is
calculated using 50 averages of the Ramsey trace. Panel c) shows Ramsey
experiment repeated for 4 hours with (right) and without (left) feedback. With
the feedback, the qubit frequency is stable during the whole duration. d)
Averaging the data in panel c) results in significant reduction in qubit
coherence if feedback is not used (left). Using the feedback, the inferred
coherence is not affected by the measurement duration (right). The coherence
time of a single Ramsey trace here is slightly lower than in a) due to smaller
value of $\tau=$500\text{\,}\mathrm{ns}$$ used in the feedback.

The qubit excited-state probability estimator
$\hat{p}_{1}=\frac{1}{N}\sum_{i=0}^{N}q_{i}$ is calculated from the
measurement record $q_{i}$ of $N=20$ repetitions of the Ramsey sequence,
providing an estimate for the frequency detuning of the qubit,
$\hat{\delta}_{\rm q}$. If in the previous measurement the qubit was measured
to be in the excited state, we virtually reset the qubit state for the current
repetition by flipping $q_{i}$ [24, 25]. The duration of a Ramsey measurement
is $T=$3.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, comprising the phase
accumulation time $\tau=$1.25\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, the
readout duration of $750\text{\,}\mathrm{ns}$, and the combined resonator
reset time and overhead from the electronics,
$1.5\text{\,}\mathrm{\SIUnitSymbolMicro s}$. Thus, a single round of frequency
estimation takes $T_{\rm N}=NT=$70\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ in
our implementation. The frequency estimation is not able to detect
fluctuations that occur faster than the repetition period of the feedback,
ultimately limiting the effective bandwidth of the noise suppression. Next we
investigate how the feedback reduces the noise spectral density of the qubit.

###  II.1 Noise spectral density

To estimate the noise power spectral density affecting the qubit frequency, we
first bias the qubit away from the sweet spot by $0.11\text{\,}$ flux quanta
$\Phi_{0}={\rm h}/2{\rm e}$, at a transition frequency $f_{\rm
q}=$4.69\text{\,}\mathrm{GHz}$$ that is sensitive to flux noise (Fig. 1a), and
repeatedly estimate its frequency $10^{5}$ times. The blue dots in Fig. 1d
show the power spectral density of the measured qubit frequency fluctuations
subject to flux noise. The measured power spectral density follows a power law
(green solid line) $S_{f_{\rm q}f_{\rm q}}(f)=A_{f_{\rm
q}}\left(\frac{$1\text{\,}\mathrm{Hz}$}{f}\right)^{\alpha}\approx$27.3\text{\times}{10}^{6}\text{\,}\mathrm{Hz}^{2}\mathrm{/}\mathrm{Hz}$\times\left(\frac{$1\text{\,}\mathrm{Hz}$}{f}\right)^{0.8}$,
but starts to deviate at frequencies above $500\text{\,}\mathrm{H}\mathrm{z}$
for $N=20$ and $\tau=$1.25\text{\,}\mathrm{\SIUnitSymbolMicro s}$$. This is
due to noise added by the finite number of samples in the estimate
$\hat{\delta}_{\rm q}$ and can be approximated as

$\delta\hat{\delta}_{\rm q}\approx\frac{1}{2\pi\tau\sqrt{N}},$ (3)

see Supplementary material for the derivation. The statistical sampling noise
is modeled as Gaussian white noise with an upper cutoff given by the duration
of the frequency estimation [26],

$S_{\rm est}(f)=\begin{cases}\frac{T}{2\pi^{2}\tau^{2}},&0\leq
f\leq\frac{1}{2NT},\\\ 0,&\text{otherwise},\end{cases}$ (4)

shown with dashed black line in Fig. 1d. The sampling noise can be suppressed
in post-processing by cross-correlating time-shifted measurement traces as
demonstrated in Ref. [26]. The results of applying this protocol are shown
with green dots in Fig. 1d. With the sampling noise suppressed, the measured
power spectrum fits well to the power law across the whole bandwidth. There is
an additional peak at $60\text{\,}\mathrm{Hz}$ corresponding to the noise from
the mains power. The sampling noise suppression is not used in the real-time
feedback signal calculation due to small amount of samples available at the
time of computation.

Figure 3: Qubit coherence at different bias points a) The qubit dephasing rate
$\Gamma_{\phi}$ is evaluated at several different bias fluxes. Qubit’s
sensitivity to flux noise increases further away from the sweet spot,
resulting in reduced dephasing times (increased dephasing rates). The error
bars show 68% confidence intervals for the median of the dephasing times
measured 60 times. The solid lines show a linear fit to the dephasing rates
with respect to the curvature of the qubit spectrum with respect to flux. The
black stars show the dephasing times extracted from a spin-echo experiment,
used as a reference.

Next, we turn on the feedback to reduce the fluctuations in the qubit
frequency. We aim to minimize the deviation of the qubit frequency $f_{\rm
q}(t)$ from the target frequency $f_{\rm d}$ by using the offset
$\hat{\delta}_{\rm q}[n]$ as the error signal in the feedback loop, see Fig.
1c for the schematic of the signal flow. Here we use $n$ to number each time
feedback is applied, sampled at times $t_{n}=nT_{\rm N}$. In practice, the
sampled error signal represents the average of the qubit frequency fluctuation
during the sampling period $T_{\rm N}$, limiting the maximum bandwidth of the
feedback to $1/(2T_{\rm N})\approx$7\text{\,}\mathrm{kHz}$$ if the time spent
on the interleaved computation step is omitted. The error signal
$\hat{\delta}_{\rm q}[n]$ is multiplied by a controllable gain $G$ and fed
into an accumulator that controls the feedback signal,
$p[n]=p[n-1]+G\delta_{\rm q}[n]$. We deliberately set $G=0.35$ to reduce the
bandwidth of the feedback frequency response $H_{\rm p}(f)$, shown with solid
black line in Fig. 1d, to be less affected by the statistical sampling noise,
see Supplementary material for additional details.

The output of the accumulator $p[n]$ is scaled and converted to an arbitrary
waveform generator voltage that drives the current responsible for creating a
magnetic flux through the qubit loop, adjusting its transition frequency. The
feedback significantly reduces the noise spectral density of the error-signal
$\hat{\delta}_{\rm q}[n]$, shown with orange dots in Fig. 1d.

After closing the feedback loop, the error signal $\hat{\delta}_{\rm q}[n]$
provides only indirect information about the actual qubit frequency
fluctuations, which are also affected by the frequency response of the
feedback (see Supplementary material for theoretical analysis). Thereby, we
employ a simulation to assess the impact of the feedback transfer function on
stabilizing the qubit frequency $f_{\rm q}(t)$. We use the fitted noise
spectral density as the starting point of the simulation (green solid line in
Fig. 1d) to generate time traces of the fluctuating qubit frequencies. Using
the same parameters as in the experiment, we simulate the estimation of the
qubit frequency and the feedback, which results in a time trace of estimated
qubit frequencies. The power spectral density of the simulated qubit frequency
estimation without the feedback is shown with solid blue line in Fig. 1d,
matching the experiment almost perfectly. The solid orange line shows the
simulated qubit frequency estimates when the feedback is turned on, again
matching well with the experiment. Finally, using the simulation we can
calculate the power spectral density of the real qubit frequency fluctuations
when the feedback is applied, shown with solid red line in Fig. 1d. The total
noise power mitigated by the feedback protocol is indicated by the shaded area
in Fig. 1d.

###  II.2 Qubit coherence

Next, we interleave the frequency estimation sequences with a Ramsey
experiment to demonstrate that the lower noise power increases the qubit
coherence time $T_{2}$, see schematic in Fig. 1b. We apply two $\pi/2$ pulses
around the $y$-axis of the Bloch sphere, and sweep the delay between the
pulses, $\tau_{\rm R}$. Fig. 2a shows the qubit excited state population
$p_{1}$ going through the Ramsey oscillations, first without the feedback, and
then with the feedback turned on. For a Gaussian-distributed noise process,
the envelope of the Ramsey oscillations decays as [27]

$\chi_{\rm R}(t)=\exp\left[-2t^{2}\pi^{2}\int_{f_{0}}^{\infty}S_{f_{\rm
q}f_{\rm q}}(f)\,\mathrm{sinc}^{2}(\pi ft)\,\mathrm{d}f\right],$ (5)

where $S_{f_{\rm q}f_{\rm q}}(f)$ is the unilateral power spectral density of
the frequency fluctuations and $f_{0}$ is the lower cutoff frequency equal to
the inverse of the total duration of the experiment. We fit the experimental
data to an oscillating function with a decay envelope corresponding to a power
spectral density of $1/f$ noise and extract the coherence time from when the
decay envelope drops below $1/\mathrm{e}$ from its original value at $t_{\rm
R}=0$. The feedback increases the coherence time from
$T_{2}=$6.2\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ to
$T_{2}=$7.7\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, or 26%. The observed
improvement in coherence time was consistent across several repetitions of the
experiment. The fitted decay envelope can be compared to a theoretical
estimate obtained by directly substituting the measured power spectra to Eq.
(5), shown with dashed lines in Fig. 2a, closely matching experimentally
observed decay envelopes. The only fit parameter is the qubit state
initialization fidelity, 92%, which scales the decay envelope amplitude. This
confirms that the decoherence of the qubit is well described by the measured
power spectral density of the frequency fluctuations.

The inferred coherence time of the qubit depends on the total duration of the
Ramsey experiment through the cutoff frequency $f_{0}$ in Eq. (5). In Fig. 2b,
the coherence time $T_{2}$ is evaluated using different numbers of averages to
calculate $p_{1}$, thereby changing the total duration of the single Ramsey
experiment and the cutoff frequency $f_{0}$. The data from the experiment is
collected only once, and is sectioned to different numbers of averages as
post-processing. Without the feedback, the coherence time gradually decreases
as the duration of the experiment increases. This is due to the increased
noise power at lower frequencies for $1/f$ noise. When the feedback is
activated, the overall coherence time is increased — similar to the experiment
shown in Fig. 2a — and remains constant independent of the cutoff frequency
$f_{0}$ due to the elimination of the low frequency noise. Moreover, the qubit
frequency remains stable during the measurement, as inferred from the
frequency of the Ramsey oscillations during the experiment, shown in inset of
Fig. 2b. The qubit frequencies in the inset are evaluated from a running
average of 50 Ramsey traces.

The stability of the qubit frequency can be maintained for hours using the
feedback protocol as demonstrated in Fig. 2c. There, Ramsey experiments are
repeated for more than four hours, either with the feedback turned off (left
panel) or on (right panel). When the measured qubit excited state populations
are overlaid (Fig. 2d), the Ramsey oscillations without the feedback are
blurred due to the constant fluctuation in the qubit frequency, whereas with
the feedback the oscillations are clearly visible. The coherence time
extracted from the average of all uncorrected (feedback off) Ramsey
experiments is only $T_{2}=$2.1\text{\,}\mathrm{\SIUnitSymbolMicro s}$$,
compared to $T_{2}=$6.9\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ with the
feedback.

Thus far, we have operated the qubit at a fixed flux bias point. We next
demonstrate that the feedback protocol improves the coherence time for a range
of flux biases. The sensitivity of the qubit to the flux noise is determined
by the curvature of its frequency spectrum with respect to flux bias, allowing
us to probe the efficiency of the feedback at various noise levels and
intrinsic coherence times. We evaluate the pure dephasing rate of the qubit at
11 different bias points, first without the feedback and then with it on. The
pure dephasing rates are extracted from the decay of the Ramsey oscillations
by subtracting the effect of the energy-relaxation rate [9]. The measured
dephasing rates (Fig. 3) are lowest close to the flux sweet spot and gradually
increase away from this spot as the qubit sensitivity to the noise increases.
In the limit where decoherence is dominated by flux noise, Eq. (5) can be used
to show that there is an (almost) linear dependence between the dephasing rate
$\Gamma_{\phi}$ and the flux sensitivity of the qubit [27, 9],

$\Gamma_{\phi}=k\left|\frac{\partial f_{\rm q}}{\partial\Phi}\right|.$ (6)

We find the coefficient $k$ from a linear fit to the data in Fig. 3a and use
the result to assess the impact of the feedback on mitigating the flux noise.
Without the feedback $k=(58\pm
1.1)$\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$, and reduces to $k=(48\pm
1.4)$\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$ when the feedback is used,
an improvement of 17%. This implies that the feedback effectively reduces the
flux noise amplitude seen by the qubit.

The ultimate goal would be to suppress the low frequency noise to the level
that the coherence time measured using a spin echo experiment — which is
insensitive to low frequency noise — would be equal or higher than the
coherence time in a feedback stabilized Ramsey experiment. In Fig. 3a, the
dephasing rate measured from an echo experiment are shown with black dots, and
the linear fit yields $k_{\rm E}=(17\pm
0.5)$\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$, which corresponds to the
$1/f$ flux-noise amplitude of $\sqrt{A_{\Phi}}=(3.3\pm
0.1)$\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$, see Supplementary
material for additional information. The main reason why the impact of the
noise on the echo experiment is lower than the on feedback-stabilized Ramsey
experiment is the limited bandwidth of our feedback implementation. While the
feedback efficiently suppresses the noise up to $1\text{\,}\mathrm{kHz}$, the
echo experiment is mostly insensitive to the noise below the inverse duration
of a single echo experiment, here approximately
$\sim$100\text{\,}\mathrm{kHz}$$.

Unlike the spin echo experiment, which is specifically designed to be
insensitive to the low frequency qubit frequency fluctuations, many quantum
algorithms or sequences of quantum gates are highly sensitive to small
deviations in the qubit frequency. While there exist several open-loop control
strategies for minimizing gate sensitivity to noise at different frequencies
[28, 29], there is always an added cost in the duration of the gate sequence
or complexity in calibration. The advantage of the feedback based
stabilization method is that no changes to the gate sequences or controls are
required.

###  II.3 Randomized benchmarking

Figure 4: Randomized benchmarking a) Randomized benchmarking of the single
qubit gates at $f_{\rm q}=$4.44\text{\,}\mathrm{GHz}$$. The orange (blue) dots
show the average of 50 realizations of the random Clifford sequences with
(without) the feedback. The shaded area shows the 68% confidence interval for
the average trajectory. The gate error is extracted from the fit to the data
(solid line). b) Randomized benchmarking is repeated 250 times with the
feedback (orange dots) and then without (blue dots). The shaded area shows the
68% confidence interval of for the fitted gate error. c) The qubit spectrum
and the operating points at which the gate errors are evaluated in d). d)
Randomized benchmarking at several different qubit bias fluxes. Black dots
show the estimated coherence limit for the gate fidelities, inferred from the
energy relaxation time $T_{1}$ of the qubit [30]. The dashed black line shows
a linear fit to the estimated coherence limits at the different operating
points.

Next, we demonstrate that the feedback protocol improves the single-qubit gate
fidelity in our device. We bias the qubit $400\text{\,}\mathrm{MHz}$ away from
the flux sweet spot so that it is highly sensitive to flux noise and perform
single-qubit randomized benchmarking [31] while stabilizing the qubit
frequency with the feedback. Fig. 4a shows the qubit excited state population
$p_{1}$ as a function of number of Clifford gates in a random sequence,
followed by a Clifford gate that would ideally bring the qubit back to the
ground state. With the feedback off (blue) and on (orange), the experiment is
repeated for 50 different random sequences, which are averaged and used to
find the average error per gate [30]. The feedback reduces the average error
per gate from $($8.5$\pm$2.1$)\times 10^{-4}$ to $($5.9$\pm$0.7$)\times
10^{-4}$, approaching the limit imposed by the energy-relaxation rate of the
qubit, see Supplementary material for additional information. While such
fidelities are commonly observed at or near the flux-insensitive point in many
devices [30, 32], it is unusual to see it so far from the sweet spot.

We attribute the reduction in the average gate error mostly to the improved
stability of the qubit frequency. Without feedback, several of the random
sequences show oscillating decay functions with the number of Clifford gates,
indicative of coherent control errors, such as the quasi-static shifts in the
qubit frequency due to the low frequency noise. In contrast, with the feedback
activated, the function monotonically decays for all random sequences,
something that is typical for incoherent errors [33, 34]. This observation is
supported by the experiment in Fig. 4b where we show the gate errors inferred
from randomized benchmarking experiments repeated 250 times with 7 different
realizations of a random gate sequence, for a total duration of
$840\text{\,}\mathrm{s}$. Without the feedback, we observe a significant drift
in the gate errors. With feedback turned on, the gate errors remain
consistently low, indicating that the feedback successfully stabilized the
drift and fluctuations in the qubit frequency. This is also manifested in
lowered uncertainty in the fitted gate error in Fig. 4a.

We observe a similar improvement on gate fidelity for all flux bias points
with the exception of the sweet spot, see Figs. 4c and d. At each flux point,
we repeated the randomized benchmarking experiment 10 times with and without
feedback, re-calibrating the qubit frequency between every repetition of the
randomized benchmarking experiment. In contrast to the coherence time $T_{2}$
— which is the highest at the sweet sport and then consistently decreases as
the flux noise sensitivity increases away from the sweet spot — the highest
gate fidelities are in fact measured at the most flux-sensitive point we
investigated. We attribute the increasing trend in the fidelities at lower
frequencies to the higher energy-relaxation time $T_{1}$ of the qubit.
Additionally, spectrally moving parasitic two-level fluctuators reduced the
best achievable single-qubit gate fidelities at certain bias points. Due to
the coupling to the qubit, the changes in the frequencies of the two-level
fluctuators caused shifts in the qubit frequency. As a result, at some bias
points there were wide variations in the measured gate fidelities between
different experiments, which are responsible for the large error bars in Fig.
4c when the feedback was not used. The variation is significantly reduced by
the feedback, as indicated by the repetition experiment shown in Fig. 4b and
further verified in Fig. 4c. This highlights that the implemented feedback
protocol can mitigate the impact of many sources of low-frequency noise, such
as frequency shifts caused by two-level fluctuators.

##  III Conclusions

In this work we have implemented a closed-loop feedback protocol to stabilize
the drift and fluctuations in the frequency of a superconducting transmon
qubit. In the probing phase we use repeated Ramsey experiments to estimate the
qubit frequency and adjust the qubit frequency to cancel the measured
frequency offset. The probing phase can be interleaved with a computational
workload such as algorithm execution.

We have demonstrated that the feedback stabilizes the qubit frequency
fluctuations even when the qubit is not operated at the noise-insensitive
operation point. This leads to reduction in the noise power observed by the
qubit, resulting in improved coherence times and improved gate fidelities. The
ability to operate qubits away from the protected bias point will help
addressing the frequency crowding problem in large quantum processors by
increasing the operable frequency band for the qubits. Additionally, the
increased operable frequency band helps to avoid spurious modes arising for
example from two-level fluctuators.

While in this work we have mostly focused on mitigating $1/f$-flux noise, the
feedback algorithm is agnostic to the source of the noise in the qubit
frequency. The proposed technique should work equally well on other low-
frequency noise sources, such as charge noise or TLS induced shifts in the
qubit frequency. Moreover, the feedback relies only on single-qubit operations
which implies that it can be extended to multi-qubit systems without
additional cost in time or complexity.

###### Acknowledgements.

This work was supported in part by the U.S. Army Research Office (ARO) Grant
W911NF-18-1-0411; by the ARO Multi-University Research Initiative
W911NF-18-1-0218; and by the Assistant Secretary of Defense for Research and
Engineering via MIT Lincoln Laboratory under Air Force Contract No.
FA8721-05-C-0002. A.K. acknowledges support from the NSF Graduate Research
Fellowship program.

## References

  * Steck  _et al._ [2004] D. A. Steck, K. Jacobs, H. Mabuchi, T. Bhattacharya, and S. Habib, Quantum feedback control of atomic motion in an optical cavity, Phys. Rev. Lett. 92, 223004 (2004).
  * Wilson  _et al._ [2015] D. J. Wilson, V. Sudhir, N. Piro, R. Schilling, A. Ghadimi, and T. J. Kippenberg, Measurement-based control of a mechanical oscillator at its thermal decoherence rate, Nature 524, 325 (2015).
  * Vijay  _et al._ [2012] R. Vijay, C. Macklin, D. Slichter, S. Weber, K. Murch, R. Naik, A. N. Korotkov, and I. Siddiqi, Stabilizing rabi oscillations in a superconducting qubit using quantum feedback, Nature 490, 77 (2012).
  * Shulman  _et al._ [2014] M. D. Shulman, S. P. Harvey, J. M. Nichol, S. D. Bartlett, A. C. Doherty, V. Umansky, and A. Yacoby, Suppressing qubit dephasing using real-time hamiltonian estimation, Nature communications 5, 1 (2014).
  * Nakajima  _et al._ [2020] T. Nakajima, A. Noiri, K. Kawasaki, J. Yoneda, P. Stano, S. Amaha, T. Otsuka, K. Takeda, M. R. Delbecq, G. Allison, A. Ludwig, A. D. Wieck, D. Loss, and S. Tarucha, Coherence of a driven electron spin qubit actively decoupled from quasi-static noise, Phys. Rev. X (2020).
  * Majumder  _et al._ [2020] S. Majumder, L. A. de Castro, and K. R. Brown, Real-time calibration with spectator qubits, npj Quantum Information 6, 1 (2020).
  * Yan  _et al._ [2016] F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, _et al._ , The flux qubit revisited to enhance coherence and reproducibility, Nature communications 7, 1 (2016).
  * Catelani  _et al._ [2011] G. Catelani, R. J. Schoelkopf, M. H. Devoret, and L. I. Glazman, Relaxation and frequency shifts induced by quasiparticles in superconducting qubits, Physical Review B 84, 064517 (2011).
  * Bylander  _et al._ [2011] J. Bylander, S. Gustavsson, F. Yan, F. Yoshihara, K. Harrabi, G. Fitch, D. G. Cory, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Noise spectroscopy through dynamical decoupling with a superconducting flux qubit, Nature Physics 7, 565 (2011).
  * Braumüller  _et al._ [2020] J. Braumüller, L. Ding, A. P. Vepsäläinen, Y. Sung, M. Kjaergaard, T. Menke, R. Winik, D. Kim, B. M. Niedzielski, A. Melville, _et al._ , Characterizing and optimizing qubit coherence based on squid geometry, Physical Review Applied 13, 054079 (2020).
  * McEwen  _et al._ [2021] M. McEwen, L. Faoro, K. Arya, A. Dunsworth, T. Huang, S. Kim, B. Burkett, A. Fowler, F. Arute, J. C. Bardin, _et al._ , Resolving catastrophic error bursts from cosmic rays in large arrays of superconducting qubits, arXiv preprint arXiv:2104.05219 (2021).
  * Wilen  _et al._ [2020] C. Wilen, S. Abdullah, N. Kurinsky, C. Stanford, L. Cardani, G. d’Imperio, C. Tomei, L. Faoro, L. Ioffe, C. Liu, _et al._ , Correlated charge noise and relaxation errors in superconducting qubits, arXiv preprint arXiv:2012.06029 (2020).
  * Vepsäläinen  _et al._ [2020] A. P. Vepsäläinen, A. H. Karamlou, J. L. Orrell, A. S. Dogra, B. Loer, F. Vasconcelos, D. K. Kim, A. J. Melville, B. M. Niedzielski, J. L. Yoder, _et al._ , Impact of ionizing radiation on superconducting qubit coherence, Nature 584, 551 (2020).
  * Koch  _et al._ [2007] J. Koch, M. Y. Terri, J. Gambetta, A. A. Houck, D. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the cooper pair box, Physical Review A 76, 042319 (2007).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369 (2020).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Barends  _et al._ [2014] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, _et al._ , Superconducting quantum circuits at the surface code threshold for fault tolerance, Nature 508, 500 (2014).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Applied Physics Reviews 6, 021318 (2019).
  * Klimov  _et al._ [2018] P. Klimov, J. Kelly, Z. Chen, M. Neeley, A. Megrant, B. Burkett, R. Barends, K. Arya, B. Chiaro, Y. Chen, _et al._ , Fluctuations of energy-relaxation times in superconducting qubits, Physical review letters 121, 090502 (2018).
  * de Graaf  _et al._ [2020] S. de Graaf, L. Faoro, L. Ioffe, S. Mahashabde, J. Burnett, T. Lindström, S. Kubatkin, A. Danilov, and A. Y. Tzalenchuk, Two-level systems in superconducting quantum devices due to trapped quasiparticles, arXiv preprint arXiv:2004.02485 (2020).
  * Burnett  _et al._ [2019] J. J. Burnett, A. Bengtsson, M. Scigliuzzo, D. Niepce, M. Kudra, P. Delsing, and J. Bylander, Decoherence benchmarking of superconducting qubits, npj Quantum Information 5, 1 (2019).
  * Schlör  _et al._ [2019] S. Schlör, J. Lisenfeld, C. Müller, A. Bilmes, A. Schneider, D. P. Pappas, A. V. Ustinov, and M. Weides, Correlating decoherence in transmon qubits: Low frequency noise by single fluctuators, Physical review letters 123, 190502 (2019).
  * Wallraff  _et al._ [2005] A. Wallraff, D. Schuster, A. Blais, L. Frunzio, J. Majer, M. Devoret, S. Girvin, and R. Schoelkopf, Approaching unit visibility for control of a superconducting qubit with dispersive readout, Physical review letters 95, 060501 (2005).
  * Rol  _et al._ [2017] M. Rol, C. Bultink, T. O’Brien, S. De Jong, L. Theis, X. Fu, F. Luthi, R. Vermeulen, J. de Sterke, A. Bruno, _et al._ , Restless tuneup of high-fidelity qubit gates, Physical Review Applied 7, 041001 (2017).
  * Werninghaus  _et al._ [2020] M. Werninghaus, D. Egger, and S. Filipp, High-speed calibration and characterization of superconducting quantum processors without qubit reset, arXiv preprint arXiv:2010.06576 (2020).
  * Yan  _et al._ [2012] F. Yan, J. Bylander, S. Gustavsson, F. Yoshihara, K. Harrabi, D. G. Cory, T. P. Orlando, Y. Nakamura, J.-S. Tsai, and W. D. Oliver, Spectroscopy of low-frequency noise and its temperature dependence in a superconducting qubit, Physical Review B 85, 174521 (2012).
  * Ithier  _et al._ [2005] G. Ithier, E. Collin, P. Joyez, P. Meeson, D. Vion, D. Esteve, F. Chiarello, A. Shnirman, Y. Makhlin, J. Schriefl, _et al._ , Decoherence in a superconducting quantum bit circuit, Physical Review B 72, 134519 (2005).
  * Sheldon  _et al._ [2016] S. Sheldon, E. Magesan, J. M. Chow, and J. M. Gambetta, Procedure for systematically tuning up cross-talk in the cross-resonance gate, Physical Review A 93, 060302 (2016).
  * Milne  _et al._ [2020] A. R. Milne, C. L. Edmunds, C. Hempel, F. Roy, S. Mavadia, and M. J. Biercuk, Phase-modulated entangling gates robust to static and time-varying errors, Physical Review Applied 13, 024022 (2020).
  * Chen [2018] Z. Chen,  _Metrology of quantum control and measurement in superconducting qubits_ , Ph.D. thesis, UC Santa Barbara (2018).
  * Knill  _et al._ [2008] E. Knill, D. Leibfried, R. Reichle, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, Randomized benchmarking of quantum gates, Physical Review A 77, 012307 (2008).
  * Sung  _et al._ [2020] Y. Sung, L. Ding, J. Braumüller, A. Vepsäläinen, B. Kannan, M. Kjaergaard, A. Greene, G. O. Samach, C. McNally, D. Kim, _et al._ , Realization of high-fidelity cz and zz-free iswap gates with a tunable coupler, arXiv preprint arXiv:2011.01261 (2020).
  * Sanders  _et al._ [2015] Y. R. Sanders, J. J. Wallman, and B. C. Sanders, Bounding quantum gate error rate based on reported average fidelity, New Journal of Physics 18, 012002 (2015).
  * Greene  _et al._ [2021] A. Greene, M. Kjaergaard, M. Schwartz, G. Samach, A. Bengtsson, M. O’Keeffe, D. Kim, M. Marvian, A. Melville, B. Niedzielski, _et al._ , Error mitigation via stabilizer measurement emulation, arXiv preprint arXiv:2102.05767 (2021).
  * Danilin  _et al._ [2018] S. Danilin, A. V. Lebedev, A. Vepsäläinen, G. B. Lesovik, G. Blatter, and G. Paraoanu, Quantum-enhanced magnetometry by phase estimation algorithms with a single artificial atom, npj Quantum Information 4, 1 (2018).
  * Oppenheim  _et al._ [1999] A. Oppenheim, R. Schafer, and J. Buck,  _Discrete-time Signal Processing_, Prentice Hall International Editions Series (Prentice Hall, 1999).
  * O’Malley  _et al._ [2015] P. O’Malley, J. Kelly, R. Barends, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, I.-C. Hoi, _et al._ , Qubit metrology of ultralow phase noise using randomized benchmarking, Physical Review Applied 3, 044009 (2015).

## Supplementary material: Improving qubit coherence using closed-loop
feedback

###  .1 Qubit frequency estimation

The qubit frequency can be estimated using a simple version of the single
qubit phase estimation algorithm where the qubit is first prepared in a
superposition state with a $\pi/2$ rotation around $Y$ axis followed by free
evolution. During the free evolution, the qubit wavefunction accumulates a
phase (in the frame rotating with the drive)
$\phi(t)=\int_{0}^{t}2\pi\delta_{\rm q}(t)\mathrm{d}t$, where $\delta_{\rm
q}(t)=f_{\rm d}-f_{\rm q}(t)$ and $f_{\rm d}$ and $f_{\rm q}(t)$ are the drive
frequency and the qubit frequency. The accumulated phase can be estimated by
measuring either $\langle X\rangle$ or $\langle Y\rangle$. During the
experiment, the qubit state evolves as

$|\psi\rangle=\exp{\left(-\frac{i\pi}{4}M\right)}\exp{\left(-i\frac{\pi}{2}\phi(\tau)Z\right)}\exp{\left(-\frac{i\pi}{4}Y\right)}|\psi_{0}\rangle,$
(S1)

where $|\psi_{0}\rangle$ is the initial state of the qubit, $\tau$ is the
free-evolution time, $M$ is the measurement operator, and
$\phi(\tau)=2\pi\tau\delta_{\rm q}$ if the detuning $\delta_{\rm q}$ is
assumed to be constant during the experiment. The probability of measuring the
excited state is given by

$p_{1}=\frac{1-\langle\psi|Z|\psi\rangle}{2}=\frac{1}{2}+\frac{1}{2}\cos(2\pi\delta_{\rm
q}\tau-\phi_{\rm m}),$ (S2)

where $\phi_{\rm m}$ is a phase factor that depends on the measurement
operator,

$\phi_{\rm m}=\begin{cases}0,&\text{if}\ M=X,\\\ \pi/2,&\text{if}\
M=Y.\end{cases}$ (S3)

Solving for $\delta_{\rm q}$ yields

$\delta_{\rm q}=\frac{\pm\arccos(2p_{1}-1)+\phi_{\rm m}+2\pi k}{2\pi\tau}.$
(S4)

By choosing $k=0$ and the negative branch $\delta_{\rm q}$ is uniquely defined
in the range $\delta_{\rm q}=[\frac{\phi_{\rm
m}-\pi}{2\pi\tau},\frac{\phi_{\rm m}}{2\pi\tau}]$, which is symmetric around
zero for $\phi_{\rm m}=\frac{\pi}{2}$.

The qubit excited state probability $\hat{p}_{1}$ can be estimated from $N$
measurements of the qubit state $q_{i}$ through
$\hat{p}_{1}=\frac{1}{N}\sum_{i=0}^{N}q_{i}$, which then gives an estimate of
the qubit frequency $\hat{\delta}_{\rm q}$. However, this frequency estimate
is noisy due to finite number of samples used to estimate the excited state
probability. The impact of the noise can be estimated by noting that
$\hat{p}_{1}$ is binomially distributed, and its standard deviation is

$\delta\hat{p}_{1}=\sqrt{\frac{p_{1}(1-p_{1})}{N}}=\frac{1}{2\sqrt{N}},\,\,\mathrm{for}\,\,p=1/2.$
(S5)

The error in the frequency estimate is then given by

$\delta\hat{\delta}_{\rm q}=\left|\frac{\partial\delta_{\rm q}}{\partial
p_{1}}\delta\hat{p}_{1}\right|\approx\left|\frac{\partial}{\partial
p_{1}}\frac{\frac{\pi}{2}+\phi_{\rm
m}+2p_{1}-1}{2\pi\tau}\delta\hat{p}_{1}\right|=\frac{1}{2\pi\tau\sqrt{N}},$
(S6)

where $\delta_{\rm q}$ has been expanded to the first order in $p_{1}$. This
implies that the estimation noise can be reduced either by increasing $\tau$
or the number of measurements used in the estimation. However, the bandwidth
of the estimation algorithm is limited by the requirement that each frequency
has to be uniquely given by $p_{1}$. This is true if $|\delta_{\rm
q}\tau|\leq\frac{1}{4}$. As a result, there is a trade-off between the
estimation bandwidth and sensitivity. Decoherence further reduces the
sensitivity, and the optimum is obtained at $\tau=T_{2}$ [35]. In practice,
$\tau$ often needs to be shorter to satisfy the bandwidth requirement. Using a
larger $N$ increases the total time taken by the frequency estimation, which
reduces the repetition rate of the feedback algorithm.

###  .2 Circuit analysis

Here we theoretically analyze the implemented feedback circuit and assess its
performance. We start by writing the block diagram for the circuit, shown in
S1a. In the diagram and the following calculation, the quantities are sampled
at times $t=nT_{\rm N}=nNT$, where $T$ is the total duration of a Ramsey
experiment and $N$ is the number of Ramsey experiments used to estimate the
qubit frequency. The target qubit frequency is controlled by $d[n]=f_{\rm d}$
which we keep fixed at all times. At every cycle of the feedback, the target
frequency $d[n]$ is compared to the estimated qubit frequency $f[n]$,
resulting in the error signal $e[n]=d[n]-(f[n]+v[n])$, where $v[n]$ is the
sampling noise. In the experiment, the error signal is given by
$e[n]=\hat{\delta}_{\rm q}[n]$, and the sampling noise
$v[n]=\delta\hat{\delta}_{\rm q}[n]$ is defined in Eq. (S6). The error signal
is multiplied by a controllable gain $G$, and is fed to an accumulator
outputting the feedback control signal $p[n-1]$ — delayed by one time step —
which then adjusts the qubit frequency.

The sampled qubit frequency $f[n]$ is estimated from the real qubit frequency
$f_{\rm q}(t)$ using $N$ Ramsey measurements as described in the main text.
Here, we approximate the sampled frequency as the average qubit frequency
during the sampling

$f[n]\approx\frac{1}{N\tau}\sum_{i=0}^{N-1}\int_{nT_{\rm N}+iT}^{nT_{\rm
N}+iT+\tau}f_{\rm q}(t)\mathrm{d}t.$ (S7)

The qubit frequency is determined by its intrinsic fluctuating value
$\tilde{f}_{\rm q}$ which is acted on by the feedback signal $p[n-1]$ so that
$f_{\rm q}(t)=\tilde{f}_{\rm q}(t)+p[n-1]$. Substituting to Eq. (S7) yields

$f[n]\approx p[n-1]+\frac{1}{N\tau}\sum_{i=0}^{N-1}\int_{nT_{\rm
N}+iT}^{nT_{\rm N}+iT+\tau}\tilde{f}_{\rm q}(t)\mathrm{d}t\equiv
p[n-1]+\tilde{f}[n],$ (S8)

forming the feedback loop.

To analyze the impact of the feedback on the signal and the sampling noise, we
start by solving the transfer function of the feedback signal $p[n]$ from

$\displaystyle e[n]$ $\displaystyle=d[n]-v[n]-\tilde{f}[n]-p[n-1],$ (S9)
$\displaystyle p[n]$ $\displaystyle=p[n-1]+y[n],$ $\displaystyle y[n]$
$\displaystyle=Ge[n],$ $\displaystyle\implies p[n]$
$\displaystyle=p[n-1]-G(-d[n]+v[n]+\tilde{f}[n]+p[n-1]).$

Figure S1: a) The diagram of the feedback circuit. b) The squared magnitude of
the frequency responses for the control and error signals $p[n]$ and $e[n]$.
c) The squared magnitude of the frequency response for the qubit frequency
during the feedback. The dashed line shows the frequency response of the qubit
to the sampling noise, and the dash-dotted line shows the frequency response
to the fluctuations in the intrinsic qubit frequency. d) The power spectral
density of the qubit frequency fluctuation with the feedback activated (red)
calculated using the frequency responses shown in c) and the experimentally
measured spectral density of the intrinsic qubit frequency fluctuation
(green). The dashed line shows the sampling noise.

Applying the $z$-transform yields

$\displaystyle P(z)$
$\displaystyle=z^{-1}P(z)-G(-D(z)+V(z)+\tilde{F}(z)+z^{-1}P(z)),$ (S10)
$\displaystyle P(z)$
$\displaystyle=\frac{G(D(z)-V(z)-\tilde{F}(z))}{1-z^{-1}+z^{-1}G},$

resulting in the transfer function

$X_{\rm p}(z)=P(z)/(D(z)-V(z)-\tilde{F}(z))=\frac{G}{1-z^{-1}+z^{-1}G}.$ (S11)

To get the frequency response $H_{\rm p}(f)$ from the transfer function in Eq.
(S11), we substitute $z=\mathrm{e}^{i2\pi fNT}$ so that $H_{\rm p}(f)=X_{\rm
p}(\mathrm{e}^{i2\pi fNT})$, which is shown in Fig. S1b for $G=0.35$. The
power spectral density of the feedback signal as a response to the frequency
fluctuations of the system is given by the squared magnitude of the frequency
response [36]

$S_{\rm pp}(f)=\left|H_{\rm p}(f)\right|^{2}S_{\rm ss}(f),$ (S12)

where $S_{\rm ss}(f)$ is the combined power spectral density of the target
signal $d[n]$, the sampled intrinsic qubit frequency $\tilde{f}[n]$ and the
sampling noise $v[n]$, $S_{\rm ss}(f)=S_{\rm
dd}(f)+S_{\rm\tilde{f}\tilde{f}}(f)+S_{\rm vv}(f)$. We assume that there is
only a minimal amount of noise in frequency of the microwave source which sets
the target frequency, implying that we can set $S_{\rm dd}\approx 0$. To
maximize the efficiency of the feedback, we limit the bandwidth of the
frequency response to the frequencies where the spectral density of the qubit
frequency fluctuations is higher than the sampling noise. Another option would
be to increase number of phase sampling experiments $N$ per frequency
estimate, which reduces the sampling noise power at the cost of a longer total
time used for feedback.

To better understand the experimentally measured spectrum in Fig. 1d of the
main text, we next calculate the frequency response of the error signal $e[n]$
from Eqs. (S9). Applying the $z$-transform and solving for $E(z)$ yields

$\displaystyle E(z)$ $\displaystyle=D(z)-V(z)-\tilde{F}(z)-z^{-1}P(z),$ (S13)
$\displaystyle P(z)$ $\displaystyle=z^{-1}P(z)+GE(z),$ $\displaystyle\implies
E(z)$
$\displaystyle=\frac{1-z^{-1}}{1-z^{-1}+Gz^{-1}}(D(z)-V(z)-\tilde{F}(z)),$
$\displaystyle\implies X_{\rm e}(z)$
$\displaystyle=\frac{1-z^{-1}}{1-z^{-1}+Gz^{-1}}.$

The frequency response $H_{\rm e}(f)=X_{\rm e}(\mathrm{e}^{i2\pi fNT})$ is
plotted in Fig. S1b. Without the feedback, the expected spectral density of
the error signal is $S_{\rm ee}^{\rm
nofeedback}(f)=S_{\rm\tilde{f}\tilde{f}}(f)+S_{\rm vv}(f)$. When the feedback
is turned on, the response is $S_{\rm ee}^{\rm
feedback}=(S_{\rm\tilde{f}\tilde{f}}(f)+S_{\rm vv}(f))|H_{\rm e}(f)|^{2}$
instead. At lower frequencies, the frequency response of the error signal goes
down, explaining the observed behavior in the spectra in Fig. 1b in the main
text. Because the sampling noise in the frequency estimation affects the qubit
frequency through the feedback, the spectral density of the error signal does
not exactly match the spectral density of the real qubit frequency.

We now analyze the impact of the feedback directly on the noise spectral
density of the qubit frequency. We first calculate the transfer function
$X_{\rm f}(z)=\frac{F(z)}{D(z)-V(z)-\tilde{F}(z)}$ of the qubit frequency
$f[n]$ right after a feedback cycle when $f[n]=p[n]+\tilde{f}[n]$. Using the
result of Eq. (S10) we get

$\displaystyle X_{\rm f}(z)$
$\displaystyle=\frac{\tilde{F}(z)}{D(z)-V(z)-\tilde{F}(z)}+X_{\rm p}(z)$ (S14)
$\displaystyle\implies X_{\rm f}(z)$ $\displaystyle=\begin{cases}X_{\rm
p}(z)-1,&{\rm for}\,V(z)=0\,\text{and}\,D(z)=0\\\ X_{\rm p}(z),&{\rm
for}\,\tilde{F}(z)=0\,\text{and}\,D(z)=0.\\\ \end{cases}$

The squared magnitude of the frequency responses of these transfer functions,
$|H_{\rm f}^{\tilde{F(z)=0}}(f)|^{2}$ and $|H_{\rm
f}^{\tilde{V(z)=0}}(f)|^{2}$, are plotted in Fig. S1c. The spectral density of
the qubit frequency fluctuations is then $S_{\rm ff}(f)=|H_{\rm
f}^{\tilde{F}(z)=0}(f)|^{2}S_{\rm vv}(f)+|H_{\rm
f}^{\tilde{V}(z)=0}(f)|^{2}S_{\rm\tilde{f}\tilde{f}}(f)$, assuming we can
independently sum the contributions from the sampling noise and the qubit
frequency fluctuations, shown in Fig. S1d for the same parameters as used in
Fig. 1d in the main text. Quite remarkably, this analysis very closely
reproduces the simulation of the qubit frequency fluctuation spectral density
shown with red in Fig. 1d.

###  .3 Spectral density of flux noise

The measured coefficients $k_{\rm E/R}$ can be connected to the spectral
density of $1/f$ flux noise, $S_{\Phi\Phi}(f)=A_{\Phi}/|2\pi f|$, as

$\sqrt{A_{\rm\Phi}}=\frac{k_{\rm E/R}}{2\pi\sqrt{\eta_{\rm E/R}}},$ (S15)

where the labels E or R refer either to echo or Ramsey sequence, and
$\eta_{\rm E/R}$ is a scaling parameter which depends on the bandwidth of the
noise to which the used pulse sequence is sensitive, $\eta_{\rm E}={\rm ln}2$,
and $\eta_{\rm R}={\rm ln}\frac{f_{\rm u}}{2\pi f_{\rm l}}\approx 13$ [27].
Here, we take the upper cutoff frequency $f_{\rm
u}=$100\text{\,}\mathrm{kHz}$\sim 1/T_{2}$ and the lower cutoff frequency
$f_{\rm l}=$40\text{\,}\mathrm{mHz}$$ corresponding to the inverse of the
total duration of the Ramsey experiment. By substituting the measured values
of $k$ to the above formula, we infer
$\sqrt{A_{\phi}}=$2.8\text{\,}\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$
without the feedback. With the feedback activated, the flux noise spectral
density no longer follows $1/f$ law due to the suppression of low frequency
noise, but for comparison we can still calculate an effective flux noise
amplitude using Eq. (S15), which yields
$\sqrt{A_{\phi}}=$2.4\text{\,}\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$.
When measured using the echo sequence we get
$\sqrt{A_{\phi}}=$3.3\text{\,}\mathrm{\SIUnitSymbolMicro}\Phi_{\mathrm{0}}$$.
Ideally, the flux noise amplitude inferred from the Ramsey experiment without
the feedback should yield the same value as the flux noise amplitude
calculated from the echo experiment, but we attribute the discrepancy to the
deviation of the noise spectrum from the ideal $1/f$ spectrum, as observed
earlier in Fig. 1b. The measured flux noise amplitude is consistent with the
earlier systematic study we performed for flux noise amplitudes of SQUID loops
of different sizes [10].

###  .4 Experimental setup and the hardware implementation

Figure S2: The experimental setup.

Fig. S2 shows the hardware setup used in the experiment. The qubit $XY$
control and readout pulses pulses are created by a Keysight PXI arbitrary
waveform generator. The in-phase and quadrature pulses are up-converted to the
qubit and the resonator frequencies using a R&S SGMA RF source, which has a
built-in IQ-mixer used for single-sideband modulation. The signals are
attenuated to mitigate thermal noise from room-temperature, and sent to the
qubit installed in a Leiden Cryogenics dilution refrigerator. After
interacting with the readout resonator, the readout signal is amplified by a
traveling wave parametric amplifier (TWPA), which is driven by a pump tone.
The TWPA is followed by a chain of conventional amplifiers and the signal is
digitized at room-temperature using a Keysight PXI digitizer.

The qubit flux bias is controlled by two channels of another Keysight PXI AWG.
The first channel is used to set the operating point for the qubit, whereas
the second channel is controlled by the feedback loop. The feedback channel is
heavily attenuated in order to use the full voltage bandwidth of the AWG to
reduce the discretization noise.

Figure S3: The block diagram of the two FPGAs performing the feedback. The
FGPA in the digitizer discriminates the state, and forwards the information to
the FPGA which controls the AWG.

The feedback loop is formed by the readout digitizer and the qubit flux
control AWG. The on-board FPGA of the digitizer is used to digitally
demodulate and integrate the readout signal and then use a threshold
discriminator to assign the qubit state $q_{i}$ to either $0$ or $1$, see Fig.
S3 for the block diagram of the FPGA operation. The discriminated states are
transferred to the AWG for the calculation of the feedback signal.

First we perform the virtual qubit state reset by flipping the discriminated
state $q_{i}$ if the previous discriminated state was excited state. This is
equivalent to $q_{i}\rightarrow{\rm XOR}(q_{i},q_{i-1})$. In the next step,
the $q_{i}$ are used to calculate the estimate for the qubit state
$s_{i}=N\hat{p}_{1}=\sum_{j=i-N+1}^{i}q_{j}=s_{i-1}-q_{i-N+1}+q_{i}$. The
value of the buffer sum $s$ is used as an index for a lookup table which
contains the binary representation of the voltages corresponding to each of
the possible $N+1$ values for $s_{i}$. In the last step, the voltage values
are fed to an accumulator, which value is updated every $N_{\rm S}$ steps. We
set $N_{\rm S}=N$ for the measurement of the qubit frequency fluctuation
spectral density, and $N_{\rm S}=N+1$ for the interleaved operation.

###  .5 Coherence limit in randomized benchmarking

The limit imposed on the gate fidelity by decoherence has been studied for
example in [37, 30]. The error probability per gate due to decoherence depends
on the gate time and the coherence times as

$\epsilon=\frac{t_{\rm gate}}{3T_{1}}+\frac{t_{\rm gate}}{3T_{\phi
1}}+\frac{1}{3}\left(\frac{t_{\rm gate}}{T_{\phi 2}}\right)^{2},$ (S16)

where $T_{1}$ is the energy-relaxation rate of the qubit, $T_{\phi 1}$ is the
exponential part of the pure dephasing time, and $T_{\phi 2}$ is the Gaussian
part of the pure dephasing time. In the experiment, we have $t_{\rm
gate}=$40\text{\,}\mathrm{ns}$$, which consists of $30\text{\,}\mathrm{ns}$
long cosine-shaped envelope and $10\text{\,}\mathrm{ns}$ delay between the
pulses. The $T_{1}$ times of our qubit were between
$26\text{\,}\mathrm{\SIUnitSymbolMicro s}$ and
$40\text{\,}\mathrm{\SIUnitSymbolMicro s}$, depending on the bias point used
in the experiment. Using Eq. (S16), we calculate $\epsilon$ at all the 11 bias
points used in the experiment, yielding values in the range of
$3\text{\times}{10}^{-4}\text{\,}5\text{\times}{10}^{-4}\text{\,}$, which are
shown with black dots in Fig. 4c in the main text. The gate errors due to the
Gaussian dephasing $\frac{1}{3}\left(\frac{t_{\rm gate}}{T_{\phi
2}}\right)^{2}<0.33\times 10^{-4}$ contribute to the total error rate
significantly less than the energy-relaxation rate, even at the most sensitive
bias point where $T_{\rm\phi 2}\approx$4\text{\,}\mathrm{\SIUnitSymbolMicro
s}$$. Feedback improved the dephasing time to $T_{\rm\phi
2}\approx$5\text{\,}\mathrm{\SIUnitSymbolMicro s}$$, which according to the
estimate would reduce the Gaussian dephasing induced error rate to $\sim
0.2\times 10^{-4}$. However, the measured decrease in gate errors was
significantly higher, from $($8.5$\pm$2.1$)\times 10^{-4}$ to
$($5.9$\pm$0.7$)\times 10^{-4}$. We attribute the higher efficiency of the
feedback to its ability to correct slow frequency fluctuations of the qubit,
not captured by a single $T_{2}$ experiment.