# Fractional damping induces resonant behavior in the Duffing oscillator

Mattia Coccolo  Nonlinear Dynamics, Chaos and Complex Systems Group,
Departamento de Física, Universidad Rey Juan Carlos, Tulipán s/n, 28933
Móstoles (Madrid), Spain  Jesús M. Seoane  Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Tulipán s/n, 28933 Móstoles (Madrid), Spain  Miguel A.F. Sanjuán  Nonlinear
Dynamics, Chaos and Complex Systems Group, Departamento de Física, Universidad
Rey Juan Carlos, Tulipán s/n, 28933 Móstoles (Madrid), Spain

###### Abstract

The interaction between the fractional order parameter and the damping
parameter can play a relevant role for introducing different dynamical
behaviors in a physical system. Here, we study the Duffing oscillator with a
fractional damping term. Our findings show that for certain values of the
fractional order parameter, the damping parameter, and the forcing amplitude
high oscillations amplitude can be induced. This phenomenon is due to the
appearance of a resonance in the Duffing oscillator only when the damping term
is fractional.

##  I Introduction

Nonlinear oscillators are dynamical systems that exhibit oscillatory behaviors
but do not follow a linear relationship between the restoring force and the
displacement. These systems often arise in various scientific and engineering
fields, and their study has profound implications in understanding complex
phenomena and designing efficient devices. When considering damping in
nonlinear oscillators, a commonly used approach is to introduce a fractional
derivative in the damping term. The use of fractional derivatives in damping
terms provides a powerful tool to capture the memory effects and non-local
behavior observed in certain physical systems. Unlike the traditional integer-
order derivatives that describe instantaneous rates of change, fractional
derivatives involve a non-integer order and account for past history in the
damping process. This allows us a more accurate representation of damping in
systems with memory and hysteresis as occurs in some properties of complex
materials. The presence of a fractional derivative damping term in nonlinear
oscillators can lead to novel behaviors such as sub-diffusive or super-
diffusive damping, where the decay of oscillations is slower or faster than in
traditional damping scenarios. Moreover, the inclusion of the fractional
derivatives can alter the stability regions, bifurcation boundaries, and
resonance responses of the oscillatory system, introducing rich and intricate
dynamics. The study of nonlinear oscillators with fractional derivative in the
damping term finds applications in disciplines such as mechanical engineering,
physics, robotics, control systems, and even social sciences Boroviec ;
Dafermos ; Lu ; Chellaboina ; De ; Elliott ; Horr ; Ding ; Li ; Strogatz .
Understanding the behavior and control of these oscillators is crucial for
optimizing energy dissipation, enhancing stability, and designing systems with
specific response characteristics.

On the other hand, nonlinear resonances Sanjuan  is a fascinating phenomenon
that occurs when a system, subjected to an external force or perturbation,
exhibits a response that is not proportional to the input. Unlike linear
resonance, where the response is directly proportional to the forcing
frequency, nonlinear resonances arise due to the intricate interplay between
multiple frequencies and energy transfer mechanisms within the system, as in
our case the fractional order parameter, the damping parameter and the forcing
amplitude. In nonlinear resonances, the system’s response can display a wide
range of complex behaviors, including the generation of harmonics, sub-
harmonics, chaotic oscillations, and bifurcations. As stated before, the
fractional derivatives are the key to explore and modeling materials with
memory effects and fractal geometry or diffusion processes, among others. In
fact, fractional derivatives can play a crucial role on them. Therefore, it is
relevant to analyze when the fractional order parameter can trigger complex
behaviors and nonlinear resonance phenomena.

Here, we analyze the Duffing oscillator with a fractional derivative in its
damping term as in Ref. Coccolo_fr2 . Specifically, we focus our study in the
numerical analysis of the fractional order parameter effects on the
oscillations amplitude of the system. We have found that the fractional order
parameter can induce resonance phenomena. Firstly, we analyze, for smaller
forcing amplitude values, beyond a certain threshold, the appearance of
resonance peaks that are triggered by the high damping parameter values and
specific values of the fractional order parameter. These high amplitude
oscillations are related with a complex dynamics of the oscillator in the
phase space. Then, we find a region of resonances for higher forcing amplitude
values. Consequently, we characterize and analyze the appearance of resonance
phenomena induced by the effect of the three main ingredients: high forcing
amplitude and damping parameter values, and some specific values of the
fractional order parameter. It is relevant to emphasize that these complex
behaviors and resonance phenomena are an effect of the fractional order
parameter, and they only appear when a fractional derivative damping term
exists.

The organization of this paper is as follows. In Sec. II, we introduce the
fractional Duffing oscillator and its dynamics. In Sec. III, we numerically
present the effects on the oscillations amplitude of the fractional order
parameter in which a resonant-like behavior is uncovered. The impact of the
forcing amplitudes due to the fractional order parameter is analyzed in Sec.
IV, where we numerically characterize the birth of this resonance-like
behavior. Conclusions and a discussion of the main results of this paper are
presented in Sec. V.

##  II Model Description

We study a Duffing oscillator with a fractional derivative of fractional order
parameter $\alpha$ instead of the first derivative, i.e., $\dot{x}(t)\to
D^{\alpha}x(t)=d^{\alpha}x(t)/dt^{\alpha}$ as a damping term. Then, the
Duffing oscillator with a fractional damping term is given by

$\frac{d^{2}x}{dt^{2}}+\mu\frac{d^{\alpha}x}{dt^{\alpha}}-x+x^{3}=F\cos{\omega
t},$ (1)

where $\mu$ is the damping parameter, $F$ and $\omega$ are the forcing
amplitude and forcing frequency, respectively, and $\alpha$ is the fractional
order parameter. In order to integrate the system, we use the Grünwald-
Letnikov fractional derivative, as in previous works Coccolo_fr ; Coccolo_fr2
. Consequently, the system is composed by a set of three fractional
differential equations that yield

