# Selective nuclear-spin interaction based on a dissipatively stabilized
nitrogen-vacancy center

Jiawen Jiang and Q<EMAIL_ADDRESS>Department of Physics and
Key Laboratory of Low Dimensional Quantum Structures and Quantum Control of
Ministry of Education, Hunan Normal University, Changsha 410081, China.

###### Abstract

Current typical methods to realize nuclear-nuclear quantum gates require a
sequence of electron-nuclear quantum gates by using dynamical decoupling
techniques, which are implemented at low temperature because of short
decoherence and relaxation time of the NV spin at room temperature. This
limitation could be overcome by using periodical resets of an NV spin as a
mediator of interaction between two nuclear spins [Chen, Schwarz, and Plenio,
119, 010801 (2017)]. However, this method works under stringent coupling
strengths condition, which makes it not applicable to heteronuclear quantum
gate operations. Here we develop this scheme by using radio-frequency (RF)
fields to control different nuclear spin species. Periodical resets of the NV
center protect the nuclear spins from decoherence and relaxation of the NV
spin. RF control provides probability to have highly selective and high
fidelity quantum gates between heteronuclear spins as well as detecting
nuclear spins by using a nuclear spin sensor under ambient conditions.

##  I Introduction

Nuclear spins in material such as diamond associated with single defects
represent a promising platform for quantum information registers
Taminiau2014UniversalCA ; 2018Pulse ; Waldherr2014QuantumEC ; 2011High  and
sensing purposes Taminiau2012DetectionAC ; Kolkowitz2012SensingDN ;
Zhao2012SensingSR ; Ermakova2013DetectionOA ; Mamin2013NanoscaleNM ;
Zaiser2016EnhancingQS ; Staudacher2013NuclearMR ; Wu2016DiamondQD  due to
their long coherence times and the potentially large number of available
spins. Nuclear spins can be initialized, controlled, and read out through the
electron spin of the nitrogen vacancy (NV) center Doherty2013TheNC  driven by
optical fields and microwave radiation 2016Optomechanical ; 2014Optically .
Recent progress manifests in a number of works, such as electron-nuclear
Casanova2016NoiseResilientQC ; Bermdez2011ElectronmediatedNI ;
zimmermann2020selective ; rong2014implementation ; hegde2020efficient ;
abobeih2019atomic ; Bradley2019ATS ; Wang2017DelayedEE ;
Tratzmiller2021ParallelSN ; Casanova2017ArbitraryNG , electron-electron
degen2021entanglement  and nuclear-nuclear quantum gates Bian2017UniversalQC .
Many schemes are proposed to realize electron-nuclear quantum gate operations,
both with Bradley2019ATS ; Wang2017DelayedEE  and without
Tratzmiller2021ParallelSN ; Casanova2017ArbitraryNG  additional radio-
frequency (RF) control on the nuclear spins themselves.

Nuclear-nuclear quantum gate is implemented by a sequence of the NV-nuclear
quantum operations Taminiau2014UniversalCA ; 2018Pulse , and it is a delicate
issue to have a complete set of quantum gates on specific nuclei in samples
Bradley2019ATS . Current related experiments are operated at low temperature,
because of the relatively short lifetime of the NV center. Several outstanding
challenges are caused by the relaxation and decoherence processes of the
electron spin as these limit quantum gate fidelities on nuclear registers as
well as spectral resolution and selectivity. A scheme was proposed to have
highly selective and high fidelity quantum gates between nuclear spins under
ambient conditions Chen2017DissipativelySQ . However, it has strict
requirements of resonance conditions and is applicable for special cases,
i.e., for nuclear spins in the same species with the parallel coupling
components between the NV and nuclear spins far smaller than the vertical
components.

Here we extend the scheme in Ref. Chen2017DissipativelySQ  to a heteronuclear
case by using radio-frequency fields to control different nuclear spin species
individually. The effective substantial second-order coupling between the
nuclear spins obtained through a MW driving NV center which is periodically
reinitialized by a dissipative process Chen2017DissipativelySQ . Similarly,
the periodical reinitialization of the NV center decouple it from the dynamics
and its effect on system is an effective weak dissipation process. Thus high
selectivity and fidelity of nuclear-spin quantum gate for different nuclear
spin species could be possible even at ambient condition, which is an
extension for quantum computation and simulation applications by using
different nuclear species controlled by NV centers. Additionally, one can use
a nuclear spin as a quantum sensor to detect nuclear spins in another species
as well as analysis of complex spin structures.

This paper is structured as follows. We start with the derivation of effective
Hamiltonian of the system by using Schrieffer-Wolf
transformationKessler2012GeneralizedSF ; Bravyi2011SchriefferWolffTF , which
demonstrates the validity of the indirect interaction between nuclear spins by
the application of suitably tuned RF fields. We take two relevant species,
i.e., carbon-13 and silicon-29 spins as examples to show our protocol
efficiency. Coherent evolution between a carbon-13 and a silicon-29 spin is
limited by the NV lifetime, this limitation is overcome through periodical
resets of NV spin. The feasibility of high selectivity and fidelity nuclear-
nuclear quantum gates by using RF controls is discussed as well as the sensing
application. Finally, we compare our approach to the previous scheme in Ref.
Chen2017DissipativelySQ .

##  II Indirect interaction between nuclear spins

We consider two nuclear spins in different species are coupled to an electron
spin of single NV center, and microwave (MW) field and radio-frequency (RF)
field are used for external control over the electron and nuclear spins as
well as for achieving selective internuclear interactions, see Fig. 1(a). The
magnetic field $B$ is applied along the NV axis (the $\hat{z}$ axis), which is
large enough to split the degenerate states of $|m_{s}=|\pm 1\rangle$. The
Hamiltonian of the whole system is given by

$\displaystyle H$ $\displaystyle=$ $\displaystyle
DS_{z}^{2}+\gamma_{e}BS_{z}+\sum_{i}\gamma_{n_{i}}BI_{i}^{z}$
$\displaystyle+S_{z}\sum_{i}\overrightarrow{A_{i}}\cdot\overrightarrow{I_{i}}+H_{w}+H_{r}.$

Note that we hereafter omit the Dirac constant $\hbar$ for simplicity. Here
$D=(2\pi)2.87$ GHz denotes zero-field splitting of the electronic ground
state, $\gamma_{e}$ and $\gamma_{n_{i}}$ are the gyromagenetic ratio of the
electron spin and nuclear spins, respectively. The interaction between the NV
center and the $i$th nucleus is mediated by the hyperfine vector
$\overrightarrow{A_{i}}$, $A_{i}=(a_{\parallel i},a_{\perp i})$ with
$a_{\parallel i}$ and $a_{\perp i}$ denotes the related coupling parallel and
perpendicular components to the nuclear spin quantization axes $a_{\parallel
i}=\overrightarrow{A}_{i}\cdot\hat{z}$ and $a_{\perp
i}=\sqrt{|\overrightarrow{A}_{i}|^{2}-a^{2}_{\parallel i}}$. Hamiltonians
$H_{w}$ and $H_{r}$ describe the action of the MW driving and RF fields,
respectively, which reads

$\displaystyle H_{w}$ $\displaystyle=$ $\displaystyle\sqrt{2}\Omega\cos\omega
tS_{x},$ (2) $\displaystyle H_{r}$ $\displaystyle=$
$\displaystyle\sum_{i}2\Omega_{rfi}\cos\omega_{rfi}tI_{i}^{x},$ (3)

where $\Omega$ ($\Omega_{rfi}$) is the Rabi frequency of the MW (RF) driving
field with the corresponding frequency $\omega$ ($\omega_{rfi}$).

The MW driving field is applied to be on resonance of the transition
$|-1\rangle\leftrightarrow|0\rangle$ as $\omega=D-\gamma_{e}B$, we can rewrite
the effective Hamiltonian of the NV center as $H_{NV}=\Omega\sigma_{z}$. Then
in the rotating frame with $H_{0}=(D-\gamma_{e}B)|-1\rangle\langle-1|$, the
Hamiltonian could be rewritten as

$\displaystyle H_{t}$ $\displaystyle=$
$\displaystyle\Omega\sigma_{z}+\sum_{i}(\gamma_{ni}B+\frac{a_{\parallel
i}}{2})I^{z}_{i}$ (4) $\displaystyle+a_{\parallel
i}\sigma_{x}I^{z}_{i}+a_{\perp i}\sigma_{x}I^{x}_{i}+H_{r},$

where $\sigma_{z}=\frac{1}{2}(|+\rangle_{e}\langle+|-|-\rangle_{e}\langle-|)$
with $|+\rangle_{e}=\frac{1}{\sqrt{2}}(|0\rangle+|-1\rangle)$ and
$|-\rangle_{e}=\frac{1}{\sqrt{2}}(|0\rangle-|-1\rangle)$.

Figure 1: (a) The NV center mediates the coupling between a silicon-29 spin
with $\gamma_{n_{1}}B=(2\pi)4$ MHz and a carbon-13 spin
$\gamma_{n_{2}}B=(2\pi)5.06$ MHz to achieve a quantum gate, while itself is
decoupled from the dynamics. RF fields are applied to control two nuclear spin
species individually, when we use continual MW field to drive the NV spin. (b)
Population evolutions are shown by considering the ideal case in which no
dissipation is included and the nuclear spins is initialized in state
$|+_{1}-_{2}\rangle$. The Rabi frequencies of RF fields are given by
$\Omega_{rf1}=\Omega_{rf2}=(2\pi)1$ kHz, microwave driving field
$\Omega=(2\pi)400$ kHz, the parallel coupling components $[a_{\parallel
1},a_{\parallel 2}]=(2\pi)[9,11]$ kHz and the effective coupling coefficient
$g_{e}=(2\pi)0.12$ kHz. These exact numerical simulations by using Eq. (5) fit
well to the theoretical derivation of the effective (eff) dynamics under
Hamiltonian Eq. (7). (c) By considering the practical case, we use the same
parameters in (b) and the NV dissipation is included as $T_{1\rho}=200$ $\mu
s$, exact numerical simulation of the master equation Eq. (10) shows that
there is no coherent evolution between the nuclear spins due to the NV life
time limitation. (d) The periodical NV resets makes it possible to extend the
coherent evolution of the heteronuclear spins well beyond the NV spin
lifetime. The exact numerical calculation is based on Eq. (17), when we
simulate the effective master equation Eq. (14) of nuclear spins. The
parameters are as the same as in (b) except that the NV center is reinitialied
to state $|-\rangle_{e}$ every $t_{re}=20$ $\mu s$. (e) Process infidelity for
different detuning and errors of the MW driving of the NV center.

Two weak RF fields are applied to individually control two different nuclear
spins with $\omega_{rfi}=\gamma_{n_{i}}B+\frac{a_{\parallel i}}{2}$. Working
in the rotating frame with
$H^{\prime}_{0}=\sum_{i}(\gamma_{ni}B+\frac{a_{\parallel i}}{2})I^{z}_{i}$, by
assuming $a_{\parallel i},a_{\perp i}\ll\Omega\ll\omega_{rfi}$ and
$a_{\parallel i}$ is comparable or larger than $a_{\perp i}$, and the total
Hamiltonian of the system can be approximated as

$\displaystyle H^{\prime}_{t}$ $\displaystyle=$
$\displaystyle\Omega\sigma_{z}+\sum_{i}\Omega_{rfi}I^{x}_{i}+a_{\parallel
i}\sigma_{x}I^{z}_{i}.$ (5)

Consider the Hamiltonian is of the form $H^{\prime}_{t}=H^{\prime}_{NV}+V$,
where $H^{\prime}_{NV}=\Omega\sigma_{z}$ denotes the Hamiltonian of the MW
driving and $V=\sum_{i}\Omega_{rfi}I^{x}_{i}+a_{\parallel
i}\sigma_{x}I^{z}_{i}$ is the weak perturbation with $a_{\parallel
i},\Omega_{rfi}\ll\Omega$. By using the Schrieffer-Wolff transformation
Kessler2012GeneralizedSF ; Bravyi2011SchriefferWolffTF  in condensed matter,
the second order expansion due to perturbation terms $V$ can be obtained as

$\displaystyle\langle\alpha|H_{e}|\beta\rangle_{e}=\frac{1}{2}\sum_{i}(\frac{\langle\alpha|V|i\rangle\langle
i|V|\beta\rangle_{e}}{E_{\alpha}-E_{i}}-\frac{\langle\alpha|V|i\rangle\langle
i|V|\beta\rangle_{e}}{E_{i}-E_{\beta}}),$ (6)

