# Effective Hamiltonians in the Quantum Rabi Problem

P. Gartner1 V. Moldoveanu1 1 National Institute of Materials Physics, PO Box
MG-7, Bucharest-Magurele, Romania

(August 27, 2024)

###### Abstract

We revisit the theoretical description of the ultrastrong light-matter
interaction in terms of exactly solvable effective Hamiltonians. A
perturbative approach based on polaronic and spin-dependent squeezing
transformations provides an effective Hamiltonian for the quantum Rabi model
up to the second order in the expansion parameter. The model consistently
includes both rotating and counter-rotating terms, going therefore beyond the
rotating wave approximation. Analytical and numerical results show that the
proposed Hamiltonian performs better than the Bloch-Siegert model when
calculating operator averages (e.g. the mean photon number and number of
excitations). This improvement is due to a refined calculation of the dressed
states within the present model. Regarding the frequency shift induced by the
qubit-photon interaction, we find a different sign from the Bloch-Siegert
value. This influences the eigenstates structure in a non-trivial way and
ensures the correct calculation of the number of excitations associated to a
given dressed state. As a consistency check, we show that the exactly solvable
independent boson model is reproduced as a special limit case of the
perturbative Hamiltonian.

##  I Introduction

The physics of strongly coupled qubit-photon hybrid systems has become a
steadily developing field Kockum ; RMP  which challenges the existing
theoretical methods. In particular, as the field-matter coupling strength
increases, the Jaynes-Cummings model JCM , which neglects the so-called off-
resonant transitions within the rotating wave approximation (RWA), becomes
less accurate.

More general models taking fully into acount the interaction between a two-
level system (e.g. a superconducting qubit or a quantum dot) and a boson mode
(cavity photon or transmission line resonator) are defined by the class of
Hamiltonians

$H=\frac{1}{2}\varepsilon\sigma_{z}+\omega
a^{\dagger}a+g_{1}(a^{\dagger}\sigma+a\sigma^{\dagger})+g_{2}(a^{\dagger}\sigma^{\dagger}+a\sigma)\,,$
(1)

where the quasi-spin notation is used for the qubit and $a,a^{\dagger}$ are
bosonic operators of the single-mode field. The corresponding energies are
$\varepsilon$ and $\omega$ respectively, and $g_{1}$ and $g_{2}$ are the
coupling strengths for the resonant (rotating) and off-resonant (counter-
rotating) processes. We take $\hbar=1$ throughout. The Jaynes-Cummings (JC)
Hamiltonian is then recovered from Eq.(1) by taking $g_{2}=0$, while the
quantum Rabi model (QRM) Meystre ; Walls-Milburn  corresponds to
$g_{1}=g_{2}=g$:

$H_{R}=\frac{1}{2}\varepsilon\sigma_{z}+\omega
a^{\dagger}a+g\sigma_{x}(a^{\dagger}+a)\,.$ (2)

Generally, the problems defined by Eq.(1) are not exactly soluble, one notable
exception being the quantum Rabi Hamiltonian whose spectrum was calculated
Braak1 ; Braak2 ; Wang  as zeroes of recursively defined new special functions
or in terms of confluent Heun functions Zhong . The mathematics is quite
involved both for the eigenvalues and for eigenvectors. In view of this
situation, in the literature one prefers to devise effective Hamiltonians RMP
to obtain approximate eigenvalues and eigenvectors, in suitably restricted
domains of the parameter space defined by $\\{\varepsilon,\omega,g\\}$. Also,
simple, analytic effective descriptions are strongly preferred when including
additional interactions to describe optical, transport Cirio ; DMG  or
relaxation processes DeLiberato ; Beaudoin ; Settineri .

A useful test for these approximations consists in inspecting limit cases in
which exact solutions are known for the full Hamiltonian. One such case is the
JC model mentioned above, in which the RWA brings the Hamiltonian to a
$2\times 2$ block-diagonal form. A second one is the independent boson model
(IBM) Mahan  which is extensively used for quantum optical and transport
calculations involving carrier interaction with photons, phonons Krummheuer
or vibrons Mitra . In the present case the IBM limit is obtained by taking
$\varepsilon=0$ in the QRM, Eq. (2), and the corresponding Hamiltonian is
diagonalized by a displacement of the bosonic oscillator centers brought in by
the Lang-Firsov, or polaronic unitary transform.

In deriving approximate effective Hamiltonians for the quantum Rabi problem
one usually employs a limited expansion of various unitary transforms in terms
of a suitable small parameter. An example is provided by the dispersive
Hamiltonian Zueco ; Blais  obtained by a unitary transform aiming at
diagonalizing the qubit-boson interaction, expanded to second order in the
parameter $\lambda=g/|\varepsilon-\omega|$, in the weak coupling limit
$\lambda\ll 1$. The resulting effective Hamiltonian reproduces the JC spectrum
to second order in $\lambda$, i.e. within the domain of parameters in which it
was derived. Moreover, it recovers the exact result in the IBM limit.

Another popular effective Hamiltonian is the Bloch-Siegert (BS) one Klimov-
Chumakov ; Rossatto . It treats the QRM problem up to second order in the
parameter $\lambda=g/(\varepsilon+\omega)$ and brings it to a renormalized
form of the JC Hamiltonian. Two succesive unitary transforms are used in the
process and higher order off-resonant terms are discarded as rapidly
oscillating, in the spirit of the RWA. Supposedly, the validity domain of this
approach is roughly $g/\omega\lesssim 0.3$ and corresponds to the so-called
perturbative ultra strong coupling regime (pUSC) RMP . The perturbative
criteria for the BS Hamiltonian was recently used Rossatto  to establish a
spectral classification of the field-matter coupling regimes in the QRM. On
the experimental side the Bloch-Siegert shift was recently observed in
strongly coupled qubit-oscillator systems Diaz-BS . An alternative improvement
of the JC model relies on the Van Vleck perturbation theory to the second
order in the qubit-field coupling Hausinger .

Surprisingly, a simple analysis shows that the IBM limit is not recovered by
the BS Hamiltonian, despite the fact that taking $\varepsilon=0$ does not
preclude the expansion parameter, now equal to $g/\omega$, to remain small.
Still, the BS model is tacitly accepted as the main candidate to describe the
ultra-strong coupling regime.

In view of this situation, in what follows we show how to obtain an effective
Hamiltonian for the QRM, in the same perturbative conditions as for the BS
model, i.e. up to the second order in $\lambda$, having the same simple,
renormalized JC structure and which additionally meets the test of recovering
the exact IBM result for $\varepsilon=0$. Moreover, it will become obvious
that the procedure is purely perturbative, i.e. it does not have to rely on
discarding off-resonant terms.

As the qubit energy $\varepsilon$ approaches resonance the spectra of the two
effective models become more similar. Nevertheless, we show that important
differences remain in the dressed states and are inherited by the expectation
values provided by these states. Especially the excitation number turns out to
be sensitive to these differences.

The rest of the paper is organized as follows. The perturbative effective
Hamiltonian is derived in Section II. A detailed comparison of this
perturbative Hamiltonian with the Bloch-Siegert model is presented and
discussed in Section III. Their spectra, average photon number, population
inversion and number of excitations are obtained and compared both
analytically and numerically in subsections III A, III B, and III C. We also
compare our result with the so called generalized rotating wave Hamiltonian
(GRWA) obtained by Irish Irish  in the deep strong coupling case (subsection
III D). This regime is defined by large values of the coupling,
$g/\omega\approx 1$ and beyond, but the numerical evidence Yu ; Xie ;
Zhang2016  shows that this approach gives good results at lower couplings too.
Therefore comparison against it is relevant. Besides, the GRWA Hamiltonian
meets the IBM limit test. Section IV is left for conclusions.

##  II An IBM-compatible effective Hamiltonian

We consider unitary transforms of the form $e^{S}$, with $S^{\dagger}=-S$ and
their expansion in iterated commutators with their generators $S$. Then any
operator $A$ transforms as

$A_{S}=e^{S}A\,e^{-S}=A+\left[S,A\right]+\frac{1}{2!}\left[S,\left[S,A\right]\right]+\dots$
(3)

Now let us rewrite the Rabi Hamiltonian of Eq.(2) as

$H_{R}=\frac{1}{2}\varepsilon\sigma_{z}+H_{IBM}\,,$ (4)

in which we separate the free qubit contribution and the independent boson
model part. The linear part in $H_{IBM}$ can be absorbed in the quadratic term
by a shift of the bosonic operators, produced by the generator

$S=\lambda(a^{\dagger}-a)\,\sigma_{x}\,$ (5)

with the value of $\lambda$ to be chosen later. Using

$[S,a^{\dagger}]=-\lambda\,\sigma_{x}\,,\qquad[S,a]=-\lambda\,\sigma_{x}\,,$
(6)

it is clear that the canonical transformation amounts to a simple shift

$a^{\dagger}_{S}=a^{\dagger}-\lambda\,\sigma_{x}\,,\qquad
a_{S}=a-\lambda\,\sigma_{x}\,.$ (7)

As a result, the transformed Hamiltonian becomes

$H_{IBM,S}=\omega
a^{\dagger}a+(g-\lambda\omega)(a^{\dagger}+a)\sigma_{x}+(\lambda^{2}\omega-2\lambda
g)\,.$ (8)

Choosing $\lambda=g/\omega$ removes the interaction term and diagonalizes the
IBM Hamiltonian.

$H_{IBM,S}=\omega a^{\dagger}a-g^{2}/\omega\,.$ (9)

Therefore the IBM problem alone is exactly solvable. The solution consists of
the oscillator keeping its frequency $\omega$, but shifted by $g/\omega$, the
shift direction being decided by the eigenstate of $\sigma_{x}$

The presence of the free qubit part in Eq.(4) spoils this simple scenario.
Both the frequency and the shift are modified, and we show below that in these
circumstances a different choice for the value of $\lambda$ is more
appropriate.

Obviously, $\sigma_{x}$ is left unchanged by the IBM transform. In contrast,
the change in $\sigma_{z}$ is rather complex. One can easily check that:

$\displaystyle[S,\sigma_{z}]=$
$\displaystyle-2\,\lambda(a^{\dagger}-a)\,i\sigma_{y}\,,$
$\displaystyle[S,i\sigma_{y}]=$
$\displaystyle-2\,\lambda(a^{\dagger}-a)\,\sigma_{z}\,.$ (10)

This means that successive commutators generate alternatively the $\sigma_{z}$
and $i\sigma_{y}$ operators, and each step brings in an additional factor
$-2\,\lambda(a^{\dagger}-a)$. The $\sigma_{z}$ operator collects the even
terms, which add up to a hyperbolic cosine function, while the $i\sigma_{y}$
terms generate a hyperbolic sine function. The net result is

$\displaystyle H_{R,S}=\omega
a^{\dagger}a+(g-\lambda\omega)(a^{\dagger}+a)\sigma_{x}+(\lambda^{2}\omega-2\lambda
g)$
$\displaystyle+\,\frac{1}{2}\varepsilon\left\\{\cosh[2\,\lambda(a^{\dagger}-a)]\,\sigma_{z}-\sinh[2\,\lambda(a^{\dagger}-a)]\,i\sigma_{y}\right\\}\,.$
(11)

Up to now no approximation is involved. A simpler Hamiltonian, valid to the
second order in $\lambda$, is obtained by expanding $\cosh(x)\approx
1+x^{2}/2$ and $\sinh(x)\approx x$. Taking into account that
$\sigma_{x}=\sigma^{\dagger}+\sigma$ and
$i\sigma_{y}=\sigma^{\dagger}-\sigma$, this expansion leads to an effective
Hamiltonian $H^{\prime}$ of the form

$\displaystyle H_{R,S}\approx
H^{\prime}=\frac{1}{2}\varepsilon\sigma_{z}+\omega
a^{\dagger}a+(\lambda^{2}\omega-2\lambda g)$
$\displaystyle+(g-\lambda\omega+\lambda\varepsilon)(a^{\dagger}\sigma+a\sigma^{\dagger})+(g-\lambda\omega-\lambda\varepsilon)(a^{\dagger}\sigma^{\dagger}+a\sigma)$
$\displaystyle+\lambda^{2}\varepsilon(a^{\dagger}-a)^{2}\sigma_{z}\,.$ (12)