$\left\\{\begin{array}[]{l}D^{\alpha}x=y,\vskip 6.0pt plus 2.0pt minus
2.0pt\\\ D^{1-\alpha}y=z,\vskip 6.0pt plus 2.0pt minus 2.0pt\\\ Dz=\mu
y+F\cos{(\omega t)}+x-x^{3},\end{array}\right.$ (2)

where $z$ is an auxiliary component coming from the transformation into a
fractional order system. The integration step used is $\pi/2000$, and the
maximum time value is $t=300$ that is large enough for the system to reach the
steady state. We have carried out numerical simulations with several initial
conditions and, although the results are numerically different, the general
trends of the curves are qualitatively similar. So, we carry on the study with
the initial condition $(x_{0},y_{0},z_{0})=(0.2,0.3,0.0)$, as in Ref. Syta ;
Coccolo_fr2 . We expect that our conclusions are of general validity and not
specific for the considered boundary conditions. In Fig. 1, we show some
examples of trajectories of the system and the impact of the $\alpha$ value on
its dynamics, when the $\mu$ parameter varies. In fact, according to Fig.
1(a), where $\mu=0.6$, and Fig. 1(b), where $\mu=0.9$, we have to change the
$\alpha$ parameter from $\alpha=0.3$ to $\alpha=0.1$, respectively, if we want
to observe similar oscillations amplitudes. In Fig. 1(c), we compare
trajectories for $\mu=0.9$ and $\alpha=0.3$ (in blue) and for $\mu=0.9$ and
$\alpha=0.1$ (in red). Having traced more trajectories than those shown, we
have seen that the phenomenon is complex and worthy of study.

Figure 1: In this figure, we have plotted the asymptotic behaviors of the
system trajectory in phase space for (a) $\mu=0.6$ and $\alpha=0.3$, (b)
$\mu=0.9$ and $\alpha=0.1$ and (c) $\mu=0.9$ and $\alpha=0.1$ (in red) and
$\alpha=0.3$ (in blue). We fix the forcing amplitude $F=1.5$ for all cases.
Note the key role of the $\alpha$ and $\mu$ parameters in the amplitude of the
oscillations.

##  III Fractional order parameter effects on the oscillations amplitude

As we previously mentioned, some preliminary computations about the resonance
induced by the fractional derivative in the Duffing damping term have been
done in Coccolo_fr2 . The authors have disclosed that for a high value of the
damping parameter $\mu=0.8$ and for distinctive values of the fractional order
parameter $\alpha$, some interesting peaks pop up when the forcing amplitude
is large enough, as shown in Fig. 2. Also, for higher $\alpha$ values we have
seen that the oscillations amplitude do not show those kind of peaks and that
the curve stays smooth under an oscillations amplitude the of value $0.5$. We
can see that two particular peaks stand out in the amplitude and the
$Q-$factor plots, Fig. 2(a) and Fig. 2(c), respectively. However, this fact
was just mentioned, though not analyzed in detail. As stated before, the
$Q-$factor provides an idea of how much the signal is amplified by a certain
parameter, in our case $\alpha$. We calculate it by computing the sine and
cosine components:

$B_{s}=\frac{2}{nT}\int^{nT}_{0}{x(t)\sin{\omega t}dt}\quad\text{and}\quad
B_{c}=\frac{2}{nT}\int^{nT}_{0}{x(t)\cos{\omega t}dt},$ (3)

where $T=2\pi/\omega$ and $n$ is an integer. Then, we can find the dependence
on $\alpha$ with

$Q=\frac{\sqrt{B^{2}_{s}+B^{2}_{c}}}{\alpha}.$ (4)

Those peaks are related with fractional order parameter values that suddenly
give birth to aperiodic orbits, as we can see in Fig. 2(b). They are
resonance-like peaks for the $\alpha$ parameter, as confirmed by the
$Q-$factor, but it is worth of study the conditions for them to appear and the
parameter values for which the gain in amplitude is higher. Our objective now
is to analyze in depth this resonant behavior. For this purpose, we plot the
gradient of the amplitude in the parameter set $\alpha-\mu$ for two values of
the forcing amplitude $F=0.01$ and $F=0.1$ in Fig. 3. In the first case, Fig.
3(a), we can see that the gradient in the parameter set is smooth, the
oscillations amplitude is not that high and the higher oscillations amplitude
zone is very easy to spot. The only interesting thing is the intermediate
oscillations amplitude smudge for $\mu\gtrsim 0.6$. Even more interesting are
the other two panels where Figs. 3(c) is a zoom of Fig. 3(b). When we reach
the forcing amplitude value $F=0.1$, the parameter set becomes more
complicated and the high and low oscillations amplitude intermingle.

From values $\mu\gtrsim 0.4$, the structure that can be spotted in Fig. 2(a)
and that ends at $\alpha\approx 0.15$, starts to take form with lots of the
resonance peaks, in which the oscillations span both potential wells. In fact,
in Fig. 4 we show the oscillations amplitude plots for different $\mu$ values
as a section of the $3-D$ figure. There, we can see that for $\mu=0.35$ there
is no sign of the structure, while it is visible for $\mu=0.55$. Clearly, for
higher $\mu$ values it becomes even more visible. In the $3-D$ figure, it is
possible to discern the built-up of such structure and how the higher
oscillations amplitude peaks diminish when $\mu$ grows. This structure is the
result of the higher oscillations amplitude being attenuated by the
dissipation term, so that the peaks are the consequence of the fractional
order parameter feat interacting with the forcing amplitude. Hence, it seems
that they are fractional-induced amplitude peaks that behave like a resonance
when only the fractional order parameter varies. In fact, if we compute the
$Q-$factor for the parameter $\mu$ the peaks do not show up in the plot.

In Fig. 4(e), we can see that the trajectories become bounded to one well when
the dissipation is too strong, $\mu=1$, for the fractional order parameter to
maintain the interwell oscillations. Figure 3(c) has been plotted to have a
better understanding of the higher amplitude oscillations zone. In particular
of the cross like patch centered at $\mu=0.6$. Here, we can find a region of
small oscillations amplitude embedded inside a high oscillations amplitude
region. It is precisely there, where the already mentioned structure, well
recognisable in Fig. 2(a), starts to take form.

Figure 2: We plot here (a) the oscillations amplitude, (b) the asymptotic
behaviors diagram and (c) the $Q-$factor for $F=0.1,\mu=0.8$. The two
resonance peaks can be recognise in the panels. Figure 3: The figure shows the
oscillations amplitude for $F=0.01$ (a) and $F=0.1$ (b). Panel (c) is a zoom
of panel (b). Here it is possible to appreciate how the difference in the
amplitude forcing complicates the topology of the parameter set.  Figure 4: In
panel (a), we show the $3D$ version of Fig. 3(b). The other panels display $4$
slices of the panel (a) along the $\alpha$ axis and for (b) $\mu=0.35$, (c)
$\mu=0.55$, (d) $\mu=0.76$, and (e) $\mu=1$. Here, the formation of the
central structure is clear. In fact, in panel (c) it is possible to see how
the high oscillations amplitude of panel (b) decreases and start to form the
central structure of panels (d) and (e). Figure 5: We show the oscillations
amplitude in the parameter set $F-\alpha$, for different $\mu$ values. Here,
it appears along the panels the formation of a well located region of high
oscillations amplitude near $F=1.5$ Figure 6:  The figure shows four slices of
Fig. 5(b) overlapped in panel (a) and the oscillations frequencies calculated
with the Fast Fourier transform in panel (b). In the first panel, the curves
related to fractional order values, $\alpha=0.1$ and $\alpha=0.3$, that cross
the high amplitude oscillations region are recognisable. The other two, after
they grow at $F\sim 0.5$, become smooth. The frequencies in panel (b) for
$F\gtrsim 0.9$ tend to the forcing frequency $\omega=1$, except for the ones
related with the resonance peak.  Figure 7:  We show the behavior of the
system for fixed values of $\alpha$ and $\mu$. Comparing panel (a) and (b), we
can observe that the system behaviors at $F=1.5$ are different. For
$\alpha=0.1$, in panel (e), the oscillations are clearly not periodic, while
for $\alpha=0.3$, in panel (f), the oscillations are periodic of period $3$.
In the other two panels, the higher oscillations amplitude enlargement is
lacking.

##  IV The forcing amplitude effect on the oscillations amplitude

To carry on our analysis, we have decided to study the effect of the forcing
amplitude $F$ on the oscillations amplitude since it will provide more easily
clear light on the resonance-like phenomenon. The main result of this
manuscript is shown in Fig. 5, where we have depicted the gradient parameter
set $F-\alpha$ of the oscillations amplitude for different values of $\mu$. As
expected, the higher the damping term the lower the oscillations amplitude.
Interestingly, a high amplitude oscillations region starts to appear when
$\mu>0.5$, in particular it becomes more evident at $\mu\approx 0.6$. This
region is centered around $F=1.5$ and spans all or part of the parameter set,
depending on the damping parameter. In fact, we can see there that all along
the panels of Fig. 5 that region shrinks along the $\alpha$ axis together with
the oscillations amplitude. That means that for a high enough damping
parameter value, the high oscillations amplitude is confined around a certain
forcing amplitude value, that is, $F=1.5$ and it is induced by a certain range
of fractional order parameter values. These results corroborate what was
observed in Fig. 1.

Besides, the region shrinks spanning smaller values of the fractional order
parameter when the dissipation grows. This phenomenon can be explained as in
the previous section. The interaction among high dissipation, certain values
of the fractional order parameter and high forcing amplitude triggers a
resonant behavior. The interesting fact is that the gain in amplitude is
confined to certain $\alpha$ values, so we can define the phenomenon as a
resonance induced by the fractional order parameter. The difference with
$F=0.1$ is that now the forcing amplitude is large enough to sustain the high
oscillations amplitude inside a large region of the parameter set and not just
for some particular values of $\alpha$.

For a better understanding and the characterization of the phenomenon, we have
deepened its analysis with the other figures below. We want to stress out that
all the figures from now on are related with Fig. 5 and they provide more
specific details of the features of the resonance phenomenon. Therefore, we
start by showing in Fig. 6(a) the oscillations amplitude versus the forcing
amplitude $F$ for $\mu=0.6$. Every curve is associated with a fractional order
parameter value as indicated there. It is possible to see that the amplitude
curves corresponding with the two values $\alpha=0.1$ and $\alpha=0.3$, for
this damping parameter, cross the high oscillations amplitude region of Fig.
5(b) and present the corresponding peaks, near $F=1.5$. On the other hand, the
curves related with the other two $\alpha$ values are smooth after a first
jump, at $\alpha\approx 0.5$. In Fig. 6(b), we show the frequencies of the
oscillations, calculated with the Fast Fourier transform. It is possible to
see that from $F=0.8$ the calculated frequencies tend naturally to the forcing
frequency $\omega=1$, but not near the forcing amplitude for which the
resonance peaks are produced for the $\alpha$ values that cross the high
amplitude region of Fig. 5(b).

In order to have a deeper insight on the dynamics of the system, we plot Fig.
7 in which we show the maxima-minima diagram for $\mu=0.6$ and for $4$ values
of the $\alpha$ parameter, along with two orbits in the case of $F=1.5$ and
$\alpha=0.1$ (Fig. 7(e)) and $\alpha=0.3$ (Figs. 7(f)). These fractional order
parameter values are the ones that cross the high oscillations amplitude in
Fig. 5(b). The other two values, see Figs. 7(c) and (d), do not cross that
region as indicated by the maxima-minima diagram smoothness. One consideration
can be made, we can see that there is an important difference between the two
plotted orbits, as seen in Figs. 7(e) and (f). The first one is not periodic,
while the second one is periodic of period $3$. It seems that when the
fractional order parameter grows, the large amplitude aperiodic trajectories
disappear becoming periodic as if the fractional order parameter contributed
to introduce order into chaos.

Now, we check the robustness of the resonance induced by the fractional order
parameter when the forcing frequency $\omega$ changes. For that purpose, we
plot the $Q-$factor in Fig. 8, which shows the oscillations amplitude and the
maxima-minima diagram in function of the forcing frequency and for $F=1.5$ and
$\mu=0.6$. In Figs. 8(a) and (b) it is possible to notice some peaks in the
$Q-$factor curve that assure us that they are related with the fractional
order parameter induced resonance. This resonance phenomenon needs an
interaction among high forcing amplitude, high damping parameter and small
values of the fractional order parameter. In Fig. 8(a), we have highlighted
three peaks at $\alpha=0.54,\alpha=1$ and $\alpha=2.01$ that will be discussed
later when we analyze higher $\mu$ values. The peaks are corroborated by the
maxima-minima diagram, as seen in Fig. 8(c). In this last figure, it is
interesting to notice that the higher peaks are related with aperiodic
trajectories.

Figure 8: We show the $Q-$factor (a), the oscillations amplitude (b) and the
maxima-minima diagram in function of the forcing frequency $\omega$ for the
following parameters: $\mu=0.6$, $\alpha=0.1$ and $F=1.5$. In the first panel,
three peaks are emphasized with three vertical lines in order to compare them
with similar peaks in other forthcoming figures. The peaks in the panels are
the result of the equilibrium reached by the high oscillations suppression due
to the damping term, the effect of the forcing amplitude and the action of the
fractional order parameter. Figure 9: Panel (a) shows four slices of Fig. 5(e)
overlapped. Panel (b) shows the oscillations frequencies calculated with the
fast Fourier transform. In panel (a), it is possible to see that the curve for
$\alpha=0.1$ is the only one showing the peak due to the high oscillations
amplitude region of Fig. 7. Also here, the oscillations frequencies tend to
the forcing frequency $\omega=1$, except for the ones related with the
oscillations amplitude peak.  Figure 10: The panels show the behavior of the
system for fixed $\alpha$ and $\mu$. Analyzing panel (a), we can see that the
system behavior at $F=1.5$ is periodic of period $3$, as shown in panel (e),
unlike in Fig. 7(e). In the other three panels, the higher oscillations
amplitude enlargement is lacking. In fact, in panel (f) we compare the orbits
for $\alpha=0.1$ and $\alpha=0.3$ and we can appreciate the red trajectory
broadening with respect to the blue orbit, due to the resonance effect. Figure
11: The panels show for $\mu=0.9$, $\alpha=0.1$ and $F=1.5$ (a) the
$Q-$factor, (b) the oscillations amplitude and (c) the maxima-minima diagram
in function of the forcing frequency $\omega$. Panels (d), (e) and (f) show
respectively the same plots for $\mu=1$. In panels (a) and (d) the already
commented three peaks are emphasized with three vertical lines in order to
compare them with similar peaks in Fig. 8. Those peaks can be seen more clear,
and others, that can be seen in Fig. 8, are disappearing as an effect of the
dissipation growth. Also they are centered to certain values of the forcing
frequency $\omega=0.5,\omega=1$ and $\omega=2$.

Now, we analyze the situation in which the damping parameter has a large value
and, therefore, we are in the overdamped regime. In Fig. 9, we can find a
similar plot as Fig. 6 for the damping parameter value $\mu=0.9$. Here, in
Fig. 9(a), we can see that only for $\alpha=0.1$ there is a resonance peak, as
expected from Fig. 5(e). The computed oscillations frequencies, Fig. 9(b),
behave as in the previous case. This result is confirmed by Fig. 10, where we
can see that the resonance amplitude can be found only for $\alpha=0.1$, Fig.
10(a), while for the others fractional values the maxima-minima diagram is
smooth for large forcing amplitudes, Figs. 10(b)-(d). Interestingly, the high
oscillations amplitude is related with periodic trajectories, as we can see in
Fig. 10(e), differently from what was happening in Fig. 7(e). We also compare
the trajectories for $F=1.5$ at $\alpha=0.1$ and at $\alpha=0.3$ in Fig. 10(f)
and we can appreciate the gain in amplitude, due to the resonance. Then, we
plot Fig. 11, by varying the forcing frequency $\omega$. Here, we can spot the
three peaks already mentioned before, as depicted in Fig. 8(a). Again, the
interaction between high damping parameter, forcing and fractional order
parameter induces the resonance. In this figure, we can see that the
$Q-$factor curves, Fig. 11(a) and Fig. 11(d), become smoother and that the
three peaks at $\omega=0.5,\omega=1$ and $\omega=2$ are more recognizable. Of
course, those peaks are related with higher oscillations amplitude, as shown
in Fig. 11(b) and Fig. 11(e). Then, in Fig. 11(c) and Fig. 11(f), we can
appreciate that the peaks match regions of non periodic dynamics of the
system. Contrarily to the $\mu=0.6$ situation, in the $\mu=0.9$ and $\mu=1$
cases, the resonance peaks are centered to specific values of the forcing
frequency $\omega$, related among them as $\omega_{0}=1$,
$\omega_{1/2}=\omega_{0}/2$ and $\omega_{2}=2\omega_{0}$. So, we obtained the
resonance frequency and its harmonic and subharmonic. Also, in Fig. 11(c) and
Fig. 11(f), we can see that the dissipation keeps doing its job by bringing
more order into chaos.

Finally, we present the oscillations amplitude gradient for $F=1$ in the
parameter set $\mu-\alpha$, Fig. 12(a) and for $F=2$ in Fig. 12(b). As
expected, in the both cases, the high oscillations amplitudes are common in
all the parameter set, although the higher ones are concentrated at small
$\mu$ and $\alpha$ values, like in Fig. 3(b). Then, in the first case,
$F=1.5$, the lower amplitude oscillations can be found in the exact opposite
corner of the figure and scattered all around its upper center. Otherwise, in
the second case, $F=2$, as the damping parameter reaches the value $\mu=0.6$,
the oscillations amplitude diminishes abruptly and for the parameter values
for which the resonance is triggered the oscillations amplitude becomes
smaller than in the $F=1.5$ case. This outcome is corroborated by Fig. 12(c),
where the sign of the difference between the oscillations amplitude at $F=2$
and $F=1.5$ is shown. The yellow points indicate a positive difference
$A_{F=2}>A_{F=1.5}$, the blue ones a negative difference $A_{F=1.5}>A_{F=2}$.
The points in this plot are clearly intermingled, but the blue points generate
a more defined region for higher $\mu$ values. This figure shows well that,
although the forcing amplitude is bigger in the $F=2$ case, the interaction
between the damping parameter and the fractional order parameter can trigger
trajectories with a higher oscillations amplitude, giving birth to a resonance
for $F=1.5$. Then, in the other panels, we show three slices of the two
gradient plots, shown in Fig. 12(a) and Fig. 12(b). Here, we can see that, for
small damping parameter values ($\mu=0.18$) the oscillations amplitude for
$F=2$ is generally comparable with the oscillations amplitude for $F=1.5$.
When the damping parameter grows, $\mu=0.5$ the oscillations amplitude for
$F=1.5$ are a little bit higher than the $F=2$ case, but not for all the
$\alpha$ values. Then, when the damping parameter grows further, $\mu=0.8$,
the oscillations amplitude for $F=1.5$ is definitively higher. Figures 12(d-f)
show the tendency of the oscillations amplitude for both forcing amplitudes.
All of this corroborate the results shown in Fig. 12(c). Again, at a high
damping parameter value the $F=1.5$ case oscillates at higher amplitudes due
to the resonance induced by the interaction among the high damping parameter,
the fractional order parameter and the forcing amplitude.

Figure 12:  The figure shows (a) the gradient of the oscillations amplitude in
the parameter set $\alpha-\mu$ for $F=1$, (b) the same gradient plot for
$F=2$, and (c) the sign of the difference $A_{F=2}-A_{F=1.5}$. In panel (c)
the yellow points show the positive sign of the difference, while the blue
points represent the negative sign. The panels (d-f) show the oscillations
amplitude at different values of the damping parameter $\mu$ and for the
forcing amplitudes: $F=2$ in blue and $F=1.5$ in red. We can see that as the
damping parameter grows, the oscillations amplitude for $F=1.5$ becomes higher
than in the case of $F=2$.

##  V Conclusions and discussion

We have studied the fractional Duffing oscillator with a damping fractional
term, using the Grünwald-Letnikov integrator. For the parameter values used,
we have find that the fractional order parameter can induce a resonance. When
the forcing is relatively small, at $F=0.1$, the peaks related with the
resonance appear for particular values of the fractional order parameter,
values that change in function of the damping parameter. For larger values of
the forcing amplitude, say $F=1.5$, the higher oscillations amplitude due to
the resonance form a region in the parameter set ($F-\alpha$). In both cases
the phenomenon is clear, the high oscillations amplitude, that the oscillator
exhibits at low damping parameter values, are attenuated at higher dissipation
but the interaction between a particular forcing amplitude and certain low
fractional order parameter values maintain them, generating a resonance. We
think that this phenomenon can be relevant to better understand the relation
between the damping parameter and the fractional order parameter. Besides, it
shows the richness in the dynamical behavior due to the presence of a
fractional term in the nonlinear oscillator. This last point can easily be
observed in the overdamped case, where the dynamics of the system becomes
still complex. Finally, we expect that this study can be useful to predict the
birth of complex or resonant behaviors in viscoelastic and porous materials or
diffusion processes, that can be only modeled by using fractional derivatives.

##  VI Acknowledgments

This work has been supported by the Spanish State Research Agency (AEI) and
the European Regional Development Fund (ERDF, EU) under Project No.
PID2019-105554GB-I00 (MCIN/AEI/10.13039/501100011033).

## References

  * (1) Borowiec M, Litak G, Syta A. Vibration of the Duffing oscillator: effect of fractional damping. Shock Vib. 2007; 14: 29-36.
  * (2) Dafermos C. Dissipation in materials with memory. In: Lodge AS, Renardy M, Nohel JA, editors. Viscoelasticity and Rheology. Proceedings of a Symposium Conducted by the Mathematics Research Center, Wisconsin–Madison: Academic Press; 1985, p. 221-234.
  * (3) Lu Z, Wang Z, Zhou Y, Lu X. Nonlinear dissipative devices in structural vibration control: A review. J. Sound Vib. 2018; 423: 18-49.
  * (4) Chellaboina V, Haddad WM. Exponentially dissipative nonlinear dynamical systems: a nonlinear extension of strict positive realness. In: Jeltsema D, Scherpen J, editors. Proceedings of the 2000 American Control Conference, Chicago: IEEE Xplore; 2000, p. 3123-3127.
  * (5) De S, Kunal K, Aluru NR. Nonlinear intrinsic dissipation in single layer MoS 2 resonators. RSC Adv. 2017; 7: 6403-6410.
  * (6) Elliott SJ, Tehrani MG, Langley RS. Nonlinear damping and quasi-linear modelling. Philos. Trans. R. Soc. A 2015; 373: 20140402.
  * (7) Horr AM, Schmidt LC. A fractional-spectral method for vibration of damped space structures. Eng. Struct. 1996; 18: 947-956.
  * (8) Ding Y, Liu X, Chen P, Luo X, Luo Y. Fractional-Order Impedance Control for Robot Manipulator. Fractal fract. 2022; 6: 684.
  * (9) Li X, Liu Z, Li J, Tisdell C. Existence and controllability for nonlinear fractional control systems with damping in Hilbert spaces. Acta Math. Sci. 2019; 39; 229-242.
  * (10) Strogatz SH. Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering. Boca Raton: CRC press; 2018.
  * (11) Rajasekar S, Sanjuán MAF. Nonlinear resonances. Switzerland: Springer International Publishing; 2016.
  * (12) Coccolo M, Seoane JM, Lenci S, Sanjuán MAF. Fractional damping effects on the transient dynamics of the Duffing oscillator. Commun. Nonlinear Sci. Numer. Simul. 2023; 117: 106959.
  * (13) Ortiz A, Yang J, Coccolo M, Seoane JM, Sanjuán MAF. Fractional damping enhances chaos in the nonlinear Helmholtz oscillator. Nonlinear Dyn. 2020; 102: 2323-2337.
  * (14) Syta A, Litak G, Lenci S, Scheffler M. Chaotic vibrations of the duffing system with fractional damping. Chaos 2014; 24: 013107.


# Can second-order perturbation theory accurately predict electron density of
open-shell molecules? The importance of self-consistency

Lan Nguyen Tran<EMAIL_ADDRESS>HCMC Institute of Physics, National
Institute of Applied Mechanics and Informatics, Vietnam Academy of Science and
Technology, Ho Chi Minh City 700000, Vietnam

###### Abstract

Electron density distribution plays an essential role in predicting molecular
properties. It is also a simple observable from which machine-learning models
for molecular electronic structure can be derived. In the present work, we
present the performance of the one-body Møller-Plesset second-order
perturbation (OBMP2) theory that we have recently developed. In OBMP2, an
effective one-body Hamiltonian including dynamic correlation at the MP2 level
is derived using the canonical transformation followed by the cumulant
approximation. We evaluate electron density and related properties of three
groups of open-shell systems: atoms and their ions, main-group radicals, and
halogen dimmers. We find that OBMP2 outperforms standard MP2 and density
functional theory in all cases considered here, and its accuracy is comparable
to coupled-cluster singles and doubles (CCSD), a higher-level method. OBMP2 is
thus believed to be an effective method for predicting accurate electron
density of open-shell molecules.

##  I Introduction

Molecular properties such as electrostatic moments (e.g., charges, dipole,
quadrupoles), electrostatic potentials, and electrostatic interaction energies
can be derived from electron density distribution. Information needed to
characterize chemical systems can be thus encoded in electron density Robert
G. Parr (1994); Fabrizio  _et al._ (2019). Particularly, electron density
distribution and its derivatives are among the simplest and generally useful
tools to investigate non-covalent compounds Fabrizio  _et al._ (2019); Kumar
_et al._ (2021); Boto  _et al._ (2020). Recently, several works have developed
density-learning models whose parameters are trained by fitting to reference
electron densities obtained from ab initio quantum chemical calculations
Fabrizio  _et al._ (2019); Cuevas-Zuviría and Pacios (2020). To achieve a good
model, one needs to prepare many high-quality training data. Low-scaling
methods able to generate accurate electron density are thus highly desirable.

Møller-Plesset Second-order perturbation theory (MP2) is the simplest
correlated wavefunction method beyond Hartree-Fock (HF). In MP2, the
correlated density matrix is evaluated using the double-excitation MP2
amplitude that relies on HF molecular orbitals and orbital energies. In many
cases, particularly for open-shell systems where unrestricted Hartree-Fock
(UHF) solutions severely suffer from spin-contamination, correlated MP2
density is not good Kossmann and Neese (2010); Tran (2021). One can relax MP2
density by solving $Z$-vector equation for the orbital response. However, as
shown previously, the singularity of the inverse of the orbital Hessian matrix
causes the violation of $N$-representability of relaxed density, leading to
unphysical discontinuities Kurlancheek and Head-Gordon (2009).

Orbital-optimized MP2 (OOMP2) and its spin-scaled variants have been developed
activelyLochan and Head-Gordon (2007); Neese  _et al._ (2009); Bozkaya  _et
al._ (2011); Bozkaya and Sherrill (2013); Bozkaya (2014a, b); Stück and Head-
Gordon (2013); Sharada  _et al._ (2015); Razban  _et al._ (2017); Lee and
Head-Gordon (2018). In OOMP2, orbitals are optimized by minimizing the
Hylleraas functional or its generalization. OOMP2 can eliminate the spin
contamination present in UHF for open-shell molecules Bozkaya (2014b);
Kossmann and Neese (2010). OOMP2 has been shown to outperform standard MP2 for
many energy-related properties. However, to our knowledge, we are only aware
of a few studies analyzing OOMP2 density-related properties Kurlancheek and
Head-Gordon (2009); Kossmann and Neese (2010); Sandhoefer  _et al._ (2013). In
the DFT community, the accuracy of electron density is a critical criterion
for evaluating if a functional is good Medvedev  _et al._ (2017). A lot of
functionals that give accurate energy-related properties have been found to
have poor electron density Medvedev  _et al._ (2017). Recent studies have
shown that double-hybrid functionals that are based on a mixing of a hybrid
functional and a perturbative second-order correlation part (MP2) can yield an
accurate prediction of both density- and energy-related properties Su  _et
al._ (2018); Hait and Head-Gordon (2018a). Note that the formal scaling of
double-hybrid functionals is the same as MP2. Double-hybrid functionals with
relaxed density are also ill-behaved for stretched bonds, giving unphysical
density-related properties Peverati and Head-Gordon (2013); Hait and Head-
Gordon (2018a, b). To resolve these issues, one can optimize orbitals of
double-hybrid functionals like OOMP2 Peverati and Head-Gordon (2013); Sancho-
Garcia  _et al._ (2016).

Recently we have developed a new self-consistent perturbation theory named
one-body MP2 (OBMP2)Tran and Yanai (2013); Tran (2021). The central idea of
OBMP2 is the use of canonical transformationYanai and Chan (2006, 2007); Chan
and Yanai (2007); Neuscamman  _et al._ (2009, 2010a, 2010b) followed by the
cumulant approximationKutzelnigg and Mukherjee (1997); Mazziotti (1998a, b);
Kutzelnigg and Mukherjee (1999) to derive an effective one-body Hamiltonian.
The resulting OBMP2 Hamiltonian, as a correlated Fock operator, is a sum of
the standard Fock operator and a one-body correlation MP2 potential. Molecular
orbitals and orbital energies are relaxed in the presence of correlation by
diagonalizing correlated Fock matrix. The double-excitation MP2 amplitudes are
then updated using those new molecular orbitals and orbital energies,
resulting in a self-consistency. Very recently, we have shown that the self-
consistency of OBMP2 can resolve issues caused by the non-iterative nature of
standard MP2 calculations for open-shell systems Tran (2021).

We explore the OBMP2 performance to predict electron density and related
properties for open-shell systems in the present work. We consider three
groups: atoms and their ions, main-row radicals, and halogen-bond dimers. We
find that OBMP2 can outperform standard MP2 and is comparable to CCSD. It is
important to emphasize that the computational scaling of OBMP2 is similar to
that of MP2 ($N^{5}$) that is lower than CCSD one ($N^{6}$). Therefore, OBMP2
is highly promising for accurately predicting the properties of large
molecules.

##  II Theory

Let us recap the OBMP2 theory whose formulation details are presented in Refs.
7 and 26. The OBMP2 approach was derived through the canonical transformation
developed by Yanai and his coworkersYanai and Chan (2006, 2007); Chan and
Yanai (2007); Neuscamman  _et al._ (2009, 2010a, 2010b). In this approach, an
effective Hamiltonian that includes dynamic correlation effects is achieved by
a similarity transformation of the molecular Hamiltonian $\hat{H}$ using a
unitary operator $e^{\hat{A}}$ :

$\displaystyle\hat{\bar{H}}=e^{\hat{A}^{\dagger}}\hat{H}e^{\hat{A}}$ (1)

with the anti-Hermitian excited operator $\hat{A}=-\hat{A}^{\dagger}$, and the
molecular Hamiltonian in spin orbitals as

$\displaystyle\hat{H}=h^{p\sigma}_{q\sigma}\hat{a}_{p\sigma}^{q\sigma}+\tfrac{1}{2}g^{p\sigma
r\sigma^{\prime}}_{q\sigma s\sigma^{\prime}}\hat{a}_{p\sigma
r\sigma^{\prime}}^{q\sigma s\sigma^{\prime}}$ (2)

where $\left\\{p,q,r,\ldots\right\\}$ indices refer to general spatial
orbitals and $\left\\{\sigma,\sigma^{\prime}\right\\}$ refer to spin indices.
In OBMP2, the cluster operator $\hat{A}$ is modeled such that including only
double excitation.

$\displaystyle\hat{A}=\hat{A}_{\text{D}}=\tfrac{1}{2}T_{i\sigma
j\sigma^{\prime}}^{a\sigma b\sigma^{\prime}}(\hat{a}_{i\sigma
j\sigma^{\prime}}^{a\sigma b\sigma^{\prime}}-\hat{a}_{a\sigma
b\sigma^{\prime}}^{i\sigma j\sigma^{\prime}})\,,$ (3)

with the MP2 amplitude

$\displaystyle T_{i\sigma j\sigma^{\prime}}^{a\sigma
b\sigma^{\prime}}=\frac{g_{i\sigma j\sigma^{\prime}}^{a\sigma
b\sigma^{\prime}}}{\epsilon_{i\sigma}+\epsilon_{j\sigma^{\prime}}-\epsilon_{a\sigma}-\epsilon_{b\sigma^{\prime}}}\,,$
(4)

where $\left\\{i,j,k,\ldots\right\\}$ indices refer to occupied spatial
orbitals and $\left\\{a,b,c,\ldots\right\\}$ indices refer to virtual spatial
orbitals. $\epsilon_{i\sigma}$ is the orbital energy of the spin-orbital
$i\sigma$.

The OBMP2 Hamiltonian is defined as

$\displaystyle\hat{H}_{\text{OBMP2}}=\hat{H}_{\text{HF}}+\left[\hat{H},\hat{A}_{\text{D}}\right]_{1}+\tfrac{1}{2}\left[\left[\hat{F},\hat{A}_{\text{D}}\right],\hat{A}_{\text{D}}\right]_{1}.$
(5)

with

$\displaystyle\hat{H}_{\text{HF}}=\hat{F}+C=f^{p\sigma}_{q\sigma}\hat{a}_{p\sigma}^{q\sigma}+\,\,C$
(6)

Where $\hat{F}$ is the Fock operator and $C$ is a constant. In Eq.5,
commutators with the subscription 1, $[\ldots]_{1}$, involve one-body
operators and constants that are reduced from many-body operators using the
cumulant approximationKutzelnigg and Mukherjee (1997); Mazziotti (1998a, b);
Kutzelnigg and Mukherjee (1999). We eventually derive at the OBMP2 Hamiltonian
as follows

$\displaystyle\hat{H}_{\text{OBMP2}}=$
$\displaystyle\,\,\hat{H}_{\text{HF}}+\hat{V}_{\text{OBMP2}}$ (7)

where $\hat{V}_{\text{OBMP2}}$ is a correlated potential composing of one-body
operators. The working expression are given in Refs. 7 and 26. We rewrite
$\hat{H}_{\text{OBMP2}}$ (Eqs. 5 and 7) in a similar form to Eq. 6 for
$\hat{H}_{\text{HF}}$ as follows:

$\displaystyle\hat{H}_{\text{OBMP2}}=$ $\displaystyle\hat{\bar{F}}+\bar{C}$
(8)

with $\hat{\bar{F}}=\bar{f}^{p\sigma}_{q\sigma}\hat{a}_{p\sigma}^{q\sigma}$.
$\bar{f}^{p\sigma}_{q\sigma}$ is so-called correlated Fock matrix and written
as

$\displaystyle\bar{f}^{p\sigma}_{q\sigma}$
$\displaystyle=f^{p\sigma}_{q\sigma}+v^{p\sigma}_{q\sigma}.$ (9)

The perturbation matrix $v^{p\sigma}_{q\sigma}$ serves as the correlation
potential altering the uncorrelated HF picture. We update the MO coefficients
and energies by diagonalizing the matrix $\bar{f}^{p}_{q}$, leading to orbital
relaxation in the presence of dynamic correlation effects. The OBMP2 method is
implemented within a local version of PySCFSun  _et al._ (2018). At the
convergence, the OBMP2 electron density is evaluated using the $T_{2}$
amplitude (Eq. 4) as in standard MP2.

##  III Results

We first consider the atomic electron density that has been often used to
examine the accuracy of DFT functionalsMedvedev  _et al._ (2017). We measured
the root-mean-square deviation (RMSD) relative to the CCSD electron density.
For comparison, we also performed calculations of standard methods, including
HF, MP2, and three DFT functional (B3LYP, HSE06, MN15-L) using pySCFSun  _et
al._ (2018). To see the effect of the orbital relaxation in OBMP2, we present
MP2 results with the unrelaxed density that corresponds to the first OBMP2
iteration.

Table 1:  Root mean square deviation (RMSD), normalized by $10^{-3}$ a.u, in
atomic electron density.

Atoms | HF | MP2a | OBMP2 | HSE06 | B3LYP | MN15-L
---|---|---|---|---|---|---
Al | 1.019 | 0.842 | 0.180 | 2.623 | 2.986 | 3.502
B | 0.890 | 0.699 | 0.427 | 0.444 | 8.959 | 10.119
B+ | 1.019 | 0.782 | 0.600 | 6.250 | 6.639 | 6.558
B- | 0.827 | 0.665 | 0.341 | 5.998 | 11.379 | 13.979
C | 0.877 | 0.681 | 0.413 | 14.655 | 8.927 | 12.771
C+ | 1.157 | 0.889 | 0.625 | 17.833 | 13.499 | 10.542
C- | 0.717 | 0.553 | 0.122 | 0.355 | 0.770 | 3.032
Cl | 1.553 | 1.108 | 0.094 | 8.204 | 10.806 | 11.529
F | 1.159 | 0.782 | 0.107 | 5.396 | 22.571 | 24.975
N | 0.462 | 0.336 | 0.100 | 0.442 | 0.922 | 3.345
N- | 2.899 | 2.810 | 0.182 | 4.861 | 9.812 | 11.627
N+ | 1.055 | 0.806 | 0.537 | 21.990 | 21.671 | 23.387
O | 0.833 | 0.596 | 0.117 | 22.532 | 22.569 | 22.464
O- | 1.571 | 1.056 | 0.295 | 19.904 | 17.967 | 9.971
O+ | 0.382 | 0.260 | 0.080 | 0.534 | 1.059 | 3.464
P | 1.359 | 1.061 | 0.079 | 1.135 | 1.530 | 3.088
S | 1.497 | 1.126 | 0.091 | 6.159 | 9.685 | 8.851
Si | 1.282 | 1.033 | 0.193 | 5.896 | 6.957 | 7.237
MADb | 1.142 | 0.894 | 0.255 | 8.067 | 9.928 | 10.580
MAXb | 2.899 | 2.810 | 0.625 | 22.532 | 22.571 | 24.975

a The MP2 density matrix is not relaxed.
b Mean absolute deviation (MAD) and maximum absolute deviation (MAX) relative
to CCSD.

Results are summarized in Table 1. All three DFT functionals have notoriously
huge errors for atomic systems considered in this work. Although HF performs
better than three DFT functionals, its errors are still quite large,
especially for N-. While MP2 can overall improve electron densities, it is
still unable to give a satisfying error for N-. When the electron density,
evaluated directly using the $T_{2}$ amplitude (Eq. 4), is relaxed via self-
consistency in OBMP2, the errors are several times smaller than MP2 errors. In
particular, as shown in Figure 1, the error for N- rapidly decreases with the
OBMP2 iteration, indicating the importance of self-consistency.

Figure 1: 3D surfaces of the difference between OBMP2 and CCSD densities
during the OBMP2 iteration for the N- system. The starting point of the self-
consistency (iter = 0) correspond to standard MP2 without relaxation. The
isosurface is 0.002.

We now analyze the Fukui functions for the removal of an electron $f^{+}(\bf
r)$ and the addition of an electron $f^{-}(\bf r)$ that are defined as
follows.

$\displaystyle f^{+}(\bf r)=\rho^{+}({\bf r})-\rho(\bf r)$ $\displaystyle
f^{-}(\bf r)=\rho({\bf r})-\rho^{-}(\bf r)$

Fukui functions$f^{+}(\bf r)$ and $f^{-}(\bf r)$ are the important
regioselectivity indicators for chemical species prone to accepting and
donating electrons in chemical reactions, respectivelyGalván  _et al._ (1986).
The reaction will take place where the Fukui functions can be found to have a
large value. 3D surfaces of errors relative to the CCSD reference of MP2 and
OBMP2 results are presented in Figure 2. Overall, except for the Fukui
function $f^{+}(r)$ of C and N, OBMP2 errors are much smaller than standard
MP2 errors and nearly vanish for $f^{+}(r)$ of O and $f^{-}(r)$ of N. The
Fukui functions and their generalization have been widely used to predict
molecular structureYañez  _et al._ (2021), we thus expect that OBMP2 will be
an affordable and accurate approach to such an interesting topic.

Figure 2: 3D surfaces of errors in the Fukui functions $f^{+}(r)$ (left
column) and $f^{-}(r)$ (right column) relative to the CCSD reference of four
atoms B, C, N, O. The isosurface is 0.002. Table 2:  Dipole moments of main-
group radicals at equilibrium geometry in the aug-cc-pCVTZ basis set.

Radicals |  | HF | MP2 |  |  |  | DSD- |  | OOMP2 |  | OBMP2 |  |  CCSDa |  |  CCSD(T)a
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
| unrelaxed |  | relaxed |  |  |  |  | BPEP86 |  |  |  |
AlH2 |  | 0.418 | 0.407 |  | 0.434 |  |  |  |  | 0.407 |  | 0.431 |  | 0.430 |  | 0.4069 |  | 0.399
BeH |  | 0.271 | 0.263 |  | 0.263 |  |  |  |  | 0.247 |  | 0.260 |  | 0.262 |  | 0.2338 |  | 0.230
BH2 |  | 0.495 | 0.497 |  | 0.498 |  |  |  |  | 0.514 |  | 0.497 |  | 0.497 |  | 0.5006 |  | 0.500
BN |  | 2.674 | 2.674 |  | 2.211 |  |  |  |  | 2.211 |  | 1.902 |  | 2.034 |  | 2.1593 |  | 2.000
BO |  | 2.904 | 2.738 |  | 2.468 |  |  |  |  | 2.374 |  | 2.135 |  | 2.337 |  | 2.4028 |  | 2.293
BS |  | 1.365 | 1.239 |  | 1.123 |  |  |  |  | 0.876 |  | 0.669 |  | 0.680 |  | 0.8451 |  | 0.730
CF |  | 0.430 | 0.508 |  | 0.735 |  |  |  |  | 0.744 |  | 0.843 |  | 0.849 |  | 0.6275 |  | 0.692
CH |  | 1.583 | 1.567 |  | 1.506 |  |  |  |  | 1.494 |  | 1.495 |  | 1.504 |  | 1.4423 |  | 1.416
CN |  | 2.181 | 2.133 |  | 2.305 |  |  |  |  | 1.663 |  | 0.896 |  | 1.360 |  | 1.5683 |  | 1.410
FCO |  | 1.283 | 1.166 |  | 0.779 |  |  |  |  | 0.782 |  | 0.653 |  | 0.685 |  | 0.8222 |  | 0.757
H2O-Al |  | 4.741 | 4.707 |  | 4.371 |  |  |  |  | 4.386 |  | 4.322 |  | 4.382 |  | 4.4076 |  | 4.355
H2O-F |  | 1.951 | 1.912 |  | 1.935 |  |  |  |  | 2.136 |  | 2.437 |  | 2.119 |  | 2.1011 |  | 2.173
LiN |  | 7.357 | 7.321 |  | 7.050 |  |  |  |  | 7.050 |  | 7.017 |  | 7.087 |  | 7.1024 |  | 7.020
NF |  | 0.325 | 0.247 |  | 0.000 |  |  |  |  | 0.041 |  | 0.183 |  | 0.078 |  | 0.007 |  | 0.069
NO |  | 0.249 | 0.145 |  | 0.031 |  |  |  |  | 0.190 |  | 0.308 |  | 0.195 |  | 0.1136 |  | 0.137
NS |  | 1.538 | 1.520 |  | 2.076 |  |  |  |  | 1.993 |  | 1.868 |  | 1.978 |  | 1.8406 |  | 1.839
OF |  | 0.420 | 0.356 |  | 0.225 |  |  |  |  | 0.095 |  | 0.212 |  | 0.016 |  | 0.0653 |  | 0.016
PO |  | 2.572 | 2.364 |  | 1.864 |  |  |  |  | 1.916 |  | 1.677 |  | 1.754 |  | 2.0908 |  | 1.935
PS |  | 0.746 | 0.695 |  | 0.464 |  |  |  |  | 0.620 |  | 0.431 |  | 0.601 |  | 0.7487 |  | 0.633
SF |  | 1.119 | 1.024 |  | 0.833 |  |  |  |  | 0.794 |  | 0.710 |  | 0.722 |  | 0.8862 |  | 0.808
SiH |  | 0.242 | 0.237 |  | 0.103 |  |  |  |  | 0.130 |  | 0.099 |  | 0.114 |  | 0.1206 |  | 0.102
MADb |  | 0.330 | 0.273 |  | 0.148 |  |  |  |  | 0.065 |  | 0.124 |  | 0.058 |  | 0.068 |  |
MAXb |  | 0.772 | 0.723 |  | 0.895 |  |  |  |  | 0.254 |  | 0.514 |  | 0.181 |  | 0.160 |  |

a Taken from Ref.22.
b Mean absolute deviation (MAD) and maximum absolute deviation (MAX) relative
to CCSD(T).

We now consider the dipole moment that has been widely used as the simplest
observable to examine the accuracy of new DFT functionals in predicting
electron densities. For example, Hait and Head-Gordon Hait and Head-Gordon
(2018a) have recently carried out an extensive assessment of the performance
of different DFT functionals belonging to different rungs of DFT Jacob’s
ladder. In the current work, we adopted 20 open-shell radicals from their work
and employed the aug-cc-pCVTZ basis set. For comparison, we performed
calculations for standard MP2 and a double-hybrid functional DSD-PBEP86 using
ORCA Neese  _et al._ (2020), and orbital-optimized MP2 (OOMP2) using Psi4Smith
_et al._ (2020). We adapted coupled-cluster results, CCSD and CCSD(T), from
Ref. 22 and use CCSD(T) as the reference.

All results are summarized in Table 2. CCSD yields the result closest to the
CCSD(T) reference. We can see that MP2 with unrelaxed density is slightly
better than HF. While the density relaxation can help to reduce MP2 errors,
they are still quite large. Surprisingly, OOMP2 errors are not much better
than standard MP2 with relaxed density, whereas the double hybrid functional
DSD-PBEP86 gives errors twice smaller than those of OOMP2. Overall, among
MP2-based approaches, OBMP2 yields the smallest errors and is comparable to
CCSD.

Figure 3: The change of dipole moments as stretching bond length for BN (left)
and CN (right). DHF stands for double-hybrid functional DSD-PBEP86. Standard
MP2 and DSD-PBEP86 densities are relaxed by solving the $Z-$vector equation
for orbital response.

As previously discussed by some authors, standard MP2 and double-hybrid
functionals with relaxed density give unphysical jumps of dipole moments in
the bond-stretched regime. This happens due to the singularity of the inverse
of the orbital Hessian matrix used for evaluating the orbital response of MP2
density matrixKurlancheek and Head-Gordon (2009). Head-Gordon and some other
authors have shown that orbital optimization is crucial to avoid unphysical
jumps of dipole moment Hait and Head-Gordon (2018b). Thus, it is interesting
to examine how OBMP2 dipole moments behave as stretching molecules. In Figure
3, we plot the change of dipole moments of BN and CN as functions of bond
length. For both radicals, it is clear that MP2 and DSD-PBEP86 dipole moments
with relaxed densities exhibit unphysical jumps. OOMP2 and OBMP2, on the
contrary, have continuous curves.

Table 3:  Atomic Mulliken charge of halogen-bond systems.

Systems |  |  | HF |  | MP2 |  |  OBMP2 |  |  HSE06 |  |  B3LYP |  |  MN15-L |  |  CCSD
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
FBr-Br | F |  | $-0.5540$ |  | $-0.5408$ |  | $-0.4573$ |  | $-0.4767$ |  | $-0.4639$ |  | $-0.4780$ |  | $-0.4654$
| Br1 |  | $\,\,\,\,0.5068$ |  | $\,\,\,\,0.4921$ |  | $\,\,\,\,0.4036$ |  | $\,\,\,\,0.5009$ |  | $\,\,\,\,0.4589$ |  | $\,\,\,\,0.4643$ |  | $\,\,\,\,0.4165$
| Br2 |  | $\,\,\,\,0.0472$ |  | $\,\,\,\,0.0486$ |  | $\,\,\,\,0.0537$ |  | $-0.0242$ |  | $\,\,\,\,0.0050$ |  | $\,\,\,\,0.0138$ |  | $\,\,\,\,0.0489$
FBr-Cl | F |  | $-0.5375$ |  | $-0.5245$ |  | $-0.4431$ |  | $-0.4608$ |  | $-0.4485$ |  | $-0.4597$ |  | $-0.4503$
| Br |  | $\,\,\,\,0.5250$ |  | $\,\,\,\,0.5114$ |  | $\,\,\,\,0.4284$ |  | $\,\,\,\,0.5403$ |  | $\,\,\,\,0.4917$ |  | $\,\,\,\,0.5194$ |  | $\,\,\,\,0.4377$
| Cl |  | $\,\,\,\,0.0125$ |  | $\,\,\,\,0.0131$ |  | $\,\,\,\,0.0147$ |  | $-0.0794$ |  | $-0.0432$ |  | $-0.0597$ |  | $\,\,\,\,0.0126$
FBr-OH | F |  | $-0.5556$ |  | $-0.5428$ |  | $-0.4607$ |  | $-0.4834$ |  | $-0.4698$ |  | $-0.4838$ |  | $-0.4682$
| Br |  | $\,\,\,\,0.5441$ |  | $\,\,\,\,0.5306$ |  | $\,\,\,\,0.4416$ |  | $\,\,\,\,0.5083$ |  | $\,\,\,\,0.4621$ |  | $\,\,\,\,0.4261$ |  | $\,\,\,\,0.4523$
| O |  | $-0.2903$ |  | $-0.2830$ |  | $-0.2696$ |  | $-0.2754$ |  | $-0.2555$ |  | $-0.3659$ |  | $-0.2694$
| H |  | $\,\,\,\,0.3018$ |  | $\,\,\,\,0.2952$ |  | $\,\,\,\,0.2887$ |  | $\,\,\,\,0.2506$ |  | $\,\,\,\,0.2632$ |  | $\,\,\,\,0.4236$ |  | $\,\,\,\,0.2854$
FCl-Br | F |  | $-0.4074$ |  | $-0.3956$ |  | $-0.3305$ |  | $-0.3618$ |  | $-0.3438$ |  | $-0.3566$ |  | $-0.3363$
| Cl |  | $\,\,\,\,0.3834$ |  | $\,\,\,\,0.3707$ |  | $\,\,\,\,0.2981$ |  | $\,\,\,\,0.3360$ |  | $\,\,\,\,0.3079$ |  | $\,\,\,\,0.3262$ |  | $\,\,\,\,0.3076$
| Br |  | $\,\,\,\,0.0240$ |  | $\,\,\,\,0.0249$ |  | $\,\,\,\,0.0323$ |  | $\,\,\,\,0.0258$ |  | $\,\,\,\,0.0359$ |  | $\,\,\,\,0.0304$ |  | $\,\,\,\,0.0287$
FCl-Cl | F |  | $-0.3931$ |  | $-0.3814$ |  | $-0.3179$ |  | $-0.3472$ |  | $-0.3291$ |  | $-0.3405$ |  | $-0.3236$
| Cl1 |  | $\,\,\,\,0.3857$ |  | $\,\,\,\,0.3738$ |  | $\,\,\,\,0.3066$ |  | $\,\,\,\,0.3463$ |  | $\,\,\,\,0.3156$ |  | $\,\,\,\,0.3234$ |  | $\,\,\,\,0.3141$
| Cl2 |  | $\,\,\,\,0.0074$ |  | $\,\,\,\,0.0077$ |  | $\,\,\,\,0.0113$ |  | $\,\,\,\,0.0009$ |  | $\,\,\,\,0.0135$ |  | $\,\,\,\,0.0172$ |  | $\,\,\,\,0.0095$
FCl-OH | F |  | $-0.4092$ |  | $-0.3975$ |  | $-0.3333$ |  | $-0.3627$ |  | $-0.3440$ |  | $-0.3547$ |  | $-0.3388$
| Cl |  | $\,\,\,\,0.4080$ |  | $\,\,\,\,0.3960$ |  | $\,\,\,\,0.3232$ |  | $\,\,\,\,0.3472$ |  | $\,\,\,\,0.3223$ |  | $\,\,\,\,0.3102$ |  | $\,\,\,\,0.3317$
| O |  | $-0.2877$ |  | $-0.2807$ |  | $-0.2667$ |  | $-0.2300$ |  | $-0.2309$ |  | $-0.3704$ |  | $-0.2668$
| H |  | $\,\,\,\,0.2890$ |  | $\,\,\,\,0.2822$ |  | $\,\,\,\,0.2768$ |  | $\,\,\,\,0.2456$ |  | $\,\,\,\,0.2526$ |  | $\,\,\,\,0.4149$ |  | $\,\,\,\,0.2738$
MADa |  |  | 0.0525 |  |  0.0433 |  |  0.0059 |  |  0.0356 |  |  0.0172 |  |  0.0446 |  |
MAXa |  |  | 0.0918 |  |  0.0783 |  |  0.0129 |  |  0.1025 |  |  0.0557 |  |  0.1411 |  |

aMean absolute deviation (MAD) and maximum absolute deviation (MAX) relative
to CCSD.

The last test set we consider is of halogen-bonding dimmers. In these systems,
the charge transfer between monomers importantly contributes to the bonding
mechanism Cavallo  _et al._ (2016). Here, we consider six open-shell dimmers
of FX–Y (X = Br, Cl, and Y = Br, Cl, OH) whose geometries were adopted from
Ref. 43. We use the aug-cc-pVTZ basis set and CCSD as a reference. WE did not
relax the density matrix of standard MP2 to show the importance of orbital
optimization in OBMP2. We also carried out three DFT functional calculations
for comparison, including B3LYP, HSE06, and MN15-L. All these calculations
were done using pySCFSun  _et al._ (2018).

Figure 4: Errors relative to the CCSD reference in dipole moments of halogen-
bond systems evaluated using some wavefunction and DFT methods.

Table 3 shows atomic Mulliken charges of the six dimmers. We can see that MP2
with unrelaxed density is slightly better than HF. Among three DFT
functionals, B3LYP performs best, and its errors are even smaller than
standard MP2. Thanks to the orbital optimization via self-consistency, OBMP2
dramatically improves upon MP2 results. Compared with other methods considered
here, its errors relative to the CCSD reference are the smallest. Figure 4
depicts errors in dipole moments relative to the CCSD reference. It is not
surprising that the DFT performance is system-dependent. While their errors
are small for HCl–Br and HCl–OH, they are large for other systems,
particularly HBr–Cl. HF and standard MP2 performance are quite stable, but
their errors are still large. OBMP2 outperforms standard MP2 and DFT and
yields the smallest errors with a maximum of 2.3% for HBr–Cl.

##  IV Conclusion

In summary, we have presented the performance of our recently-developed OBMP2
method to predict electron density and related properties of open-shell
systems. We have considered three groups: atoms and their ions, main-group
radicals, and halogen-bonding dimmers. We compared the OBMP2 performance
against standard methods, including HF, MP2, DFT, and coupled clusters. We
have found that OBMP2 yields result closest to CCSD in all cases considered in
this work. It is worth reminding that the formal scaling of OBMP2 is the same
as standard MP2 ($N^{5}$), which is lower than CCSD ($N^{6}$). OBMP2 is thus
believed to be an effective method for generating accurate electron densities.
Further speeding up OBMP2 by employing local and density-fitting
approximations is in progress.

## Acknowledgments

This work is supported by the Vietnam Academy of Science and Technology (VAST)
through the VAST Program for Young Researchers under the grant number
DLTE00.02/22-23. We performed all calculations on the High-Performance
Computing system at the Center for Informatics and Computing (CIC), VAST.

## References

  * Robert G. Parr (1994) Yang Weitao Robert G. Parr,  _Density functional theory of atoms and molecules_ (Oxford University Press, 1994).
  * Fabrizio  _et al._ (2019) Alberto Fabrizio, Andrea Grisafi, Benjamin Meyer, Michele Ceriotti, and Clemence Corminboeuf, “Electron density learning of non-covalent systems,” Chem. Sci. 10, 9424 (2019).
  * Kumar  _et al._ (2021) Nandan Kumar, Soumen Saha, and G Narahari Sastry, “Towards developing a criterion to characterize non-covalent bonds: a quantum mechanical study,” Phys. Chem. Chem. Phys. 23, 8478 (2021).
  * Boto  _et al._ (2020) Roberto A Boto, Francesca Peccati, Rubén Laplaza, Chaoyu Quan, Alessandra Carbone, Jean-Philip Piquemal, Yvon Maday, and Julia Contreras-García, “Nciplot4: Fast, robust, and quantitative analysis of noncovalent interactions,” J. Chem. Theory Comput. 16, 4150 (2020).
  * Cuevas-Zuviría and Pacios (2020) Bruno Cuevas-Zuviría and Luis F Pacios, “Analytical model of electron density and its machine learning inference,” J. Chem. Inf. Model. 60, 3831 (2020).
  * Kossmann and Neese (2010) Simone Kossmann and Frank Neese, “Correlated ab initio spin densities for larger molecules: orbital-optimized spin-component-scaled mp2 method,” J. Phys. Chem. A 114, 11768 (2010).
  * Tran (2021) Lan Nguyen Tran, “Improving perturbation theory for open-shell molecules via self-consistency,” J. Phys. Chem. A 125, 9242 (2021).
  * Kurlancheek and Head-Gordon (2009) Westin Kurlancheek and Martin Head-Gordon, “Violations of $n$-representability from spin-unrestricted orbitals in møller–plesset perturbation theory and related double-hybrid density functional theory,” Mol. Phys. 107, 1223 (2009).
  * Lochan and Head-Gordon (2007) Rohini C Lochan and Martin Head-Gordon, “Orbital-optimized opposite-spin scaled second-order correlation: An economical method to improve the description of open-shell molecules,” J. Chem. Phys. 126, 164101 (2007).
  * Neese  _et al._ (2009) Frank Neese, Tobias Schwabe, Simone Kossmann, Birgitta Schirmer, and Stefan Grimme, “Assessment of orbital-optimized, spin-component scaled second-order many-body perturbation theory for thermochemistry and kinetics,” J. Chem. Theory Comput. 5, 3060 (2009).
  * Bozkaya  _et al._ (2011) Uğur Bozkaya, Justin M Turney, Yukio Yamaguchi, Henry F Schaefer III, and C David Sherrill, “Quadratically convergent algorithm for orbital optimization in the orbital-optimized coupled-cluster doubles method and in orbital-optimized second-order møller-plesset perturbation theory,” J. Chem. Phys. 135, 104103 (2011).
  * Bozkaya and Sherrill (2013) Uğur Bozkaya and C David Sherrill, “Analytic energy gradients for the orbital-optimized second-order møller–plesset perturbation theory,” J. Chem. Phys. 138, 184103 (2013).
  * Bozkaya (2014a) Ugur Bozkaya, “Orbital-optimized second-order perturbation theory with density-fitting and cholesky decomposition approximations: An efficient implementation,” J. Chem. Theory Comput. 10, 2371 (2014a).
  * Bozkaya (2014b) Ugur Bozkaya, “Analytic energy gradients and spin multiplicities for orbital-optimized second-order perturbation theory with density-fitting approximation: an efficient implementation,” J. Chem. Theory Comput. 10, 4389–4399 (2014b).
  * Stück and Head-Gordon (2013) David Stück and Martin Head-Gordon, “Regularized orbital-optimized second-order perturbation theory,” J. Chem. Phys. 139, 244109 (2013).
  * Sharada  _et al._ (2015) Shaama Mallikarjun Sharada, David Stück, Eric J Sundstrom, Alexis T Bell, and Martin Head-Gordon, “Wavefunction stability analysis without analytical electronic hessians: application to orbital-optimised second-order møller–plesset theory and vv10-containing density functionals,” Mol. Phys. 113, 1802 (2015).
  * Razban  _et al._ (2017) Rostam M Razban, David Stück, and Martin Head-Gordon, “Addressing first derivative discontinuities in orbital-optimised opposite-spin scaled second-order perturbation theory with regularisation,” Mol. Phys. 115, 2102 (2017).
  * Lee and Head-Gordon (2018) Joonho Lee and Martin Head-Gordon, “Regularized orbital-optimized second-order møller–plesset perturbation theory: A reliable fifth-order-scaling electron correlation model with orbital energy dependent regularizers,” J. Chem. Theory Comput. 14, 5203 (2018).
  * Sandhoefer  _et al._ (2013) Barbara Sandhoefer, Simone Kossmann, and Frank Neese, “Derivation and assessment of relativistic hyperfine-coupling tensors on the basis of orbital-optimized second-order møller–plesset perturbation theory and the second-order douglas–kroll–hess transformation,” J. Chem. Phys. 138, 104102 (2013).
  * Medvedev  _et al._ (2017) Michael G Medvedev, Ivan S Bushmarinov, Jianwei Sun, John P Perdew, and Konstantin A Lyssenko, “Density functional theory is straying from the path toward the exact functional,” Science 355, 49 (2017).
  * Su  _et al._ (2018) Neil Qiang Su, Zhenyu Zhu, and Xin Xu, “Doubly hybrid density functionals that correctly describe both density and energy for atoms,” Proc. Natl. Acad. Sci. 115, 2287 (2018).
  * Hait and Head-Gordon (2018a) Diptarka Hait and Martin Head-Gordon, “How accurate is density functional theory at predicting dipole moments? an assessment using a new database of 200 benchmark values,” J. Chem. Theory Comput. 14, 1969 (2018a).
  * Peverati and Head-Gordon (2013) Roberto Peverati and Martin Head-Gordon, “Orbital optimized double-hybrid density functionals,” J. Chem. Phys. 139, 024110 (2013).
  * Hait and Head-Gordon (2018b) Diptarka Hait and Martin Head-Gordon, “Communication: xdh double hybrid functionals can be qualitatively incorrect for non-equilibrium geometries: Dipole moment inversion and barriers to radical-radical association using xyg3 and xygj-os,” J. Chem. Phys. 148, 171102 (2018b).
  * Sancho-Garcia  _et al._ (2016) Juan-Carlos Sancho-Garcia, AJ Pérez-Jiménez, Marika Savarese, E Brémond, and Carlo Adamo, “Importance of orbital optimization for double-hybrid density functionals: Application of the oo-pbe-qidh model for closed-and open-shell systems,” J. Phys. Chem. A 120, 1756–1762 (2016).
  * Tran and Yanai (2013) Lan Nguyen Tran and Takeshi Yanai, “Correlated one-body potential from second-order møller-plesset perturbation theory: Alternative to orbital-optimized mp2 method,” J. Chem. Phys. 138, 224108 (2013).
  * Yanai and Chan (2006) Takeshi Yanai and Garnet Kin-Lic Chan, “Canonical transformation theory for multireference problems,” J. Chem. Phys. 124, 194106 (2006).
  * Yanai and Chan (2007) Takeshi Yanai and Garnet Kin-Lic Chan, “Canonical transformation theory from extended normal ordering,” J. Chem. Phys. 127, 104107 (2007).
  * Chan and Yanai (2007) Garnet Kin-Lic Chan and Takeshi Yanai, “Canonical Transformation Theory for Dynamic Correlations in Multireference Problems,” in  _Adv. Chem. Phys., Volume 134_ , edited by D. A. Mazziotti (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2007) reduced-de ed., Chap. 13, p. 343.
  * Neuscamman  _et al._ (2009) Eric Neuscamman, Takeshi Yanai, and Garnet Kin-Lic Chan, “Quadratic canonical transformation theory and higher order density matrices,” J. Chem. Phys. 130, 124102 (2009).
  * Neuscamman  _et al._ (2010a) Eric Neuscamman, Takeshi Yanai, and Garnet Kin-Lic Chan, “Strongly contracted canonical transformation theory,” J. Chem. Phys. 132, 024106 (2010a).
  * Neuscamman  _et al._ (2010b) Eric Neuscamman, Takeshi Yanai, and Garnet Kin-Lic Chan, “A review of canonical transformation theory,” Int. Rev. Phys. Chem. 29, 231 (2010b).
  * Kutzelnigg and Mukherjee (1997) Werner Kutzelnigg and Debashis Mukherjee, “Normal order and extended wick theorem for a multiconfiguration reference wave function,” J. Chem. Phys. 107, 432 (1997).
  * Mazziotti (1998a) David A Mazziotti, “Contracted schrödinger equation: Determining quantum energies and two-particle density matrices without wave functions,” Phys. Rev. A 57, 4219 (1998a).
  * Mazziotti (1998b) David A Mazziotti, “Approximate solution for electron correlation through the use of schwinger probes,” Chem. Phys. Lett. 289, 419 (1998b).
  * Kutzelnigg and Mukherjee (1999) Werner Kutzelnigg and Debashis Mukherjee, “Cumulant expansion of the reduced density matrices,” J. Chem. Phys. 110, 2800 (1999).
  * Sun  _et al._ (2018) Qiming Sun, Timothy C Berkelbach, Nick S Blunt, George H Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D McClain, Elvira R Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan, “Pyscf: the python-based simulations of chemistry framework,” WIREs: Comput. Mol. Sci. 8, e1340 (2018).
  * Galván  _et al._ (1986) Marcelo Galván, José L Gázquez, and Alberto Vela, “Fukui function: Spin-density and chemical reactivity,” J. Chem. Phys. 85, 2337 (1986).
  * Yañez  _et al._ (2021) Osvaldo Yañez, Rodrigo Báez-Grez, Diego Inostroza, Ricardo Pino-Rios, Walter A Rabanal-León, Julia Contreras-García, Carlos Cardenas, and William Tiznado, “Kick–fukui: A fukui function-guided method for molecular structure prediction,” . Chem. Inf. Model. 61, 3955 (2021).
  * Neese  _et al._ (2020) Frank Neese, Frank Wennmohs, Ute Becker, and Christoph Riplinger, “The orca quantum chemistry program package,” J. Chem. Phys. 152, 224108 (2020).
  * Smith  _et al._ (2020) Daniel GA Smith, Lori A Burns, Andrew C Simmonett, Robert M Parrish, Matthew C Schieber, Raimondas Galvelis, Peter Kraus, Holger Kruse, Roberto Di Remigio, Asem Alenaizan, _et al._ , “Psi4 1.4: Open-source software for high-throughput quantum chemistry,” J. Chem. Phys. 152, 184108 (2020).
  * Cavallo  _et al._ (2016) Gabriella Cavallo, Pierangelo Metrangolo, Roberto Milani, Tullio Pilati, Arri Priimagi, Giuseppe Resnati, and Giancarlo Terraneo, “The halogen bond,” Chem. Rev. 116, 2478 (2016).
  * Bandyopadhyay and Seikh (2020) Prasanta Bandyopadhyay and Md Motin Seikh, “Components of the interaction energy of the odd-electron halogen bond: an ab initio study,” Phys. Chem. Chem. Phys. 22, 15389 (2020).


# Towards Compute-Optimal Transfer Learning

Massimo Caccia
Alexandre Galashov
Arthur Douillard22footnotemark: 2
Amal Rannen-Triki22footnotemark: 2
Dushyant Rao22footnotemark: 2
Michela Paganini22footnotemark: 2
Laurent Charlin
Marc’aurelio Ranzato22footnotemark: 2
Razvan Pascanu22footnotemark: 2
Current affiliation: Mila - Quebec AI Institute. Work done while interning at
DeepMind. email<EMAIL_ADDRESS>LondonMila, HEC Montréal,
Canada CIFAR AI Chair

###### Abstract

The field of transfer learning is undergoing a significant shift with the
introduction of large pretrained models which have demonstrated strong
adaptability to a variety of downstream tasks. However, the high computational
and memory requirements to finetune or use these models can be a hindrance to
their widespread use. In this study, we present a solution to this issue by
proposing a simple yet effective way to trade computational efficiency for
_asymptotic_ performance which we define as the performance a learning
algorithm achieves as compute tends to infinity. Specifically, we argue that
zero-shot structured pruning of pretrained models allows them to increase
compute efficiency with minimal reduction in performance. We evaluate our
method on the Nevis’22 continual learning benchmark that offers a diverse set
of transfer scenarios. Our results show that pruning convolutional filters of
pretrained models can lead to more than 20% performance improvement in low
computational regimes.

##  1 Asymptotic Performance versus Computational Efficiency

(a)

(b)

(c)

Figure 1:  Discover the ideal model architecture and size that best fits your
target data and training compute budget. The efficiency of different model
architectures/sizes varies in terms of asymptotic performance vs compute
efficiency. The left plot shows the goal of finding a model that optimally
fits your compute budget, but a pretrained foundation model is a recommended
approach. The middle plot highlights that the computational complexity of the
foundation model may exceed your budget. The right plot presents our solution:
zero-shot structured pruning, which increases computational efficiency at the
cost of reduced asymptotic performance, allowing the foundation model to fit
your compute budget.

An additional important aspect is that unstructured sparsity while in theory
should save compute, it rarely can be exploited currently, on traditional
hardware and due to the lack of low-level efficient sparse operations. To
alleviate this, we focus on structured sparsity that can be readily lead to
reduction in computation. We show that depending on the computational budget,
there exists a sparsity level that will lead to a better performance/compute
tradeoff then finetuning the original dense model.

Structured pruning is a popular technique to reduce the size of deep neural
networks and improve computational efficiency. One of the most commonly used
structured pruning approaches is channel pruning (He et al., 2017), which
reduces the number of channels in a layer. The computational gains from
channel pruning are substantial, as the FLOPs decrease at the square of the
pruning rate. However, as the representations become smaller, the model has
less capacity for learning.

An alternative approach to channel pruning is convolution filter pruning (Li
et al., 2016), which prunes at a finer granularity. With convolution filter
pruning, the representation size is preserved, but more aggressive pruning is
required to achieve the same computational speedup. This is because there are
fewer activations to connect between layers in channel pruning. Therefore, for
a fixed sparsity value, convolutional pruning enjoys more representations but
channel more weights.

Experiments will help us understand the relative merits of channel pruning and
convolution filter pruning. For a visual comparison, see Appendix A.

##  2 Related Work

The closest work to ours is Chen et al. (2021), which examines the existence
of sparse subnetworks in _dense_ pretrained computer vision models using the
lottery ticket hypothesis. With unstructured iterative magnitude pruning (IMP)
they find sparse subnetworks that retain the transfer performance of full
models. Our work extends this by using structured pruning for computational
speedup, zero-shot pruning to save compute, and empirically examining the
tradeoff between performance and compute.

In the field of transferring sparse pretrained models, several works have
explored the effect of different unstructured pruning techniques during
pretraining on transfer performance (Mehta, 2019; Morcos et al., 2019;
Paganini & Forde, 2020b; Sabatelli et al., 2020; Sun et al., 2022; Liu et al.,
2022). Another line of work focuses on pruning channels (He et al., 2017) and
convolution filters (Li et al., 2016). Additionally, structured IMP has been
studied in Chen et al. (2022); Paganini & Forde (2020a); Rachwan et al.
(2022).

Another work that is relevant to our paper is Frantar & Alistarh (2023), which
structurally sparsifies a pretrained GPT without much loss in performance, We
refer to Hoefler et al. (2021) for a thourought investigation of sparse neural
networks.

##  3 Experiments

###  3.1 Experimental Details

In order to study methods in both compute and performance space, we conduct
experiments using the new Nevis’22 benchmark (Bornschein et al., 2022). It
consists of a stream of over 100 visual classification tasks, sorted
chronologically and extracted from papers sampled from computer vision
proceedings spanning the last three decades. The stream reflects the research
community’s understanding of meaningful tasks at any given point in time, and
while limited to classification, it includes a diverse range of tasks
including OCR, texture analysis, crowd counting, scene recognition, and more.
The diversity of tasks is also reflected in the range of dataset sizes, which
span over four orders of magnitude.

For our experiments, we focus on the SHORT stream version of Nevis, which
consists of 25 tasks, including 16 object recognition tasks and 8 non-object
tasks (OCR, scenes, faces, texture, medical). To evaluate approaches that do
not start from scratch, Nevis’22 prescribes ImageNet pretraining. Accordingly,
we have removed ImageNet from the SHORT stream. The non-object tasks are
further considered out-of-distribution (OoD) with respect to ImageNet
pretraining.

The backbone of our model is a ResNet50 architecture, taken from the model
zoo, a platform of pretrained deep learning models. The ResNet is composed of
four ResNet blocks which we label 0 to 3. For the majority of our experiments,
we perform pruning locally, i.e., we prune the smaller weights of individual
layers. We prune such that the resulting model achieves a desired fraction of
the initial required FLOPs required for the forward and backward passes of the
model. For weight and convolutional pruning, a sparsity of X% means that X% of
the weights or convolution filters have been removed, resulting in only 1-X%
of the initial FLOPs. For channel pruning, we instead remove $1-\sqrt{1-X\%}$
because the FLOPs decrease at the square of the pruning rate

###  3.2 Zero-shot structured pruning of dense pretrained models efficiently
trades-off performance for compute

In our experiments, we first evaluate the performance of structured pruning in
the transfer learning scenario. The results, as shown in Figures 2(a) and
2(b), indicate that convolutional filter pruning performs similarly to
unstructured weight pruning, where the latter does not lead to computational
efficiency gains on accelerators. Furthermore, compared to channel pruning,
convolutional filter pruning is the most efficient way to prune the dense
pretrained model. This is especially true for transfer learning tasks that are
closer to the pretraining distribution, i.e. objects. It is possible that for
target distributions that are closer to the pretraining distribution, more of
the learned features are transferable, and therefore, it is necessary to keep
all of them.

Based on these results, we choose to stick with convolutional pruning for the
remaining experiments. We also evaluate the performance of convolutional
pruning in the continual learning scenario (see Figure 2(c) and 2(b)) and find
that it provides an efficient way to trade off performance for compute.

(a) Transfer learning – Object

(b) Transfer learning – Non-object

(c) Continual learning – Object

(d) Continual learning – Non-bject

Figure 2: Zero-shot structure pruning of dense pretrained models efficiently
trades-off performance for compute. In the first set of experiments, we
compare different pruning levels on the same pretrained model in transfer
learning (top plots). Unstructured Weight pruning achieves the desired
tradeoff, but doesn’t offer any computational gains on accelerators. Channel
pruning is not so efficient in the tradeoff. Alluringly, convolutional pruning
achieves almost the same tradeoff as weight pruning. In the second set of
experiments (bottom plots), we compare compare convolutional pruning with
finetuning the full model in the continual learning setting. Convolutional
pruning is again an effective solution.

###  3.3 Zero-shot structured pruning is better than other approaches

There are different approaches to reducing the computational complexity of a
pretrained model in transfer learning. The first set of baselines involves
finetuning only the last layers of the pretrained model. This approach saves
computation by only performing backpropagation on the last layers, which is
twice as expensive as the forward pass. The notation ft_block_Xto3 means
finetuning the ResNet blocks X to 3 (the last one), and visual support can be
found in the Appendix B.

The results are presented in Figure 3. The proposed method outperforms the
compared baselines in transfer learning tasks across different datasets,
especially in the object recognition domain. This is because the baselines
cannot modify the initial feature representation, which may require more
adaptation as the target data distribution deviates from the pre-training
distribution.

However, the pruning approach does not achieve superior performance over the
baselines in the continual learning scenario on object recognition datasets.
This outcome is expected, as the ImageNet pre-training already provides the
model with informative representations, leaving limited scope for knowledge
transfer between tasks. In this case, it may be more effective to freeze the
early layers to avoid catastrophic forgetting.

The second set of baselines involves selecting only the first layers and
training everything. This approach can modify the earlier representations of
the data but does not have the same depth as other approaches, and thus incurs
a loss of expressivity. The notation block_0toX_ft_all means picking the first
X blocks and finetuning everything, with visual support available in the
Appendix B.

The results are again displayed in Figure 3. The performance of the
block_0toX_ft_all baselines is inferior to that of the ft_block_Xto3. The
proposed zero-shot structured pruning approach demonstrates favorable results
across all scenarios. Although the baselines can modify the initial feature
representations, the loss in expressivity may render them less effective.

The final baseline utilizes 1x1 convolutional adapters (Rebuffi et al., 2017).
At first glance, the use of adapters may appear as a suboptimal trade-off
between performance and compute, as they do not reduce the computational
demands of the forward and backward passes. However, some literature suggests
that this approach may not be computationally efficient . Our experiments
confirm that adapters (denoted as adp1x1) are not an effective solution for
the studied use case, which is unsurprising.

(a) Transfer learning – Object

(b) Transfer learning – Non-object

(c) Continual learning – Object

(d) Continual learning – Non-bject

Figure 3: Our results demonstrate that the zero-shot structured pruning
approach is more effective in balancing the trade-off between performance and
compute compared to other methods. In both object and non-object recognition
tasks and for transfer and continual learning scenarios, the zero-shot
structured pruning approach exhibits superior efficiency. The pruning approach
benefits from both full modification of the feature representations and the
expressivity of the pretrained model, whereas the baselines are limited to one
or the other.

#### Ablations

We have conducted further empirical analysis, namely that of the relationship
between pruning rate and the trade-off between performance and compute (refer
to Appendix C). Additionally, we found that the choice of pruning signal is
crucial to the performance of the method, with pruning based on the magnitude
of the weights being more effective than random pruning (see Appendix D).

##  4 Future Work & Discussion

The study of scaling laws, which typically explores the relationship between
computation, data, and model size, is a topic of ongoing research. The
development of scaling laws for sparsity and transfer would have the potential
to predict the ideal level of sparsity for a specific transfer task and
computational budget. Another avenue of research is adaptive pruning, where
the pruning process is dynamically adjusted during training via the target
loss. Further investigation is also needed to better understand the effects of
pretraining on sparsity, including the number of pretraining steps, supervised
versus self-supervised training, and model size.

Comparing the results of structured pruning with other methods such as
distillation is crucial in evaluating its effectiveness and limitations.
Another promising area for future investigation is the use of zero-shot
structured pruning in the context of pretrained transformers for natural
language processing. In this case, channel pruning is analogous to pruning the
heads or activations in the transformer layer. Finally, we could consider
pruning foundation models in reinforcement learning for computationally-
efficient transfer reinforcement learning (Zhu et al., 2020) and continual
reinforcement learning (Caccia et al., 2022).

## References

  * Bornschein et al. (2022) Jorg Bornschein, Alexandre Galashov, Ross Hemsley, Amal Rannen-Triki, Yutian Chen, Arslan Chaudhry, Xu Owen He, Arthur Douillard, Massimo Caccia, Qixuang Feng, et al.  Nevis’22: A stream of 100 tasks sampled from 30 years of computer vision research.  _arXiv preprint arXiv:2211.11747_ , 2022.
  * Caccia et al. (2022) Massimo Caccia, Jonas Mueller, Taesup Kim, Laurent Charlin, and Rasool Fakoor.  Task-agnostic continual reinforcement learning: In praise of a simple baseline.  _arXiv preprint arXiv:2205.14495_ , 2022.  URL https://arxiv.org/abs/2205.14495.
  * Chen et al. (2021) Tianlong Chen, Jonathan Frankle, Shiyu Chang, Sijia Liu, Yang Zhang, Michael Carbin, and Zhangyang Wang.  The lottery tickets hypothesis for supervised and self-supervised pre-training in computer vision models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pp. 16306–16316, 2021.
  * Chen et al. (2022) Tianlong Chen, Xuxi Chen, Xiaolong Ma, Yanzhi Wang, and Zhangyang Wang.  Coarsening the granularity: Towards structurally sparse lottery tickets.  In _International Conference on Machine Learning_ , pp. 3025–3039. PMLR, 2022.
  * Frantar & Alistarh (2023) Elias Frantar and Dan Alistarh.  Massive language models can be accurately pruned in one-shot.  _arXiv preprint arXiv:2301.00774_ , 2023.
  * He et al. (2017) Yihui He, Xiangyu Zhang, and Jian Sun.  Channel pruning for accelerating very deep neural networks.  In _Proceedings of the IEEE international conference on computer vision_ , pp. 1389–1397, 2017.
  * Hoefler et al. (2021) Torsten Hoefler, Dan Alistarh, Tal Ben-Nun, Nikoli Dryden, and Alexandra Peste.  Sparsity in deep learning: Pruning and growth for efficient inference and training in neural networks.  _The Journal of Machine Learning Research_ , 22(1):10882–11005, 2021.
  * Li et al. (2016) Hao Li, Asim Kadav, Igor Durdanovic, Hanan Samet, and Hans Peter Graf.  Pruning filters for efficient convnets.  _arXiv preprint arXiv:1608.08710_ , 2016.
  * Liu et al. (2022) Yuanxin Liu, Fandong Meng, Zheng Lin, Peng Fu, Yanan Cao, Weiping Wang, and Jie Zhou.  Learning to win lottery tickets in BERT transfer via task-agnostic mask training.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pp. 5840–5857, Seattle, United States, July 2022. Association for Computational Linguistics.  doi: 10.18653/v1/2022.naacl-main.428.  URL https://aclanthology.org/2022.naacl-main.428.
  * Mehta (2019) Rahul Mehta.  Sparse transfer learning via winning lottery tickets.  _arXiv preprint arXiv:1905.07785_ , 2019.
  * Morcos et al. (2019) Ari Morcos, Haonan Yu, Michela Paganini, and Yuandong Tian.  One ticket to win them all: generalizing lottery ticket initializations across datasets and optimizers.  _Advances in neural information processing systems_ , 32, 2019.
  * Paganini & Forde (2020a) Michela Paganini and Jessica Forde.  On iterative neural network pruning, reinitialization, and the similarity of masks.  _arXiv preprint arXiv:2001.05050_ , 2020a.
  * Paganini & Forde (2020b) Michela Paganini and Jessica Zosa Forde.  Bespoke vs. pr$\backslash$^ et-$\backslash$a-porter lottery tickets: Exploiting mask similarity for trainable sub-network finding.  _arXiv preprint arXiv:2007.04091_ , 2020b.
  * Rachwan et al. (2022) John Rachwan, Daniel Zügner, Bertrand Charpentier, Simon Geisler, Morgane Ayle, and Stephan Günnemann.  Winning the lottery ahead of time: Efficient early network pruning.  In _International Conference on Machine Learning_ , pp. 18293–18309. PMLR, 2022.
  * Rebuffi et al. (2017) Sylvestre-Alvise Rebuffi, Hakan Bilen, and Andrea Vedaldi.  Learning multiple visual domains with residual adapters.  _Advances in neural information processing systems_ , 30, 2017.
  * Sabatelli et al. (2020) Matthia Sabatelli, Mike Kestemont, and Pierre Geurts.  On the transferability of winning tickets in non-natural image datasets.  _arXiv preprint arXiv:2005.05232_ , 2020.
  * Sun et al. (2022) Xing Sun, Ali Hassani, Zhangyang Wang, Gao Huang, and Humphrey Shi.  Disparse: Disentangled sparsification for multitask model compression.  _2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pp. 12372–12382, 2022.
  * Zhu et al. (2020) Zhuangdi Zhu, Kaixiang Lin, Anil K Jain, and Jiayu Zhou.  Transfer learning in deep reinforcement learning: A survey.  _arXiv preprint arXiv:2009.07888_ , 2020.

##  Appendix A Visual support for pruning type

(a) pretrained dense layer

(b) channel pruning

(c) unstructured convolution filter pruning

(d) structured convolution filter pruning

Figure 4:  Visualizing the Impact of Structured Pruning on Convolution Layers.
In this illustration, circles represent input/output feature maps ($H\times
W$), cubes symbolize convolution kernels ($K\times K$), and the $\otimes$
operator denotes the convolution. The rows and columns of the convolution
layer index the output and input channels, respectively. The pruning rate is
set at 50%. The channel pruning approach (b) removes groups of filters
responsible for a specific channel, reducing the output representation size
and FLOPs by half. In contrast, unstructured convolution filter pruning (c)
removes a fraction of filters for each channel while maintaining the output
representation size and FLOPs. Structured convolution filter pruning (d)
involves restructuring the layer, resulting in a halving of FLOPs, but
necessitating the broadcasting of inputs, increasing memory usage. However,
the increased memory usage is not prohibitive, as it only occurs locally at
the layer level and does not accumulate from one layer to the next.

##  Appendix B Visual support for parameter efficient methods

(a) ft_block_Xto3 baselines

(b) block_0toX_ft_all baselines

Figure 5: Visual support for parameter efficient methods. The pretrained
ResNet is formed of four blocks. Green blocks are trainable, blue are frozen
and white are removed.

##  Appendix C Pruning rate controls the asymptotic performance vs compute
efficiency trade-off

(a) Transfer Learning - Object

(b) Transfer Learning - Non-object

Figure 6: Pruning rate modulates the asymptotic performance vs computational
efficiency trade-off

##  Appendix D Pruning signal matters

(a) Transfer Learning - Object

(b) Transfer Learning - Non-object

Figure 7: Pruning signal matters Randomly pruning the network is not as
effective as pruning based on weight magnitude.