# Disorder information from conductance: a quantum inverse problem

S. Mukim1, F. P. Amorim3, A. R. Rocha3, R. B. Muniz4, C. Lewenkopf4 and M. S.
Ferreira1,2 1School of Physics, Trinity College Dublin, Dublin 2, Ireland
2Centre for Research on Adaptive Nanostructures and Nanodevices (CRANN) &
Advanced Materials and Bioengineering Research (AMBER) Centre, Trinity College
Dublin, Dublin 2, Ireland  3 Instituto de Física Teórica, São Paulo State
University, 01140-070, São Paulo, Brazil  4 Instituto de Física, Universidade
Federal Fluminense, Brazil

###### Abstract

It is straightforward to calculate the conductance of a quantum device once
all its scattering centers are fully specified. However, to do this in
reverse, i.e., to find information about the composition of scatterers in a
device from its conductance, is an elusive task. This is particularly more
challenging in the presence of disorder. Here we propose a procedure in which
valuable compositional information can be extracted from the seemingly noisy
spectral conductance of a two-terminal disordered quantum device. In
particular, we put forward an inversion methodology that can identify the
nature and respective concentration of randomly-distributed impurities by
analyzing energy-dependent conductance fingerprints. Results are shown for
graphene nanoribbons as a case in point using both tight-binding and density
functional theory simulations, indicating that this inversion technique is
general, robust and can be employed to extract structural and compositional
information of disordered mesoscopic devices from standard conductance
measurements.

††preprint: APS/123-QED

Structures whose dimensions are comparable to or smaller than the electronic
mean free path display transport features not associated with the classical
ohmic behaviour [1]. These are quantum features found by solving the
Schrödinger equation once the system Hamiltonian is known. Indeed, it is
straightforward to describe how current flows in a quantum device by
calculating its conductance once all scattering centres are specified.
However, to do this in reverse, i.e., to find information about scatterers in
a quantum device by simply looking at its conductance is rather challenging.
This is a kind of Inverse Problem (IP) which consists of obtaining from a set
of observations the causal factors that generated them in the first place. IP
are intrinsic parts of numerous visualization tools [2, 3, 4, 5] but are not
as common in the quantum realm, and even less so in the presence of disorder.

Disorder makes the description of the impurity potential by inverse scattering
methods quite a daunting task. Multiple scattering depends on the scatterers’
locations and therefore likely to affect the electronic dynamics in seemingly
unpredictable ways. To make matters worse theory shows that quantum
interference in chaotic [6] and/or diffusive systems [7] gives rise to
fluctuations whose statistical properties are universal, i.e. system
independent, indicating that standard IP methods are neither practical nor
useful in such situations.

Here we give a different twist to quantum IP approaches and demonstrate that,
instead of detrimental, disorder may be actually beneficial to extracting
information about scattering centres in a quantum device. In particular, we
focus on the energy-dependent conductance of a quantum system, hereafter
referred to as the conductance spectrum, which will serve as the only input of
the inversion procedure described here. This is a quantity normally obtained
by standard experimental setups of a gated two-terminal device but may also be
found by calculation once the underlying Hamiltonian is fully specified. Here
we introduce our inversion methodology by using the latter as a proxy for the
former, i.e, calculated conductance spectra representing their experimental
equivalent. The advantage of using calculated input functions is that we can
refer back to the Hamiltonian that generated them in the first place, making
it possible to assess the success of the inversion procedure.

Let us start by defining the system to be used throughout the manuscript. It
consists of two electrodes separated by a scattering region of length $L$ and
width $W$, as illustrated in Fig. 1(a), a rather typical setup of electronic
transport. The distinction between the leads and the central region is that
the latter contains impurities as scattering centres [8]. In the linear
response regime, the Landauer conductance reads ${\cal G}=2e^{2}/h\int
dE(-\partial f/\partial E)\Gamma(E)$, where $f(E)$ is the Fermi distribution
and $\Gamma(E)$ is the dimensionless conductance (or transmission) given by
[9, 10]

$\Gamma(E)={\rm tr}\big{[}{\bf G}^{r}(E){\bm{\Gamma}}_{R}(E){\bf
G}^{a}(E){\bm{\Gamma}}_{L}(E)\big{]}\,$ (1)

Here ${\bf G}^{r}$ (${\bf G}^{a}=[{\bf G}^{r}]^{\dagger})$ is the full
retarded (advanced) Green function and ${\bm{\Gamma}}_{L}$
(${\bm{\Gamma}}_{R}$) is the line width function accounting for the injection
and lifetime of states in the left (right) contact. For simplicity, we
consider the zero-temperature limit ${\cal G}=2e^{2}/h\,\Gamma(E)$. Thermal
corrections will be shown to have little effect on our procedure.

Figure 1: (a) Schematic diagram of the system under consideration. Two semi-
infinite leads of width $W$ are separated by a region of length $L$ containing
a concentration $n$ of impurities represented by red dots. (b) Calculated
conductance for a system made of a graphene nanoribbon of length $L=100$ unit
cells as a function of the chemical potential (in units of $t$). The dashed
line corresponds to the conductance of the pristine system, i.e., in the
absence of impurities. The solid line shows the conductance for a disordered
configuration containing $N=42$ substitutional impurities ($\epsilon=0.5t$),
which accounts for a concentration of $n=3\%$.

The expression for $\Gamma$ in Eq. (1) is model-independent, i.e., once the
Hamiltonian is known one can find the corresponding Green function and obtain
the energy dependent conductance of the system [9]. We focus on systems with
relatively simple electronic structures, namely graphene nanoribbons (GNR).
Although not an essential requirement, it helps to illustrate the methodology
since GNR are well described by the tight-binding model [11, 12]. In this case
the nearest-neighbour hopping $t$ and the on-site energies fully define the
electronic structure of the nanoribbon. Fig. 1(b) shows the energy-dependent
conductance for an armchair-edged GNR of width $W=3a$, $a=2.7$Å being the
graphene lattice parameter. The dashed line is the conductance spectrum
$\Gamma_{0}$ for the pristine GNR. Results are shown for positive energies
knowing that $\Gamma_{0}(E)=\Gamma_{0}(-E)$.

We now introduce $N$ substitutional impurities. It is convenient to express
the impurity number also as the percentage concentration $n$ defined as
$n=100\times(N/N_{\rm tot})$, where $N_{\rm tot}$ is the total number of sites
in the central region. Both $N$ and $n$ will be used interchangeably. The
scattering strength of the impurities is characterized by the contrast between
their on-site potential $\epsilon$ relative to that of the host, chosen to be
zero. The solid line of Fig. 1(b) shows the conductance for the GNR with
$N=42$ impurities, $\epsilon=0.5\,t$ and $L=100$ unit cells ($n=3\%$).
Impurity locations were randomly selected but kept fixed in the underlying
Hamiltonian that generated the conductance $\Gamma(E)$ of Fig. 1(b).

Being the conductance sensitive to the locations of scattering centres, it is
difficult to devise an inversion tool capable of spatially mapping all
impurities from the conductance spectrum information alone. The brute-force
method of comparing the input conductance $\Gamma(E)$ with those of every
possible disorder configuration is not viable because the number of
combinations is too large for any practical situation. Machine-learning
strategies are currently being attempted to overcome this combinatorial hurdle
[13, 14, 15, 16, 17, 18, 19, 20] but spatial mapping of quantum devices
through inversion remains challenging.

Nevertheless, given the conductance spectrum $\Gamma(E)$ of a device, one
might ask whether it is possible to find the exact number of scattering
centres in it. This may not reveal the position of every impurity but it is a
valuable piece of information and a lot more feasible to obtain. While the
brute-force approach remains impractical, we must account for as many disorder
realizations as possible. We define the configurationally averaged (CA)
conductance $\langle\Gamma\rangle$ by summing over $M$ realizations with the
same number $N$ of impurities (or concentration $n$), i.e.,

$\langle\Gamma(E,n)\rangle=\frac{1}{M}\sum_{j=1}^{M}\Gamma_{j}(E)\,\,,$ (2)

where $j$ labels the different configurations. See SM [21] for a discussion on
the suitable choices for $M$.

The deviation between an arbitrary conductance result and its CA counterpart
is given by

$\Delta\Gamma(E,n)=\Gamma(E)-\langle\Gamma(E,n)\rangle$ (3)

We reiterate that $\Gamma(E)$ acts as the input conductance spectrum of a
single realization and represents the immutable conductance of the device
under investigation. The CA conductance spectrum, on the other hand, reflects
the contribution of very many configurations and depends, in addition, on the
impurity concentration. By treating $n$ as a variable parameter, we can look
for minimization trends in $\Delta\Gamma(E,n)$ that might indicate the real
concentration in the device. Unfortunately, when plotted as a function of $n$
in Fig. 2(a), the deviation $\Delta\Gamma$ for a fixed energy $E$ is
featureless with wide error bars that result from repeating the calculation
1000 times.

Figure 2: (a) Absolute value of the conductance deviation
$|\Delta\Gamma(E,n)|$ as a function of the impurity concentration $n$ (in
percentage) for a fixed energy ($E=0.42t$). (b) Misfit function $\chi$ (in
arbitrary units) as a function of $n$. The vertical (red) dashed line on the
lower part of the panel indicates the real number of impurities used to
generate the sample conductance, which coincides with the minimum of
$\chi(n)$. The (black) dotted line in the upper part of the panel is the
approximate concentration $n^{*}$. Integration limits were ${\cal E}_{-}=0.5t$
and ${\cal E}_{+}=1.5t$ and the solid (blue) line is simply a guide to the
eyes.

However, much cleaner trends are seen when $\Delta\Gamma$ is used in the form
of a functional that measures how good a match $\Gamma(E)$ and
$\langle\Gamma(E,n)\rangle$ are. This quantity is the misfit function
$\chi(n)$ defined as