Now it is clear that an appropriate choice for $\lambda$ is

$\lambda=\frac{g}{\varepsilon+\omega}\,$ (13)

which cancels out the off-resonant term linear in the bosonic operators. This
is the same small parameter as the one used in the Bloch-Siegert formalism
Klimov-Chumakov . The effective Hamiltonian has now the expression

$\displaystyle H^{\prime}$
$\displaystyle=\frac{1}{2}\varepsilon\sigma_{z}+\omega
a^{\dagger}a-(2\,\lambda g-\lambda^{2}\omega)$
$\displaystyle+2\lambda\varepsilon(a^{\dagger}\sigma+a\sigma^{\dagger})+\lambda^{2}\varepsilon(a^{\dagger}-a)^{2}\sigma_{z}\,.$
(14)

Finally, one has to handle the last, quadratic term. To this end one considers
a second unitary transform, with the generator

$S^{\prime}=\frac{1}{2}\eta(a^{\dagger\,2}-a^{2})\sigma_{z}\,.$ (15)

This $\sigma_{z}$-dependent squeezing transform Comm-Zhang  leads to

$[S^{\prime},a^{\dagger}]=-\eta\,a\,\sigma_{z}\,,\qquad[S^{\prime},a]=-\eta\,a^{\dagger}\sigma_{z}\,,$
(16)

in which again succesive commutators generate alternatively the same two
operators. One gets a Bogolyubov-type transform

$\displaystyle a^{\dagger}_{S^{\prime}}$
$\displaystyle=\cosh\eta\,a^{\dagger}-\sinh\eta\,a\,\sigma_{z}\,,$
$\displaystyle a_{S^{\prime}}$
$\displaystyle=-\sinh\eta\,a^{\dagger}\,\sigma_{z}+\cosh\eta\,a\,.$ (17)

Different Hamiltonian pieces change accordingly, for instance the photon
number operator

$\displaystyle a_{S^{\prime}}^{\dagger}a_{S^{\prime}}=\cosh
2\eta\,a^{\dagger}a+\sinh^{2}\eta$ $\displaystyle-\frac{1}{2}\sinh
2\eta(a^{\dagger}-a)^{2}\sigma_{z}-\frac{1}{2}\sinh
2\eta(2a^{\dagger}a+1)\sigma_{z}\,,$ (18)

and

$\displaystyle(a_{S^{\prime}}^{\dagger}-a_{S^{\prime}})$
$\displaystyle=(\cosh\eta+\sinh\eta\,\sigma_{z})(a^{\dagger}-a)\,,$
$\displaystyle(a_{S^{\prime}}^{\dagger}-a_{S^{\prime}})^{2}\sigma_{z}$
$\displaystyle=[\cosh 2\eta\,\sigma_{z}+\sinh 2\eta](a^{\dagger}-a)^{2}\,.$
(19)

From Eqs.(18),(19) above the quantities in $H^{\prime}_{S^{\prime}}$
containing $(a^{\dagger}-a)^{2}$ are identified as:

$\displaystyle X_{z}$ $\displaystyle:=$
$\displaystyle\sigma_{z}\left(\lambda^{2}\varepsilon\cosh
2\eta-\frac{\omega}{2}\sinh 2\eta\right)(a^{\dagger}-a)^{2}\,,$ (20)
$\displaystyle X_{0}$ $\displaystyle:=$
$\displaystyle\lambda^{2}\varepsilon\sinh 2\eta(a^{\dagger}-a)^{2}\,.$ (21)

The terms contained in $X_{z}$ cancel out by choosing $\eta$ so that

$\frac{1}{2}\sinh 2\eta=\frac{\lambda^{2}\varepsilon}{\omega}\cosh 2\eta\,.$
(22)

The condition makes $\eta$ of second order in $\lambda$ and, neglecting higher
contributions, one has

$\eta=\frac{\lambda^{2}\varepsilon}{\omega}\,.$ (23)

In this case the term $X_{0}$ in Eq.(21) is of fourth order in $\lambda$ and
is discarded.

The terms in the Hamiltonian Eq.(14) not considered up to now are the free
qubit energy, which commutes with $S^{\prime}$, and the JC interaction
$2\lambda\,\varepsilon(a^{\dagger}\sigma+a\sigma^{\dagger})$, which is already
a first order quantity, and therefore its transformation by
$S^{\prime}\sim\lambda^{2}$ brings in changes of at least third order.
Therefore both are kept unchanged.

Finally, expanding systematically the hyperbolic functions to second order in
$\lambda$ (i.e. first order in $\eta$) and spelling out the expressions of
$\lambda$ and $\eta$ in terms of the model parameters, leaves us with
$H^{\prime}_{S^{\prime}}\approx H_{P}$ where the perturbative IBM-compatible
effective Hamiltonian is given by

$\displaystyle H_{P}=$
$\displaystyle\frac{\varepsilon}{2}\left(1-\frac{2g^{2}}{(\varepsilon+\omega)^{2}}\right)\sigma_{z}+\left(\omega-\frac{2g^{2}\varepsilon}{(\varepsilon+\omega)^{2}}\,\sigma_{z}\right)\,a^{\dagger}a$
$\displaystyle+2\,\frac{g\,\varepsilon}{\varepsilon+\omega}(a^{\dagger}\sigma+a\,\sigma^{\dagger})-\frac{g^{2}}{\varepsilon+\omega}\,\frac{2\,\varepsilon+\omega}{\varepsilon+\omega}\,.$
(24)

This expression of the effective Hamiltonian $H_{P}$ is the main result of the
paper.

It is now obvious that for $\varepsilon=0$, $H_{P}$ recovers exactly the IBM
result, Eq. (9) and, moreover, this holds for all values of the coupling
constant $g$.

##  III Comparing effective Hamiltonians

The Bloch-Siegert and the perturbative effective Hamiltonian given by Eq. (24)
have both the structure of a renormalized JC model

$H=\frac{1}{2}(\varepsilon+\Omega)\sigma_{z}+(\omega+\Omega\sigma_{z})a^{\dagger}a+\gamma(a^{\dagger}\sigma+a\sigma^{\dagger})-\delta\,.$
(25)

The difference comes from the parameters $\Omega,\gamma$ and $\delta$. They
are given by Rossatto

$\Omega=\frac{g^{2}}{\omega+\varepsilon},\quad\gamma=g,\quad\delta=\frac{1}{2}\,\frac{g^{2}}{\omega+\varepsilon}\,$
(26)

for $H=H_{BS}$, and by

$\Omega=-\frac{g^{2}}{\omega+\varepsilon}\,\frac{2\varepsilon}{\omega+\varepsilon},\quad\gamma=g\,\frac{2\varepsilon}{\omega+\varepsilon},\quad\delta=\frac{g^{2}}{\omega+\varepsilon}\,\frac{\omega+2\varepsilon}{\omega+\varepsilon}\,$
(27)

for $H=H_{P}$.

The JC structure of Eq.(25) reduced the problem to $2\times 2$ blocks, indexed
by the photon number $n$. In the corresponding basis
$\left\\{\ket{n,\uparrow},\ket{n+1,\downarrow}\right\\}$ the $n$-th block has
the form

$\displaystyle{\hat{h}}_{n}$ $\displaystyle=$
$\displaystyle\alpha_{n}\sigma_{0}+\frac{1}{2}\Delta_{n}\sigma_{z}+\gamma\sqrt{n+1}\sigma_{x}$
(28) $\displaystyle=$
$\displaystyle\alpha_{n}\sigma_{0}+R_{n}(\cos\theta_{n}\sigma_{z}+\sin\theta_{n}\sigma_{x})\,,$

where $\sigma_{0}$ is the $2\times 2$ unit matrix and we introduced

$\displaystyle\alpha_{n}$ $\displaystyle=$
$\displaystyle\omega\,(n+1/2)-\left(\frac{\Omega}{2}+\delta\right)$ (29)
$\displaystyle R_{n}$ $\displaystyle=$
$\displaystyle\sqrt{\Delta_{n}^{2}/4+\gamma^{2}(n+1)}\,,$ (30)

with $\Delta_{n}=\varepsilon-\omega+2\,\Omega\,(n+1)$ acting as an effective,
$n$-dependent detuning. The angle $\theta_{n}$ is defined by the relations
$\cos\theta_{n}=\Delta_{n}/(2R_{n})$ and
$\sin\theta_{n}=\gamma\sqrt{n+1}/R_{n}$. Note that $\Omega/2+\delta$ turns out
to take the same value $g^{2}/(\omega+\varepsilon)$ in both effective
Hamiltonians considered, therefore the middle points of their spectra
coincide.

It is known that the spectrum of ${\hat{h}}_{n}$ is $E_{n,\pm}=\alpha_{n}\pm
R_{n}$ with the eigenvectors (i.e. dressed states) given by the linear
combination

$\ket{\varphi_{n,\pm}}:=u_{\pm}\ket{n,\uparrow}+v_{\pm}\ket{n+1,\downarrow}$
(31)

where $(u_{+},v_{+})=(\cos\theta_{n}/2,\sin\theta_{n}/2)$ for the higher state
$\ket{\varphi_{n,+}}$ and $(u_{-},v_{-})=(\sin\theta_{n}/2,-\cos\theta_{n}/2)$
for the lower state $\ket{\varphi_{n,-}}$. To complete the picture we also
calculated give the energy of the ground state $\ket{0,\downarrow}$ as
$E_{0}=-\varepsilon/2-(\Omega/2+\delta)$.

By examining Eqs.(26) and (27), it is clear that the differences between the
two models are more pronounced for small values of $\varepsilon$. As already
mentioned, at $\varepsilon=0$, $H_{P}$ recovers the exact IBM result, while
$H_{BS}$ does not.

On the other hand at resonance, $\varepsilon=\omega$, the parameters of the
two models coincide, except for the sign of $\Omega$. The latter has no
bearing on the spectra, which therefore are the same for the two models, but
it modifies the structure of the eigenstates. Indeed, at resonance the sign of
$\Omega$ carries over to a sign change in $\cos\theta_{n}$ which in turn leads
to exchanging $\cos\theta_{n}/2$ and $\sin\theta_{n}/2$. In other words, even
though at resonance the spectra of $H_{BS}$ and $H_{P}$ coincide, because of
the sign difference in $\Delta_{n}$ the attribution of the eigenvalues is
flipped: for example, the weight of the two basis vectors $\ket{n,\uparrow}$
and $\ket{n+1,\downarrow}$ in the upper state $\ket{\varphi_{n,+}}$ of the BS
Hamiltonian is the same as for the lower state $\ket{\varphi_{n,-}}$ of the
perturbative Hamiltonian. Further consequences of this fact are analyzed in
subsection III.3.

To resume the situation, the two models predict different spectra and states,
as suggested by these two examples. Intermediate values for
$\varepsilon\in[0,\omega]$ are expected to lead to similar conclusions. Let us
also remark that $\Omega$ is interpreted as a bosonic frequency shift which
depends on the qubit state, and detecting the former gives information about
the latter RMP . Therefore the sign difference means that this way of reading
the qubit state works in opposite directions in the predictions of the two
models.

Still, these arguments are merely heuristic. Eigenstates are not directly
comparable, $H_{P}$ and $H_{BS}$ having been obtained by different unitary
transforms, i.e. the two effective Hamiltonians work in different bases. The
relevant comparison is on basis-independent quantities, to which we turn
below.

###  III.1 Energy spectra

Figure 1: Two spectral branches of the Rabi (R), Bloch-Siegert (BS) and
perturbative (P) Hamiltonians as a function of the matter-photon coupling $g$
at different values of the qubit energy: (a) $\varepsilon=0$, (b)
$\varepsilon=0.33$ and (c) $\varepsilon=0.66$. The spectrum of the QRM is
calculated numerically. $g$ is restricted to the perturbative USC regime.
Panel (a) also includes the branches $E_{3,+}$ and $E_{5,-}$ (see the text).