in which $\alpha,\beta,i=\\{+,-\\}$, $H_{e}$ is defined as the effective
Hamiltonian of the system and $E_{k}$ ($k=\\{+,-\\}$) is the energy
corresponding to the state $|k\rangle_{e}$. By simple calculations, we find
$\langle+|H_{e}|+\rangle_{e}=g_{e}I^{z}_{1}I^{z}_{2}$ in which
$|i\rangle_{e}=|-\rangle_{e}$ is the channel state of virtual electron spin
flip, when
$\langle-|H_{e}|-\rangle_{e}=-\langle+|H_{e}|+\rangle_{e}=g_{e}I^{z}_{1}I^{z}_{2}$
with the channel state $|i\rangle_{e}=|+\rangle_{e}$, here
$g_{e}\approx\frac{a_{\parallel 1}a_{\parallel 2}}{2\Omega}$. Thus we
adiabatically eliminate the fast electronic degrees of freedom from the slow
nuclear dynamics by a Schrieffer-Wolff (SW) transformation, and derive the
effective Hamiltonian

$\displaystyle
H_{e}\approx\sum_{i=1,2}\Omega_{rfi}I^{x}_{i}+g_{e}I^{z}_{1}I^{z}_{2}\otimes(|+\rangle_{e}\langle+|-|-\rangle_{e}\langle-|).$
(7)

To investigate the dynamics of the whole system described by density matrix
$\rho$, the system is governed by the master equation

$\displaystyle\frac{d\rho}{dt}=-i[H^{\prime}_{t},\rho]+L_{e}\rho
L_{e}^{\dagger}-\frac{1}{2}(L_{e}L_{e}^{\dagger}\rho+\rho
L_{e}L_{e}^{\dagger})$ (8)

in which $L_{e}=\sqrt{\gamma_{e}}|-\rangle_{e}\langle+|$ with
$\gamma_{e}=1/T_{1\rho}$. We apply this method and evaluate the performance of
the gate between the first spin (carbon-13) with $\gamma_{n_{1}}B=(2\pi)4$ MHz
and the second spin (silicon-29) $\gamma_{n_{2}}B=(2\pi)5.06$ MHz coupled to
an NV center with $[a_{\parallel 1},a_{\parallel 2}]=(2\pi)[9,11]$ kHz, [see
Fig. 1(a)]. Here the driving fields are given as RF Rabi frequencies
$\Omega_{rf1}=\Omega_{rf2}=(2\pi)1$ kHz, and the MW Rabi frequency
$\Omega=(2\pi)400$ kHz. The NV center is initialized to state $|-\rangle_{e}$,
when the nuclear spins are initialized to $|+_{1}-_{2}\rangle$ with
$|+_{i}\rangle=\frac{1}{2}(|\uparrow_{i}\rangle+|\downarrow_{i}\rangle)$.
Perfect state transfer between nuclear spins silicon-29 and carbon-13
$|+_{1}-_{2}\rangle\rightarrow|-_{1}+_{2}\rangle$ is shown in Fig. 1(b), in
which the population evolutions by using effective second order Hamiltonian
Eq. (7) fit well with the exact total Hamiltonian of the system Eq. (5). The
effective Hamiltonian Eq. (7) includes a three-body interaction and there are
two channels that can mediate internuclear interaction via virtual electron
spin flips through the NV microwave dressed states $|+\rangle_{e}$ and
$|-\rangle_{e}$. However, the two channels can be mixed up if the nuclear-
nuclear interaction time is longer than the NV relaxation time. As shown in
Fig. 1(c), in the case when the life time of NV center $T_{1\rho}=200$ $\mu
s$, there is no coherent evolution of the nuclear spins. Therefore, the system
is limited by the NV relaxation at room temperature.

Figure 2: Nuclear-nuclear quantum gate of two nuclear spins mediated by an NV
center in diamond. The basic logic-operations contain the initialization, gate
operation based on second-order coupling and readout of nuclear spins. The
electron spin works as a quantum bus and is stabilized in a quasi-steady state
by periodically resetting in $|-\rangle_{e}$, which is shown in the yellow
dashed box. Here P and M stand for the electron (nuclear) spin polarization,
$Y_{\theta}$ presents operation with rotation angles $\theta$ about the $Y$
axes and $N$ is the total number of resets.

##  III The effect of electron spin resets

In our scheme, in order to overcome the limitation of NV relaxation, we intend
to periodical reset the electron spin to state $|-\rangle_{e}$ within its
relaxation time $T_{1\rho}$. It is reasonable to consider the steady-state
polarization of the NV spin in a quasi-steady state as

$\displaystyle\langle 2\sigma_{z}\rangle=e^{-t_{re}/T_{1\rho}}=p.$ (9)

The three-body interaction in the effective Hamiltonian Eq. (7) can be
estimated as
$g_{e}I^{z}_{1}I^{z}_{2}\otimes(|+\rangle_{e}\langle+|-|-\rangle_{e}\langle-|)\approx
g_{e}I^{z}_{1}I^{z}_{2}\langle 2\sigma_{z}\rangle=pg_{e}I^{z}_{1}I^{z}_{2}$.
Therefore, the NV spin is not involved in the nuclear-nuclear interaction
Hamiltonian. Resets also bring another dissipation item in the master equation
and it is reasonable to have the effective relaxation of the NV spin as
$\gamma_{r}=1/T_{1\rho}+1/t_{re}$. Thus, we can rewrite the master equation of
the system as

$\displaystyle\frac{d\rho}{dt}=-i[H^{\prime}_{t},\rho]+L_{r}\rho
L_{r}^{\dagger}-\frac{1}{2}(L_{r}^{\dagger}L_{r}\rho+\rho
L_{r}^{\dagger}L_{r})$ (10)

in which $L_{r}=\sqrt{\gamma_{r}}|-\rangle_{e}\langle+|$. The non-Hermitian
Hamiltonian of the quantum jump formalism Reiter2012EffectiveOF  of the NV
spin is given by

$\displaystyle H_{H}=\Omega\sigma_{z}-\frac{1}{2}L_{r}^{\dagger}L_{r}$ (11)

According to the second expansion of SW transformation, we have energy
difference between $|+\rangle_{e}$ and $|-\rangle_{e}$ as a complex energy
$E_{+}-E_{-}=\Omega-\frac{i}{2}\gamma_{r}$, one can derive effective
Hamiltonian

$\displaystyle
H^{\prime}_{e}\approx\sum_{i}\Omega_{rfi}I^{x}_{i}+g^{\prime}_{e}I^{z}_{1}I^{z}_{2}\otimes(|+\rangle_{e}\langle+|-|-\rangle_{e}\langle-|),$
(12) $\displaystyle L_{n}\approx\sum_{k}\frac{\sqrt{\gamma_{r}}a_{\parallel
k}}{\Omega-\frac{i}{2}\gamma_{r}}I_{k}^{z}\otimes(|+\rangle_{e}\langle+|-|-\rangle_{e}\langle-|),$
(13)

in which the effective dissipation item is given by
$\langle\alpha|L_{e}|\beta\rangle=\sum_{i}\frac{\langle\alpha|L_{r}|i\rangle\langle
i|V|\beta\rangle}{E_{\alpha}-E_{i}}$, and the effective coupling is

$g^{\prime}_{e}\approx\frac{\Omega a_{\parallel 1}a_{\parallel
2}}{2(\Omega^{2}+\frac{\gamma_{r}^{2}}{4})}.$

By considering the electron spin is periodically reset in state
$|-\rangle_{e}$ with the definition $\langle 2\sigma_{z}\rangle=p$, we have
the effective master equation of the nuclear spins as

$\displaystyle\frac{d}{dt}\rho_{n}$ $\displaystyle=$
$\displaystyle-i[H_{N},\rho_{n}]+L_{N}\rho_{n}L_{N}^{\dagger}$ (14)
$\displaystyle-\frac{1}{2}(L_{N}^{\dagger}L_{N}\rho_{n}+\rho_{n}L_{N}^{\dagger}L_{N}),$

with the effective Hamiltonian and Lindblad operators

$\displaystyle H_{N}$ $\displaystyle\approx$
$\displaystyle\sum_{i}\Omega_{rfi}I^{x}_{i}+pg^{\prime}_{e}I^{z}_{1}I^{z}_{2},$
(15) $\displaystyle L_{N}$ $\displaystyle\approx$
$\displaystyle\sum_{k}\frac{\sqrt{p\gamma_{r}}a_{\parallel
k}}{\Omega-\frac{i}{2}\gamma_{r}}I_{z}^{k}.$ (16)

The effective dissipation rate is given by

$\gamma_{N}=\sum_{i,j}\frac{p\gamma_{r}a_{\parallel i}a_{\parallel
j}}{\Omega^{2}+\frac{\gamma_{r}^{2}}{4}}.$

Therefore, in principle, if resonant condition $\Omega_{rf1}=\Omega_{rf2}$ and
$pg^{\prime}_{e}\gg\gamma_{N}$ are satisfied, perfect coherent state
transition between the nuclear spins can be possible.

Figure 3: Numerical simulations are show with the same parameters in the Fig.
1, i.e., parallel couplings are given as $[a_{\parallel 1},a_{\parallel
2}]=(2\pi)[9,11]$ kHz. (a) The detection bandwidth of the swept RF frequency
based on different MW Rabi frequencies. The RF field is applied to the
carbon-13 sensor resonantly when the frequency of the other RF field is swept.
The evolution time is fixed as $T=[8.8,4.4]$ ms for $\Omega=(2\pi)[400,200]$
kHz. (b) The infidelity of quantum gate between carbon-13 and silicon-29 spins
when there is another silicon-29 spin. The infidelity varies with the detuning
$\delta_{3}$, $\delta_{3}=(a_{\parallel 3}-a_{\parallel 2})/2$.

For comparison, it is necessary to employ exact numerical simulations as well
to show the efficiency of our theory. Exact state evolutions of the system are
simulated as follows. The NV center is reinitialized periodically to the state
$|-\rangle_{e}$ of the dressed state basis, namely,
$\rho(Nt_{re})\rightarrow[\operatorname{Tr}_{e}\rho(Nt_{re})]\otimes|-\rangle_{e}\langle-|$,
where $\operatorname{Tr}_{e}$ denotes the partial trace over the electron spin
and N is an integer and reset of the NV spin every $t_{re}$ introduces an
effective interaction time in each cycle. The density matrix of the system
evolves according to

$\displaystyle\rho_{n}\rightarrow\cdot\cdot\cdot
U_{t}\operatorname{Tr}_{e}[U_{t}(\rho_{n}\otimes|-\rangle_{e}\langle-|)U^{\dagger}_{t}]\otimes|-\rangle_{e}\langle-|U^{\dagger}_{t},$
(17)

in which $U_{t}$ is the time evolution operator according master equation Eq.
(10). Both of the state evolutions by using exact numerical simulation and
effective second-order master equation based on Eq. (17) and Eq. (14) are
shown in Fig. 1(d), and the importance for periodic reinitializing of electron
spin is illustrated.

We consider the same case of a silicon-29 and a carbon-13 spin coupled to an
NV spin which is presented in Section II except that the NV center is
reinitialized to state $|-\rangle_{e}$ every $t_{re}=20$ $\mu s$. The resonant
condition is matched as $\Omega_{rf1}=\Omega_{rf2}=(2\pi)1$ kHz, both of the
simulations of exact numerical calculations and effective master equation show
perfect state transfer between nuclear spins silicon-29 and carbon-13
$|+_{1}-_{2}\rangle\rightarrow|-_{1}+_{2}\rangle$, [see Fig. 1(d)]. The state
evolutions based on effective master equation are slightly different from the
exact simulations. The possible reason is high orders of SW transformation
expansions are not included.

In the absence of periodic reinitialization of the electron spin, channel
mixing up leads to no coherent evolution of the nuclear spins [see Fig. 1(c)],
while periodic reinitialization of the NV center to state $|-\rangle_{e}$
provides a priority channel during the operations and coherent evolution
between the nuclear spins is possible to extend beyond the life time of the NV
spin. High fidelity of the nuclear-nuclear quantum gate in one-step is
possible ($>$0.99), which could be significantly higher than the fidelity
($<$0.66) of nuclear-nuclear gates achieved so far with NV centers by using 4
electron-nuclear spin quantum gates (if each fidelity $<$0.90
Taminiau2014UniversalCA ). Additionally, because of the second-order coupling,
the scheme is not sensitive to the detunings from the resonance and Rabi
frequency errors of the MW driving of the NV center [see Fig. 1(e)]. The
process fidelity of two nuclear spins is given by the overlap of the states
corresponding to the implemented evolution
$|\psi\rangle=U^{zz}\otimes\sum_{i,j}|i_{1}\rangle|j_{2}\rangle/4$ (i,j=+,-)
to the target states.