$\chi(n)=\int_{{\cal E}_{-}}^{{\cal
E}_{+}}dE\,\left[\,\Gamma(E)-\langle\Gamma(E,n)\rangle\,\right]^{2}\,\,,$ (4)

where ${\cal E}_{-}$ and ${\cal E}_{+}$ establish the energy window over which
the integration takes place. Fig. 2(b) shows $\chi$ as a function of $n$ and
displays a more distinctive trend with smaller error bars. The plot indicates
that there is a sweet spot in impurity concentration for which the integrated
deviation is minimal. Remarkably, this agrees with the actual number of
impurities used in the calculation of $\Gamma(E)$, shown as a vertical (red)
dashed line in the lower part of Fig. 2(b). Such a coincidence suggests that
it might be possible to use $\chi(n)$ as an inversion tool to find the number
of impurities in a quantum device from simple conductance measurements. Note
that no prior knowledge of the actual number of impurities was necessary to
identify the misfit-function minimum.

To demonstrate that the agreement between the minimizing concentration of
$\chi$ and the actual impurity concentration is not a fortuitous coincidence,
we proceed to write the CA conductance [10] to linear order in $n$

$\langle\Gamma(E,n)\rangle=\Gamma_{0}(E)-\beta(E)\,n,$ (5)

where $\beta(E)$ is the derivative of $\langle\Gamma(E,n)\rangle$ with respect
to $n$ evaluated at $n=0$. The misfit function $\chi(n)$ will naturally
develop a minimum at $n^{*}=-B/A$, where

$A=\int_{{\cal E}_{-}}^{{\cal E}_{+}}\\!\\!dE\beta^{2}(E)\;\
\mbox{and}\;B=\int_{{\cal E}_{-}}^{{\cal
E}_{+}}\\!\\!dE\beta(E)\big{[}\Gamma(E)-\Gamma_{0}(E)\big{]}.$ (6)

The vertical (black) dotted line shown on the upper part of Fig. 2(b)
indicates the value of $n^{*}$. Note that both the dotted and dashed lines are
aligned, i.e., $n^{*}$ coincides not only with the actual concentration $n$
but also with the impurity concentration that minimizes $\chi(n)$. In fact,
$n^{*}$ is an excellent approximation up to 3% impurity concentration and
provides a simple yet accurate way of identifying the misfit function minimum.
111An alternative definition of $\beta(E)$, shown in the SM [21], may extend
the validity of the approximation to up to 7%.

Besides the impurity concentration, other degrees of freedom can be added to
the IP in question. Let us now consider an arbitrary sample conductance but
this time assume that nothing is known about the impurity, i.e., neither its
concentration $n$ nor its on-site potential $\epsilon$. In this case, the
input $\Gamma(E)$ is simply treated as the conductance of a system with an
unknown number of unspecified impurities. To calculate the misfit function we
must compute the CA conductance which will now depend on $\epsilon$ as well as
$n$. A 2D contour plot of $\chi$ as a function of these two quantities is
shown in Fig. 3(a).

Figure 3: 2D contour plots of the logarithm of the misfit function. Dashed
lines in the plots indicate the values of the respective quantities used in
the underlying Hamiltonians. (a): $\chi$ as a function of $n$ and $\epsilon$.
(b): $\chi$ as a function of $n_{a}$ and $n_{b}$.

A distinctive minimum is seen which, once again, coincides with the exact
values used to generate the input conductance, shown as dashed lines.
Consequently, both the type and concentration of scattering centres inside a
quantum device can be identified through its energy-dependent conductance
fingerprints. Furthermore, still using two degrees of freedom, the IP can also
be implemented in the case of two types of impurities with unknown
concentrations. In other words, two impurities described by known on-site
potentials $\epsilon_{A}$ and $\epsilon_{B}$ are randomly dispersed with
respective concentrations $n_{A}$ and $n_{B}$, which are unknown. Writing the
CA conductance as a function of both concentrations leads to the misfit
function being plotted now as a function of the same quantities in Fig. 3(b).
The dashed lines stand for the actual values $n_{A}$ and $n_{B}$. They
accurately match the concentration values that minimize the misfit function.
This indicates that it is possible to increase the number of degrees of
freedom in the inversion procedure, offering variety and versatility in how we
wish to interrogate the system. While such an increase may lead to the
appearance of more minima in the misfit function, it remains a straightforward
numerical task to identify them all.

Figure 4: (a) The (red) solid line is the conductance spectrum $\Gamma(E)$
obtained from a specific disorder realization of $N=10$ nitrogen impurities,
i.e., $n=0.7\%$. The (blue) dashed line represents the CA conductance
($M=2000$) calculated for the concentration that minimizes the misfit function
$\chi$, shown in (b) as a function of $n$. The misfit function was generated
with an energy window defined by ${\cal E}_{-}=0$ and ${\cal E}_{+}=2\,{\rm
eV}$. Dashed line indicates the exact concentration used to generate the input
conductance.

To show that this method is indeed model-independent we have also tested it in
the case of a system fully described by Density Functional Theory (DFT). In
particular, DFT-based conductance calculations were performed for armchair GNR
of sizes $L=100$ unit cells and $W=3a$ containing a specific spatial
distribution of $N=10$ nitrogen atoms as substitutional impurities, which
corresponds to $n=0.7\%$. See SM [21] for details of the DFT calculations. The
solid line in Fig. 4(a) shows the conductance spectrum $\Gamma(E)$ of the
system obtained within DFT, which will now serve as the input conductance for
the inversion procedure. Regarding the CA conductance, there are two options
on how $\langle\Gamma(E,n)\rangle$ can be obtained. First option is to repeat
the steps taken to generate $\Gamma(E)$ of Fig. 4(a) over several different
disorder configurations. Bearing in mind that the impurity concentration $n$
must be kept as a variable parameter, the averaging procedure defined in Eq.
(2) requires $M$ different configurations for every single value of $n$, which
indicates how computationally demanding this task might become if carried out
entirely within DFT. The alternative option is to use the tight-binding (TB)
model to carry out the CA calculation. In this case the TB model provides a
fast averaging strategy without necessarily compromising in accuracy.

Here, we have selected the latter option and made use of DFT conductance
spectra to identify suitable TB parameters that describe the nitrogen
impurities in the GNR. See SM [21] for details, including a brief discussion
on the choices available to extract TB parameters from DFT calculations [23,
24, 25, 26] and their implementation in large scale transport calculations
[27, 28, 29]. The CA conductance $\langle\Gamma(E,n)\rangle$ is then generated
within the TB model employing $M=2000$ realizations, leading subsequently to
the misfit function $\chi(n)$ of Fig. 4(b). Once again, $\chi(n)$ displays a
minimum at exactly the same nitrogen concentration used to generated the input
conductance $\Gamma(E)$. Another indication of success can be seen by plotting
the CA conductance evaluated at the minimizing concentration of $\chi(n)$.
Shown as a dashed line in Fig. 4(a), it has all the key features of
$\Gamma(E)$ even though both curves were calculated independently. Despite the
computational complexity of obtaining the misfit function entirely within DFT,
this has also been tested and, reassuringly, we find exactly the same answer
(see SM [21] for details).

All input spectra in this manuscript were calculated from known Hamiltonians
because it is then straightforward to test the success of the inversion
method. Ultimately, inversions must be performed based on experimental
conductance data of systems for which we do not have the full Hamiltonian.
This calls for an inversion tool that is reliable, general and robust. The
fact that our inversion strategy works for systems whose electronic structures
are described by a simple TB model as well as by DFT calculations is
indicative of the generality and robustness of this approach.

That disorder is beneficial to this inversion procedure is made evident by the
distinction between the two panels of Fig. 2. Viewed at a fixed energy, the
deviation between the input and the CA conductance spectra reveals very little
about the system. In contrast, the misfit function entails a lot of
information. The efficiency of the procedure relies on an ergodic hypothesis:
conductance fluctuations of a single sample versus energy are related to
sample to sample fluctuations at a fixed energy. More precisely, the ergodic
hypothesis assumes that a running average over a continuous parameter upon
which the conductance depends is equivalent to sampling different impurity
configurations. Here the continuous parameter is the energy, but the concept
can be extended to a range of other quantities, e.g., magnetic field, gate
voltage, etc. A thorough mathematical discussion of this issue is beyond the
scope of this manuscript, but it is worth mentioning that the ergodic
hypothesis can be proven exactly for certain models of disordered systems[30,
31, 32]. That explains why the energy integration induces a distinctive
minimum in $\chi(n)$ since it is analogous to considering a wider universe of
disorder configurations in the CA procedure. The success of the inversion
procedure will therefore depend on how wide the integration range is. See SM
[21] for details on the suitable choices for integration limits.

Another feature of the conductance fluctuations of quantum devices [7, 33] is
that the correlation function $C(\delta X)=\langle\Gamma(X+\delta
X/2)\Gamma(X-\delta X/2)\rangle-\langle\Gamma(X)\rangle^{2}$ is also universal
[31]. The latter is characterized by a correlation length $\delta X_{\rm cor}$
that is system dependent. Therefore, the method works best when the
integration interval $\Delta X$ is such that $\Delta X/\delta X_{\rm cor}\gg
1$. In our case, where we consider variations of the system chemical
potential, $\delta E_{\rm cor}$ scales with the mean level density times the
transmission [34], which makes $\delta E_{\rm cor}$ much larger than typical
experimental temperatures, justifying our earlier claim that temperature has
little effect on our inversion procedure. This also explains why the energy
window used in the misfit-function definition does not have to be very large.