At resonance, $\varepsilon=\omega$, all parameters in Eq.(28) are exactly the
same for the two cases, except for sign of $\Delta_{n}$, and thus the spectra
are identical. On the other hand, in the limit $\varepsilon=0$, $H_{P}$
becomes exact. Therefore, it is expected that for the whole subresonant range,
$\varepsilon\leqslant\omega$, our approach is an improvement over the BS
efective Hamiltonian.

To illustrate this fact we compare in Fig. 1 the spectra of the two effective
models (BS and perturbative), and those produced by numerical diagonalization
of the full Rabi Hamiltonian. For the latter a cutoff at $N_{{\rm ph}}=30$
photons is found to be sufficient. We use $\omega=1$ units throughout, since
the problem is homogeneous of degree $1$ in the parameters, i.e. it is defined
essentially by only two dimensionless parameters $\varepsilon/\omega$ and
$g/\omega$.

We selected the pair of spectral branches $E_{4,\pm}$ for three values of the
qubit energy $\varepsilon$. However, in the IBM limit (i.e for
$\varepsilon=0$, see Fig. 1(a)) two additional branches are added, to
illustrate that the coinciding spectra of the Rabi and perturbative
Hamiltonians are doubly degenerate, while the Bloch-Siegert spectrum displays
a degeneracy lifting as the qubit-photon coupling strenght increases. Note
that for the sake of clarity the spectral branches of $H_{P}$ are represented
by red circles.

To understand these differences one has to recall that in the limit
$\varepsilon=0$, $H_{P}$ becomes exact, while the BS spectra are correct only
to $\lambda^{2}$ accuracy, as it is easy to check on the expression for
$E_{n,\pm}$. With increasing $\varepsilon$ the spectral difference between the
models diminishes, see Fig. 1(b),(c), as predicted by the analytic expression
of $H_{P}$ and $H_{BS}$ until, for $\varepsilon=1$ they coincide (not shown).
In the process the perturbative values remain closer to the exact Rabi result
than the BS ones.

###  III.2 Operator averages. Photon population

When calculating the expectation value of an operator on a given eigenstate of
the effective Hamiltonian, one has to take into acccount the unitary
transforms that lead to it. The operators should be ’rotated’ by the same
procedure as the Hamiltonian.

In the case of $H_{P}$, we use again the two generators
$S=\lambda(a^{\dagger}-a)\sigma_{x}$ and
$S^{\prime}=\frac{1}{2}\eta(a^{\dagger\,2}-a^{2})\sigma_{z}$, with the
parameters given in Eqs.(13) and (23). The observable we are interested in is
the photon number $\hat{n}=a^{\dagger}a$. From Eq.(7) we obtain by the first
transform

$a^{\dagger}_{S}a_{S}=a^{\dagger}a-\lambda(a^{\dagger}+a)\sigma_{x}+\lambda^{2}\,.$
(32)

The second step involves $S^{\prime}$ and is apparently more complicated but
it gets simpler by limiting the expansion to second order in $\lambda$ and
taking into account that $\eta\sim\lambda^{2}$. For instance, it is clear that
the second term in Eq.(32) remains unchanged. Indeed, changes are at least of
first power of $\eta$ and with the prefactor $\lambda$ they are third order
quantities.

For the first term $a^{\dagger}a$ one uses Eq.(18) expanded to second order in
$\lambda$, with the net result that after both transforms one has

$\hat{n}_{S,S^{\prime}}\approx\hat{n}-\eta(a^{\dagger\,2}+a^{2})\sigma_{z}-\lambda(a^{\dagger}+a)\sigma_{x}+\lambda^{2}\,.$
(33)

This is the operator to be averaged on the eigenstates of the form given in
Eq.(31). Things turn out to be further simplified by noting that operators
like $a^{2},a\sigma$ and their hermitic conjugates do not contribute to the
expectation value, since they do not connect vectors belonging to the same
two-dimensional subspace.

To conclude, one has

$\displaystyle\left<\hat{n}\right>_{P,\pm}=$
$\displaystyle\bra{\varphi_{n,\pm}}[a^{\dagger}a-\lambda(a^{\dagger}\sigma+a\sigma^{\dagger})+\lambda^{2}]\ket{\varphi_{n,\pm}}$
$\displaystyle=$ $\displaystyle u^{2}_{\pm}n+v^{2}_{\pm}(n+1)-2\lambda
u_{\pm}v_{\pm}\sqrt{n+1}+\lambda^{2}\,,$ (34)

where the index $P$ recalls that the average is calculated w.r.t. the dressed
states of $H_{P}$. Obviously, $\left<\hat{n}\right>_{P,\pm}$ depends on the
specific two-dimensional subspace $n$ as well but we drop this additional
index for the simplicity of writing. Replacing $u_{\pm}$ and $v_{\pm}$ leads
to the values

$\displaystyle\left<\hat{n}\right>_{P,\pm}=\,n+\frac{1}{2}(1\mp\cos\theta_{n})\mp\lambda\sin\theta_{n}\sqrt{n+1}+\lambda^{2}\,.$
(35)

A negative sign of $\Omega$ means a lower $\cos\theta_{n}$ value, and this
encourages a higher photon population in the upper state. For the ground state
one obtains $\left<{\hat{n}}\right>_{P,0}=\lambda^{2}$.

The prediction of the Bloch-Siegert Hamiltonian for the same quantity is
evaluated by the same procedure. Some differences arise though, beginning with
the first generator which now reads
$\tilde{S}=\lambda(a^{\dagger}\sigma^{\dagger}-a\sigma)$ with the result that
now

$a_{\tilde{S}}^{\dagger}a_{\tilde{S}}=a^{\dagger}a-\lambda(a^{\dagger}\sigma^{\dagger}+a\sigma)+\lambda^{2}\sigma\sigma^{\dagger}-\lambda^{2}a^{\dagger}a\sigma_{z}\,.$
(36)

The second transform has the same form
$\tilde{S}^{\prime}=\frac{1}{2}\eta(a^{\dagger\,2}-a^{2})\sigma_{z}$ with
$\eta\sim\lambda^{2}$.

As before, the terms containing $\lambda$ remain unchanged up to the second
order, and the first term $a^{\dagger}a$ is modified by a combination of
$a^{2}$ and $a^{\dagger\,2}$, which do not contribute to the expectation
value. The new feature is that no contribution comes from the second term in
Eq.(36) above. Therefore one is lead to

$\displaystyle\left<\hat{n}\right>_{BS,\pm}=$
$\displaystyle\bra{\varphi_{n,\pm}}[a^{\dagger}a+\lambda^{2}\sigma\sigma^{\dagger}-\lambda^{2}a^{\dagger}a\sigma_{z}]\ket{\varphi_{n,\pm}}$
$\displaystyle=$
$\displaystyle\,u^{2}_{\pm}(n-\lambda^{2}n)+v^{2}_{\pm}(n+1+\lambda^{2}(n+2))\,.$
(37)

As a result one has

$\displaystyle\left<\hat{n}\right>_{BS,\pm}=\,n+\frac{1}{2}+\lambda^{2}\mp\left(\frac{1}{2}+\lambda^{2}(n+1)\right)\cos\theta_{n}\,.$
(38)

For the ground state one finds $\left<\hat{n}\right>_{BS,0}=\lambda^{2}$.

Figure 2: The mean photon number $N$ associated to the eigenstates
$|\varphi_{4,\pm}\rangle$ of the Rabi (R), Bloch-Siegert (BS) and perturbative
(P) Hamiltonians as a function of the matter-photon coupling $g$ at different
values of the qubit energy: (a) $\varepsilon=0$, (b) $\varepsilon=0.33$ and
(c) $\varepsilon=0.66$.

The numerical results presented in Fig. 2 complement our analysis. We
calculated the mean photon number $N$ corresponding to the eigenstates of the
fourth two-dimensional Jaynes-Cummings subspace for the exact and the
effective Hamiltonians. Again, we find that in the IBM limit (see Fig. 2(a))
the perturbative Hamiltonian reproduces the exact expectation value of the
Rabi model (again, we used red circles to highlight the results of the
perturbative model). Note that for $\varepsilon=0$ one has $S^{\prime}=0$ too.
Thus one is left with only the first transform Eq.(32), which is exact. Thus
not only the Hamitonian, but the ’rotated’ form of the $\hat{n}$ operator is
exact. In contrast, the results obtained from the Bloch-Siegert model show
noticeable differences starting around the qubit-photon coupling strenght
$g\sim 0.15$. In particular, the mean photon number corresponding to
$|\varphi_{4,+}\rangle$ decreases as $g$ increases. These differences are even
more pronounced for higher two-dimensional subspaces. As in the case of the
spectrum, the photon number expectation values given by the two models become
closer as $\varepsilon$ increases, with the perturbative ones remaining nearer
to the exact Rabi values (see Figs. 2(b),(c)).

###  III.3 Operator averages. Population inversion and number of excitations

More clearcut differences exist between the two models than suggested by the
above results, especially at resonance. We show this by performing the same
comparison of expectation values, this time for the population inversion
operator $\sigma_{z}$ and for the number of excitations
$\hat{N}_{x}=\hat{n}+1/2(1+\sigma_{z})$.

For $H_{P}$ the change of $\sigma_{z}$ under the first unitary transform is
obtained as in the derivation of Eqs.(10),(11).

$\displaystyle\sigma_{z,S}=$
$\displaystyle\,\cosh[2\,\lambda(a^{\dagger}-a)]\,\sigma_{z}-\sinh[2\,\lambda(a^{\dagger}-a)]\,i\sigma_{y}$
$\displaystyle\approx$
$\displaystyle\,\sigma_{z}+2\lambda^{2}(a^{\dagger}-a)^{2}\sigma_{z}-2\lambda(a^{\dagger}-a)i\sigma_{y}\,.$
(39)

The expression remains unchanged by the second transform, which commutes with
$\sigma_{z}$. The terms containing $\lambda$ are left unchanged, as discussed
in the previous cases. Also, as before, the operators $a^{2},a^{\dagger\,2}$
and the off-resonant ones $a^{\dagger}\sigma^{\dagger},a\sigma$ do not
contribute to the expectation value. We are left with

$\sigma_{z,S,S^{\prime}}=\,\sigma_{z}+2\lambda(a^{\dagger}\sigma+a\sigma^{\dagger})-2\lambda^{2}(2a^{\dagger}a+1)\sigma_{z}\,,$
(40)

which by straightforward calculations leads to

$\displaystyle\left<\sigma_{z}\right>_{P,\pm}$
$\displaystyle=(u^{2}_{\pm}-v^{2}_{\pm})[1-2\lambda^{2}(2n+1)]$
$\displaystyle+\,4\lambda u_{\pm}v_{\pm}\sqrt{n+1}+4\lambda^{2}v^{2}\,.$ (41)

Specifically for the two states one obtains

$\displaystyle\left<\sigma_{z}\right>_{P,\pm}=$
$\displaystyle\,\pm\cos\theta_{n}[1-4\lambda^{2}(n+1)]$ $\displaystyle\pm$
$\displaystyle 2\lambda\sin\theta_{n}\sqrt{n+1}+2\lambda^{2}\,.$ (42)

A negative value for $\cos\theta_{n}$ hints to a higher population inversion
on the lower state. For the ground state one finds
$\left<\sigma_{z}\right>_{P,0}=-1+2\lambda^{2}$.

Finally, we address the same problem for $H_{BS}$. The transform of
$\sigma_{z}$ under $\tilde{S}=\lambda(a^{\dagger}\sigma^{\dagger}-a\sigma)$ is
complicated, but for the second order in $\lambda$ the quantities
$[\tilde{S},\sigma_{z}]$ and $[\tilde{S},[\tilde{S},\sigma_{z}]]$ are
sufficient. It is easy to check that