##  IV Quantum gate implementation

In the Fig. 2, we include all the related operations of the nuclear ZZ gate
implementation for experiments. The initialization consists of the electron
(nuclear) spin polarization P (M). The electron spin of NV center can be
optically initialized and read out by using laser illumination. Here, P is
obtained by the optical pumping cycle available for NV centers
Jelezko2002SingleSS ; Jelezko2004ObservationOC , whereas M is based on the
techniques developed for the nuclear single-shot measurement
Neumann2010SingleShotRO , followed by the electron state-dependent
fluorescence Jelezko2002SingleSS ; Jelezko2004ObservationOC . Gate operations
includes 3 steps: (i) Polarize the carbon-13 and silicon-29 nuclear spins by
using the NV center, which leads to the initial state
$|-\rangle_{e}\otimes|+_{1}-_{2}\rangle$ of the system. (ii) Implement the
gate operation, which is governed by the effective master equation of Eq.
(14), $U^{ZZ}(t)=\exp(-iH_{N}t)$. In this process, we periodically reset the
NV center and we control the heteronuclear spins via RF fields individually.
(iii) Use the NV center to read out the states of nuclear spins. Here we
suppose that the carbon-13 and silicon-29 spin can be well controlled by a
SWAP gate between NV and nuclear spin; i.e., it can be polarized and its
position can be detected precisely.

Assuming that the RF field is applied to the carbon-13 spin resonantly and
$\delta_{2}=\gamma_{n2}B+\frac{a_{\parallel 2}}{2}-\omega_{rf2}$ is a detuning
from the resonance of RF driving to the silicon-29 spin, the effective
Hamiltonian is given by

$\displaystyle
H^{\prime}_{N}\approx\bar{\Omega}_{rf2}I^{x^{\prime}}_{2}+\Omega_{rf1}I^{x}_{1}+pg^{\prime}_{e}I^{z}_{1}I^{z}_{2},$
(18)

in which $\bar{\Omega}_{rf2}=\sqrt{\delta_{2}^{2}+\Omega_{rf2}^{2}}$. The
initial state of the system is $|+_{1}-_{2}\rangle$, the population ($P_{+}$)
of the initial state $|+_{1}\rangle$ of carbon-13 spin is given by

$\displaystyle
P_{+}=1-\frac{(pg^{\prime}_{e}\cos\theta)^{2}\sin^{2}[t\sqrt{(pg^{\prime}_{e}\cos\theta)^{2}+\Delta^{2}}/2]}{(pg^{\prime}_{e}\cos\theta)^{2}+\Delta^{2}},$
(19)

$\cos\theta=\frac{\Omega_{rf2}}{\sqrt{\delta_{2}^{2}+\Omega_{rf2}^{2}}}$ and
$\Delta=\bar{\Omega}_{rf2}-\Omega_{rf1}$. The frequency of the RF driving to
silicon-29 spin is swept to show the detection bandwidth of our method [see
Fig. 3]. The perfect population transfer could be possible if
$|\bar{\Omega}_{rf2}-\Omega_{rf1}|\ll pg^{\prime}_{e}\cos\theta$. Namely, the
dip position with $\omega_{rf2}=\gamma_{n2}B+\frac{a_{\parallel
2}}{2}=(2\pi)5065.5$kHz indicates $\delta_{2}=0$ and
$\Omega_{rf2}=\Omega_{rf1}$. The dip height depends on the longitudinal
dipolar coupling with $P_{+}=\cos^{2}(pg^{\prime}_{e}T/2)$ with $T$ the total
evolution time. Increasing the Rabi frequency of the MW driving induces a
smaller effective coupling between the nuclear spins, decreasing the frequency
bandwidth [see Fig. 3], which also gives a tunable frequency filter. We
investigate the case when there is a third silicon-29 nuclear spin with
coupling $a_{\parallel_{3}}$ to NV center and show the efficiency of the
selectivity of our scheme. In Fig. 3(b), we initialize the nuclei in state
$|+_{1}-_{2}-_{3}\rangle$. We calculate the infidelity of nuclear ZZ gate
between two nuclear spins affected by the third silicon-29 spin. The first
carbon-13 and second silicon spins are coupled to an NV center with
$[a_{\parallel 1},a_{\parallel 2}]=(2\pi)[11,9]$ kHz. The coupling of the
third spin matches $\delta_{3}=(a_{\parallel 2}-a_{\parallel 3})/2>(2\pi)0.5$
kHz gives high fidelity of the gate operation of target nuclear spins
($>0.95$). Thus, our method can address two heteronuclear spins by using the
RF fields individually, and implement a near perfect quantum gate at room
temperature with high selectivity.

##  V Sensing application

Another important application of our scheme is to detect nuclear spins in
another species outside of the diamond. For example, one can use a carbon-13
spin which is well controlled by the NV center to detect hydrogen-1 spins
outside of the diamond. Thus carbon-13 spin is initialized in state
$|+_{1}\rangle$ and hydrogen-1 spins are in maximally mixed states, one can
detect the signal by measuring the probability of that the sensor carbon-13
spin remains in the state $|+_{1}\rangle$. The related operations are shown in
Fig. 2, except of the initialization and readout of hydrogen-1 spins. To
evaluate the effectiveness of our sensing application, we initialize the
carbon-13 spin to $|+_{1}\rangle$ state to detect the other three nuclear
(hydrogen-1) spins, see Fig. 4. Based on Eq. (19), we can adjust the parallel
components and control frequency of RF to match the resonance condition for
target (hydrogen-1) nuclear spins. Once the resonance condition is achieved,
the signal dips mark the presence of the nuclear spins. Similar to the case in
the previous scheme Chen2017DissipativelySQ , the frequency resolution is not
limited by the NV relaxation $T_{1\rho}$, it requires
$|\bar{\Omega}_{rfi}-\Omega_{rf1}|>pg^{\prime}_{e}\cos\theta$ which is limited
by the target spin decoherence time $T_{2}=200$ ms, the sensitivity per unit
time of our scheme is proportional to
$pg^{\prime}_{e}/\sqrt{1/\Gamma_{N}}\sim(a_{\parallel i}/4)/\sqrt{T_{1\rho}}.$

Figure 4: Assuming a carbon-13 and three hydrogen-1 spins coupled to the NV
center with parallel couplings $[a_{\parallel 1},a_{\parallel 2},a_{\parallel
3},a_{\parallel 4}]=(2\pi)[11,4,9,11]$ kHz. The carbon-13 spin is the sensor
with the initial state $|+_{1}\rangle$ and the other three hydrogen-1 are the
targets which are in maximum mixed states. The population of the initial state
$|+_{1}\rangle$ of the carbon-13 sensor is calculated by exact simulation of
Eq. (17) with evolution time $T=8$ ms for $\Omega=(2\pi)400$ kHz and the
target spins decoherence time $T_{2}=200$ ms.

##  VI Conclusion

Before summarizing, we would like to compare our new method with the previous
scheme Chen2017DissipativelySQ . Both of schemes realize the effective
coupling of target nuclear spins by periodical resetting of electron spin, and
the impact of NV decoherence and relaxation processes on these nuclei are
suppressed. The difference mainly comes from that we employ two weak RF fields
to control the heteronuclear spins individually. In the previous approach,
nuclear spins in the same species are considered, coarse tuning of the
direction of the magnetic field in advance is necessary for matching the
resonant condition. Our method allows for effective interaction of nuclear
spins in different species which loosen the requirement of the stringent
resonant condition of the previous scheme. The signal accumulation and
resonant condition depend on longitudinal dipolar coupling component
$a_{\parallel_{i}}$, which leads to no information of the transverse coupling
component $a_{\perp_{i}}$ to be detected in this new scheme. Therefore, less
information is detected when one uses this new method for sensing nuclear
spins near the NV center. But it also bring advantage that the quantum gate
implementation is not limited by the transverse coupling component.

In this work we presented an extension of dissipatively stabilized NV center
to a case, in which one achieves a tunable second-order effective coupling
between distant nuclear spins in different species. We employ RF fields and MW
driving to match the resonant condition, and control electron and nuclear
spins individually. The coupling is mediated by a dissipatively decoupled
electron spin of a NV center, when the impact of NV decoherence and relaxation
processes on these nuclei are suppressed. Thus the effective indirect
interaction enables selectively initialization and coherent control of nuclear
spins as well as analysis of complex spin structures at ambient condition.

_Acknowledgements —_ The authors thank Natural Science Foundation of Hunan
Province, China (2019JJ10002), Hunan Provincial Hundred People Plan (2019),
Huxiang High-level Talent Gathering Project (2019RS1043).

## References

  * (1) T. Taminiau, J. Cramer, T. van der Sar, V. Dobrovitski, and R. Hanson, Nature nanotechnology 9 3, 171 (2014).
  * (2) J. Zhang, S. S. Hegde, and D. Suter, Physical Review A 98, 042302 (2018).
  * (3) G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, et al., Nature 506, 204 (2014).
  * (4) L. Robledo, L. Childress, H. Bernien, B. Hensen, P. Alke-made, and R. Hanson, Nature 477, 574 (2011).
  * (5) T. Taminiau, J. Wagenaar, T. van der Sar, F. Jelezko, V. Dobrovitski, and R. Hanson, Physical review letters 109 13, 137602 (2012).
  * (6) S. Kolkowitz, Q. Unterreithmeier, S. Bennett, and M. Lukin, Physical review letters 109 13, 137601 (2012).
  * (7) N. Zhao, J. Honert, B. Schmid, M. Klas, J. Isoya, M. Markham, D. Twitchen, F. Jelezko, R.-B. Liu, H. Fedder, et al., Nature nanotechnology 7 10, 657 (2012).
  * (8) A. Ermakova, G. Pramanik, J.-M. Cai, G. Algara-Siller, U. Kaiser, T. Weil, Y.-K. Tzeng, H. C. Chang, L. McGuinness, M. Plenio, et al., Nano letters 13 7, 3305 (2013).
  * (9) H. J. Mamin, M. Kim, M. Sherwood, C. Rettner, K. Ohno, D. Awschalom, and D. Rugar, Science 339, 557 (2013).
  * (10) S. Zaiser, T. Rendler, I. Jakobi, T. Wolf, S.-Y. Lee, S.Wagner, V. Bergholm, T. Schulte-Herbrüggen, P. Neumann, and J. Wrachtrup, Nature Communications 7 (2016).
  * (11) T. Staudacher, F. Shi, S. Pezzagna, J. Meijer, J. Du, C. A. Meriles, F. Reinhard, and J. Wrachtrup, Science 339, 561 (2013).
  * (12) Y. Wu, F. Jelezko, M. Plenio, and T. Weil, Angewandte Chemie 55 23, 6586 (2016).
  * (13) M. Doherty, N. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. Hollenberg, Physics Reports 528, 1 (2013).
  * (14) D. A. Golter, T. Oo, M. Amezcua, K. A. Stewart, and H. Wang, Physical review letters 1116, 143602 (2016).
  * (15) D. A. Golter and H. Wang, Physical review letters 112, 116403 (2014).
  * (16) J. Casanova, Z.-Y. Wang, and M. Plenio, Physical review letters 117 13, 130502 (2016).
  * (17) A. Bermúdez, F. Jelezko, M. Plenio, and A. Retzker, Physical review letters 107 15, 150503 (2011).
  * (18) J. Zimmermann, P. London, Y. Yirmiyahu, F. Jelezko, A. Blank, and D. Gershoni, Physical Review B 102, 245408 (2020).
  * (19) X. Rong, J. Geng, Z. Wang, Q. Zhang, C. Ju, F. Shi, C.K. Duan, and J. Du, Physical review letters 112, 050503 (2014).
  * (20) S. S. Hegde, J. Zhang, and D. Suter, Physical review letters 124, 220501 (2020).
  * (21) M. Abobeih, J. Randall, C. Bradley, H. Bartling, M. Bakker, M. Degen, M. Markham, D. Twitchen, and T. Taminiau, Nature 576, 411 (2019).
  * (22) C. Bradley, J. Randall, M. Abobeih, R. Berrevoets, M. Degen, M. Bakker, M. Markham, D. Twitchen, and T. Taminiau, Physical Review X 9, 031045 (2019).
  * (23) Z.-Y. Wang, J. Casanova, and M. Plenio, Nature Communications 8 (2017).
  * (24) B. Tratzmiller, J. Haase, Z.-Y. Wang, and M. Plenio, Physical Review A 103, 012607 (2021).
  * (25) J. Casanova, Z.-Y. Wang, and M. Plenio, Physical Review A 96, 032314 (2017).
  * (26) M. Degen, S. Loenen, H. Bartling, C. Bradley, A. Meinsma, M. Markham, D. Twitchen, and T. Taminiau, Nature Communications 12, 1 (2021).
  * (27) J. Bian, M. Jiang, J. Cui, X. Liu, B. Chen, Y. Ji, B. Zhang, J. Blanchard, X. Peng, and J. Du, Physical Review A 95, 052342 (2017).
  * (28) Q. Chen, I. Schwarz, and M. Plenio, Physical review letters 119 1, 010801 (2017).
  * (29) E. Kessler, Physical Review A 86, 012126 (2012).
  * (30) S. Bravyi, D. DiVincenzo, and D. Loss, Annals of Physics 326, 2793 (2011).
  * (31) F. Reiter and A. Sørensen, Physical Review A 85, 032111 (2012).
  * (32) F. Jelezko, I. Popa, A. Gruber, C. Tietz, J. Wrachtrup, A. Nizovtsev, and S. Kilin, Applied Physics Letters 81, 2160 (2002).
  * (33) F. Jelezko, T. Gaebel, I. Popa, A. Gruber, and J. Wrachtrup, Physical review letters 92 7, 076401 (2004).
  * (34) P. Neumann, J. Beck, M. Steiner, F. Rempp, H. Fedder, P. Hemmer, J. Wrachtrup, and F. Jelezko, Science 329, 542 (2010).
  * (35) A. Laraoui, F. Dolde, C. Burk, F. Reinhard, J. Wrachtrup, and C. A. Meriles, Nature Communications 4, 1651 (2013).
  * (36) M. Pfender, N. Aslam, H. Sumiya, S. Onoda, P. Neumann, J. Isoya, C. A. Meriles, and J. Wrachtrup, Nature Communications 8 (2017).
  * (37) J. Cai, F. Jelezko, M. Plenio, and A. Retzker, New Journal of Physics 15, 013020 (2013).