In summary, the inversion procedure presented here provides a mechanism
capable of identifying the composition of scattering centres in a quantum
device by simply looking at the energy dependence of the two-terminal device
conductance. Assuming that the impurity type within the device is known, the
procedure establishes the exact impurity concentration in it. Alternatively,
if no information is known a priori about the scatterers, the inversion
identifies their scattering strength and respective concentration. Finally,
with a mixture of two different types of impurities we are able to establish
the fractional concentration of each component of the device. Despite being
presented with GNR, the technique is not material-specific and performs
remarkably well in the ballistic, diffusive and at the onset of localised
transport regimes. The method is based on the notion that conductance
fluctuations carry little system-specific information, the average conductance
depends smoothly on the variables of interest and on a standard ergodic
hypothesis [30], requiring only an effective (single-particle) description of
the system. 222 The applicability of the method for (strongly) interacting
systems, i.e. beyond the mean field approximation, is not clear yet. This
small set of conditions renders a very robust and versatile methodology that
can extract structural and compositional information of quantum devices from
standard transport measurements.

###### Acknowledgements.

The authors acknowledge support from FAPESP (Grant # FAPESP 2017/02317-2,
2016/01343-7, and 2017/10292-0), CNPq (Grants 308801/2015-6 and
312716/2018-4), FAPERJ (Grant E-26/202.882/2018), and the ICTP-Simons
Foundation Associate Scheme. The authors also thank the National Laboratory
for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of
the SDumont supercomputer, which have contributed to the results reported
here. Discussions with S. R. Power and M. Kucukbas are greatly acknowledged.

## References

  * Datta [2005] S. Datta,  _Quantum Transport: Atom to Transistor_ (Cambridge University Press, 2005).
  * Bertero and Piana [2006] M. Bertero and M. Piana, Inverse problems in biomedical imaging: modeling and methods of solution, in  _Complex Systems in Biomedicine_ , edited by A. Quarteroni, L. Formaggia, and A. Veneziani (Springer Milan, Milano, 2006) pp. 1–33.
  * Virieux  _et al._ [2017] J. Virieux, A. Asnaashari, R. Brossier, L. Metivier, A. Ribodetti, and W. Zhou, An introduction to full waveform inversion, in  _Encyclopedia of Exploration Geophysics_ (2017) p. 1.
  * Tromp  _et al._ [2008] J. Tromp, D. Komatitsch, and Q. Liu, Spectral-element and adjoint methods in seismology, Communications in Computational Physics 3, 1 (2008).
  * Dias and Vieira Neto [2004] E. Dias and H. Vieira Neto, A novel approach to environment mapping using sonar sensors and inverse problems, in  _Lecture Notes in Computer Science_ , Vol. 9287, edited by C. Dixon and K. Tuyls (Springer, 2004).
  * Lewenkopf and Weidenmüller [1991] C. H. Lewenkopf and H. A. Weidenmüller, Stochastic versus semiclassical approach to quantum chaotic scattering, Ann. Phys. 212, 53 (1991).
  * Lee and Stone [1985] P. A. Lee and A. D. Stone, Universal conductance fluctuations in metals, Phys. Rev. Lett. 55, 1622 (1985).
  * Tapasztó  _et al._ [2008] L. Tapasztó, G. Dobrik, P. Lambin, and L. P. Biro, Tailoring the atomic structure of graphene nanoribbons by scanning tunnelling microscope lithography, Nat. Nanotechnol. 3, 397 (2008).
  * Meir and Wingreen [1992] Y. Meir and N. S. Wingreen, Landauer formula for the current through an interacting electron region, Phys. Rev. Lett. 68, 2512 (1992).
  * Datta [1995] S. Datta,  _Electronic Transport in Mesoscopic Systems_ (Cambridge University Press, Cambridge, 1995).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Das Sarma  _et al._ [2011] S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Rev. Mod. Phys. 83, 407 (2011).
  * Lopez-Bezanilla and von Lilienfeld [2014] A. Lopez-Bezanilla and O. A. von Lilienfeld, Modeling electronic quantum transport with machine learning, Phys. Rev. B 89, 235411 (2014).
  * Hansen  _et al._ [2013] K. Hansen, G. Montavon, F. Biegler, S. Fazli, M. Rupp, M. Scheffler, O. A. Von Lilienfeld, A. Tkatchenko, and K.-R. Müller, Assessment and validation of machine learning methods for predicting molecular atomization energies, J. Theor. Comput. Chem. 9, 3404 (2013).
  * Himmetoglu [2016] B. Himmetoglu, Tree based machine learning framework for predicting ground state energies of molecules, J. Chem. Phys. 145, 134101 (2016).
  * Li  _et al._ [2018] H. Li, C. Collins, M. Tanha, G. J. Gordon, and D. J. Yaron, A density functional tight binding layer for deep learning of chemical hamiltonians, J. Chem. Theory Comput. 14, 5764 (2018).
  * Panosetti  _et al._ [2019] C. Panosetti, A. Engelmann, L. Nemec, K. Reuter, and J. T. Margraf, Learning to use the force: Fitting repulsive potentials in density-functional tight-binding with gaussian process regression (2019).
  * Xia and Kais [2018] R. Xia and S. Kais, Quantum machine learning for electronic structure calculations, Nat. Commun. 9, 4195 (2018).
  * Dral [2020] P. O. Dral, Quantum chemistry in the age of machine learning, J. Phys. Chem. Lett. 11, 2336 (2020).
  * Carrasquilla and Melko [2017] J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nat. Phys. 13, 431 (2017).
  * [21] See Supplemental Material at [URL will be inserted by publisher] for more details.
  * Note [1] An alternative definition of $\beta(E)$, shown in the SM [21], may extend the validity of the approximation to up to 7%.
  * Gresch  _et al._ [2017] D. Gresch, G. Autes, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, Z2pack: Numerical implementation of hybrid wannier centers for identifying topological materials, Phys. Rev. B 95, 075146 (2017).
  * Pizzi  _et al._ [2020] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications, J. Phys.: Condens. Matter 32, 165902 (2020).
  * Nardelli  _et al._ [2018] M. B. Nardelli, F. T. Cerasoli, M. Costa, S. Curtarolo, R. De Gennaro, M. Fornari, L. Liyanage, A. R. Supka, and H. Wang, PAOFLOW: A utility to construct and operate on ab initio hamiltonians from the projections of electronic wavefunctions on atomic orbital bases, including characterization of topological materials, Comput. Mater. Sci. 143, 462 (2018).
  * Vanderbilt [2018] D. Vanderbilt, The pythtb package, in  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018) pp. 327–362.
  * Groth  _et al._ [2014] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, Kwant: a software package for quantum transport, New J. Phys. 16, 063065 (2014).
  * Lima  _et al._ [2018] L. Lima, A. Dusko, and C. Lewenkopf, Efficient method for computing the electronic transport properties of a multiterminal system, Phys. Rev. B 97, 165405 (2018).
  * João  _et al._ [2020] S. M. João, M. Anđelković, L. Covaci, T. G. Rappoport, J. M. Lopes, and A. Ferreira, Kite: high-performance accurate modelling of electronic structure and response functions of large molecules, disordered crystals and heterostructures, R. Soc. Open Sci. 7, 191809 (2020).
  * French  _et al._ [1978] J. French, P. Mello, and A. Pandey, Ergodic behavior in the statistical theory of nuclear reactions, Phys. Lett. B 80, 17 (1978).
  * Brézin and Zee [1994] E. Brézin and A. Zee, Correlation functions in disordered systems, Phys. Rev. E 49, 2588 (1994).
  * Guhr  _et al._ [1998] T. Guhr, A. Müller-Groeling, and H. A. Weidenmüller, Random-matrix theories in quantum physics: common concepts, Phys. Rep. 299, 189 (1998).
  * Lee  _et al._ [1987] P. A. Lee, A. D. Stone, and H. Fukuyama, Universal conductance fluctuations in metals: Effects of finite temperature, interactions, and magnetic field, Phys. Rev. B 35, 1039 (1987).
  * Alves and Lewenkopf [2002] E. R. P. Alves and C. H. Lewenkopf, Conductance fluctuations and weak localization in chaotic quantum dots, Phys. Rev. Lett. 88, 256805 (2002).
  * Note [2] The applicability of the method for (strongly) interacting systems, i.e. beyond the mean field approximation, is not clear yet.


# Towards Extremely Fast Bilevel Optimization with Self-governed Convergence
Guarantees

Risheng Liu${}^{\ast 1}\quad$ Xuan Liu${}^{1}\quad$ Wei Yao${}^{2}\quad$
Shangzhi Zeng${}^{3}\quad$ Jin Zhang2,4
1International School of Information Science & Engineering, DUT
2Department of Mathematics, SUSTech 3Department of Mathematics and Statistics,
UVIC
4National Center for Applied Mathematics Shenzhen

###### Abstract

Gradient methods have become mainstream techniques for Bi-Level Optimization
(BLO) in learning and vision fields. The validity of existing works heavily
relies on solving a series of approximation subproblems with extraordinarily
high accuracy. Unfortunately, to achieve the approximation accuracy requires
executing a large quantity of time-consuming iterations and computational
burden is naturally caused. This paper is thus devoted to address this
critical computational issue. In particular, we propose a single-level
formulation to uniformly understand existing explicit and implicit Gradient-
based BLOs (GBLOs). This together with our designed counter-example can
clearly illustrate the fundamental numerical and theoretical issues of GBLOs
and their naive accelerations. By introducing the dual multipliers as a new
variable, we then establish Bilevel Alternating Gradient with Dual Correction
(BAGDC), a general framework, which significantly accelerates different
categories of existing methods by taking specific settings. A striking feature
of our convergence result is that, compared to those original unaccelerated
GBLO versions, the fast BAGDC admits a unified non-asymptotic convergence
theory towards stationarity. A variety of numerical experiments have also been
conducted to demonstrate the superiority of the proposed algorithmic
framework.

##  1 Introduction

Bi-Level Optimization (BLO) has recently attracted growing interests due to
its wide applications in the field of deep learning [1], especially for
hyperparameter optimization [2, 3, 4], meta learning [5, 6, 7], neural
architecture search [8, 9, 10], adversarial learning [11], and reinforcement
learning [12], etc. BLO tackles nested optimization structures appearing in
applications. In the last decade, BLO has emerged as a prevailing optimization
technique for modern machine learning tasks with an underlying hierarchy.