$\displaystyle[\tilde{S},\sigma_{z}]=$
$\displaystyle\,-2\lambda(a^{\dagger}\sigma^{\dagger}+a\sigma)\,,$
$\displaystyle[\tilde{S},a^{\dagger}\sigma^{\dagger}+a\sigma]=$
$\displaystyle\,2\lambda(a^{\dagger}a\sigma_{z}-\sigma\sigma^{\dagger})\,.$
(43)

This leads to the result

$\sigma_{z,\tilde{S}}=\sigma_{z}-2\lambda(a^{\dagger}\sigma^{\dagger}+a\sigma)-2\lambda^{2}(a^{\dagger}a\sigma_{z}-\sigma\sigma^{\dagger})\,,$
(44)

which remains unchanged under the second unitary transform. Thus

$\displaystyle\left<\sigma_{z}\right>_{BS,\pm}=$
$\displaystyle\bra{\varphi_{n,\pm}}[\sigma_{z}-2\lambda^{2}(a^{\dagger}a\sigma_{z}-\sigma\sigma^{\dagger})]\ket{\varphi_{n,\pm}}$
$\displaystyle=$
$\displaystyle(u_{\pm}^{2}-v_{\pm}^{2})(1-2\lambda^{2}n)+4\lambda^{2}v_{\pm}^{2}\,.$
(45)

Specifically, for the two states one has

$\displaystyle\left<\sigma_{z}\right>_{BS,\pm}=$
$\displaystyle\,\pm\cos\theta_{n}\,(1-2\lambda^{2}n)+2\lambda^{2}(1\mp\cos\theta_{n})\,$
$\displaystyle=$
$\displaystyle\,\pm\cos\theta_{n}\,[1-2\lambda^{2}(n+1)]+2\lambda^{2}\,.$ (46)

For the ground state one finds
$\left<\sigma_{z}\right>_{BS,0}=-1+2\lambda^{2}$.

The average of the operator $\hat{N}_{x}$ is now easily done for both
effective Hamiltonians by assembling Eqs.(35), (38),(42) and (46). The
surprising result is that the outcome is the same for both effective models,
and reads

$N_{x}^{\pm}=\langle\varphi_{n,\pm}|{\hat{N}}_{x}|\varphi_{n,\pm}\rangle=n+1+2\lambda^{2}\mp
2\lambda^{2}(n+1)\cos\theta_{n}\,.$ (47)

Figure 3: Two branches of the number of excitations $N_{x}^{\pm}$ associated
to the pair of dressed states $|\varphi_{4,\pm}\rangle$ as a function of the
qubit-photon coupling strenght $g$ in the resonant regime $\varepsilon=1$. (a)
$N_{x}^{-}$, the expectation value of the excitation number for the lower
spectral branch; (b) $N_{x}^{+}$, the same for the upper spectral branch. The
obvious similarity between the results of the Rabi and the perturbative models
confirms the correct assignment of $N_{x}^{\pm}$ to the upper/lower spectral
branches.

This is not a mere coincidence, and relies essentially on the fact that the
generators of the respective unitary transforms differ by
$\bar{S}:=S-\tilde{S}=\lambda(a^{\dagger}\sigma-a\,\sigma^{\dagger})$, which
commutes with $\hat{N}_{x}$ Comm-Nx .

Again, the interesting situation is provided by the resonant case. There one
can see analytically that both models predict exactly the same results for a
given spectral subspace. The curves for $N_{x}^{\pm}$ are superimposed since,
except for the sign, the cosine values coincide. But because of this sign
difference, inherited from the opposite signs of $\Omega_{P}$ and
$\Omega_{BS}$, the assignment of $N_{x}$ expectation values to the energy
branches is switched, $N_{x,BS}^{\pm}=N_{x,P}^{\mp}$. More precisely, in the
BS model $N_{x}^{-}$, corresponding to the lower branch, takes the higher
value, while $N_{x}^{+}$, stemming from the upper branch, is smaller. In the
perturbative case, it is the other way around. This is an important
difference, expressing the distinct structure of the eigenstates, and which
cannot gradually disappear as $\lambda$ gets smaller. To decide which
attribution is the correct one we rely on comparison with the result of the
full Rabi Hamiltonian.

The outcome is shown in Fig. 3, in which the values of the two models are
plotted, along with the numerical calculation using the dressed states
provided by diagonalization of the Rabi Hamiltonian. The left panel
corresponds to the eigenstates $\ket{\varphi_{4,-}}$ and the right one to
$\ket{\varphi_{4,+}}$. One can check by inspection of Figs. 3 (a) and (b) that
$N_{x,P}^{\pm}=N_{x,BS}^{\mp}$, and that comparison with the Rabi solutions
decides in favor of the perturbative results.

Figure 4: The number of excitations $N_{x}$ associated to the pair of dressed
states in the fourth spectral subspace as a function of the field-matter
coupling strength $g$, for the perturbative (P), Bloch-Siegert (BS) and Rabi
(R) Hamiltonians. (a) The IBM-limit $\varepsilon=0$, (b) $\varepsilon=0.33$,
(c) $\varepsilon=0.66$.

Finally we compare in Fig. 4 the number of excitations for the same values of
the qubit energy considered in previous discussions. In the IBM-limit
$N_{x,P}^{\pm}$ follow qualitatively the exact curves but also deviate from
them as $g$ increases. (Let us stress that for $\varepsilon=0$ the exact
result is analytically available in terms of Laguerre polynomials Mahan , but
here we stick to Eq.(47), in which the operators are rotated only to second
order). On the other hand, the excitation numbers of the BS model are quite
far from the $N_{x,R}^{\pm}$ for all values of $\varepsilon$. In particular,
at $\varepsilon=0.66$ the branches of $N_{x,BS}^{\pm}$ exhibit an incorrect
crossing point around $g=0.235$. Moreover, the BS attribution of the
excitation number to a given spectral branch is still flipped and the
inequality $N_{x,BS}^{+}<N_{x,BS}^{-}$ still holds. On the contrary, the
branches $N_{x,P}^{\pm}$ follow closely the exact ones in the sub-resonant
regime too.

###  III.4 The generalized RWA effective Hamiltonian

A successful effective Hamiltonian suitable in the deep strong coupling (DSC)
regime defined by $g/\omega\geqslant 1\,$ RMP  was proposed by Irish Irish .
This so called generalized RWA (GRWA) Hamiltonian follows from $H_{R}$ (see
Eq.(2)) by the IBM transform. The free qubit part of the Hamiltonian then
couples all pairs of shifted oscillator states.

To be specific, the model uses a single canonical transform, generated by
$S=\lambda(a^{\dagger}-a)\sigma_{x}$, this time with $\lambda=g/\omega$. As a
result one obtains, as for Eq.(11)

$\displaystyle H_{R,S}=\omega a^{\dagger}a-g^{2}/\omega$
$\displaystyle+\,\frac{1}{2}\varepsilon\left\\{\cosh[2\,\lambda(a^{\dagger}-a)]\,\sigma_{z}-\sinh[2\,\lambda(a^{\dagger}-a)]\,i\sigma_{y}\right\\}\,.$
(48)

It is lengthy but rather straighforward to expand the hyperbolic functions as
the even and odd part of the exponential (see Xie et al. Xie  but beware of
typos)

$\displaystyle e^{2\,\lambda(a^{\dagger}-a)}$
$\displaystyle=f_{0}(\hat{n},\lambda)+$
$\displaystyle\sum_{q>0}\left[a^{\dagger\,q}f_{q}(\hat{n},\lambda)+(-1)^{q}f_{q}(\hat{n},\lambda)a^{q}\right]\,,$
(49)

where the functions of $\hat{n}$ are defined, as usual, by their action on the
elements of the Fock basis
$f_{q}(\hat{n},\lambda)\ket{n}=f_{q}(n,\lambda)\ket{n}$. They are related to
the generalized Laguerre polynomials $L_{n}^{(q)}(x)$ by

$f_{q}(n,\lambda)=e^{-2\lambda^{2}}(2\lambda)^{q}\,\frac{n!}{(n+q)!}\,L^{(q)}_{n}(4\lambda^{2})\,.$
(50)

This translates into the following expression for $H_{R,S}$

$\displaystyle H_{R,S}=\,$
$\displaystyle\omega\,a^{\dagger}a-g^{2}/\omega+\frac{\varepsilon}{2}f_{0}(\hat{n},\lambda)\sigma_{z}$
$\displaystyle+\frac{\varepsilon}{2}\sum_{q>0,\rm{even}}\,\left[a^{\dagger\,q}f_{q}(\hat{n},\lambda)+f_{q}(\hat{n},\lambda)a^{q}\right]\\}\,\sigma_{z}$
$\displaystyle+\frac{\varepsilon}{2}\sum_{q>0,\rm{odd}}\,\left[a^{\dagger\,q}f_{q}(\hat{n},\lambda)-f_{q}(\hat{n},\lambda)a^{q}\right](\sigma-\sigma^{\dagger})\,.$
(51)

No approximation is involved so far. The $f_{0}$ term is diagonal, while terms
with higher $q$-values are further away from the diagonal in the Fock basis.

The idea of the GRWA is to keep, besides the diagonal terms, only the off-
diagonal ones with $q=1$ and which, moreover, are also retained in the RWA, as
not counter-rotating. More precisely, one has

$\displaystyle H_{{\rm GRWA}}=$
$\displaystyle\,\omega\,a^{\dagger}a-g^{2}/\omega+\,\frac{\varepsilon}{2}\,f_{0}(\hat{n},\lambda)\,\sigma_{z}$
$\displaystyle+$ $\displaystyle\,\frac{\varepsilon}{2}\,[a^{\dagger}\sigma
f_{1}(\hat{n},\lambda)+f_{1}(\hat{n},\lambda)\,a\sigma^{\dagger}]\,,$ (52)

which reduces the Hamiltonian to the familiar JC structure, renormalized by
$\hat{n}$-dependent factors. Thus the problem becomes $2\times 2$ block-
diagonal and therefore solvable.

Obviously, this approach is not a systematic expansion to $\lambda^{2}$.
Indeed, $f_{q}\sim\lambda^{q}$ and therefore terms with $f_{2}$ should have
been kept as second order contribution. Instead they were discarded as off-
resonant, as prescribed by the RWA argument.

Figure 5: Several spectral branches of the Rabi (R), GRWA and perturbative (P)
Hamiltonians as a function of the matter-photon coupling $g$ in the (a)
resonant regime $\varepsilon=1$ and (b) supra-resonant regime $\varepsilon=2$.
$g$ is restricted to the perturbative USC regime.

The success of the model in the DSC regime relies on the fact that all terms
arising from the transformed qubit energy become small at large $g$, due to
the exponential prefactor. Suprisingly, its predictions work quite well at
lower couplings too. This means that at least in certain circumstances the
off-resonant terms are negligible indeed. On the contrary, in situations like
$\varepsilon\approx 2\omega$ this is not the case any more. Terms like
$a^{\dagger\,2}\sigma$ now describe resonant transitions in the system and
thus are no longer negligible. In contrast, these processes are included in
the present perturbative approach by its second unitary transform
$S^{\prime}$.

The situation is illustrated in Fig. 5 where we compare a part of the low-
energy spectra of the Rabi, GRWA and perturbative Hamiltonians in the resonant
($\varepsilon=\omega$) and supra-resonant ($\varepsilon>\omega$) regimes. The
eigenvalues of the GRWA Hamiltonian Irish  are given analytically by Eq.(52).
Clearly, the eigenvalues of both GRWA and perturbative Hamiltonians reasonably
follow the Rabi spectrum in the resonant case along the whole range of the
pUSC regime (see Fig. 5(a)). However, Fig. 5(b) reveals that in the supra-
resonant regime the GRWA spectrum is no longer accurate in the pUSC regime
while the errors of the perturbative model are quite small.

##  IV Conclusions

The effective Hamiltonian given by Eq.(24) was derived from the quantum Rabi
model systematically to second order in the small parameter
$\lambda=g/(\varepsilon+\omega)$. No other approximations (like discarding
terms as off-resonant) are involved. In this sense the proposed model is
purely perturbative. The expansion parameter $\lambda$ is the same as in the
Bloch-Siegert theory. In contrast to the latter, the present formalism
recovers exactly the IBM Hamiltonian in the limit $\varepsilon\to 0$. The
proposed model allows analytical calculations of various operator averages
which are used for a detailed comparison with the predictions of the BS
approach.