11institutetext:  Moscow Institute of Physics and Technology (National
Research University), Dolgoprudny, Russia
11email<EMAIL_ADDRESS>
# On the Computational Efficiency of Catalyst Accelerated Coordinate
Descent††thanks: D.A. Pasechnyuk’s research was supported by the A.M.
Raigorodsky Scholarship in the field of optimization and RFBR grant
19-31-51001 (Scientific mentoring). The work of V.V. Matyukhin was supported
by the Ministry of Science and Higher Education of the Russian Federation
(state assignment) No. 075-00337-20-03, project number 0714-2020-0005.

Dmitry Pasechnyuk 11 0000-0002-1208-1659 Vladislav Matyukhin 11
0000-0003-4961-1577

###### Abstract

This article is devoted to one particular case of using universal accelerated
proximal envelopes to obtain computationally efficient accelerated versions of
methods used to solve various optimization problem setups. We propose a
proximally accelerated coordinate descent method that achieves the efficient
algorithmic complexity of iteration and allows taking advantage of the data
sparseness. It was considered an example of applying the proposed approach to
optimizing a SoftMax-like function, for which the described method allowing
weaken the dependence of the computational complexity on the dimension $n$ in
$\mathcal{O}(\sqrt{n})$ times and, in practice, demonstrates a faster
convergence in comparison with standard methods. As an example of applying the
proposed approach, it was shown a variant of obtaining on its basis some
efficient methods for optimizing Markov Decision Processes (MDP) in a minimax
formulation with a Nesterov smoothed target functional.

###### Keywords:

Proximal accelerated method Catalyst Accelerated coordinate descent method
SoftMax Markov Decision Processes.

##  1 Introduction

One of the most important theoretical results in convex optimization was the
development of accelerated optimization methods [26]. At the initial stage of
implementation of this concept, many accelerated algorithms for different
problem setups were proposed. But each such case required special
consideration of the possibility of acceleration. Therefore, the proposed
designs were significantly different and did not allowing assume a way to
generalize them. A significant step towards the development of a universal
scheme for accelerating optimization methods was the work in which an
algorithm called Catalyst proposed, based on the idea of [29, 30] and allowing
to accelerate other optimization methods, using them for the sequential
solution of several Moreau-Yosida regularized auxiliary problems [18, 19].
Following these ideas, many variants of the applications of this method and
its modifications [14, 17, 28] were proposed. Among the most recent results
until the time of writing this paper, the generalizations of the discussed
approach to tensor methods [5, 23, 10, 11] were also described. The
corresponding representation of the accelerated proximal envelope to the
authors’ knowledge is the most general of those described in the literature,
and therefore, this work is focused primarily on the methods proposed in the
works [10, 11].

The main motivation of this work is to describe the possibilities of the
practical application of universal accelerated proximal envelopes for
constructing computationally and oracle efficient optimization methods. Let us
consider the classical coordinate descent method [4], the iteration of which
for the convex function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ has the form:

$x_{k+1}^{i}=x_{k}^{i}-\eta\nabla_{i}f(x_{k}),\quad
i\sim\mathcal{U}\\{1,...,n\\},\;\;\eta>0.$

One of the many applications of this method is the optimization of
functionals, where the calculation of the one component of the gradient is
significantly more efficient than the calculation of the full gradient vector
of these functionals (in particular, many problems in the case of sparse
formulations satisfy this condition). However, the oracle complexity of this
method, provided that the method stops when the $\varepsilon$-small residual
by the function value is reached, is
$\mathcal{O}\left(n\frac{\overline{L}R^{2}}{\varepsilon}\right)$, where
$R^{2}=\|x_{0}-x_{*}\|_{2}^{2}$, $\overline{L}=\frac{1}{n}\sum_{i=1}^{n}L_{i}$
is the average of the Lipschitz constants of the gradient components;
moreover, this estimate is not optimal for the class of convex problems. Let
us now consider the accelerated coordinate descent method proposed by Yu.E.
Nesterov [27], the oracle complexity of this method corresponds to the optimal
estimate:
$\mathcal{O}\left(n\sqrt{\frac{\widetilde{L}R^{2}}{\varepsilon}}\right)$,
where $\sqrt{\widetilde{L}}=\frac{1}{n}\sum_{i=1}^{n}\sqrt{L_{i}}$ is the mean
of square roots of the Lipschitz constants of the gradient components. At the
same time, the situation changes drastically when the algorithmic complexity
of the method considered: even if the computation of one component of the
gradient has the complexity $\mathcal{O}(s)$, $s\ll n$, the complexity of
iterating the accelerated coordinate descent method will be $\mathcal{O}(n)$,
unlike the standard method, the iteration complexity of which is
$\mathcal{O}(s)$, this means, essentially, that the degree of the sparseness
of the problem when using the accelerated coordinate descent method does not
significantly affect the complexity of the algorithm, and besides, the
complexity in this case quadratically depends on the dimension of the problem:
together, this to some extent devalues the use of the coordinate descent
method in this case. Thus, an interesting problem is the construction of an
accelerated coordinate descent method, the iteration complexity of which, as
in the standard version of the method, is $\mathcal{O}(s)$, this is possible
due to the application of the universal accelerated proximal envelope
‘‘Accelerated Meta-algorithm’’ [11]. Note that the approach to accelerating
the cordinate descent method considered in this work is, in fact, one of the
cases of applying the technique described in [14]. Thus, the content of the
paper is devoted to a more detailed analysis of this particular approach,
including implementation features and possible applications.

This article consists of an introduction, conclusion and the main Section 2.
It describes the theoretical results on the convergence and algorithmic
complexity of the coordinate descent method, accelerated by using the
‘‘Accelerated Meta-algorithm’’ envelope (Section 2.1). Using the example of
the SoftMax-like function optimization problem, it was experimentally tested
method’s effectiveness with relation to its working time, there were described
the possibilities of its computationally efficient implementation, also it was
carried out a comparison with standard methods (Section 2.2). Further, as an
example of applying the proposed approach, it was provided a method for
optimizing Markov Decision Processes in a minimax formulation, based on
applying the method introduced in this paper to the Nesterov smoothed target
functional. The proposed approach obtains estimates close to that for several
efficient and practical methods for optimizing the discounted MDP and matches
the best estimates for the averaged MDP problem (Section 2.4).

##  2 Accelerated Meta-algorithm and coordinate descent

###  2.1 Theoretical guarantees

Let us consider the following optimization problem of the function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$:

$\min_{x\in\mathbb{R}^{n}}f(x),$

subject to:

  1. 1.

$f$ is differentiable on $\mathbb{R}^{n}$;

  2. 2.

$f$ is convex on $\mathbb{R}^{n}$;

  3. 3.

$\nabla_{i}f$ is component-wise Lipschitz continuous, i.e. $\forall
x\in\mathbb{R}^{n}$ and $u\in\mathbb{R},$ $\exists\leavevmode\nobreak\
L_{i}\leavevmode\nobreak\ \in\leavevmode\nobreak\ \mathbb{R}\;(i=1,\ldots,n)$,
such that

$\left|\nabla_{i}f\left(x+ue_{i}\right)-\nabla_{i}f(x)\right|\leq L_{i}|u|,$

where $e_{i}$ is the $i$-th unit basis vector, $i\in\\{1,...,n\\}$;

  4. 4.

$\nabla f$ is $L$-Lipschitz contiouous.

Let us turn to the content of the work [11], where a general version of the
‘‘Accelerated Meta-algorithm’’ for solving convex optimization problems for
composite functionals of the form $F(x)=f(x)+g(x)$ was proposed. For the
considered formulation of the problem, such generality not required; it is
sufficient to apply a special case of the described scheme for $p=1$, $f\equiv
0$ (using the designations of the corresponding work), in which the described
envelope takes the form of an accelerated proximal method. This method is
listed as Algorithm 1.

Input: $H>0$, $x_{0}\in\mathbb{R}^{n}$;

$\lambda\leftarrow\nicefrac{{1}}{{2H}}$;

$A_{0}\leftarrow 0$; $v_{0}\leftarrow x_{0}$;

for  _k = 0, …, $\widetilde{N}-1$_ do

$\displaystyle a_{k+1}\leftarrow\frac{\lambda+\sqrt{\lambda^{2}+4\lambda
A_{k}}}{2}$;

$A_{k+1}\leftarrow A_{k}+a_{k+1}$;

$\displaystyle\widetilde{x}_{k}\leftarrow\frac{A_{k}v_{k}+a_{k+1}x_{k}}{A_{k+1}}$;

By running the method $\mathcal{M}$,

find the solution of the following auxiliary problem

with an accuracy $\varepsilon$ by the argument:

$\displaystyle
v_{k+1}\in\text{Arg}^{\varepsilon}\min_{y\in\mathbb{R}^{n}}\left\\{f(y)+\frac{H}{2}\|y-\widetilde{x}_{k}\|^{2}_{2}\right\\}$;

$\displaystyle x_{k+1}\leftarrow x_{k}-a_{k+1}\nabla f(v_{k+1})$;

end for

return $v_{\widetilde{N}}$;

Algorithm 1 Accelerated Meta-algorithm for First-order Method $\mathcal{M}$

Before formulating any results on the convergence of the proposed accelerated
coordinate descent method described below, it is necessary to start with a
detailed consideration of the process of solving the auxiliary problems, where
its analytical solution is available only in rare cases. Therefore, one should
apply numerical methods to find its approximate solution, and that is
inaccurate. The solving process of the auxiliary problem will be until the
following stop condition is satisfied ([16], Appendix B):

$\left\|\nabla\left\\{F(y_{\star}):=f(y_{\star})+\frac{H}{2}\|y_{\star}-\widetilde{x}_{k}\|^{2}_{2}\right\\}\right\|_{2}\leq\frac{H}{2}\|y_{\star}-\widetilde{x}_{k}\|_{2},$
(1)

where $y_{\star}$ is an approximate solution of the auxiliary problem,
returned by internal method $\mathcal{M}$. Due to the $\|\nabla
F(y_{*})\|_{2}=0$ (where $y_{*}$ denotes an exact solution of the considered
problem), and due to the $(L+H)$-Lipschitz continuity of $\nabla F$, we have
got:

$\|\nabla F(y_{\star})\|_{2}\leq(L+H)\|y_{\star}-y_{*}\|_{2}.$ (2)

Writing out the triangle inequality:
$\|\widetilde{x}_{k}-y_{*}\|_{2}-\|y_{\star}-y_{*}\|_{2}\leq\|y_{\star}-\widetilde{x}_{k}\|_{2}$,
and using together the inequalities (1), (2), we have got the final view of
stop condition:

$\|y_{\star}-y_{*}\|_{2}\leq\frac{H}{3H+2L}\|\widetilde{x}_{k}-y_{*}\|_{2}.$
(3)

Input: $y_{0}\in\mathbb{R}^{n}$;

$Z\leftarrow\sum_{i=1}^{n}(H+L_{i})$;

$p_{i}\leftarrow(H+L_{i})/Z,\quad i\in\\{1,...,n\\}$;

Discrete probability distribution $\pi$

with probabilities $p_{i}$;

for  _k = 0, …, $N-1$_ do

$i\sim\pi\\{1,...,n\\}$;

$y_{k+1}\leftarrow y_{k}$;

$\displaystyle y_{k+1}^{i}=y_{k}^{i}-\frac{1}{H+L_{i}}\nabla_{i}F(y_{k})$;

end for

return $y_{N}$;

Algorithm 2 Coordinate descent method

Essentially this implies that the required argument accuracy of solving the
auxiliary problem does not depend on the accuracy of the main problem. That
makes it possible to simplify the obtaining of further results.

Let us now consider the main method used for solving auxiliary problems: the
content of the coordinate descent method [25] (in the particular case, when
$\gamma=1$) is listed as Algorihtm 2. For this method, in the case of the
considered auxiliary problems, the following result holds:

###### Theorem 2.1

([4], theorem 6.8) Let $F$ be $H$-strongly convex function with respect to
$\|\cdot\|_{2}$. Then for the sequence $\\{y_{k}\\}_{k=1}^{N}$ generated by
the described coordinate descent algorithm 2, it holds the following
inequality:

$\displaystyle\mathbb{E}[F(y_{N})]-F(y_{*})\leq\left(1-\frac{1}{\kappa}\right)^{N}(F(y_{0})-F(y_{*})),$
(4)
$\displaystyle\text{where}\quad\kappa=\frac{H}{Z},\;\;Z=\sum_{i=1}^{n}(H+L_{i}),$
(5)

where $\mathbb{E}[\cdot]$ denotes the mathematical expectation of the
specified random variable with respect to the randomness of methods trajectory
induced by a random choice of components $i$ at each iteration.

Using this result, formulate the following statement on the number of
iterations of the coordinate descent sufficient to satisfy the stop condition
(3).

###### Corollary 1

The expectation $\mathbb{E}[y_{N}]$ of the point resulting from the coordinate
descent method (Algorithm 2) satisfies the condition (3) if the following
inequality on iterations number holds:

$\displaystyle N\geq
N(\widetilde{\varepsilon})=\Bigg{\lceil}\frac{Z}{H}\ln{\left\\{\left(1+\frac{L}{H}\right)\left(3+\frac{2L}{H}\right)^{2}\right\\}}\Bigg{\rceil},$
(6)
$\displaystyle\text{where}\quad\widetilde{\varepsilon}=\frac{H}{2}\left(\frac{H}{3H+2L}\right)^{2}\|y_{0}-y_{*}\|_{2}^{2}.$
(7)

###### Proof

is in Appendix 0.A.

Now that the question of the required accuracy and oracle complexity of
solving the auxiliary problem using the proposed coordinate descent method is
clarified, we can proceed to the results on the convergence of the Accelerated
Meta-algorithm. For the used stop condition (3) of the method that solves the
auxiliary problem, the following result on the convergence of the Accelerated
Meta-algorithm holds.

###### Theorem 2.2

([11], theorem 1) For $H>0$ and the sequence
$\\{v_{k}\\}_{k=1}^{\widetilde{N}}$ generated by the Accelerated Meta-
algorithm with some non-stochastic internal method, it holds the following
inequality:

$f(v_{\widetilde{N}})-f(x_{*})\leq\frac{48}{5}\frac{H\|x_{0}-x_{*}\|_{2}^{2}}{\widetilde{N}^{2}}.$
(8)

Based on the last statement, one can formulate a theorem on the convergence of
the Accelerated Meta-algorithm in the case of using the stochastic method and,
in particular, coordinate gradient descent method.

###### Theorem 2.3

For $H>0$ and some $0<\delta<1$, the point $v_{\widetilde{N}}$ resulting from
the Accelerated Meta-algorithm uses coordinate descent method to solve the
auxiliary problem, solving it $N_{\delta}$ iterations, satisfies the condition

$Pr(f(v_{\widetilde{N}})-f(x_{*})<\varepsilon)\geq 1-\delta,$

where $Pr(\cdot)$ denotes the probability of the specified event, if

$\displaystyle\widetilde{N}\geq\Bigg{\lceil}\frac{4\sqrt{15}}{5}\sqrt{\frac{H\|x_{0}-x_{*}\|_{2}^{2}}{\varepsilon}}\;\Bigg{\rceil},$
(9) $\displaystyle N_{\delta}\geq
N\left(\frac{\widetilde{\varepsilon}\delta}{\widetilde{N}}\right)=\Bigg{\lceil}\frac{Z}{H}\ln{\left\\{\frac{\widetilde{N}}{\delta}\left(1+\frac{L}{H}\right)\left(3+\frac{2L}{H}\right)^{2}\right\\}}\Bigg{\rceil}.$
(10)

###### Proof

The corollary 1 presents an estimate of the number of iterations sufficient to
satisfy the following condition for the expected value of the function at the
resulting point of the method:

$\mathbb{E}[F(y_{N(\widetilde{\varepsilon})})]-F(y_{*})\leq\widetilde{\varepsilon}.$

Use the Markov inequality and obtain the formulation of this condition in
terms of the estimate of the probability of large deviations [1]: deliberately
choose the admissible value of the probability of non-fulfillment of the
stated condition, so that $0<\delta/\widetilde{N}<1$, where $\widetilde{N}$
expressed from (8); then

$Pr\left(F\left(y_{N(\widetilde{\varepsilon}\delta/\widetilde{N})}\right)-F(y_{*})\geq\widetilde{\varepsilon}\right)\leq\frac{\delta}{\widetilde{N}}\cdot\frac{\mathbb{E}\left[F\left(y_{N(\widetilde{\varepsilon}\delta/\widetilde{N})}\right)\right]-F(y_{*})}{\widetilde{\varepsilon}\cdot\delta/\widetilde{N}}=\frac{\delta}{\widetilde{N}}.$

Since the probability that the obtained solution of some separately taken
auxiliary problem will not satisfy the stated condition is equal to
$\delta/\widetilde{N}$, it means that the probability that for $\widetilde{N}$
iterations of the Accelerated Meta-algorithm the condition will not be
satisfied although if for one of the problems, there is
$\widetilde{N}\cdot\delta/\widetilde{N}=\delta$, whence the proved statement
follows.

Further, combining the estimates given in Theorem 2.3, we can obtain an
asymptotic estimate for the total number of iterations of the coordinate
descent method, sufficient to obtain a solution of the considered optimization
problem with a certain specified accuracy, as well as an estimate for the
optimal $H$ value:

###### Corollary 2

In order to the point $v_{\widetilde{N}}$, which is the result of the
Accelerated Meta-algorithm, to satisfy the condition

$Pr(f(v_{\widetilde{N}})-f(x_{*})<\varepsilon)\geq 1-\delta,$

it is sufficient to perform a total of

$\hat{N}\geq\widetilde{N}\cdot
N_{\delta}=\mathcal{O}\left(\frac{Z\|x_{0}-x_{*}\|_{2}}{\sqrt{H}}\cdot\frac{1}{\varepsilon^{1/2}}\log\left\\{\frac{1}{\varepsilon^{1/2}\delta}\right\\}\right)$
(11)

iterations of coordinate descent method to solve the auxiliary problem. In
this case, the optimal value of the regularization parameter $H$ of the
auxiliary problem should be chosen as $H\simeq\frac{1}{n}\sum_{i=1}^{n}L_{i}$
($\simeq$ denotes equality up to a small factor of the $\log$ order).

###### Proof

The expression for $\hat{N}$ can be obtained by the direct substitution of one
of the estimates given in (9) into another, and their subsequent
multiplication. If we exclude from consideration a small factor of order
$\log(L/H)$, the constant in the estimate will depend on $H$ as:

$\sqrt{H}\cdot\frac{Z/n}{H}=\sqrt{H}\left(1+\frac{\frac{1}{n}\sum_{i=1}^{n}L_{i}}{H}\right).$

By minimizing the presented expression by $H$, we get the specified result.

A similar statement can be formulated for the expectation of the number of
total iterations of the coordinate descent method, without resorting to
estimates of the probabilities of large deviations, following the reasoning
scheme proposed in [9]:

###### Theorem 2.4

The expectation $\mathbb{E}[\hat{N}]$ of the sufficient total number of the
iterations of the coordinate descent method, to obtain a point
$v_{\widetilde{N}}$ satisfying the following condition:

$f(v_{\widetilde{N}})-f(x_{*})<\varepsilon$

can be bounded as follows:

$\mathbb{E}[\hat{N}]\leq\widetilde{N}\cdot(N(\widetilde{\varepsilon})+1)=\mathcal{O}\left(\sqrt{\frac{\overline{L}\|x_{0}-x_{*}\|^{2}_{2}}{\varepsilon}}\right),\quad\text{
where}\quad\overline{L}=Z/n=\frac{1}{n}\sum_{i=1}^{n}L_{i}.$

###### Proof

is in Appendix 0.B.

As we can see, the proposed scheme of reasoning allows us to reduce the
logarithmic factor in estimating the number of iterations of the method.
However, this result is less constructive than the one presented in Corollary
2. Indeed, in the above reasoning, we operated with the number of iterations
$N$, after which the stopping condition for the internal method is satisfied,
but in the program implementation, stopping immediately after fulfilment of
this condition is not possible, if only because, which is impossible the
verification of this criterion due to the natural lack of information about
$y_{*}$. So the last result is more relevant from the point of view of
evaluating the theoretical effectiveness of the method, while when considering
specific practical cases, one should apply the estimate (11).

Let us now consider in more detail the issue of the algorithmic complexity of
the proposed accelerated coordinate gradient descent method.

###### Theorem 2.5

Let the complexity of computing one component of the gradient of $f$ is
$\mathcal{O}(s)$. Then the algorithmic complexity of the Accelerated Meta-
algorithm with coordinate descent as internal method is

$T=\mathcal{O}\left(s\cdot
n\cdot\sqrt{\frac{\overline{L}\|x_{0}-x_{*}\|_{2}^{2}}{\varepsilon}}\log\left\\{\frac{1}{\varepsilon^{1/2}\delta}\right\\}\right).$

###### Proof

is in Appendix 0.C.

Note also that the memory complexity of the method is $\mathcal{O}(n)$, as
well as the complexity of the preliminary calculations (for the coordinate
descent method, there is no need to perform them again for every iteration).

Let us compare the estimates obtained for proposed approach (Catalyst CDM)
with estimates of other methods that can be used to solve problems in the
described setting: Fast Gradient Method (FGM), classical Coordinate Descent
Method (CDM) and Accelerated Coordinate Descent Method in the version of Yu.E.
Nesterov (ACDM). The estimates are shown in Table 1, below. As can be seen
from the above asymptotic estimates of the computational complexity, the
proposed method allows to achieve a convergence rate that is not inferior to
other methods with respect to the nature of the dependence on the dimension of
the problem $n$ and the required accuracy $\varepsilon$, at a certain price
for this in the form of a logarithmic factor assessment. Note, in addition,
that despite the significant similarity of estimates, between the two most
efficient methods in the table (FGM and Catalyst CDM) there is also a
difference in the constants characterizing the smoothness of the function,
i.e. $L$ in FGM and $\overline{L}$ in Catalyst CDM, thus the behavior of the
considered method for various problems directly depends on the character of
its componentwise smoothness.

Table 1: Comparison of the effectiveness of methods Algorithm | Iteration complexity | Comp. complexity | Source
---|---|---|---
FGM | $\mathcal{O}\left(s\cdot n\right)$ | $\displaystyle\mathcal{O}\left(s\cdot n\cdot\frac{1}{\varepsilon^{1/2}}\right)$ | [26]
CDM | $\mathcal{O}\left(s\right)$ | $\displaystyle\mathcal{O}\left(s\cdot n\cdot\frac{1}{\varepsilon}\right)$ | [4]
ACDM | $\mathcal{O}\left(n\right)$ | $\displaystyle\mathcal{O}\left(n^{2}\cdot\frac{1}{\varepsilon^{1/2}}\right)$ | [27]
Catalyst CDM | $\mathcal{O}\left(s\right)$ | $\displaystyle\widetilde{\mathcal{O}}\left(s\cdot n\cdot\frac{1}{\varepsilon^{1/2}}\right)$ | this paper