Recently, Gradient-based BLOs (GBLOs) have become mainstream techniques for
BLO in learning and vision fields. In this type, an approximate version of BLO
can be solved by taking into explicit account the optimization dynamics for
the Lower-Level (LL) objective. The issue of the quality of this approximation
has been well studied in [2, 13, 14, 15]. A major advantage of GBLOs is that
an approximate hypergradient can be directly calculated by automatic
differentiation based on the trajectory of the LL variable. However, GBLOs
theoretically require a large number of LL optimization steps to achieve the
desired accuracy [13]. Unfortunately, to achieve the accuracy requires
executing a large quantity of time-consuming iterations. Thus, computational
and memory burden are naturally caused. In practice, one simple heuristic to
accelerate is to only perform a small number of iterations for the LL problem,
see, e.g., neural architecture search [8, 16], adversarial learning [17] and
reinforcement learning [11]. However, the theoretical properties are still
unclear, although competitive empirical performance has been reported.

On the other hand, to remove the need for differentiating through the LL
optimization path, especially when the LL dynamic system iterates many times,
a type of implicit gradient-based BLO methods (called I-GBLOs for simplicity)
is employed in [18, 6, 19] for hyperparameter optimization and meta-learning.
In fact, I-GBLOs first leverages the implicit differentiation to derive an
analytical expression of hypergradient, and then solves inexactly the LL
problem up to a tolerance and performs an approximate matrix inversion
accurately. Hence, I-GBLOs requires executing a large quantity of iterations
and computational burden is naturally caused.

###  1.1 Our Motivations and Contributions

As mentioned above, to access the hypergradient accurately, the number of
iteration steps in both GBLOs and I-GBLOs theoretically goes to infinity [13,
18, 6]. Hence, computational burden is naturally caused. One simple heuristic
to accelerate is to only performs one-step iteration for the LL problem, see,
e.g., DARTS [8] and its applications for various tasks, cf. e.g., neural
architecture search [8, 16], adversarial learning [17], reinforcement learning
[11], hyperparameter optimization [20] , meta learning [21] and image
processing [22]. Although competitive empirical performance has been reported
widely, the corresponding theoretical properties are still unclear. Recently,
a broad collection of algorithms have been proposed to solve (stochastic) BLOs
with using the one-step LL problem iteration idea, see, e.g., [23, 24, 25].
However, a large quantity of time-consuming matrix multiplications or Hessian-
vector products is required there. Moreover, the strong convexity of the LL
objective is required for obtaining the convergence results of all these
proposed algorithms.

In this paper, we first construct an interesting example in Section 2.1
explicitly to indicate that existing approaches including GBLOs and I-GBLOs
with adopted naive acceleration may lead to incorrect solutions. Therefore, an
important question which motives this work arises. Can we design a extremely
fast bilevel acceleration strategy which enjoy favorably both extremely low
computational cost and reliable theoretical guarantee?

In response to the above question, we first replace the LL problem with its
first-order optimality condition and the resulting problem is a single-level
optimization problem with equality constraints.. Thus we can introduce
naturally the dual multipliers as a new variable. Our first insight is to well
exploit the bilevel structure, keeping the update of the LL variable by
gradient descent of the LL objective and its variants. Our second insight is
that a prospective point where the hypergradient of BLO vanishes happens to be
a KKT solution of the above single-level optimization problem. Thereby, to
guarantee the non-asymptotic convergence, we update cautiously the Upper-Level
(UL) and dual multipliers variables with the LL variable, alternatively
according to the KKT condition. Compared with naive accelerations, the newly
involved multipliers variable can be regarded as a dual correction. Therefore,
we establish Bilevel Alternating Gradient with Dual Correction (BAGDC), an
extremely fast BLO algorithm, to ease both computational and memory burden.
Furthermore, compared to those original unaccelerated GBLO and I-GBLO
versions, we analyze the non-asymptotic convergence behavior of BAGDC towards
self-governed KKT stationary solutions.

The main contributions of this paper can be summarized as follows:

  * •

We first demonstrate fundamental theoretical issues arising from the naive
single-inner-step acceleration which has been widely-witnessed in many real-
world learning and vision applications, by explicitly exhibiting the
trajectory of GBLOs on our designed counter-example.

  * •

By reformulating BLO as a single-level optimization with equality constraints
and introducing dual multipliers from the viewpoint of KKT condition, we then
establish BAGDC, an alternating gradient method, to successfully avoid solving
a series of time-consuming nested optimization subproblems, and hence
remarkably faster than a variety of existing GBLOs and I-GBLOs.

  * •

An important feature of our theoretical result relies on a series of new non-
asymptotic convergence guarantees of BAGDC towards KKT stationary points.
Compared to exiting works which are either lack of convergence towards
stationarity or convergence rate, our result is evidently stronger and general
enough to cover BLOs with multiple LL minimizers.

##  2 Bilevel Alternating Gradient with Dual Correction

In this work, we consider the BLO problem in the following form:

$\mathop{\min}\limits_{\mathbf{x},\mathbf{y}}F(\mathbf{x},\mathbf{y}),\quad\mathrm{s.t.}\quad\mathbf{y}\in\mathcal{S}(\mathbf{x}):=\arg\min_{\mathbf{y}\in\mathbb{R}^{m}}f(\mathbf{x},\mathbf{y}),$
(1)

where $\mathbf{x}\in\mathbb{R}^{n}$, $\mathbf{y}\in\mathbb{R}^{m}$ are UL and
LL variables respectively. Here the UL and LL objectives are smooth functions
with Lipschitz continuous first and second order derivatives.111This
assumption is common in the BLO literature, see, e.g., [26, 15, 24, 23]. A
formal description of the assumption is stated in the Supplemental Material.

###  2.1 Fundamental Issues and Counter-Example

BLOs are usually considered as nested structural optimization problems. Thus,
the existing BLO methods including GBLOs and I-GBLOs always update the UL and
LL variables in a nested manner. This roughly causes the computational and
theoretical issues for the existing methods. To design an extremely fast BLO
algorithm, we aim to update the UL and LL variables by alternating gradient,
cf. the well-known Gradient Descent-Ascent (GDA) for minimax optimization
problems. But alternating gradient for BLOs is challenging since naive
accelerations of the existing approaches may lead to incorrect solutions. Let
us take a widely used acceleration technique to exemplify this issue. Indeed,
if one only performs one-step LL iteration in GBLOs, the resulting algorithm,
called Naive One-Step Acceleration (NOSA for simplicity of notation), reads as

$\displaystyle\mathbf{y}_{k+1}=$
$\displaystyle\mathbf{y}_{k}-\beta\nabla_{\mathbf{y}}f(\mathbf{x}_{k},\mathbf{y}_{k}),$
(2) $\displaystyle\mathbf{x}_{k+1}=$
$\displaystyle\mathbf{x}_{k}-\alpha\left(\nabla_{\mathbf{x}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})-\beta\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\mathbf{x}_{k},\mathbf{y}_{k})\nabla_{\mathbf{y}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})\right),$

where $\alpha,\beta$ are stepsizes. Next, we design an interesting counter-
example (Example 1 below) to illustrates that NOSA in Eq. (2) may lead to
incorrect solutions.

######  Example 1.

(Counter-Example) Consider the following BLO problem:

$\displaystyle\min_{\mathbf{x}\in\mathbb{R}^{n}}\frac{1}{2}\|\mathbf{x}-\mathbf{z}_{0}\|^{2}+\frac{1}{2}\mathbf{y}^{*}(\mathbf{x})^{\top}\mathbf{A}\mathbf{y}^{*}(\mathbf{x})\quad\mathrm{s.t.}\quad\mathbf{y}^{*}(\mathbf{x})=\arg\min_{\mathbf{y}\in\mathbb{R}^{n}}f(\mathbf{x},\mathbf{y})=\frac{1}{2}\mathbf{y}^{\top}\mathbf{A}\mathbf{y}-\mathbf{x}^{\top}\mathbf{y},$
(3)

where $\mathbf{x}\in\mathbb{R}^{n}$, $\mathbf{y}\in\mathbb{R}^{n}$,
$\mathbf{A}\in\mathbb{S}^{n\times n}$ is a positive-definite symmetric matrix
and $\mathbf{z}_{0}\neq 0$ is a given point in $\mathbb{R}^{n}$. Note that the
LL objective is strongly convex w.r.t. the LL variable and
$\mathbf{x},\mathbf{y}$ have the same dimension. Assuming that
$\left(\mathbf{x}_{k},\mathbf{y}_{k}\right)$ in Eq. (2) converges to
$(\bar{\mathbf{x}},\bar{\mathbf{y}})$. Then
$\mathbf{A}\bar{\mathbf{y}}-\bar{\mathbf{x}}=0$ and
$(\bar{\mathbf{x}}-\mathbf{z}_{0})+\beta\mathbf{A}\bar{\mathbf{y}}=0$. This
implies that $\bar{\mathbf{x}}=\mathbf{z}_{0}/(1+\beta)$. Now we show that
$\bar{\mathbf{x}}$ is not a stationary point of
$\varphi(\mathbf{x}):=F(\mathbf{x},\mathbf{y}^{*}(\mathbf{x}))$ if
$\mathbf{z}_{0}$ is not a eigenvector of $\mathbf{A}$. Indeed, a simple
computation gives
$\nabla\varphi(\bar{\mathbf{x}})=(\bar{\mathbf{x}}-\mathbf{z}_{0})+\mathbf{A}^{-1}\bar{\mathbf{x}}=\left(\mathbf{A}^{-1}\mathbf{z}_{0}-\beta\mathbf{z}_{0}\right)/(1+\beta)\neq
0$.