Even in the cases when the spectra of two models are close, we find that the
main difference from the BS results is seen in the structure of the dressed
eigenstates. The source of this dissimilarity is traced back to the different
sign of the frequency shift induced by the qubit-field coupling. This is
furher reflected in the expectation values of relevant observables, like the
excitation number, and their attribution to the appropriate eigenstates.

The energy spectrum and the dressed state structure associated to the full
quantum Rabi model are well approximated by the IBM-compatible effective
Hamiltonian in the resonant and off-resonant perturbative USC regime. Our
analysis shows that recovering the independent-boson model in the limit
$\varepsilon\to 0$ is a valuable test for effective Hamiltonian approaches to
the quantum Rabi model.

###### Acknowledgements.

V.M. acknowledges the financial support from the Romanian Core Program PN19-03
(contract no. 21 N/08.02.2019).

## References

  * (1) A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta and F. Nori, Nature Reviews Physics 1, 19 (2019).
  * (2) P. Forn-Díaz, L. Lamata, E. Rico, J. Kono, and E. Solano, Rev. Mod. Phys. 91, 025005 (2019).
  * (3) E.T. Jaynes, F.W. Cummings, Proc. IEEE. 51, 89 (1963).
  * (4) P. Meystre, M. Sargent, Elements of Quantum Optics, Springer (2007).
  * (5) D. F. Walls, Gerald J. Milburn, Quantum Optics, 2nd Edition, Springer (2008).
  * (6) D. Braak, Phys. Rev. Lett. 107, 100401 (2011).
  * (7) F. A. Wolf, M. Kollar, and D. Braak, Phys. Rev. A 85, 053817 (2012).
  * (8) Qing-Wei Wang, Yu-Liang Liu, J. Phys. A:Math. Theor. 46, 435303 (2013).
  * (9) H. Zhong, Q. Xie, M. T. Batchelor, C. Lee, J. Phys. A:Math. Theor. 46, 415302 (2013).
  * (10) M. Cirio, S. De Liberato, N. Lambert, and F. Nori, Phys. Rev. Lett. 116, 113601 (2016).
  * (11) I. V. Dinu, V. Moldoveanu, and P. Gartner, Phys. Rev. B 97, 195442 (2018).
  * (12) S. De Liberato, D. Gerace, I. Carusotto, and C. Ciuti, Phys. Rev. A 80, 053810 (2009).
  * (13) F. Beaudoin, J. M. Gambetta, and A. Blais, Phys. Rev. A 84, 043832 (2011).
  * (14) A. Settineri, V. Macrí, A. Ridolfo, O. Di Stefano, A. F. Kockum, F. Nori, and S. Savasta, Phys. Rev. A 98, 053834 (2018).
  * (15) G. D. Mahan, Many-Particle Physics, Plenum Press (1990).
  * (16) B. Krummheuer, V. M. Axt, and T. Kuhn Phys. Rev. B 65, 195313 (2002).
  * (17) A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
  * (18) D. Zueco, G. M. Reuther, S. Kohler, and P. Hänggi, Phys. Rev. A 80, 033846 (2009).
  * (19) A. Blais, R-S Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 69, 062320 (2004).
  * (20) A. B. Klimov, S. M. Chumakov, A Group-Theoretical Approach to Quantum Optics: Models of Atom-Field Interactions, Wiley (2009).
  * (21) D. Z. Rossatto, C. J. Villas-Bôas, M. Sanz, and E. Solano, Phys. Rev. A 96, 013849 (2017).
  * (22) P. Forn-Díaz, J. Lisenfeld, D. Marcos, J. J. García-Ripoll, E. Solano, C. J. P. M. Harmans, and J. E. Mooij, Phys. Rev. Lett. 105, 237001 (2010).
  * (23) J. Hausinger, M. Grifoni, New Journal of Physics 10, 115015 (2008).
  * (24) Other technicalities play a role too: the commutator $[\tilde{S},\bar{S}]$ generates only off-diagonal terms, which do not influence the $\hat{N}_{x}$ expectation value.
  * (25) E. K. Irish, Phys. Rev. Lett. 99, 173601 (2007).
  * (26) L. Yu, S. Zhu, Q. Liang, G. Chen, and S. Jia, Phys. Rev. A 86, 015803 (2012).
  * (27) Weijie Xie et al. J. Phys. A: Math. Theor. 53, 095302 (2020).
  * (28) Y-Y. Zhang, Phys. Rev. A 94, 063824 (2016).
  * (29) Y-Y. Zhang, X-Y. Chen, Phys. Rev. A 96, 063821 (2017).
  * (30) Let us stress that the pure squeezing transform used in Refs.Zhang2017 ; Xie  fails to cancel the quadratic operators.


# Compositional Estimation of Lipschitz Constants
for Deep Neural Networks

Yuezhu Xu and S. Sivaranjani∗ *Corresponding author.The authors are with the
School of Industrial Engineering, Purdue University, West Lafayette, IN, USA.
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Lipschitz constant plays a crucial role in certifying the robustness of
neural networks to input perturbations and adversarial attacks, as well as the
stability and safety of systems with neural network controllers. Therefore,
estimation of tight bounds on the Lipschitz constant of neural networks is a
well-studied topic. However, typical approaches involve solving a large matrix
verification problem, the computational cost of which grows significantly for
deeper networks. In this letter, we provide a compositional approach to
estimate Lipschitz constants for deep feedforward neural networks by obtaining
an exact decomposition of the large matrix verification problem into smaller
sub-problems. We further obtain a closed-form solution that applies to most
common neural network activation functions, which will enable rapid robustness
and stability certificates for neural networks deployed in online control
settings. Finally, we demonstrate through numerical experiments that our
approach provides a steep reduction in computation time while yielding
Lipschitz bounds that are very close to those achieved by state-of-the-art
approaches.

##  I Introduction

Neural networks are increasingly common in control to capture the nonlinear
behavior of dynamical systems and to design sophisticated controllers [1][2],
with several applications including but not limited to autonomous vehicles
[3], robotics [4], demand forecasting in smart grids [5], and modeling in HVAC
systems [6]. In such applications, it is crucial to guarantee the safety and
robustness of the system when utilizing neural networks in control loops. In
the context of learning, it is also desirable to guarantee robustness against
adversarial attacks [7][8] and overfitting to enhance generalization [9].

The Lipschitz constant, which quantifies how a neural network’s output varies
in response to changes in its inputs, is a crucial measure for certifying
robustness and safety. Mathematically, a smaller Lipschitz constant indicates
greater robustness to input perturbations. In fact, regularizing neural
networks to achieve a desired Lipschitz constant is a common practice to
enhance robustness during training [10]-[12]. In control applications, the
Lipschitz constant can be used to provide stability guarantees [13][14], and
guarantee safety-critical properties [15][16], especially when neural networks
are used in perception-control [17].

However, it is NP-hard to calculate the exact Lipschitz constant for
complicated nonlinear functions like neural networks. Therefore, in recent
years, efforts have been made in estimating tight upper bounds for the
Lipschitz constant of feedforward neural networks and other architectures such
as convolutional neural networks (CNNs), either by formulating the problem as
a polynomial optimization [18] or by upper bounding the Lipschitz constant via
quadratic constraints and semidefinite programming (SDP) [19]-[21].

While such approaches are successful in providing tight Lipschitz bounds, the
computational cost explodes as the number of layers increases. One common
strategy to address this issue involves ignoring some coupling constraints
among the neurons to reduce the number of decision variables, yielding a more
scalable algorithm at the expense of estimation accuracy [19]. Another
strategy is to exploit the sparsity of the SDP using graph-theoretic
approaches in order to decompose it into smaller linear matrix inequalities
(LMI) [20][22]. Along similar lines, [21] employs a dissipativity-based method
to derive layer-wise LMIs. Yet, in practice, this requires pre-specifying the
dissipativity indices of each layer, which may require impractical manual
tuning, and may restrict the feasibility of the SDP. On the other hand,
treating these indices as decision variables in the SDP drastically increases
computational complexity, especially for deep networks. Finally, for deep
neural networks, the computation time, even for state-of-the-art algorithms,
is too large to be practical in control applications, where neural networks
need to be verified during real-time operation.

In this letter, we overcome these issues by developing a compositional method
to provide tight Lipschitz estimates for deep feedforward neural networks.
Specifically, the key contributions of this letter are as follows. First, we
obtain a series of sub-problems that are exactly equivalent to the SDP
verification problem in [13][19], wherein the sub-problems are only as large
as the size of each layer, rather than the entire neural network. Then, we
develop a compositional algorithm wherein these sub-problems can be solved
layer-by-layer for the neural network. Second, in contrast to works like [21],
we provide an approach to determine the best auxiliary parameters in sub-
problems to obtain a tight Lipschitz estimate. In fact, we show that these
sub-problems admit exact closed-form solutions for most commonly used neural
network activation functions. Finally, we demonstrate through experiments that
our approach provides a significant reduction in computation time for both
deeper and wider networks, while providing estimates that are very close to
those obtained using state-of-the-art algorithms [19]. This advantage is
particularly pronounced when we consider deeper neural networks. Overall, the
scalability of our approach will enable rapid robustness and stability
certificates for neural networks deployed in online control settings.

The letter is organized as follows. We introduce the problem formulation in
Section II, and present our methodology in Section III. In Section IV, we
propose a compositional algorithm based on the methodology in Section III, and
derive an exact solution that applies to most common neural network activation
functions. Finally, in Section V, we benchmark the performance of our
algorithm with respect to a state-of-the-art Lipschitz constant estimation
approach, and provide directions for future work in Section VI. The proofs of
all results are presented in the Appendix.

Notation: We denote the sets of positive real numbers and $n$-dimensional real
vectors by $\mathbb{R}_{++}$ and $\mathbb{R}^{n}$, respectively. We define
$\mathbb{Z}_{N}=\\{1,\ldots,N\\}$, where $N$ is a natural number excluding
zero. Given a matrix $A\in\mathbb{R}^{m\times n}$,
$A^{T}\in\mathbb{R}^{n\times m}$ represents its transpose. A symmetric
positive-definite matrix $P\in\mathbb{R}^{n\times n}$ is represented as $P>0$
(and as $P\geq 0$, if it is positive semi-definite). Similarly, a symmetric
negative definite matrix $P\in\mathbb{R}^{n\times n}$ is represented as $P<0$
(and as $P\leq 0$, if it is negative semi-definite). The standard identity
matrix is denoted by $I$, with dimensions clear from the context. $\|x\|_{2}$
represents the 2-norm of vector $x$ and $\|A\|_{2}$ means the matrix norm of
$A$ induced by the vector 2-norm. We denote by $eig_{min}(A)$ and
$eig_{max}(A)$ the smallest and the largest eigenvalue, respectively, for
symmetric matrix $A\in\mathbb{R}^{n\times n}$, where $A=A^{T}$.

##  II Problem Formulation and Background

###  II-A Feed-forward Neural Network and Lipschitz Bound

Consider a feedforward neural network (FNN) of $l$ layers with input
$z\in\mathbb{R}^{d_{0}}$ and output $y\in\mathbb{R}^{d_{l}}$ defined as
$y=f(z)$. The function $f$ is recursively formulated with layers
$\mathbf{L}_{i},i\in\mathbb{Z}_{l}$, whose mapping is defined as

$\begin{split}&\mathbf{L}_{i}:\,z^{(i)}=\phi(v^{(i)})\quad\forall
i\in\mathbb{Z}_{l-1},\\\ &\mathbf{L}_{l}:\,y=f(z)=z^{(l)}=v^{(l)},\end{split}$
(1)