###  2.2 Numerical experiments

This section describes the character of the practical behavior of the proposed
method by the example of the following optimization problem for the SoftMax-
like function:

$\min\limits_{x\in\mathbb{R}^{n}}\leavevmode\nobreak\ \leavevmode\nobreak\
\\{f(x)=\gamma\ln\left(\sum_{j=1}^{m}\exp\left(\frac{\left[Ax\right]_{j}}{\gamma}\right)\right)-\langle
b,x\rangle\\},$ (12)

where $b\in\mathbb{R}^{n}$, $A\in\mathbb{R}^{m\times
n},\gamma\in\mathbb{R}_{+}$. Problems of this kind are essential for many
applications, in particular, they arise in entropy-linear programming problems
as a dual problem [8, 12], in particular in optimal transport problem, is also
a smoothed approximation of the $\max$ function (which gave the function the
name SoftMax) and, accordingly, the norm $\|\cdot\|_{\infty}$, which may be
needed in some formulations of the PageRank problem or for solving systems of
linear equations. Moreover, in all the described problems, an important
special case is the sparse setting, in which the matrix $A$ is sparse, that
is, the average number of nonzero entries in the row $A_{j}$ does not exceed
some $s\ll n$ (it will also be convenient to assume the possibility one of the
strings $A_{j}$ is completely non-sparse).

Let us formulate the possessed properties by the function $f$ [13]:

  1. 1.

$f$ is differentiable;

  2. 2.

$\nabla f$ satisfies the Lipschitz condition with the constant
$L=\max_{j=1,...,m}\|A_{j}\|_{2}^{2}$;

  3. 3.

$\nabla_{i}f$ satisfy the component-wise Lipschitz condition with the
constants $L_{i}=\max_{j=1,...,m}|A_{ji}|$.

Let us write the expression for the $i$-th component of the gradient of the
function $f$:

$\nabla_{i}f(x)=\frac{\sum_{j=1}^{m}{A_{ji}\exp\left(\left[Ax\right]_{j}\right)}}{\sum_{j=1}^{m}{\exp\left(\left[Ax\right]_{j}\right)}}.$

As we can see, the naive calculation of this expression can take time
comparable to the calculation of the whole gradient and it will significantly
affect the computational complexity, and hence the working time of the method.
However, at the same time, many terms in this expression can be recalculated
either infrequently or in a component-wise manner, and used as members of
methods additional sequences when performing a step of method, so that the
complexity of the iteration will remain efficient, and the use of the
coordinate descent methods will be justified. For the convenience of
describing the computational used methods, we write the step of the coordinate
descent algorithm in such form:

$y_{k+1}=y_{k}+\eta e_{i},$

where $\eta$ is the step size, multiplied by the corresponding gradient
component, $e_{i}$ is the $i$-th unit basis vector.

In order to the performed experiments:

  1. 1.

We will store a sequence of values
$\left\\{\exp\left(\left[Ay_{k}\right]_{j}\right)\right\\}_{j=1}^{m}$, used to
calculate the sum in the numerator. Updating these values after executing a
method step takes algorithmic complexity $\mathcal{O}(s)$, due to the
$Ay_{k+1}=Ay_{k}+\eta A_{i}$, and at the same time $A_{i}$ has at most $s$
nonzero components, which means that it will be necessary to calculate only
less than or equal to $s$ correcting factors and multiply the corresponding
values from the sequence by them.

  2. 2.

From the first point, it can be understood that the multiplication of sparse
vectors should be performed in $\mathcal{O}(s)$, considering only nonzero
components. In terms of program implementation, this means the need to use a
sparse representation for cached values and for rows of the matrix $A$, that
is, storing only index-value pairs for all nonzero elements. Then, obviously,
the complexity of arithmetic operations for such vectors will be proportional
to the complexity of a loop with elementary arithmetic operations, the number
of iterations of which is equal to the number of nonzero elements (in the
python programming language, for example, this storage format is implemented
in the method scipy.sparse.csr_matrix [31]).

  3. 3.

Similarly, we will store the value
$\sum_{j=1}^{m}{\exp\left(\left[Ay_{k}\right]_{j}\right)}$, which is the
denominator of the presented expression. Its updating is carried out with the
same complexity as updating a sequence (by calculating the sum of nonzero
terms added to each value from the sequence).

  4. 4.

Since evaluating the specified expression requires evaluating exponent values,
type overflow can occur. To solve this problem, the standard technique is exp-
normalize trick [3]. However, to use it, one should also store the value
$\max_{j=1,...,m}\left[Ay_{k}\right]_{j}$. At the same time, there is no need
to maintain exactly this value, or, otherwise, its approximation to keep the
exponents values small, so this value can be recalculated much rarely: for
example, once in $m$ iterations (in this case, the amortized complexity will
also be equal to $\mathcal{O}(s)$).

So, in the further reasoning, one can assume that the iteration of the
coordinate descent algorithm for solving the corresponding auxiliary problem
has amortized complexity $\mathcal{O}(s)$.

Further, let us consider in more detail the question of the values of the
smoothness constants of this functional. We can write down asymptotic formulas
for $L$ and $\overline{L}=\frac{1}{n}\sum_{i=1}^{n}L_{i}$:

$L=\max_{j=1,...,m}\|A_{j}\|_{2}^{2}=\mathcal{O}(n),\quad\overline{L}=\frac{1}{n}\sum_{i=1}^{n}\max_{j=1,...,m}|A_{ji}|=\mathcal{O}(1).$

Using these estimates, let us refine the computational complexity of the FGM
and Catalyst CDM (CCDM) methods as applied to this problem:

$T_{FGM}=\mathcal{O}\left(s\cdot
n^{3/2}\cdot\frac{1}{\varepsilon^{1/2}}\right),\quad
T_{CCDM}=\widetilde{\mathcal{O}}\left(s\cdot
n\cdot\frac{1}{\varepsilon^{1/2}}\right).$

Thus, in theory, the application of the Catalyst CDM method for solving this
problem allows, in comparison with FGM, to reduce the factor of order
$\mathcal{O}(\sqrt{n})$ in the asymptotic estimate of the computational
complexity. In practice, this means that it is reasonable to apply the
proposed method to problems of large dimensions.

Let us now compare the performance of the proposed approach (Catalyst CDM)
with a number of alternative approaches: Gradient Method (GM), Fast Gradient
Method (FGM), Coordinate Descent Method (CDM) and Accelerated Coordinate
Descent Method (ACDM), by the example of the problem (12) with an artificially
generated matrix $A$ in two different ways. Fig. 2 and 2 present plots of the
convergence of the methods under consideration: in x-axis we show the running
time of the methods in seconds, and in the y-axis we show the residual of the
function on a logarithmic scale ($f_{*}$ calculated by searching for the
corresponding point $x_{*}$ using the FGM method, tuned for an accuracy that
is obviously much higher than that possible to achieve at the selected time
interval).

Figure 1: Convergence of methods for the SoftMax problem (12) with a uniformly
sparse random matrix.

Figure 2: Convergence of methods for the SoftMax problem (12) with
heterogeneously sparse matrix.

In Fig. 2, the case is presented for which all elements of the matrix $A$ are
i.i.d. random variables from the discrete uniform distribution
$A_{ji}\in\mathcal{U}\\{0,1\\}$, the number of nonzero elements is $s\approx
0.2m$, and the parameter $\gamma=0.6$ (as well as in the second case). In this
setting, the proposed method demonstrates faster convergence compared to all
methods under consideration, except the FGM. At the same time, in the setting
shown in Fig. 2, in which the number of nonzero elements, in comparison with
the first case, is increased to $s\approx 0.75m$, and the matrix is generated
heterogeneously in accordance with the rule: $0.9m$ rows with $0.1n$ nonzero
elements and $0.1m$ rows with $0.9n$ nonzero elements, and also one row of the
matrix is completely nonsparse, the proposed method (Accelerated Meta-
algorithm with coordinate descent as internal method) converges faster than
FGM. This is explained by the fact that in this case $L=n$, but $\overline{L}$
is still quite small and, as a result, the constant in the proposed method has
a noticeably smaller effect on the computational complexity than in the case
of FGM. From the results of the experiment, it can also be noted that the
character of its componentwise smoothness affects the efficiency of the
proposed method much more significantly than the sparseness of the problem.

###  2.3 Review of the results of related works

In view of the conceptual simplicity of the scheme of the universal
accelerated proximal envelope, especially in some of its special cases, such
as those described in [18], several variants of its program implementation
were recently proposed, and also experiments were presented that represent the
efficiency of these envelopes in the case of some typical problems. In order
to capture the picture of situation in its entirety and to indicate the
influence of subtle factors of implementation and adjustment of methods on
their effectiveness, it makes sense to analyze some of the remarks and results
of papers that consider ideologically similar problems in a slightly modified
software environment.

Let us refer to the recent review [9], and specifically to the Catalyst clause
of Section 5.6, offering remarks that are useful for the implementation of
Catalyst-like methods, to the class of which the considered in this work
special case of the ‘‘Accelerated Meta-algorithm’’ envelope belongs. One of
the assumptions in this section is that efficient implementations of Catalyst-
type methods (in particular, the implementation from [22]) rely in solving
auxiliary problems on the stop conditions, which are of an absolute nature, as
opposed to relative (multiplicative) criteria, classically used to analyze the
effectiveness of these methods.

The proposed in [22] approach is to calculate the duality gap, it is based on
the construction described in Appendix D.2 of [20]. Following it, based on the
information about the function at the current point, it is possible to obtain
an upper bound for the absolute residual by function. The proposed solution is
really elegant and possible, however, it has a number of drawbacks that do not
allow us to call this solution quite convenient:

  1. 1.

the expression for the duality gap is computationally complex and implies the
calculation of several matrix-vector products for potentially large
dimensions;

  2. 2.

the duality gap in the described version can take on infinite values;

  3. 3.

the proposed relaxed version of the duality gap is even more computationally
difficult, since requires additional implementation of the quadratic knapsack
solver.

At the same time, in this work there is a result that allows us to change the
view on the problem of choosing the accuracy of solving the auxiliary problem
(note that this problem, generally speaking, in the case of many other
approaches remains open). Namely, in Corollary 1 it was proved that to obtain
a solution to an auxiliary problem that satisfies a certain condition
sufficient for its use within the outer shell, it is sufficient to perform a
fixed number of iterations (depending only on the characteristics of the
problem being solved). Moreover, this number (even in a certain strict sense)
is small (if we discard the proportionality of the type $\sim n$, which is
inevitable when using component-wise methods). Thus, following this scheme of
proof, one can completely get away from the need to choose $\varepsilon_{tol}$
to the possibility of choosing the parameter $N$ is the number of iterations
of the auxiliary method, without worrying about changing this parameter with
an increase in the number of external iterations of the method. The practical
benefit of this remark is that for many applied problems it is sometimes
difficult to compute the constants characterizing the problem ($L$,
$\overline{L}$), especially when optimizing functions for which the analytical
expression is not known in advance. In this case, the $\varepsilon_{tol}$
parameter (as well as $N$) becomes the hyperparameter of the algorithm, and
its selection becomes the task of the engineer applying the method. However,
the search space for $N$ is discrete and essentially small, in contrast to the
continually large set of $\varepsilon_{tol}$ values variants, and both with
manual selection or selection using techniques such as grid search, the
problem of selecting $N$ turns out to be much simpler.

Let us consider the results of experiments presented in the previously
mentioned work [22]. They are remarkable in that they assess the acceleration
effect using the Catalyst envelope methods for optimization problems for sum-
like functions with different regularizers, that is, statements close to
machine learning problems, using some practical model examples. The general
conclusion from the experiments carried out in this article is as follows:
accelerating SVRG using the Catalyst envelope has no effect in the case when
the optimal problem regularization parameter is $\approx 1/n$. It is important
to note, however, that in the same situation, the MISO [21] method can be
accelerated. In the work itself, for a different setting, an explanation is
given for this fact: MISO better handle sparse matrices (no need to code lazy
update strategies, which can be painful to implement). At the same time, as
you can see from the figures, the Catalyst loss to directly accelerated SVRG
is small. Hence, we can conclude that in the case when the implementation of
lazy operations for the problem under consideration is not so onerous (the
SoftMax case, considered in detail in this article, is just acceptable,
besides the regularization parameter in this setting is equal to zero), the
use of Catalyst remains quite justified and efficient.

###  2.4 Application to Optimization of Markov Decision Processes

We denote an MDP instance by a tuple
$\mathcal{MDP}:=(\mathcal{S},\mathcal{A},\mathcal{P},r,\gamma)$ with
components defined as follows:

  1. 1.