By using the same counter-example 1 and a similar argument, one can easily
check that the same result also holds for other naive one-step accelerations
of I-GBLOs, such as Naive One-Step Accelerations of Conjugate Gradient (CG)
and Neumann series (NS) methods. For more details, we refer the reader to the
Supplemental Material.

###  2.2 A Single-Level Reformulation and Dual Correction

To figure out the reason why NOSA in Eq. (2) does not work well, we first
transform BLO into a single-level optimization problem with equality
constraint:

$\min_{\mathbf{x},\mathbf{y}}\
F(\mathbf{x},\mathbf{y})\quad\mathrm{s.t.}\quad\nabla_{\mathbf{y}}f(\mathbf{x},\mathbf{y})=0.$
(4)

Clearly, BLO in Problem (1) and Problem (4) are equivalent when the LL
objective $f(\mathbf{x},\mathbf{y})$ is convex w.r.t. $\mathbf{y}$ for any
fixed $\mathbf{x}$. Denote
$\mathcal{L}(\mathbf{x},\mathbf{y},\mathbf{v})=F(\mathbf{x},\mathbf{y})-\mathbf{v}^{T}\nabla_{\mathbf{y}}f(\mathbf{x},\mathbf{y})$,
where $\mathbf{v}\in\mathbb{R}^{m}$ is the dual multipliers. We say that
$(\mathbf{x}^{*},\mathbf{y}^{*})$ is a KKT point of Problem (4) if there
exists a dual multipliers $\mathbf{v}^{*}\in\mathbb{R}^{m}$ such that the
following KKT condition holds:

$\displaystyle\nabla\mathcal{L}(\mathbf{x}^{*},\mathbf{y}^{*},\mathbf{v}^{*})=\left(\begin{array}[]{ccc}\nabla_{\mathbf{x}}F(\mathbf{x}^{*},\mathbf{y}^{*})-\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\mathbf{x}^{*},\mathbf{y}^{*})\mathbf{v}^{*}\\\
\nabla_{\mathbf{y}}F(\mathbf{x}^{*},\mathbf{y}^{*})-\nabla_{\mathbf{y}\mathbf{y}}^{2}f(\mathbf{x}^{*},\mathbf{y}^{*})\mathbf{v}^{*}\\\
-\nabla_{\mathbf{y}}f(\mathbf{x}^{*},\mathbf{y}^{*})\end{array}\right)=0.$ (8)

To better describe what we are going to do with the KKT condition, we briefly
provide a strong connection between the KKT condition and the hypergradient
using in the existing BLO methods. In fact, the first part (called
$\mathbf{x}$-counterpart) of the KKT condition is in consistence with the
hypergradient, which is equal to
$\nabla_{\mathbf{x}}F(\mathbf{x}^{*},\mathbf{y}^{*})-\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\mathbf{x}^{*},\mathbf{y}^{*})\mathbf{v}^{*}$
when both the second part (called $\mathbf{y}$-counterpart) and the third one
(called feasible condition) of the KKT condition hold, if the Hessian
$\nabla_{\mathbf{y}\mathbf{y}}^{2}f(\mathbf{x}^{*},\mathbf{y}^{*})$ is
nonsingular. The latter condition is common in the literature. But the
nonsingularity of
$\nabla_{\mathbf{y}\mathbf{y}}^{2}f(\mathbf{x}^{*},\mathbf{y}^{*})$ is not
necessary for the KKT condition. Therefore, the KKT condition is more
applicable to general BLOs, e.g., BLOs with multiple LL minimal points.

Now we can uniformly understand the existing explicit and implicit Gradient-
based BLOs. From the viewpoint of non-asymptotic convergence analysis, roughly
speaking, the existing BLO methods can be divided into three steps: (1) given
UL variable, update the LL variable by solving the LL optimization problem up
to tolerance, which satisfies prescribed dynamic accuracy requirements. From
the view point of the KKT condition, this step is to guarantee the feasible
condition; (2) given UL and LL variable, update the dual multipliers by
approximating the correct dual multipliers accurately. There are mainly two
kinds of techniques, i.e., automatic differentiation based on the trajectory
of LL variable in RHG [2] and solving a linear system based on the Conjugate
Gradient (CG) or Neumann Series (NS) [6, 19, 14]. The correct dual multipliers
here is the one satisfying the $\mathbf{y}$-counterpart of the KKT condition;
(3) update the UL variable by the accurate approximate hypergradient. As we
discussed in the previous paragraph, this corresponds to the
$\mathbf{x}$-counterpart of the KKT condition. From the above analysis, one
can clearly see the reason why NOSA in Eq. (2) does not work well. In fact,
NOSA did not update the dual multipliers properly. The hidden update
$\mathbf{v}_{k+1}=\beta\nabla_{\mathbf{y}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})$
in NOSA can not guarantee the $\mathbf{y}$-counterpart of the KKT condition.

The above issue motivate us to consider the following question: how to update
the dual multipliers if the update of the LL variable by gradient descent of
the LL objective is reserved? This question is important and nontrivial since
the dual multipliers are always updated by the equality constraint in the
standard methods of solving nonlinear optimization with equality constraints,
see, e.g., [27, Chapter17]. The reason why we prefer updating the LL variable
by
$\mathbf{y}_{k+1}=\mathbf{y}_{k}-\beta\nabla_{\mathbf{y}}f(\mathbf{x}_{k},\mathbf{y}_{k})$
and its variants is that the equality constraint of Problem (4) comes from the
LL minimization problem of BLO. Hence, to well exploit the bilevel structure,
it is reasonable for us to keep the update of the LL variable by gradient
descent. This can also guarantee the feasible condition at the limit point. We
also keep the update of the UL variable by
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\alpha\left(\nabla_{\mathbf{x}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})-\left[\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\mathbf{x}_{k},\mathbf{y}_{k})\right]\mathbf{v}_{k+1}\right)$
and its variants to ensure the $\mathbf{x}$-counterpart of the KKT condition.
Hence, to guarantee the convergence of algorithm towards self-governed KKT
stationary solutions, the update of the multipliers variable should be to
ensure the $\mathbf{y}$-counterpart of the KKT condition.

When the LL problem is merely convex, BLO can be with multiple LL minimizers.
Then updating of the LL variable by gradient descent of the LL objective may
not necessarily generate a sequence towards minimizer of UL objective over LL
problem solutions. Inspired by BDA [13], we use aggregation parameters to
aggregate both UL and LL objectives. Specifically, we define

$\psi_{\mu}(\mathbf{x},\mathbf{y}):=\mu\lambda
F(\mathbf{x},\mathbf{y})+(1-\mu)f(\mathbf{x},\mathbf{y}),$ (9)

where $\mu,\lambda$ are aggregation parameters. We then replace the LL
objective $f$ by $\psi_{\mu_{k}}$ and let the aggregation parameter $\mu_{k}$
goes to $0$ as the iteration goes on. When the LL problem is strongly convex
w.r.t. $\mathbf{y}$, we just take the aggregation parameter $\mu=0$ in the
following algorithm. Now, we propose the following Bilevel Alternating
Gradient with Dual Correction (BAGDC) in Algorithm 1 for BLOs.

Algorithm 1 Bilevel Alternating Gradient with Dual Correction (BAGDC)

Initialize: $\mathbf{x}_{0},\mathbf{y}_{0},\mathbf{v}_{0}$, stepsizes
$\alpha_{k},\beta_{k},\eta_{k}$, aggregation parameters $\mu_{k},\lambda$;
Let: $\psi_{\mu_{k}}(\mathbf{x},\mathbf{y}):=\mu_{k}\lambda
F(\mathbf{x},\mathbf{y})+(1-\mu_{k})f(\mathbf{x},\mathbf{y})$;

1:for $k=0,1,\dots,K-1$ do

2: update
$\mathbf{y}_{k+1}=\mathbf{y}_{k}-\beta_{k}\nabla_{\mathbf{y}}\psi_{\mu_{k}}(\mathbf{x}_{k},\mathbf{y}_{k})$;

3: update
$\mathbf{v}_{k+1}=\mathbf{v}_{k}+\eta_{k}\left(\nabla_{\mathbf{y}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})-\left[\nabla_{\mathbf{y}\mathbf{y}}^{2}\psi_{\mu_{k}}(\mathbf{x}_{k},\mathbf{y}_{k+1})\right]\mathbf{v}_{k}\right);$

4: update
$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\alpha_{k}\left(\nabla_{\mathbf{x}}F(\mathbf{x}_{k},\mathbf{y}_{k+1})-\left[\nabla_{\mathbf{x}\mathbf{y}}^{2}\psi_{\mu_{k}}(\mathbf{x}_{k},\mathbf{y}_{k})\right]\mathbf{v}_{k+1}\right).$

5:end for

To help understand BAGDC in Algorithm 1, the following points should be
emphasized.

  * •

In BAGDC, to handle a general BLO, e.g., BLOs with multiple LL minimizers, we
use aggregation parameters $\mu_{k},\lambda$ to aggregate the LL objective $f$
and a positive multiple of the UL objective $\lambda F$. An explicitly
strategy for choosing stepsizes $\alpha_{k},\beta_{k},\eta_{k}$ and
aggregation parameters $\mu_{k},\lambda$ is presented in Section 3 based on
theoretical convergence analysis. Specially, when the LL objective is strongly
convex w.r.t. $\mathbf{y}$, the aggregation parameter $\mu$ is taken to be
$0$. In this case, the stepsizes $\alpha,\beta,\eta$ can be taken to be
constant stepsizes satisfying certain appropriate conditions given in Section
3.

  * •

Compared with naive accelerations, e.g., NOSA in Eq. (2), the newly involved
multipliers variable $\mathbf{v}$ can be regarded as a dual correction. Its
appearance turns bilevel alternating gradient into a practical acceleration
strategy which enjoy favorably both extremely simpler implementation and
theoretical convergence guarantee. For more details we refer the reader to the
theoretical investigation in Section 3 and experiments in Section 4.

##  3 Theoretical Investigation