where $z^{(0)}=z$, $v^{(i)}=W_{i}z^{(i-1)}+b_{i}$ with $W_{i}$ and $b_{i}$
representing the weight and bias for layer $\mathbf{L}_{i}$ respectively, and
$\phi:\mathbb{R}^{d_{i}}\to\mathbb{R}^{d_{i}}$ is a nonlinear activation
function that acts element-wise on its argument. The last layer
$\mathbf{L}_{l}$ is termed the output layer of the network. We denote the
number of neurons in layer $\mathbf{L}_{i}$ by $d_{i}$, $i\in\mathbb{Z}_{l}$.

###### Assumption 1

We assume that the function $f:\mathbb{R}^{d_{0}}\to\mathbb{R}^{d_{l}}$ is
Lipschitz continuous on $\mathcal{Z}\subseteq\mathbb{R}^{n}$. In other words,
there exists a constant $L_{0}\in\mathbb{R}_{++}$ such that

$\|f(z_{1})-f(z_{2})\|_{2}\leq L_{0}\|z_{1}-z_{2}\|_{2},\forall
z_{1},z_{2}\in\mathcal{Z}.$ (2)

The smallest positive $L_{0}$ satisfying (2) is termed the Lipschitz constant
of the function $f$.

The Lipschitz constant is widely used to quantify the robustness of neural
networks. Typically, a smaller Lipschitz constant indicates that the neural
network is more robust to input perturbations. However, it is not easy to
compute the Lipschitz constant for this complex nonlinear function $f$, and no
closed-form solution is possible. The common approach is to find a tight upper
bound [19][21], which can be used to provide robustness certificates on
downstream tasks that utilize the neural network. However, existing methods to
estimate the Lipschitz bound involve verification of the definiteness of a
large matrix that scales poorly to deep neural networks. Our goal is to
provide a scalable compositional approach to estimate the Lipschitz bound by
decomposing the large matrix verification problem into a series of layer-by-
layer verification problems with each involving matrices that are only as
large as the size of the layer.

##  III Methodology

In this section, we will first formulate the problem of estimating the
Lipschitz bound as a semidefinite program (SDP), which verifies the
definiteness of a large matrix. We then provide an exact decomposition of this
problem into layer-by-layer sub-problems that can be solved compositionally.

###  III-A Lipschitz Certificates for Feed-forward Neural Networks

We begin by deriving sufficient conditions for the neural network (1) to have
Lipschitz constant $L_{0}$. An important property used here is the assumption
of slope-restrictedness on the activation function, stated below.

######  Assumption 2 (Slope-restrictedness)

For the neural network defined in (1), the activation function $\phi$ is
slope-restricted in $[\alpha,\beta]$, $\alpha<\beta$ in the sense that
$\forall v_{1},v_{2}\in\mathbb{R}^{n}$, we have
$\alpha(v_{1}-v_{2})\leq\phi(v_{1})-\phi(v_{2})\leq\beta(v_{1}-v_{2})$
element-wise. Consequently, we have the following:

$\begin{bmatrix}v_{1}-v_{2}\\\
\phi(v_{1})-\phi(v_{2})\end{bmatrix}^{T}\begin{bmatrix}pI&-mI\\\
-mI&I\end{bmatrix}\begin{bmatrix}v_{1}-v_{2}\\\
\phi(v_{1})-\phi(v_{2})\end{bmatrix}\leq 0,$ (3)

where $p=\alpha\beta$ and $m=\frac{\alpha+\beta}{2}$.

###### Remark 1

The slope-restrictedness assumption is reasonable, since it holds for all
commonly used activation functions. For example, Assumption 2 holds with
$\alpha=0$, $\beta=1$, and $p=0$ for ReLU, sigmoid, tanh, exponential linear
functions. For the Leaky ReLU function $\phi(x)=\max(ax,x)$ ($a>0$), we have
$\alpha=\min(a,1)$, $\beta=\max(a,1)$, and $p>0$ [23, Proposition 2].

With this assumption, we can now provide a Lipschitz certificate, which is
used to verify if a constant $L\in\mathbb{R}_{++}$ is a valid Lipschitz bound
for the FNN (1).

$\displaystyle
P_{l}=\begin{bmatrix}cI+\lambda_{1}pW_{1}^{T}W_{1}&-\lambda_{1}mW_{1}^{T}&0&...&&0\\\
-\lambda_{1}mW_{1}&\lambda_{1}I+\lambda_{2}pW_{2}^{T}W_{2}&-\lambda_{2}mW_{2}^{T}&...&&0\\\
0&-\lambda_{2}mW_{2}&\lambda_{2}I+\lambda_{3}pW_{3}^{T}W_{3}&...&&0\\\
&&\vdots&&&\\\
0&...&-\lambda_{l-2}mW_{l-2}&\lambda_{l-2}I+\lambda_{l-1}pW_{l-1}^{T}W_{l-1}&-\lambda_{l-1}mW_{l-1}^{T}&0\\\
0&0&...&-\lambda_{l-1}mW_{l-1}&\lambda_{l-1}I+\lambda_{l}W_{l}^{T}W_{l}&-\lambda_{l}W_{l}^{T}\\\
0&0&...&0&-\lambda_{l}W_{l}&\lambda_{l}I\end{bmatrix}$ (4)

######  Theorem 1 (Lipschitz certificates)

Let $P_{l}$ be defined as in (4) and
$e_{i}^{l}\in\mathbb{R}^{\sum_{j=0}^{l}d_{j}}$ be defined as a vector with the
$(\sum_{j=0}^{l-1}d_{j}+i)-th$ entry being 1 and the remaining entries being
0. For the FNN (1), if there exist $c,L,\lambda_{i}\in\mathbb{R}_{++}$,
$i\in\mathbb{Z}_{l}$ such that

$P_{l}-\frac{1}{L^{2}}\sum_{i=1}^{d_{l}}e_{i}^{l}(e_{i}^{l})^{T}>0$ (5)

holds, then
$\|z_{2}^{(l)}-z_{1}^{(l)}\|_{2}\leq\sqrt{c}L\|z_{2}^{(0)}-z_{1}^{(0)}\|_{2}$.
In other words, the Lipschitz constant $L_{0}$ with respect to the input and
the output layer $\mathbf{L}_{l}$ is upper bounded by $\sqrt{c}L$.

###### Remark 2

The large matrix inequality in (5) corresponds exactly to LipSDP-Layer method
in [19] and is the fastest method among the three algorithms proposed in [19].
We will show in Section V that our compositional method outperforms the
LipSDP-Layer method in computation time while providing comparable estimates
for the Lipschitz constant.

Theorem 1 has two applications. First, we can directly utilize the
certificates in Theorem 1 to verify if some given $L_{0}\in\mathbb{R}_{++}$
upper bounds the Lipschitz constant. Second, with given constant $c$, we can
derive a tight estimate for the Lipschitz constant by formulating an SDP:

$\displaystyle H^{*}=\operatorname*{arg\,min}_{L}\quad-H$ (6) s.t.
$\displaystyle P_{l}>0,\quad H>0,\quad\lambda_{i}>0\quad i\in\mathbb{Z}_{l},$

where $H\triangleq\frac{1}{L^{2}}$. We denote $H^{*}$ to be the optimal
solution of (6); correspondingly, $L^{*}=1/\sqrt{H^{*}}$. Then, the best
estimate for $L_{0}$ is $L_{0}^{*}=\sqrt{c}L^{*}=\sqrt{c/{H^{*}}}$.

The computational cost of solving (6) explodes as the FNN (1) grows deeper.
This motivates us to provide a compositional approach to overcome this issue.

###  III-B Compositional Certification and Tight Lipschitz Estimate

In order to avoid directly solving the large matrix inequality (5), we will
leverage a sequential Cholesky decomposition similar to that developed in [24]
to derive necessary and sufficient conditions that only involve inequalities
on matrices of the same size as the number of weights in a single layer.

###### Theorem 2

The inequality (5) holds if and only if the following sequence of matrix
inequalities is satisfied:

$(M_{i})_{l}>0,\quad\forall i\in\\{0\\}\cup\mathbb{Z}_{l}$ (7)

where

$\displaystyle(M_{i})_{l}=\begin{cases}cI+\lambda_{1}pW_{1}^{T}W_{1}&\text{if
}i=0,\\\ \lambda_{i}I+\lambda_{i+1}pW_{i+1}^{T}W_{i+1}\\\
\qquad\qquad-\lambda_{i}^{2}m^{2}W_{i}(M_{i-1})_{l}^{-1}W_{i}^{T}&\text{if}\
i\in\mathbb{Z}_{l-1},\\\
\lambda_{i}I-\frac{1}{L^{2}}I-\lambda_{i}^{2}m^{2}W_{i}(M_{i-1})_{l}^{-1}W_{i}^{T}&\text{if
}i=l.\end{cases}$ (8)

Notice that in Theorem 2, all the calculations in the verification process
involve only the weights on a single layer and can be carried out in a
recursive manner. Therefore, Theorem 2 provides a faster compositional
approach to verify the Lipschitz constant of a deep FNN.

###### Remark 3

If the parameters $\lambda_{i},i\in\mathbb{Z}_{l}$ are known, this problem can
be solved as follows. Starting from $i=0$, we can verify condition (7)
sequentially for $i\in\\{0\\}\cup\mathbb{Z}_{l-1}$, since each $(M_{i})_{l}$
depends on the previous $(M_{i-1})_{l}$. Finally, if all verifications for
$i\in\\{0\\}\cup\mathbb{Z}_{l-1}$ are successful, we can optimize for the
smallest $L$ such that $(M_{l})_{l}>0$. However, if the parameters
$\lambda_{i},i\in\mathbb{Z}_{l}$ are unknown, it is challenging to compute
their optimal values to obtain a tight Lipschitz estimate. In the next
section, we will further analyze the layer-by-layer compositional structure in
Theorem 2 to propose a practical algorithm to obtain
$\lambda_{i},i\in\mathbb{Z}_{l}$ that results in a tight Lipschitz bound for
the FNN (1).

##  IV Algorithm and Practical Issues

Our ultimate goal is to find a tight Lipschitz bound $L_{0}$ with respect to
the inputs and the outputs of (1) using the fast compositional verification
result in Theorem 2. In order to develop a practical algorithm that yields the
best $\lambda_{i},i\in\mathbb{Z}_{l}$ satisfying Theorem 2, we will split the
FNN into two parts: $\mathbf{L}_{i},i\in\mathbb{Z}_{l-1}$, and
$\mathbf{L}_{l}$, and deal with them separately. Concretely, with given
constant $c\in\mathbb{R}_{++}$, we will first find the best choices for the
parameters $\lambda_{i},i\in\mathbb{Z}_{l-1}$, which will generate the
smallest Lipschitz estimate with respect to the input layer and layer
$\mathbf{L}_{l-1}$, termed $\bar{L}_{0}$. Next, in dealing with the last
layer, we will establish a relationship between $\bar{L}_{0}$ and $L_{0}$ to
obtain a tight Lipschitz bound for the FNN.

####  IV-1 Finding the best $\lambda_{i}$ for layers
$\mathbf{L}_{i},i\in\mathbb{Z}_{l-1}$

In dealing with the first $(l-1)$ layers, the sequence of matrices in (7) will
become $(M_{i})_{l-1}$, which is obtained by replacing each occurrence of the
index $l$ in (8) by $(l-1)$. Denote the solution to this problem by $\bar{L}$.
Then, similar to Theorem 1, we have the Lipschitz estimate
$\bar{L}_{0}=\sqrt{c}\bar{L}$ for the first $(l-1)$ layers. Recall that (8)
can be verified sequentially starting from the input layer and proceeding
layer-by-layer until the output layer. We observe in (8) that in the
verification problem for layer $\mathbf{L}_{i}$, $(M_{i})_{l-1}$ will involve
information from the next layer, namely $W_{i+1}$ and $\lambda_{i+1}$. In
order to develop a compositional algorithm, we first decouple the verification
conditions for layer $\mathbf{L}_{i}$ from parameters $\lambda_{i+1}$ and
weights $W_{i+1}$ as follows.

We define a sequence of matrices $X^{i},i\in\mathbb{Z}_{l-1}$ as