$\mathcal{S}$ is a finite set of states, $|\mathcal{S}|=n$;

  2. 2.

$\mathcal{A}=\bigcup_{i\in\mathcal{S}}\mathcal{A}_{i}$ is a finite set of
actions that is a collection of sets of actions $\mathcal{A}_{i}$ for states
$i\in\mathcal{S}$, $|\mathcal{A}|=m$;

  3. 3.

$\mathcal{P}$ is the collection of state-to-state transition probabilities
where
$\mathcal{P}:=\\{p_{ij}(a_{i})|i,j\in\mathcal{S},a_{i}\in\mathcal{A}_{i}\\}$;

  4. 4.

$r$ is the vector of state-action transitional rewards where
$r\in[0,1]^{\mathcal{A}}$, $r_{i,a_{i}}$ is the instant reward received when
taking action ai at state $i\in\mathcal{S}$;

  5. 5.

$\gamma$ is the discount factor of MDP, by which one down-weights the reward
in the next future step. When $\gamma\in(0,1)$, we call the instance a
discounted MDP (DMDP) and when $\gamma=1$, we call the instance an average-
reward MDP (AMDP).

Let as denote by $P\in\mathbb{R}^{\mathcal{A}\times\mathcal{S}}$ the state-
transition matrix where its $(i,a_{i})$-th row corresponds to the transition
probability from state $i\in\mathcal{S}$ where $a_{i}\in\mathcal{A}_{i}$ to
state $j$. Correspondingly we use $\hat{I}$ as the matrix with $a_{i}$-th row
corresponding to $e_{i}$, for all $i\in\mathcal{S}$,
$a_{i}\in\mathcal{A}_{i}$. Our goal is to compute a random policy which
determines which actions to take at each state. A random policy is a
collection of probability distributions $\pi:=\\{\pi_{i}\\}_{i\in\mathcal{S}}$
, where $\pi_{i}\in\Delta_{|\mathcal{A}_{i}|}$, $\pi_{i}(a_{j})$ denotes the
probability of taking $a_{j}\in\mathcal{A}_{i}$ at state $i$. One can extend
$\pi_{i}$ to the set of $\Delta_{m}$ by filling in zeros on entries
corresponding to other states $j\neq i$. Given an MDP instance
$\mathcal{MDP}=(\mathcal{S},\mathcal{A},\mathcal{P},r,\gamma)$ and an initial
distribution over states $q\in\Delta_{n}$ , we are interested in finding the
optimal $\pi^{\star}$ among all policies $\pi$ that maximizes the following
cumulative reward $\overline{v}_{\pi}$ of the MDP:

$\displaystyle\pi^{\star}:=\arg\max_{\pi}\overline{v}^{\pi},$
$\displaystyle\overline{v}^{\pi}:=\left\\{\begin{array}[]{*{3}c}&\displaystyle\mathbb{E}^{\pi}\left[\sum_{t=1}^{\infty}\gamma^{t-1}r_{i_{t},a_{t}}\big{|}i_{1}\sim
q\right]&\text{ in the case of }DMDP,\\\
&\displaystyle\lim_{T\rightarrow\infty}\frac{1}{T}\mathbb{E}^{\pi}\left[\sum_{t=1}^{T}r_{i_{t},a_{t}}\big{|}i_{1}\sim
q\right]&\text{ in the case of }AMDP.\end{array}\right.$

For AMDP, $\overline{v}^{\star}$ is the optimal average reward if and only if
there exists a vector $v^{\star}=(v_{i}^{\star})_{i\in\mathcal{S}}$ satisfying
its corresponding Bellman equation [2]:

$\overline{v}^{\star}+v_{i}^{\star}=\max_{a_{i}\in\mathcal{A}_{i}}\left\\{\sum_{j\in\mathcal{S}}p_{ij}(a_{i})v_{j}^{\star}+r_{i,a_{i}}\right\\},\forall
i\in\mathcal{S}.$

For DMDP, one can show that at optimal policy $\pi^{\star}$, each state
$i\in\mathcal{S}$ can be assigned an optimal cost-to-go value $v_{i}^{\star}$
satisfying the following Bellman equation:

$v_{i}^{\star}=\max_{a_{i}\in\mathcal{A}_{i}}\left\\{\sum_{j\in\mathcal{S}}\gamma
p_{ij}(a_{i})v_{j}^{\star}+r_{i,a_{i}}\right\\},\forall i\in\mathcal{S}.$

One can further write the above Bellman equations equivalently as the
following primal linear programming problems.

$\displaystyle\min_{\overline{v},v}\quad\overline{v}$
$\displaystyle\;\text{s.t.}\;\;\;\overline{v}\textbf{1}+(\hat{I}-P)v-r\geq 0$
(LP AMDP) $\displaystyle\min_{v}\quad(1-\gamma)q^{\top}v$
$\displaystyle\;\text{s.t.}\;\;\;(\hat{I}-\gamma P)v+r\geq 0$ (LP DMDP)

By standard linear duality, we can recast the problem formulation using the
method of Lagrangian multipliers, as bi-linear saddle-point (minimax) problem.
The equivalent minimax formulations are

$\displaystyle\min_{v\in\mathbb{R}^{n}}\left\\{F(v)=\max_{\mu\in\Delta_{m}}\left(\mu^{\top}((P-\hat{I})v+r)\right)\right\\}$
(AMDP)
$\displaystyle\min_{v\in\mathbb{R}^{n}}\left\\{F_{\gamma}(v)=\max_{\mu\in\Delta_{m}}\left((1-\gamma)q^{\top}v+\mu^{\top}((\gamma
P-\hat{I})v+r)\right)\right\\}$ (DMDP)

Then, one can apply the Nesterov smoothing technique to the presented max-type
functional, according to [24]. (The calculation is presented for the case of
DMDP, as a more general one. A smoothed version of the AMDP functional is
obtained similarly with the notation $A:=P-\hat{I}$ and $\gamma=1$):

$\displaystyle F_{\gamma}(v)$
$\displaystyle=\max_{\mu\in\Delta_{m}}\left(\underbrace{(1-\gamma)q^{\top}}_{b}v+\mu^{\top}(\underbrace{(\gamma
P-\hat{I})}_{A}v+r)\right)=\max_{\mu\in\Delta_{m}}\left(\sum_{j=1}^{m}\mu_{j}(\left[Av\right]_{j}+r_{j})\right)+\langle
b,v\rangle$
$\displaystyle\rightarrow\max_{\mu\in\Delta_{m}}\left(\sum_{j=1}^{m}\mu_{j}(\left[Av\right]_{j}+r_{j})-\sigma\sum_{j=1}^{m}\mu_{j}\ln\left(\frac{\mu_{j}}{1/m}\right)\right)+\langle
b,v\rangle=$
$\displaystyle=\sigma\ln\left(\sum_{j=1}^{m}\exp\left(\frac{\left[Av\right]_{j}+r_{j}}{\sigma}\right)\right)-\sigma\ln
m-\langle b,v\rangle=:f_{\gamma}(v),\text{ where
}\sigma:=\varepsilon/(2\ln{m}).$

The resulting problem has the form of a SoftMax function, discussed in detail
in the previous section. Taking into account the form of the matrix A, we can
calculate the average component-wise Lipschitz constant:

$P_{ji}\in[0,1],\hat{I}_{ji}\in\\{0,1\\}\quad\Longrightarrow\quad\overline{L}=\frac{1}{\sigma}\frac{1}{n}\sum_{i=1}^{n}\max_{j=1,...,m}|A_{ji}|\leq\frac{\gamma}{\sigma}=\frac{2\gamma\ln
m}{\varepsilon}$

So, one can get the following estimates, which are also given in the Tables 2,
3 (transition from the duality gap accuracy $\varepsilon$ in matrix game
setting to the $\widetilde{\varepsilon}$ accuracy to obtain the
$\widetilde{\varepsilon}$-approximate optimal policy satisfying the condition
in expectation
$\mathbb{E}\overline{v}^{\pi}\geq\overline{v}^{\star}-\widetilde{\varepsilon}$
is carried out according to the rules described in more detail in the work
[15]):

$\displaystyle\varepsilon\sim\frac{1}{2}\cdot\frac{\widetilde{\varepsilon}}{3}\quad\Longrightarrow\quad$
$\displaystyle
T_{CCDM}\;\;=\widetilde{\mathcal{O}}\left(\text{nnz}(P)\sqrt{\log
m}\cdot\widetilde{\varepsilon}^{-1}\right),$
$\displaystyle\varepsilon_{\gamma}\sim\frac{1}{2}\cdot\frac{(1-\gamma)\widetilde{\varepsilon}}{3}\quad\Longrightarrow\quad$
$\displaystyle
T_{CCDM,\gamma}=\widetilde{\mathcal{O}}\left(\gamma^{1/2}(1-\gamma)^{-1}\cdot\text{nnz}(P)\sqrt{\log
m}\cdot\widetilde{\varepsilon}^{-1}\right),$

where $\text{nnz}(P)\leq n\cdot m$ denotes the number of nonzero elements in
matrix $P$. In addition to the result corresponding to the considered method,
the Tables 2, 3 contain complexity bounds of other known approaches to solving
matrix games and MDP problems (for the sake of compactness, it is used the
notation $\text{nnz}^{\prime}(P)=\text{nnz}(P)+(m+n)\log^{3}(mn)$). It can be
seen that for the case of $\gamma=1$, the described approach allowing to
obtain one of the best among the known estimates, and in the case of
$\gamma<1$ it is close in efficiency to many modern approaches. Moreover, to
describe the method used in this article, it was enough to apply only a
special case of the universal accelerated proximal envelope for the classical
coordinate descent method. This approach is conceptually much simpler than the
other methods cited here (which, by the way, are often applicable only to very
particular settings), and allows one to obtain complexity bounds for AMDP
problem that notedly competitive with the best alternatives.

Table 2: Comparison of the effectiveness of approaches ($\gamma=1$ case) Computational complexity | Source
---|---
$\displaystyle\widetilde{\mathcal{O}}\left(\text{nnz}(P)\sqrt{\log m}\cdot\widetilde{\varepsilon}^{-1}\right)$ | this paper
$\displaystyle\widetilde{\mathcal{O}}\left(\text{nnz}(P)\sqrt{m/n}\cdot\widetilde{\varepsilon}^{-1}\right)$ | [6]
$\displaystyle\widetilde{\mathcal{O}}\left(\log^{3}(mn)\;\sqrt{\text{nnz}(P)\cdot\text{nnz}^{\prime}(P)}\cdot\widetilde{\varepsilon}^{-1}\right)$ | [7]
Table 3: Comparison of the effectiveness of approaches ($\gamma\in(0,1)$ case) Computational complexity | Source
---|---
$\displaystyle\widetilde{\mathcal{O}}\left(\gamma^{1/2}(1-\gamma)^{-1}\;\text{nnz}(P)\sqrt{\log m}\cdot\widetilde{\varepsilon}^{-1}\right)$ | this paper
$\displaystyle\widetilde{\mathcal{O}}\left(\gamma(1-\gamma)^{-1}\;\text{nnz}(P)\sqrt{m/n}\cdot\widetilde{\varepsilon}^{-1}\right)$ | [6]
$\displaystyle\widetilde{\mathcal{O}}\left(\gamma(1-\gamma)^{-1}\;\log^{3}(mn)\;\sqrt{\text{nnz}(P)\cdot\text{nnz}^{\prime}(P)}\cdot\widetilde{\varepsilon}^{-1}\right)$ | [7]
$\displaystyle\widetilde{\mathcal{O}}\left(nm\;\left(n+\left(1-\gamma\right)^{-3}\right)\cdot\log{\left(\widetilde{\varepsilon}^{-1}\right)}\right)$ | [32]

##  3 Conclusion

In this paper, we propose a version of the Coordinate Descent Method,
accelerated using the universal proximal envelope ‘‘Accelerated Meta-
algorithm’’. The performed theoretical analysis of the proposed method allows
us to assert that the dependence of its computational complexity on the
dimension of the problem and the required solution accuracy is not inferior to
other methods used to optimize convex Lipschitz smooth functions, and the
computational complexity is comparable to that of the Fast Gradient Method. At
the same time, the proposed scheme retains the properties of the classical
Coordinate Descent Method, including the possibility of using the properties
of componentwise smoothness of the function. The given numerical experiments
confirm the practical efficiency of the method, and also emphasize the
particular relevance of the proposed approach for the problem of optimizing a
SoftMax-like function that often arises in various applications. As an example
of such an application, it was considered the problem of optimizing the MDP,
and using the described approach, a method was proposed for solving the
averaged version of the MDP problem, which gives a complexity bound that
competes with the most efficient known ones.

## References

  * [1] Anikin, A., Dvurechensky, P., Gasnikov, A., Golov, A., Gornov, A., Maximov, Y., Mendel, M., Spokoiny, V.: Modern efficient numerical approaches to regularized regression problems in application to traffic demands matrix calculation from link loads. In: Proceedings of International conference ITAS-2015. Russia, Sochi (2015)
  * [2] Bertsekas, D.P.: Dynamic programming and optimal control. Athena scientific Belmont, MA (1995)
  * [3] Blanchard, P., Higham, D.J., Higham, N.J.: Accurately Computing the Log-Sum-Exp and Softmax Functions (2019)
  * [4] Bubeck, S.: Convex optimization: Algorithms and complexity. arXiv preprint arXiv:1405.4980 (2014)
  * [5] Bubeck, S., Jiang, Q., Lee, Y.T., Li, Y., Sidford, A.: Near-optimal method for highly smooth convex optimization. In: Conference on Learning Theory. pp. 492–507. PMLR (2019)
  * [6] Carmon, Y., Jin, Y., Sidford, A., Tian, K.: Variance reduction for matrix games. In: Advances in Neural Information Processing Systems. pp. 11381–11392 (2019)
  * [7] Carmon, Y., Jin, Y., Sidford, A., Tian, K.: Coordinate methods for matrix games. arXiv preprint arXiv:2009.08447 (2020)
  * [8] Chernov, A.: Direct-dual method for solving the entropy-linear programming problem. Intelligent systems. Theory and applications 20(1), 39–59 (2016)
  * [9] d’Aspremont, A., Scieur, D., Taylor, A.: Acceleration methods (2021)
  * [10] Doikov, N., Nesterov, Y.: Contracting proximal methods for smooth convex optimization. SIAM Journal on Optimization 30(4), 3146–3169 (2020)
  * [11] Dvinskikh, D., Kamzolov, D., Gasnikov, A., Dvurechensky, P., Pasechnyuk, D., Matykhin, V., Chernov, A.: Accelerated meta-algorithm for convex optimization. arXiv preprint arXiv:2004.08691 (2020)
  * [12] Gasnikov, A., Gasnikova, E., Nesterov, Y., Chernov, A.: Efficient numerical methods for entropy-linear programming problems. omputational Mathematics And Mathematical Physics 56(4), 523–534 (2016)
  * [13] Gasnikov, A.: Universal gradient descent. arXiv preprint arXiv:1711.00394 (2017)
  * [14] Ivanova, A., Pasechnyuk, D., Grishchenko, D., Shulgin, E., Gasnikov, A., Matyukhin, V.: Adaptive catalyst for smooth convex optimization. arXiv preprint arXiv:1911.11271 (2019)
  * [15] Jin, Y., Sidford, A.: Efficiently solving mdps with stochastic mirror descent. In: International Conference on Machine Learning. pp. 4890–4900. PMLR (2020)
  * [16] Kamzolov, D., Gasnikov, A., Dvurechensky, P.: On the optimal combination of tensor optimization methods. arXiv preprint arXiv:2002.01004 (2020)
  * [17] Kulunchakov, A., Mairal, J.: A generic acceleration framework for stochastic composite optimization. In: Advances in Neural Information Processing Systems. pp. 12556–12567 (2019)
  * [18] Lin, H., Mairal, J., Harchaoui, Z.: A universal catalyst for first-order optimization. Advances in neural information processing systems 28, 3384–3392 (2015)
  * [19] Lin, H., Mairal, J., Harchaoui, Z.: Catalyst acceleration for first-order convex optimization: from theory to practice. The Journal of Machine Learning Research 18(1), 7854–7907 (2017)
  * [20] Mairal, J.: Sparse coding for machine learning, image processing and computer vision. Ph.D. thesis, Cachan, Ecole normale supérieure (2010)
  * [21] Mairal, J.: Incremental majorization-minimization optimization with application to large-scale machine learning. SIAM Journal on Optimization 25(2), 829–855 (2015)
  * [22] Mairal, J.: Cyanure: An open-source toolbox for empirical risk minimization for python, c++, and soon more. arXiv preprint arXiv:1912.08165 (2019)
  * [23] Monteiro, R.D., Svaiter, B.F.: An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods. SIAM Journal on Optimization 23(2), 1092–1125 (2013)
  * [24] Nesterov, Y.: Smooth minimization of non-smooth functions. Mathematical programming 103(1), 127–152 (2005)
  * [25] Nesterov, Y.: Efficiency of coordinate descent methods on huge-scale optimization problems. SIAM Journal on Optimization 22(2), 341–362 (2012)
  * [26] Nesterov, Y.: Lectures on convex optimization, vol. 137. Springer (2018)
  * [27] Nesterov, Y., Stich, S.U.: Efficiency of the accelerated coordinate descent method on structured optimization problems. SIAM Journal on Optimization 27(1), 110–123 (2017)
  * [28] Paquette, C., Lin, H., Drusvyatskiy, D., Mairal, J., Harchaoui, Z.: Catalyst acceleration for gradient-based non-convex optimization. arXiv preprint arXiv:1703.10993 (2017)
  * [29] Parikh, N., Boyd, S.: Proximal algorithms. Foundations and Trends in optimization 1(3), 127–239 (2014)
  * [30] Rockafellar, R.T.: Monotone operators and the proximal point algorithm. SIAM journal on control and optimization 14(5), 877–898 (1976)
  * [31] SciPy.org: Python Scipy documentation: scipy.sparse.csr_matrix (2020 (accessed January 5, 2021)), https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html
  * [32] Sidford, A., Wang, M., Wu, X., Ye, Y.: Variance reduced value iteration and faster algorithms for solving markov decision processes. In: Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium on Discrete Algorithms. pp. 770–787. SIAM (2018)

## Appendices

##  Appendix 0.A Proof of the Corollary 1

Due to the $(H+L)$-Lipschitz smoothness of the function $F$ it holds the
following inequality: $\displaystyle
F(y_{0})-F(y_{*})\leq\frac{H+L}{2}\|y_{0}-y_{*}\|_{2}^{2}$. Using this
inequality together with the estimation (4), one can write out a condition of
achieving the required accuracy $\widetilde{\varepsilon}$ with respect to the
function value:
$\displaystyle\frac{H+L}{2}\left(1-\frac{1}{\kappa}\right)^{N}\|y_{0}-y_{*}\|_{2}^{2}\leq\widetilde{\varepsilon}$.
Also, the relation $1-1/\kappa\leq\exp\\{-1/\kappa\\}$ holds and therefore one
can state out:

$\displaystyle\frac{H+L}{2}\exp\\{\kappa/N\\}\|y_{0}-y_{*}\|_{2}^{2}\leq\widetilde{\varepsilon}.$

By taking the logarithm and substituting the expression for $\kappa$, we get
an expression for the number of iterations from $\widetilde{\varepsilon}$:

$N(\widetilde{\varepsilon})=\Bigg{\lceil}\frac{Z}{H}\ln{\left\\{\frac{(H+L)\|y_{0}-y_{*}\|_{2}^{2}}{2\widetilde{\varepsilon}}\right\\}}\Bigg{\rceil}.$
(13)

Due to the $H$-strongly convexity of the function $F$ it holds the following
inequality
$\displaystyle\|\overline{y}_{N}-y_{*}\|_{2}^{2}\leq\frac{2}{H}(F(\overline{y}_{N})-F(y_{*}))$,
where $\overline{y}_{N}=\mathbb{E}[y_{N}]$. The function $F$ is convex, and
therefore $F(\overline{y}_{N})\leq\mathbb{E}[F(y_{N})]$ (Jensen inequality),
and from this, together with (3), one can obtain a sufficient condition of
achieving a solution to the auxiliary problem:

$\mathbb{E}[F(y_{N})]-F(y_{*})\leq\frac{H}{2}\left(\frac{H}{3H+2L}\right)^{2}\|y_{0}-y_{*}\|_{2}^{2}.$

Substituting into the formula (13) instead of $\widetilde{\varepsilon}$ the
expression from the right-hand side of this inequality, we arrive at the
desired result.

##  Appendix 0.B Proof of the Theorem 2.4

Let us denote by $N$ the number of iterations of the internal method,
performed in a particular random case until the stop condition is satisfied.
Following the context of the proof of Corollary 1 and applying Markov
inequality, we obtain the following estimate for the probability that the
number of iterations of interest to us will exceed some given $k$:

$\displaystyle Pr(N\geq k)$ $\displaystyle\leq
Pr\left(F(y_{k})-F(y_{*})\geq\frac{H}{2}\left(\frac{H}{3H+2L}\right)^{2}\|y_{0}-y_{*}\|_{2}^{2}\right)$
$\displaystyle\leq\min\left\\{1,\;\;\frac{\mathbb{E}[F(y_{k})]-F(y_{*})}{\frac{H}{2}\left(\frac{H}{3H+2L}\right)^{2}\|y_{0}-y_{*}\|_{2}^{2}}\right\\}\quad\text{(Markov
inequality)}$
$\displaystyle\leq\min\left\\{1,\;\;\left(1+\frac{L}{H}\right)\left(3+\frac{2L}{H}\right)^{2}\exp\left\\{\frac{\kappa}{k}\right\\}\right\\}$

Further, we can estimate by the infinite series the expectation of the number
of iterations of the internal method:

$\displaystyle\mathbb{E}[\hat{N}]$ $\displaystyle=\sum_{k=1}^{\infty}Pr(N\geq
k)$
$\displaystyle\leq\int_{0}^{N(\widetilde{\varepsilon})}dk+\left(1+\frac{L}{H}\right)\left(3+\frac{2L}{H}\right)^{2}\exp\left\\{\kappa\right\\}\cdot\int_{N(\widetilde{\varepsilon})}^{\infty}\exp\left\\{-k\right\\}$
$\displaystyle=N(\widetilde{\varepsilon})+\left(1+\frac{L}{H}\right)\left(3+\frac{2L}{H}\right)^{2}\exp\left\\{\frac{\kappa}{N(\widetilde{\varepsilon})}\right\\}\leq
N(\widetilde{\varepsilon})+1.$

The number of external iterations of the method is constant, which means that
when calculating the mean of the total number of iterations, this value can be
taken out of the mean by applying for each internal iteration a single
estimate given above, whence the proved statement follows.

##  Appendix 0.C Proof of the Theorem 2.5

Let us formulate the following obvious auxiliary proposition.

###### Proposition 1

The algorithmic complexity of the accelerated meta-algorithm with coordinate
descent as an internal method is

$T=\mathcal{O}\left(\widetilde{N}(T_{out}+N_{\delta}T_{inn})\right),$

where $T_{out}$ is the amortized algorithmic complexity of computations
performed at the iteration of the accelerated meta-algorithm, $T_{inn}$ is the
amortized algorithmic complexity of the iteration of the coordinate descent
algorithm.

Let us now reformulate the estimate from the Proposition 1 as follows:

$T=\mathcal{O}\left(\hat{N}\cdot T_{iter}\right),$

where $T_{iter}$ is the amortized algorithmic complexity of an elementary
iteration of the accelerated meta-algorithm, that is, an iteration that can be
both an internal iteration of the coordinate descent algorithm and the main
iteration of the meta-algorithm. The complexity of the main iteration of the
meta-algorithm is determined, first of all, by calculating the full gradient
vector of $f$, and the complexity of this procedure (from the condition of the
theorem) is $\mathcal{O}(s\cdot n)$. At the same time, due to the
$Z=n\overline{L}$, also holds $N_{\delta}=\widetilde{\mathcal{O}}(n)$, where
symbol $\widetilde{\mathcal{O}}(\cdot)$ means the same as
$\mathcal{O}(\cdot)$, but with the possible presence of factors of
$\log(\cdot)$ order. Since the main iteration of the meta-algorithm is
performed every $N_{\delta}$ elementary iterations, where $N_{\delta}$ is
constant, this, using any of the methods of amortization analysis, trivially
yields an amortized algorithmic complexity of the main iteration of the meta-
algorithm, which is $\mathcal{O}(s)$. The complexity of the coordinate descent
algorithm’s iteration (if all operations performing ‘‘in place‘‘, instead of
copying the values of every point, it is quite acceptable for this
construction) is determined by calculating one component of the gradient, and
is also $\mathcal{O}(s)$. From this obtain $T_{iter}=\mathcal{O}(s)$. Using
the estimate (11) and substituting the optimal value of $H$, we obtain the
provided complexity of the method.