This section is devoted to the convergence analysis of our proposed algorithm
BAGDC with and without strong convexity assumption of LL problem. The UL and
LL objectives $F(\mathbf{x},\mathbf{y})$ and $f(\mathbf{x},\mathbf{y})$ are
assumed to be smooth functions with Lipschitz continuous first and second
order derivatives.222This assumption is common in the BLO literature. A formal
description of the assumption is stated in the Supplemental Material. Please
notice that all the proofs of our theoretical results given here and other
auxiliary results are stated in the Supplemental Material.

###  3.1 General Results for LL Merely Convex Case

To handle general BLOs, e.g., BLOs with multiple LL minimizers, we take the
following standard assumption as in [13].

######  Assumption 3.1.

For any fixed $\mathbf{x}$, $f(\mathbf{x},\cdot)$ is convex and
$F(\mathbf{x},\cdot)$ is $\sigma_{F}$-strongly convex.

Recall that

$\psi_{\mu}(\mathbf{x},\mathbf{y}):=\mu\lambda
F(\mathbf{x},\mathbf{y})+(1-\mu)f(\mathbf{x},\mathbf{y}),$

where $\mu,\lambda$ are aggregation parameters. Let $0<\mu\leq 1/2$ and
$\lambda>0$, for any fixed $\mathbf{x}$, as $F(\mathbf{x},\cdot)$ is assumed
to be $\sigma_{F}$-strongly convex, the aggregation function
$\psi_{\mu}(\mathbf{x},\mathbf{y})$ is $\sigma_{\psi_{\mu}}$-strongly convex
in $\mathbf{y}$ with $\sigma_{\psi_{\mu}}=\mu\lambda\sigma_{F}>0$. Hence
$\psi_{\mu}(\mathbf{x},\mathbf{y})$ has a unique minimal point for every
$\mathbf{x}$, denoted by $\mathbf{y}^{*}_{\mu}(\mathbf{x})$. Under the
smoothness assumption of $f$ and $F$, $\psi_{\mu}(\mathbf{x},\mathbf{y})$ is
also $L_{\psi_{\mu}}$-smooth in $\mathbf{y}$ for any fixed $\mathbf{x}$.

To analyze the convergence of BAGDC, we define the following Lyapunov function
for the generated sequence
$\\{(\mathbf{x}_{k},\mathbf{y}_{k},\mathbf{v}_{k})\\}$ by BAGDC 1:

$V_{k}:=F(\mathbf{x}_{k},\mathbf{y}^{*}_{\mu_{k}}(\mathbf{x}_{k}))+\frac{1}{2}\|\mathbf{y}_{k}-\mathbf{y}^{*}_{\mu_{k}}(\mathbf{x}_{k})\|^{2}+\frac{1}{2}\|\mathbf{v}_{k}-\mathbf{v}^{*}_{\mu_{k}}(\mathbf{x}_{k})\|^{2},$
(10)

where
$\mathbf{v}^{*}_{\mu}(\mathbf{x})=\left[\nabla_{\mathbf{y}\mathbf{y}}^{2}\psi_{\mu}(\mathbf{x},\mathbf{y}^{*}_{\mu}(\mathbf{x}))\right]^{-1}\nabla_{\mathbf{y}}F(\mathbf{x},\mathbf{y}^{*}_{\mu}(\mathbf{x}))$
is the correct dual multipliers. The first term quantifies the descent of the
overall UL objective functions, the second term characterizes the LL solution
errors, and the third term characterizes the multipliers errors. We can prove
the following descent property of Lyapunov function $V_{k}$ for BAGDC.

######  Proposition 3.1.

Suppose Assumption 3.1 holds. If we choose aggregation parameter
$\mu_{k}\downarrow 0$ and step sizes

$\beta_{k}\leq\frac{2}{\sigma_{\psi_{\mu_{k}}}+L_{\psi_{\mu_{k}}}},\quad\eta_{k}\leq
c_{\eta}\mu_{k}^{4}\beta_{k},\quad\alpha_{k}\leq
c_{\alpha}\mu_{k}^{7}\eta_{k},$ (11)

then when $\mu_{k}$ is sufficiently small, the sequence of
$\mathbf{x}_{k},\mathbf{y}_{k},\mathbf{v}_{k}$ generated by BAGDC satisfies

$\displaystyle V_{k+1}-V_{k}\leq$
$\displaystyle-\frac{\alpha_{k}}{2}\|\nabla\varphi_{\mu_{k}}(\mathbf{x}_{k})\|^{2}-\frac{\widehat{\alpha}_{k}}{2}\|\mathbf{x}_{k+1}-\mathbf{x}_{k}\|^{2}-\frac{\widehat{\beta}_{k}}{2}\|\mathbf{y}_{k}-\mathbf{y}_{\mu_{k}}^{*}(\mathbf{x}_{k})\|^{2}-\frac{\widehat{\eta}_{k}}{2}\|\mathbf{v}_{k}-\mathbf{v}_{\mu_{k}}^{*}(\mathbf{x}_{k})\|^{2}$
$\displaystyle+C_{1}\left(\frac{\mu_{k}-\mu_{k+1}}{\mu_{k}}\right)+C_{2}\frac{1}{\alpha_{k}}\left(\mu_{k}-\mu_{k+1}\right)^{2}+C_{3}\frac{1}{\alpha_{k}}\left(\frac{\mu_{k}}{\mu_{k+1}}\right)^{2}\left(\mu_{k}-\mu_{k+1}\right)^{2},$

where $\widehat{\alpha}_{k}$, $\widehat{\beta}_{k}$, $\widehat{\eta}_{k}$,
$c_{\eta}$, $c_{\alpha}$ $C_{1}$, $C_{2}$ and $C_{3}$ are all positive
constants, whose formal formulae are stated in the Supplemental Material.

According to this descent property, we can notice that to guarantee that
$\|\nabla\varphi_{\mu_{k}}(\mathbf{x}_{k})\|$ converges to $0$, we need to
choose aggregation parameter $\mu_{k}$ and step size $\alpha_{k}$
appropriately such that

$\sum_{k=0}^{\infty}\alpha_{k}=\infty,\
~{}\sum_{k=0}^{\infty}\frac{1}{\alpha_{k}}\left(\mu_{k}-\mu_{k+1}\right)^{2}<\infty,\
\ \mathrm{and}\
~{}\sum_{k=0}^{\infty}\frac{1}{\alpha_{k}}\left(\frac{\mu_{k}}{\mu_{k+1}}\right)^{2}\left(\mu_{k}-\mu_{k+1}\right)^{2}<\infty.$

In following, we present a simple and handy strategy for choosing aggregation
parameter $\mu_{k}$ and step sizes $\alpha_{k}$, $\beta_{k}$ and $\eta_{k}$
for guaranteeing the convergence of our proposed BAGDC.

######  Theorem 3.1.

Suppose Assumption 3.1 holds. If we choose
$\mu_{k}=\bar{\mu}\left(\frac{1}{k+1}\right)^{p}$ with $0<p<\frac{1}{11}$ and

$\alpha_{k}=\bar{\alpha}\mu_{k}^{11},\quad\beta_{k}=\bar{\beta},\quad\eta_{k}=\bar{\eta}\mu_{k}^{4},$
(12)

where $\bar{\beta}\in(0,1/L_{\psi_{\mu_{k}}})$, then when all of positive
constants $\bar{\alpha}$, $\bar{\eta}$, $\bar{\mu}$ are sufficiently small,
the sequence of $\mathbf{x}_{k},\mathbf{y}_{k}$ generated by BAGDC satisfies

$\min_{0\leq k\leq
K}\left\\{\|\nabla\varphi_{\mu_{k}}(\mathbf{x}_{k})\|^{2}\right\\}=O\left(\frac{\ln
K}{K^{1-11p}}\right)\quad\mathrm{and}\quad\min_{0\leq k\leq
K}\left\\{\|\mathbf{y}_{k}-\mathbf{y}^{*}_{\mu_{k}}(\mathbf{x}_{k})\|^{2}\right\\}=O\left(\frac{\ln
K}{K^{1-p}}\right).$

To appropriately characterize the convergence rate of BAGDC, we define the KKT
residual for $(\mathbf{x},\mathbf{y},\mathbf{v})$ as following

$\mathrm{KKT}(\mathbf{x},\mathbf{y},\mathbf{v}):=\|\nabla\mathcal{L}(\mathbf{x},\mathbf{y},\mathbf{v})\|^{2}.$

Then, $\mathrm{KKT}(\mathbf{x},\mathbf{y},\mathbf{v})=0$ if and only if
$(\mathbf{x},\mathbf{y})$ is a KKT point to Problem (4) and $\mathbf{v}$ is a
multipliers. Under the linear independent constraint qualification (LICQ), we
can have following non-asymptotic convergence rate result of BAGDC from the
perspective of KKT residual.

######  Theorem 3.2.

Suppose Assumption 3.1 holds and aggregation parameter $\mu_{k}$ and step
sizes $\alpha_{k}$, $\beta_{k}$ and $\eta_{k}$ are chosen as same as in
Theorem 3.1. Let $(\mathbf{x}_{k},\mathbf{y}_{k},\mathbf{v}_{k})$ be the
sequence of generated by BAGDC with sufficiently small $\bar{\alpha}$,
$\bar{\eta}$, $\bar{\mu}$, if $(\mathbf{x}_{k},\mathbf{y}_{k})$ is bounded and
for any limit point $(\bar{\mathbf{x}},\bar{\mathbf{y}})$ of
$(\mathbf{x}_{k},\mathbf{y}_{k})$, LICQ of Problem (4) holds at
$(\bar{\mathbf{x}},\bar{\mathbf{y}})$, i.e., rows of
$\left(\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\bar{\mathbf{x}},\bar{\mathbf{y}})^{T},\nabla_{\mathbf{y}\mathbf{y}}^{2}f(\bar{\mathbf{x}},\bar{\mathbf{y}})^{T}\right)$
are linearly independent, then we have