$\displaystyle X^{0}\triangleq cI$ (9) $\displaystyle
X^{i}\triangleq\lambda_{i}I-\lambda_{i}^{2}m^{2}W_{i}(M_{i-1})_{l-1}^{-1}W_{i}^{T},\quad
i\in\mathbb{Z}_{l-1}.$ (10)

Then, we have

$\displaystyle(M_{i})_{l-1}=\begin{cases}X^{i}+\lambda_{i+1}pW_{i+1}^{T}W_{i+1}&\text{if
}i\in\\{0\\}\cup\mathbb{Z}_{l-2},\\\ X^{i}-(1/\bar{L}^{2})I&\text{if}\
i=l-1.\end{cases}$ (11)

Accordingly, we have the following recursive relationship corresponding to the
layer $\mathbf{L}_{i},i\in\mathbb{Z}_{l-1}$:

$\displaystyle
X^{i}=\lambda_{i}I-\lambda_{i}^{2}m^{2}W_{i}(X^{i-1}+\lambda_{i}pW_{i}^{T}W_{i})^{-1}W_{i}^{T}.$
(12)

Note that $X_{i}$ no longer depends on any information from any layers
$\mathbf{L}_{j},j>i$.

###### Theorem 3

For any $c\in\mathbb{R}_{++}$,
$\lambda_{i}\in\mathbb{R}_{++},i\in\mathbb{Z}_{l-1}$, $p\geq 0$, and $\bar{L}$
that satisfies

$X^{i}>0\ \ \forall i\in\mathbb{Z}_{l-2},\quad X^{l-1}-(1/{\bar{L}^{2}})I>0,$
(13)

we have $(M_{i})_{l-1}>0$.

In other words, Theorem 3 is a sufficient condition for Theorem 2 with respect
to the first $(l-1)$ layers. Based on Theorem 3, we can provide a Lipschitz
certificate for the first $(l-1)$ layers as follows. First, given a constant
$c$, we can determine $\lambda_{1}$ by treating it as a variable that must
satisfy condition (13) in Theorem 3 with $i=1$. Similarly, at the $i$-th
stage, $i>1$, and $\lambda_{j},j\in\mathbb{Z}_{i-1}$, we can find
$\lambda_{i}$ satisfying condition (13) in Theorem 3, with the knowledge of
$W_{i}$ and $(M_{i-1})_{l-1}$ from the previous stage.

If we can somehow determine the ‘best’ $\lambda_{i},i\in\mathbb{Z}_{l-1}$, we
can then directly derive a closed-form expression for the smallest $\bar{L}$
with respect to the input layer and layer $\mathbf{L}_{l-1}$, as

$\displaystyle\bar{L}=\sqrt{1/(eig_{min}(\lambda_{l-1}I-\lambda_{l-1}^{2}m^{2}W_{l-1}(M_{l-2})_{l-1}^{-1}W_{l-1}^{T})})$
(14)

to ensure $(M_{l-1})_{l-1}=X^{l-1}-(1/\bar{L}^{2})I>0$ as in (13).

To formally define and thus set up the criterion for the ‘best’
$\lambda_{i},i\in\mathbb{Z}_{l-1}$, we need to first understand how $c$ and
the $\lambda_{i},i\in\mathbb{Z}_{l-1}$ impact the feasible region of $\bar{L}$
in the last stage. We analyze the problem recursively in a backward manner. At
the $(l-1)$-th stage, the constant $c$ and all $\lambda_{i}$,
$i\in\mathbb{Z}_{l-2}$ are already determined and all weights are known. Based
on (13), obtaining smallest estimate $\bar{L}$ is equivalent to finding the
best $\lambda_{i},i\in\mathbb{Z}_{l-1}$ that yields the largest
$eig_{min}(X^{l-1})$. When we move backwards to stage $i$, $\forall
i\in\mathbb{Z}_{l-2}$, the recursive relationship characterized in (12) is
much more complicated. However, we can still formulate a series of
optimization problems as follows to obtain the
$\lambda_{i},i\in\mathbb{Z}_{l-1}$ at each stage:

$\displaystyle\tilde{\lambda_{i}}=\operatorname*{arg\,min}_{\lambda_{i}}\quad-
eig_{min}(X^{i})$ (15) s.t. $\displaystyle X^{i}>0,\quad\lambda_{i}>0,\quad
W_{i},X^{i-1}\ \text{known}$

The intuition behind this formulation is that a larger $eig_{min}(X^{i-1})$ is
generally expected to result in a better (larger) $eig_{min}(X^{i})$.
Therefore, with fixed $c$, $\lambda_{j}$, $j\in\mathbb{Z}_{i-1}$ and known
$W_{j}$s, $j\in\mathbb{Z}_{i}$, we search for $\lambda_{i}$ that maximizes
$eig_{min}(X^{i})$, $i\in\mathbb{Z}_{l-1}$.

###### Remark 4

As we transform the centralized certificate (5) into distributed sufficient
but not necessary conditions (15), there is a natural trade-off between the
verification speed and the tightness of the Lipschitz bound. Later, in Section
V, we will experimentally demonstrate that this formulation significantly
reduces computation time while still yielding a tight bound that is only
slightly larger than the centralized estimates obtained from (5). This is
particularly useful in providing Lipschitz certificates for fairly deep FNNs.

####  IV-2 Exact Solution

Recall from Remark 1 that for almost all commonly used activation functions
(with the exception of Leaky ReLU), it holds that $p=0$. When $p=0$, we find
that the optimization problem (15) in fact admits an exact closed-form
solution as follows. When $p=0$, we have

$eig_{min}(X^{i-1})=\lambda_{i}-m^{2}\lambda_{i}^{2}eig_{max}(W_{i}(X^{i-1})^{-1}W_{i}^{T}),$

which is quadratic in $\lambda_{i}$. Therefore, the solution of (15) can be
explicitly written as

$\displaystyle\tilde{\lambda_{i}}=1/(2m^{2}eig_{max}(W_{i}(X^{i-1})^{-1}W_{i}^{T})).$
(16)

Additionally, we note that when $p=0$, from (11), we have exactly
$(M_{i})_{l-1}=X^{i},i\in\mathbb{Z}_{l-1}$. Therefore, Theorem 3 becomes
exactly equivalent to Theorem 2 for the first $l-1$ layers of FNN (1).

###### Remark 5

When $p>0$, it is hard to obtain the exact optimal solution to (15). Note that
the output layer poses an additional complication, since $p\neq 0$ for this
layer. In fact, $p=m=1$ for the output layer, since it is a linear layer with
$\alpha=\beta=1$. Therefore, we handle the output layer separately.

####  IV-3 Dealing with the last layer

We now derive the relationship between $\bar{L}_{0}$ and $L_{0}$, which are
the Lipschitz estimates for the first $(l-1)$ layers and entire neural network
including the output layer $\mathbf{L}_{i},i\in\mathbb{Z}_{l}$, respectively.
In fact, $L_{0}$ at the output layer can be computed as

$\displaystyle L_{0}=\bar{L}_{0}\|W_{l}\|_{2}=\sqrt{c}\bar{L}\|W_{l}\|_{2}.$
(17)

The first equality directly follows from the definition of the matrix norm and
the fact that $z_{2}^{(l)}-z_{1}^{(l)}=W_{l}(z_{2}^{(l-1)}-z_{1}^{(l-1)})$.
The second equality is from Theorem 1.

####  IV-4 Compositional Algorithm

We summarize the compositional algorithm to compute a Lipschitz constant for
(1) for the case where $p=0$ in Algorithm 1.

Algorithm 1 Compositional Lipschitz Constant Estimation

Input Weights $\\{W_{i}\\}_{i=1}^{l}$ from a FNN (1) with activation function
slope-restricted in $[0,1]$; constant $c\in\mathbb{R}_{++}$
Output Lipschitz estimate $L_{0}$

1: Set $X^{0}=cI$

2: for $i=1,2,...,l-1$ do

3: Obtain $\tilde{\lambda_{i}}$ from (16)

4: Obtain $X^{i}$ from (12) with with $p=0$ and
$\lambda_{i}=\tilde{\lambda}_{i}$ from Step 3

5: end for

6: Obtain $\bar{L}$ from (14)

7: Return $L_{0}=\sqrt{c}\bar{L}\|W_{l}\|_{2}$

Note that Algorithm (1) utilizes the exact solution from (16). For the case
where $p>0$, we can solve (15) numerically for all $l$ layers to obtain the
Lipschitz estimate.

##  V Experiments

In this section, we will implement our compositional Lipschitz constant
estimation procedure in Algorithm 1 on several FNNs formulated as (1), and
compare the performance to both the trivial upper bound and the state-of-the-
art Lipschitz constant estimation approach LipSDP-Layer [19]. The neural
networks on which we experiment are chosen to have $l\in\\{4,8,12,16,20\\}$
layers. The number of neurons in the hidden layers are selected to be the same
for each neural network, and are chosen as $d_{i}\in\\{4,8,16,32\\}$. We
select the constant $c$ to be 100 for all FNNs. Note that all the neural
networks are set to have inputs $z^{(0)}\in\mathbb{R}^{4}$ and outputs
$z^{(l)}\in\mathbb{R}$. All experiments are implemented locally on a 12-core
CPU with 16GB of RAM using MATLAB. The LipSDP-Layer approach is implemented
using YALMIP/Penlab [25][26], while our algorithm does not require any package
thanks to the availability of closed-form solutions as described in Section
III. We will demonstrate the performance of our algorithm in terms of the
computation time and the tightness of the Lipschitz estimates for two
experimental cases.

#### Case 1 - Varying neural network depth (number of layers)

We implement Algorithm 1 to estimate the Lipschitz constants for FNNs with
number of layers $l=4,8,12,16,20$ respectively. The number of neurons
$d_{i}=16$ for the layers $\mathbf{L}_{i},i\in\mathbb{Z}_{l-1}$. We summarize
the results on the Lipschitz bound in Fig. 1. We can see that as the number of
layers increases, the computation time for the LipSDP-Layer algorithm
explodes. In contrast, the computation time for our algorithm stays low and
grows only linearly with respect to the number of layers, indicating that our
algorithm is highly scalable. Furthermore, the Lipschitz estimate given by our
algorithm stays close to the one from LipSDP-Layer while significantly
outperforming the trivial bound (Fig. 1). The trade-off is very worthwhile in
that our algorithm leads to significant improvements in computation time with
only a slight sacrifice in the tightness of the Lipschitz estimate.

Figure 1: Case 1 \- Lipschitz constant estimates (left) and computation time
(right) with respect to depth for an FNN with 16 neurons per hidden layer.

#### Case 2 - Varying neural network width (number of neurons per layer)

We now examine the performance of our compositional approach for wider (more
hidden neurons per layer), rather than deeper neural networks (with more
layers). This corresponds to practical situations where a large neural network
is essential to capture complicated nonlinear behavior of a system. We select
a deep neural network with 20 layers. Then, we execute our experiments where
the number of neurons for the layers $\mathbf{L}_{i},i\in\mathbb{Z}_{l-1}$ are
set to $d_{i}=4,8,16,32$. We can observe from Fig. 2 that the computation time
needed for LipSDP-Layer also explodes as the number of neurons increases,
while the computation time of our method still grows linearly. Meanwhile, the
Lipschitz estimates of our method are only slightly larger than the one given
by LipSDP-Layer. Thus, we can conclude that our method also outperforms
standard centralized approaches for wider neural networks.

Figure 2: Case 2 \- Lipschitz constant estimates (left) and computation time
(right) with respect to the number of hidden neurons per layer for an FNN with
20 layers.

##  VI Conclusions and Future Work

We offer a compositional method to estimate a tight upper bound for the
Lipschitz constant on a FNN by exploiting the underlying structure of the
large but sparse verification matrix to break down the corresponding
certificates into multiple smaller sub-problems. Our approach provides tight
Lipschitz bounds with a steep reduction in computation time compared to state-
of-the-art techniques, which provides a significant advantage in applications
where neural networks are utilized as part of control loops and need to be
certified during real-time operation. Future directions of interest include
extension to less commonly used activation functions (such as LeakyReLU) and
other neural network structures such as convolutional neural networks.

##  VII Appendix

We state the proofs of all the results in the paper here.

Proof of Theorem 1: Given any two inputs $z_{1}^{(0)}$ and $z_{2}^{(0)}$, let
$z^{(i)}_{1},z^{(i)}_{2},v^{(i)}_{1},v^{(i)}_{2},i\in\mathbb{Z}_{l}$ be
calculated as in (1). We denote $\Delta z^{(i)}=z^{(i)}_{1}-z^{(i)}_{2}$,
$i\in\\{0\\}\cup\mathbb{Z}_{l}$ and $\Delta v^{(j)}=v^{(j)}_{1}-v^{(j)}_{2}$,
$j\in\mathbb{Z}_{l}$. We both left and right multiply the left hand side in
(5) by the vector $[\Delta z^{(0)},\Delta z^{(1)},...,\Delta z^{(l)}]^{T}$.
Then, we have

$\displaystyle\sum_{i=1}^{l}\left(\begin{bmatrix}\Delta z^{(i-1)}\\\ \Delta
z^{(i)}\end{bmatrix}^{T}\begin{bmatrix}\lambda_{i}pW_{i}^{T}W_{i}&-\lambda_{i}mW_{i}^{T}\\\
-\lambda_{i}mW_{i}&\lambda_{i}I\end{bmatrix}\begin{bmatrix}\Delta z^{(i-1)}\\\
\Delta z^{(i)}\end{bmatrix}\right)$ (18) $\displaystyle+c(\Delta
z^{(0)})^{2}-\frac{1}{L^{2}}(\Delta z^{(l)})^{2}>0.$

Now we show that every summand in the upper inequality is negative
semidefinite. In fact, with Assumption 2, using notation $\Delta
v^{(i)}=v^{(i)}_{1}-v^{(i)}_{2}$ and
$\Delta\phi(v^{(i)})=\phi(v_{1}^{(i)})-\phi(v_{2}^{(i)})$,
$i\in\mathbb{Z}_{l}$, we can write

$\begin{bmatrix}\Delta v^{(i)}\\\
\Delta\phi(v^{(i)})\end{bmatrix}^{T}\begin{bmatrix}pI&-mI\\\
-mI&I\end{bmatrix}\begin{bmatrix}\Delta v^{(i)}\\\
\Delta\phi(v^{(i)})\end{bmatrix}\leq 0.$ (19)

Note that $\Delta
v^{(i)}=(W_{i}z_{1}^{(i-1)}+b_{i})-(W_{i}z_{2}^{(i-1)}+b_{i})=W_{i}\Delta
z^{(i-1)}$ and $\Delta\phi(v^{i})=\Delta z^{i}$. We can express these
relationships in matrix form as

$\begin{bmatrix}\Delta v^{(i)}\\\
\Delta\phi(v^{(i)})\end{bmatrix}=\begin{bmatrix}W_{i}&0\\\
0&I\end{bmatrix}\begin{bmatrix}\Delta z^{(i-1)}\\\ \Delta
z^{(i)}\end{bmatrix}.$ (20)

Substituting (20) into (19) and multiplying both sides by
$\lambda_{i}\in\mathbb{R}_{++}$, we finally have

$\begin{bmatrix}\Delta z^{(i-1)}\\\ \Delta
z^{(i)}\end{bmatrix}^{T}\begin{bmatrix}\lambda_{i}pW_{i}^{T}W_{i}&-\lambda_{i}mW_{i}^{T}\\\
-\lambda_{i}mW_{i}&\lambda_{i}I\end{bmatrix}\begin{bmatrix}\Delta z^{(i-1)}\\\
\Delta z^{(i)}\end{bmatrix}\leq 0.$ (21)

Combining (LABEL:ineq:_left_and_right_multiply_large_matrix) and (21), we have

$c(\Delta z^{(0)})^{2}-\frac{1}{L^{2}}(\Delta z^{(l)})^{2}>0,$

which gives a upper bound $\sqrt{c}L$ for the Lipschitz constant.

Proof of Theorem 2: We first restate Lemma 2 in [24] here, which will be
crucial in proving this result. A symmetric block tri-diagonal matrix defined
as

$\begin{bmatrix}\mathcal{P}_{0}&\mathcal{R}_{1}&0&...&&0\\\
\mathcal{R}_{1}^{T}&\mathcal{P}_{1}&\mathcal{R}_{2}&...&&0\\\
0&\mathcal{R}_{2}^{T}&\mathcal{P}_{2}&\mathcal{R}_{3}&...&0\\\ &&\vdots&&&\\\
0&...&0&\mathcal{R}_{l-1}^{T}&\mathcal{P}_{l-1}&\mathcal{R}_{l}\\\
0&...&&0&\mathcal{R}_{l}^{T}&\mathcal{P}_{l}\\\ \end{bmatrix},$ (22)

is positive definite if and only if

$M_{i}>0,\quad\forall i\in\\{0\\}\cup\mathbb{Z}_{l},$ (23)

where

$\displaystyle M_{i}=\begin{cases}\mathcal{P}_{i}&\text{if }i=0,\\\
\mathcal{P}_{i}-\mathcal{R}_{i}^{T}M_{i-1}^{-1}\mathcal{R}_{i}&\text{if}\
i\in\mathbb{Z}_{l}.\\\ \end{cases}$ (24)

Now, we define

$\mathcal{R}_{i}=-\lambda_{i}mW_{i}^{T},i\in\mathbb{Z}_{l},$ (25)

and

$\displaystyle\mathcal{P}_{i}=\begin{cases}cI+\lambda_{1}pW_{1}^{T}W_{1}&\text{if
}i=0,\\\ \lambda_{i}I+\lambda_{i+1}pW_{i+1}^{T}W_{i+1}&\text{if}\
i\in\mathbb{Z}_{l-1},\\\ \lambda_{i}I&\text{if }i=l.\end{cases}$ (26)

Then, with $\mathcal{P}_{i}$ and $\mathcal{R}_{i}$ defined above, we directly
have $(M_{i})_{l}=M_{i}$, where $(M_{i})_{l}$ is defined in (8). In other
words, $(M_{i})_{l}>0,i\in\\{0\\}\cup\mathbb{Z}_{l}$ is a sufficient and
necessary condition for (5).

Proof of Theorem 3: With $\lambda_{i}\in\mathbb{R}_{++},i\in\mathbb{Z}_{l}$,
$p\geq 0$ and $W_{i}^{T}W_{i}\geq 0$, we know that
$\lambda_{i}pW_{i}^{T}W_{i}\geq 0$. Then, according to (11),
$X^{i}\geq(M_{i})_{l-1},i\in\mathbb{Z}_{l-2}$. When $i=l-1$,
$X^{l-1}=(M_{l-1})_{l-1}$. When $i=0$, with $c\in\mathbb{R}_{++}$ and
$\lambda_{1}pW_{1}^{T}W_{1}\geq 0$, $(M_{0})_{l-1}>0$ is automatically
satisfied. Gathering all these cases, (13) implies (8).

## References

  * [1] J. Sarangapani, _Neural network control of nonlinear discrete-time systems_. CRC press, 2018.
  * [2] A. Aswani, H. Gonzalez, S. S. Sastry, and C. Tomlin, “Provably safe and robust learning-based model predictive control,” _Automatica_ , vol. 49, no. 5, pp. 1216–1226, 2013.
  * [3] S. Kuutti, R. Bowden, Y. Jin, P. Barber, and S. Fallah, “A survey of deep learning applications to autonomous vehicle control,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 2, pp. 712–733, 2020.
  * [4] D. Nguyen-Tuong and J. Peters, “Model learning for robot control: a survey,” _Cognitive processing_ , vol. 12, pp. 319–340, 2011.
  * [5] M. Q. Raza and A. Khosravi, “A review on artificial intelligence based load demand forecasting techniques for smart grid and buildings,” _Renewable and Sustainable Energy Reviews_ , vol. 50, pp. 1352–1372, 2015.
  * [6] Z. Afroz, G. Shafiullah, T. Urmee, and G. Higgins, “Modeling techniques used in building hvac control systems: A review,” _Renewable and sustainable energy reviews_ , vol. 83, pp. 64–84, 2018.
  * [7] E. Wong and Z. Kolter, “Provable defenses against adversarial examples via the convex outer adversarial polytope,” in _International conference on machine learning_. PMLR, 2018, pp. 5286–5295.
  * [8] A. Raghunathan, J. Steinhardt, and P. Liang, “Certified defenses against adversarial examples,” _arXiv preprint arXiv:1801.09344_ , 2018.
  * [9] P. L. Bartlett, D. J. Foster, and M. J. Telgarsky, “Spectrally-normalized margin bounds for neural networks,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [10] H. Gouk, E. Frank, B. Pfahringer, and M. J. Cree, “Regularisation of neural networks by enforcing lipschitz continuity,” _Machine Learning_ , vol. 110, pp. 393–416, 2021.
  * [11] Y. Tsuzuku, I. Sato, and M. Sugiyama, “Lipschitz-margin training: Scalable certification of perturbation invariance for deep neural networks,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [12] P. Pauli, N. Funcke, D. Gramlich, M. A. Msalmi, and F. Allgöwer, “Neural network training under semidefinite constraints,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 2731–2736.
  * [13] P. Pauli, A. Koch, J. Berberich, P. Kohler, and F. Allgöwer, “Training robust neural networks using lipschitz bounds,” _IEEE Control Systems Letters_ , vol. 6, pp. 121–126, 2021.
  * [14] Y. Xu and S. Sivaranjani, “Learning dissipative neural dynamical systems,” _IEEE Control Systems Letters_ , vol. 7, pp. 3531–3536, 2023.
  * [15] X. Huang, D. Kroening, W. Ruan, J. Sharp, Y. Sun, E. Thamo, M. Wu, and X. Yi, “A survey of safety and trustworthiness of deep neural networks: Verification, testing, adversarial attack and defence, and interpretability,” _Computer Science Review_ , vol. 37, p. 100270, 2020.
  * [16] L. Brunke, M. Greeff, A. W. Hall, Z. Yuan, S. Zhou, J. Panerati, and A. P. Schoellig, “Safe learning in robotics: From learning-based control to safe reinforcement learning,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 5, pp. 411–444, 2022.
  * [17] S. Dean, N. Matni, B. Recht, and V. Ye, “Robust guarantees for perception-based control,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 350–360.
  * [18] F. Latorre, P. Rolland, and V. Cevher, “Lipschitz constant estimation of neural networks via sparse polynomial optimization,” _arXiv preprint arXiv:2004.08688_ , 2020.
  * [19] M. Fazlyab, A. Robey, H. Hassani, M. Morari, and G. Pappas, “Efficient and accurate estimation of lipschitz constants for deep neural networks,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [20] A. Xue, L. Lindemann, A. Robey, H. Hassani, G. J. Pappas, and R. Alur, “Chordal sparsity for lipschitz constant estimation of deep neural networks,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 3389–3396.
  * [21] P. Pauli, D. Gramlich, and F. Allgöwer, “Lipschitz constant estimation for 1d convolutional neural networks,” in _Learning for Dynamics and Control Conference_. PMLR, 2023, pp. 1321–1332.
  * [22] M. Newton and A. Papachristodoulou, “Exploiting sparsity for neural network verification,” in _Learning for Dynamics and Control_. PMLR, 2021, pp. 715–727.
  * [23] M. Fazlyab, M. Morari, and G. J. Pappas, “Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 1, pp. 1–15, 2020.
  * [24] E. Agarwal, S. Sivaranjani, V. Gupta, and P. Antsaklis, “Sequential synthesis of distributed controllers for cascade interconnected systems,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 5816–5821.
  * [25] J. Lofberg, “Yalmip: A toolbox for modeling and optimization in matlab,” in _2004 IEEE international conference on robotics and automation (IEEE Cat. No. 04CH37508)_. IEEE, 2004, pp. 284–289.
  * [26] J. Fiala, M. Kočvara, and M. Stingl, “Penlab: A matlab solver for nonlinear semidefinite optimization,” _arXiv preprint arXiv:1311.5240_ , 2013.