$\min_{0\leq k\leq
K}\left\\{\mathrm{KKT}(\mathbf{x}_{k},\mathbf{y}_{k},\mathbf{v}_{k})\right\\}=O\left(\max\left\\{\frac{\ln
K}{K^{1-11p}},\frac{1}{K^{2p}}\right\\}\right).$

Furthermore, any limit point $(\bar{\mathbf{x}},\bar{\mathbf{y}})$ of
$\left(\mathbf{x}_{k},\mathbf{y}_{k}\right)$ is a KKT point of Problem (4).

###  3.2 Special Discussion for LL Strongly Convex Case

For the special case where the LL objective $f(\mathbf{x},\cdot)$ is strongly
convex with constant $\sigma_{f}>0$ for any fixed $\mathbf{x}$, the strong
convexity assumption on UL objective $F$ is not needed, and BAGDC can be shown
to have a $O(1/K)$ non-asymptotic convergence rate with taking $\mu_{k}=0$ and
constant step sizes $\alpha_{k}$, $\beta_{k}$ and $\eta_{k}$. We use the
following assumption to replace Assumption 3.1 in this section.

######  Assumption 3.2.

For any fixed $\mathbf{x}$, $f(\mathbf{x},\cdot)$ is $\sigma_{f}$-strongly
convex.

Note that when the LL objective $f(\mathbf{x},\cdot)$ is strongly convex,
$\mathbf{y}^{*}(\mathbf{x})$ is well defined and be differentiable. Hence, it
is proper to use the gradient of
$\varphi(\mathbf{x})=F(\mathbf{x},\mathbf{y}^{*}(\mathbf{x}))$ to characterize
the optimality of the generated sequence by the algorithms.

######  Theorem 3.3.

Suppose Assumption 3.2 holds. If we choose $\mu_{k}=0$ and

$\alpha_{k}=\bar{\alpha},\quad\beta_{k}=\bar{\beta},\quad\eta_{k}=\bar{\eta},$
(13)

where $\bar{\alpha},\bar{\beta}$, and $\bar{\eta}$ are positive constants
satisfying certain conditions 333A formal description of the conditions is
stated in the Supplemental Material., then $\mathbf{x}_{k}$ generated by BAGDC
satisfies

$\min_{0\leq k\leq
K}\\{\|\nabla\varphi(\mathbf{x}_{k})\|^{2}\\}=O\left(\frac{1}{K}\right).$

##  4 Experiments

In this section, we conduct extensive experiments on both synthetic data and
real applications to verify our theoretical findings and demonstrate the
efficacy of BAGDC, especially on very high-dimensional tasks. The experiments
were done on a PC with Intel i7-9700K CPU (4.2 GHz), 32GB RAM, and an NVIDIA
RTX 2060S GPU. The algorithm was implemented using PyTorch 1.6 with CUDA 10.2.

###  4.1 Numerical Verification

In this subsection, we start by discussing the effect of dual correction under
different settings. We then compare the performance of BAGDC with other GBLOs
and I-GBLOs by using different accuracies for LL problem. The third section
demonstrates the efficacy of BAGDC compared to RHG, CG and NS. Finally, we
compare the empirical performance of BAGDC for BLOs with multiple LL
minimizers. We conduct experiments in Eq. (3) which satisfies the assumptions
in Theorem 3.3. In particular we set $\mathbf{A}=\mathbf{I}$,
$\mathbf{z}_{0}=\mathbf{e}$, where $\mathbf{e}$ is the vector whose elements
are all equal to 1. The unique solution is
$\mathbf{x}^{*}=\mathbf{y}^{*}=\mathbf{e}/2$. For the existing method, we set
the number of lower level iterations to 100.
$\mathbf{d}=\mathbf{d}(\mathbf{x},\mathbf{y},\mathbf{v})$ is the direction of
descent of $\mathbf{x}$ calculated by different methods, e.g. for NOSA,
$\mathbf{d}=\left(\nabla_{\mathbf{x}}F(\mathbf{x},\mathbf{y})-\left[\nabla_{\mathbf{x}\mathbf{y}}^{2}f(\mathbf{x},\mathbf{y})\right]\nabla_{\mathbf{y}}F(\mathbf{x},\mathbf{y})\right)$.
We define the gradient of $\mathbf{x}$ as $\nabla\varphi$ where
$\varphi(\mathbf{x})=F(\mathbf{x},\mathbf{y}^{*}(\mathbf{x}))$.

Effect of Dual Correction with Different Settings. Although we proved in
Section 3 that the addition of dual variable $\mathbf{v}$ guarantees
convergence, the update step size $\eta$ of $\mathbf{v}$ still depends on
properties of $F$ and $f$ (e.g. the Lipschitz constant and the strongly convex
coefficients). Since this is usually agnostic in practical applications, we
first verify the corrective effect of $\mathbf{v}$ for convergence for
different $\eta$ settings. In addition to different fixed $\eta$, we also
validate a strategy for the Adaptive (Adapt.) computation of $\eta$, a CG-
inspired strategy.444More details is stated in the Supplemental Material.
Figure 1 (a) shows the convergence of NOSA and BAGDC with different $\eta$. It
can be seen that all $\eta$ improve the convergence performance of non-
convergence NOSA. For a fixed $\eta$, any $\eta$ within the theoretical bound
is guaranteed to converge, although smaller $\eta$ lead to a lower speed of
convergence. On the other hand, $\eta$ outside the theoretical bound may lead
to a collapse after a while. Note that despite the convergence result of
Theorem 3.3, BAGDC is not guaranteed to be monotonically decreasing, and in
fact, some degree of oscillation is expected for some $\eta$. $\eta$ with
adaptive strategy achieved a convergence speed second only to $\eta=C_{f}$.
Since adaptive $\eta$ does not need to know the nature of $F$ and $f$, we used
adaptive $\eta$ as our strategy in all subsequent experiments. However,
considering that adaptive method may require more computational cost in
computing $\eta$, we further show in Figure 1 (b) the Average (Avg.) step time
and the overall time to converge to $\|\mathbf{d}\|\leq 10^{-5}$ for different
$\eta$ settings. As can be seen, although the adaptive method takes slightly
more time for a single step, the overall convergence time is still second only
to $\eta=C_{f}$.

|
---|---
|
(a) Convergence curves in different $\eta$ | (b) Overall/Avg. Step
Figure 1: Convergence comparison of update strategies for $\mathbf{v}$ in
different $\eta$. $C_{f}$ is the theoretical upper bound. (a) Convergence
curve with different $\eta$. (b) Overall time for convergence and Average
(Avg.) step time. The theoretical upper bound has the best convergence speed,
but the sub-optimal adaptive $\eta$ can be applied to real applications where
the theoretical upper bound is difficult to compute.

Solving LL with Different Accuracies vs BAGDC. We next show that BAGDC
effectively decouples the dependence on the accuracy of the solution to the LL
problem when applied to different BLO methods. We use the inner iterations $T$
for RHG, error $\epsilon$ for CG and sequence lengths $M$ for NS as indicators
of the accuracy of LL problem. In Figure 2, we compare the convergence of
existing methods for different LL problem solving accuracies and existing
methods with our BAGDC. We show the error of direction of descent
$\|\mathbf{d}-\nabla F\|$ to illustrate that an incorrect estimate of
$\mathbf{v}$ for the existing methods leads to the wrong direction of descent,
and $\|\mathbf{x}-\mathbf{x}^{*}\|/\|\mathbf{x}^{*}\|$ to illustrate the
effect of this on the final convergence result. It can be seen that
accelerating existing method by directly simplifying the LL iterations will
reduce the quality of the solution. By introducing $\mathbf{v}$ as a
correction in BAGDC, our method can continuously improve accuracy for these
BLOs with alternating solving UL and LL problem.

|  |
---|---|---
|  |
(a) RHG | (b) CG | (c) NS
Figure 2: The convergence curves of our BAGDC and various methods (RHG, CG and
NS) in different accuracy of LL problem, i.e., iterations $T$ for RHG, error
$\epsilon$ for CG, sequence lengths $M$ for NS. Our BAGDC method with only one
step for LL problem achieves the best convergence results.

Significant Speed Advantage of BAGDC. We then show the improvement in
computational efficiency of BAGDC, the extremely simpler implementation of
BLO. For existing methods, the main part of the computational complexity is
the large number of matrix-vector products in nested iteration. Therefore by
changing the dimension of the matrix $\mathbf{A}$, we can easily extend Eq (3)
to the high dimensional case to observe the computational efficiency. Here we
show the convergence in different dimensions on Eq. (3) with limited total
time $10^{4}$ s. Figure 3 (a) shows the convergence time of different
indicator $\|\mathbf{d}\|\leq 10^{-3}$, $\|F-F^{*}\|/\|F^{*}\|\leq 10^{-3}$
and $\|\mathbf{x}-\mathbf{x}^{*}\|/\|\mathbf{x}^{*}\|\leq 10^{-4}$. Since each
step usually contain only one matrix vector product, the alternating
optimization BAGDC saves a lot of time compared to the multi-step per
iteration of nested optimization. Thus our method requires much less time
compared to other methods of all scales. Figure 3(b) illustrates the limiting
problem dimension that existing methods can achieve when reaching the time
limit. BAGDC does not contain repeated nested matrix-vector products allowing
our method to be solved on problems of higher dimensionality.

|
---|---
(a) Convergence time |  (b) Convergence dimension
Figure 3: Computational efficiency comparison of existing BLOs (RHG, CG and
NS) and ours BAGDC. (a) The time required for different convergence metrics in
different dimensions. BAGDC reduces the time required compared to existing
methods. (b) The problem dimension when reached the time limit. BAGDC greatly
increases the dimension of the problem to reach the time limit.

Solving BLOs with multiple LL minimizers. We next verify that BAGDC can also
solve the BLOs with multiple LL minimizers as mentioned in Theorem 3.1. We use
the example without LLS given in Supplemental Material. Figure 4 compares the
convergence curves with BLO methods. As we can see, only BAGDC and BDA are
able to converge to the optimum point, the other methods only converge to an
incorrect stability point. Thanks to the much simpler iteration format
resulting from alternating optimisation, BAGDC converge much faster than BDA.

---
Figure 4: Convergence curves of existing methods under BLO beyond LLS. BAGDC
has a very high convergence speed with guaranteed correct convergence because
of the much simpler iteration.

###  4.2 Real-world Applications

In this section, we verify the performance of our BAGDC on real data and high-
dimensional problems. We show two widely used applications of BLOs, data
hyper-cleaning, and few-shot learning.

Figure 5: Comparison of the validation loss $F$ and accuracy for hyper-
cleaning. BAGDC dramatically increases the training speed.

Data Hyper-cleaning. Data hyper-cleaning is a widely used application for bi-
level optimization. Assuming that some of the labels in our dataset are
contaminated, data hyper-cleaning aims to reduce the impact of incorrect
samples by adding hyper-parameters to them. In the experiment, we follow the
settings in [13] and based on FashionMNIST datasets. To demonstrate the
advantage of our method in terms of computational efficiency, we show in the
left part of Table 1 the F1 scores and the time required for different methods
to achieve the same accuracy. BAGDC achieves the same accuracy and similar F1
scores while reducing the time requirement to an order of magnitude. Figure 5
shows the validation loss for the different methods. It can be seen that our
method has the fastest convergence rate and can still be the fastest at higher
accuracies after quickly achieving $81\%$ accuracy rate.

Few-shot Learning. We then test the application in higher dimensions to verify
the computational efficiency of our method in few-shot learning. Few-shot
learning (i.e., N-way M-shot) is a multi-task N-way classification and it aims
to learn the hyperparameter such that each task can be solved only with M
training samples. In the experiment, we follow the settings in [13] and based
on DoubleMNIST datasets [28]. In the right part of Table 1 , we shows the time
that existing methods achieve the same accuracy (95% for 5-way and 85% for-20
way). Our method achieves the same accuracy rate while significantly reducing
the time required. Therefore, our alternating optimization method can be
effectively applied to large-scale neural networks.

Table 1: Comparison of the results of existing methods for hyper-cleaning and few-shot learning. We compared the time and F1 score for hyper-cleaning when achieving almost the same accuracy ($81\%$) and the time for few-shot learning ($95\%$ for 5-way and $85\%$ for 20-way). BAGDC has the highest computational efficiency while training to the same accuracy rate. Method | Hyper-cleaning |  |  | 5-way | 20-way
---|---|---|---|---|---
Acc. | F1 score | Time |  |  | Acc. | Time | Acc. | Time
RHG | 81.07 | 87.10 | 41.47 |  |  | 95.05 | 1614.38 | 85.06 | 2645.27
BDA | 81.09 | 88.13 | 64.84 |  |  | 95.07 | 1864.90 | 85.09 | 3017.33
CG | 81.05 | 87.85 | 19.72 |  |  | 95.08 | 1056.91 | 85.00 | 2060.94
NS | 81.03 | 85.53 | 26.01 |  |  | 95.07 | 1005.26 | 85.05 | 1942.02
BAGDC | 81.00 | 87.06 | 2.92 |  |  | 95.06 | 98.16 | 85.11 | 381.21

##  5 Conclusion

This work first designed a counter-example to demonstrate that there actually
exist fundamental theoretical issues in these naive single-inner-step
acceleration for existing GBLOs. Then by reformulating BLO as a constrained
single-level optimization problem and introducing dual multipliers from the
viewpoint of KKT condition, an extremely fast gradient method, named BAGDC, is
proposed for BLOs. Our theoretical results guaranteed that BAGDC can
successfully obtain non-asymptotic convergence sequences towards self-governed
KKT stationary points. We also extended this result for the more challenging
BLO problems, i.e., with multiple LL solutions.

## References

  * Liu et al. [2021] Risheng Liu, Jiaxin Gao, Jin Zhang, Deyu Meng, and Zhouchen Lin.  Investigating bi-level optimization for learning and vision from a unified perspective: A survey and beyond.  _IEEE TPAMI_ , 2021.
  * Franceschi et al. [2017] Luca Franceschi, Michele Donini, Paolo Frasconi, and Massimiliano Pontil.  Forward and reverse gradient-based hyperparameter optimization.  In _ICML_ , 2017.
  * Okuno et al. [2018] Takayuki Okuno, Akiko Takeda, and Akihiro Kawana.  Hyperparameter learning via bilevel nonsmooth optimization.  _arXiv preprint arXiv:1806.01520_ , 2018.
  * Mackay et al. [2019] Matthew Mackay, Paul Vicol, Jonathan Lorraine, David Duvenaud, and Roger Grosse.  Self-tuning networks: Bilevel optimization of hyperparameters using structured best-response functions.  In _ICLR_ , 2019.
  * Franceschi et al. [2018] Luca Franceschi, Paolo Frasconi, Saverio Salzo, Riccardo Grazzi, and Massimiliano Pontil.  Bilevel programming for hyperparameter optimization and meta-learning.  In _ICML_ , 2018.
  * Rajeswaran et al. [2019] Aravind Rajeswaran, Chelsea Finn, Sham M Kakade, and Sergey Levine.  Meta-learning with implicit gradients.  In _NeurIPS_ , 2019.
  * Zügner and Günnemann [2018] Daniel Zügner and Stephan Günnemann.  Adversarial attacks on graph neural networks via meta learning.  In _ICLR_ , 2018.
  * Liu et al. [2018] Hanxiao Liu, Karen Simonyan, and Yiming Yang.  Darts: Differentiable architecture search.  In _ICLR_ , 2018.
  * Liang et al. [2019] Hanwen Liang, Shifeng Zhang, Jiacheng Sun, Xingqiu He, Weiran Huang, Kechen Zhuang, and Zhenguo Li.  Darts+: Improved differentiable architecture search with early stopping.  _arXiv preprint arXiv:1909.06035_ , 2019.
  * Chen et al. [2019] Xin Chen, Lingxi Xie, Jun Wu, and Qi Tian.  Progressive differentiable architecture search: Bridging the depth gap between search and evaluation.  In _ICCV_ , 2019.
  * Pfau and Vinyals [2016] David Pfau and Oriol Vinyals.  Connecting generative adversarial networks and actor-critic methods.  _arXiv preprint arXiv:1610.01945_ , 2016.
  * Yang et al. [2019] Zhuoran Yang, Yongxin Chen, Mingyi Hong, and Zhaoran Wang.  Provably global convergence of actor-critic: A case for linear quadratic regulator with ergodic cost.  In _NeurIPS_ , 2019.
  * Liu et al. [2020] Risheng Liu, Pan Mu, Xiaoming Yuan, Shangzhi Zeng, and Jin Zhang.  A generic first-order algorithmic framework for bi-level programming beyond lower-level singleton.  In _ICML_ , 2020.
  * Grazzi et al. [2020] Riccardo Grazzi, Luca Franceschi, Massimiliano Pontil, and Saverio Salzo.  On the iteration complexity of hypergradient computation.  In _ICML_ , 2020.
  * Ji et al. [2021] Kaiyi Ji, Junjie Yang, and Yingbin Liang.  Bilevel optimization: Convergence analysis and enhanced design.  In _ICML_ , 2021.
  * Elsken et al. [2020] Thomas Elsken, Benedikt Staffler, Jan Hendrik Metzen, and Frank Hutter.  Meta-learning of neural architectures for few-shot learning.  In _CVPR_ , 2020.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  _NeurIPS_ , 2014.
  * Pedregosa [2016] Fabian Pedregosa.  Hyperparameter optimization with approximate gradient.  In _ICML_ , 2016.
  * Lorraine et al. [2020] Jonathan Lorraine, Paul Vicol, and David Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  In _AISTATS_ , 2020.
  * Cao et al. [2021] Xiangyong Cao, Xueyang Fu, Danfeng Hong, Zongben Xu, and Deyu Meng.  Pancsc-net: A model-driven deep unfolding method for pansharpening.  _IEEE Transactions on Geoscience and Remote Sensing_ , 2021.
  * Finn et al. [2017] Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _ICML_ , 2017.
  * Jin et al. [2021] Dian Jin, Long Ma, Risheng Liu, and Xin Fan.  Bridging the gap between low-light scenes: Bilevel learning for fast adaptation.  In _Proceedings of the 29th ACM International Conference on Multimedia_ , pages 2401–2409, 2021.
  * Chen et al. [2022] Tianyi Chen, Yuejiao Sun, Quan Xiao, and Wotao Yin.  A single-timescale method for stochastic bilevel optimization.  In _AISTATS_ , 2022.
  * Khanduri et al. [2021] Prashant Khanduri, Siliang Zeng, Mingyi Hong, Hoi-To Wai, Zhaoran Wang, and Zhuoran Yang.  A near-optimal algorithm for stochastic bilevel optimization via double-momentum.  _NeurIPS_ , 2021.
  * Yang et al. [2021] Junjie Yang, Kaiyi Ji, and Yingbin Liang.  Provably faster algorithms for bilevel optimization.  _NeurIPS_ , 2021.
  * Ghadimi and Wang [2018] Saeed Ghadimi and Mengdi Wang.  Approximation methods for bilevel programming.  _arXiv preprint arXiv:1802.02246_ , 2018.
  * Wright et al. [1999] Stephen Wright, Jorge Nocedal, et al.  Numerical optimization.  _Springer Science_ , 35(67-68):7, 1999.
  * Sitzmann et al. [2020] Vincent Sitzmann, Eric Chan, Richard Tucker, Noah Snavely, and Gordon Wetzstein.  Metasdf: Meta-learning signed distance functions.  _NeurIPS_ , 2020.