#  Underdamped stochastic thermodynamic engines
in contact with a heat bath with arbitrary temperature profile

Olga Movilla Miangolarra  Department of Mechanical and Aerospace Engineering,
University of California, Irvine, CA 92697, USA  Rui Fu  Department of
Mechanical and Aerospace Engineering, University of California, Irvine, CA
92697, USA  Amirhossein Taghvaei  Department of Mechanical and Aerospace
Engineering, University of California, Irvine, CA 92697, USA  Yongxin Chen
School of Aerospace Engineering, Georgia Institute of Technology, Atlanta, GA
30332, USA  Tryphon T. Georgiou  Department of Mechanical and Aerospace
Engineering, University of California, Irvine, CA 92697, USA

###### Abstract

We study thermodynamic processes in contact with a heat bath that may have an
arbitrary time-varying periodic temperature profile. Within the framework of
stochastic thermodynamics, and for models of thermodynamic engines in the
idealized case of underdamped particles in the low friction regime subject to
a harmonic potential, we derive explicit bounds as well as optimal control
protocols that draw maximum power and achieve maximum efficiency at any
specified level of power.

††preprint: APS/123-QED

For more than 200 years, Carnot’s idealized concept of a heat engine [1] has
been the cornerstone of equilibrium thermodynamics [2]. It provided a standard
as well as led to the discovery of the absolute temperature scale, entropy,
and the second law of thermodynamics. The sequel, beyond equilibrium and
quasi-static operation, proved a challenging one. Statistical mechanics and
stochastic fluctuation theories have since, and for over a hundred years, been
advancing our insight into the irreversibility of nonequilibrium thermodynamic
transitions. In this endeavor, a particular promising model of far-from-
equilibrium transitions is that of Stochastic Thermodynamics (also, Stochastic
Energetics) [3, 4, 5, 6] that followed at the heels of fluctuation theories
[7, 8, 9, 10, 11].

In the present work we adopt the framework of Stochastic Thermodynamics so as
to quantify the amount of power as well as the efficiency that thermodynamic
engines can achieve when harvesting energy from a heat bath with an arbitrary
but periodically fluctuating temperature profile. Our study is expected to be
of value, in particular, to theorists and engineers that work with nano
machines and nano sensing devices at the scale of molecular and thermal
fluctuations. Specifically, in our work, the Langevin equation

$\displaystyle dX_{t}$ $\displaystyle=v_{t}dt$ (1a) $\displaystyle mdv_{t}$
$\displaystyle=\\!-\\!\nabla_{x}U(t,X_{t})dt\\!-\\!\gamma
v_{t}dt\\!+\\!\sqrt{2\gamma k_{B}T(t)}dB_{t},$ (1b)

provides a model for a molecular system that interacts with a thermal
environment. In this, $X_{t}\in\mathbb{R}$ denotes the location of a particle
(which is taken to be one-dimensional for notational simplicity),
$v_{t}\in\mathbb{R}$ represents its velocity at time $t$, $m$ its mass,
$\gamma$ the viscosity coefficient of an ambient medium, $k_{B}$ the Boltzmann
constant, $T(t)$ the time-varying temperature of the heat bath at time $t$
that a statistical ensemble of particles is in contact with, $B_{t}$ a
standard Brownian motion that represents thermal excitation from a heat bath,
and $U(t,x)$ a time-varying potential exerting a force $-\nabla_{x}U(t,x)$ on
a particle at location $x\in\mathbb{R}$. This potential represents a control
action that models interaction of our ensemble of particles with the external
world.

This Langevin model has been extensively used to describe heat engines in a
Carnot-like cycle, that cycles through contact with heat baths of different
temperature, in order to quantify maximum power that can be drawn as well as
efficiency at maximum power. Specifically, foundational work has been carried
out in the overdamped regime [12, 13, 14] where inertial effects are
neglected, as well as in the underdamped regime [15, 16] under a low friction
assumption, where inertial effects now dominate dissipation. Most importantly,
experimental realization and verification of findings have been reported in
different settings [17, 18, 19].

Energy and entropy exchange between a heat bath and a particle ensemble is not
particular to (thermodynamic) engines. Fluctuations in chemical concentrations
in conjunction with the variability of electrochemical potentials may provide
a universal source of cellular energy [20]. Indeed, a similar mechanism
appears to be at play in biological engines. In these, energy exchange appears
to be mediated by continuous processes and energy differentials [21], in
contrast to the classical Carnot cycle which alternates between heat baths of
constant temperature in an idealized engine. Thus, as contemplated by E.
Schrödinger [22, 23], the ability of life to maintain a distance from
equilibrium may be tied to mechanisms that allow drawing energy out of
naturally or self-induced fluctuating temperature or chemical concentrations.

It is partly the above circle of ideas that led us to consider the rudimentary
model of a cyclic process (periodical temperature profile of a heat bath) that
fuels a molecular engine represented by stochastically driven Langevin
dynamics. Another paradigm that may provide context is that of molecular
engines that may some day be engineered to tap onto cyclic temperature
fluctuations for work production. The work we present below is an extension of
recent studies for periodic temperature profile in the linear regime [24, 25,
26]. It addresses the nonlinear regime and, in particular, the case of
underdamped dynamics where the low friction assumption applies.

##  I Energetics and low friction approximation

Throughout the paper we will consider the Langevin dynamics (1) with a
harmonic potential

$U(t,x)=\frac{1}{2}q(t)x^{2},$

where $q(t)$ represents its intensity and constitutes our control variable.
The total energy of a particle is the sum of kinetic and potential energies,

$E_{t}=\frac{1}{2}mv_{t}^{2}+\frac{1}{2}q(t)X_{t}^{2}.$

The work and heat exchange at the level of a single particle, $d\bar{}\hskip
1.99997ptW_{t}$ and $d\bar{}\hskip 1.99997ptQ_{t}$, respectively, are [3]

$\displaystyle d\bar{}\hskip 1.99997ptW_{t}$
$\displaystyle=\frac{1}{2}\dot{q}(t)X_{t}^{2}dt,$ (2a) $\displaystyle
d\bar{}\hskip 1.99997ptQ_{t}$ $\displaystyle=-\gamma
v_{t}^{2}{\color[rgb]{0,0,0}{dt}}+\frac{\gamma k_{B}T(t)}{m}dt+\sqrt{2\gamma
k_{B}T(t)}v_{t}dB_{t}.$ (2b)

These definitions are in agreement with the first law of thermodynamics, i.e.
$dE=d\bar{}\hskip 1.99997ptQ+d\bar{}\hskip 1.99997ptW$. Taking the expectation
of (2) yields the average rate of work and heat exchange,

$\displaystyle d\bar{}\hskip 1.99997pt\mathcal{W}_{t}$
$\displaystyle=\frac{1}{2}\dot{q}(t)\Sigma_{x}(t)dt,$ (3a) $\displaystyle
d\bar{}\hskip 1.99997pt\mathcal{Q}_{t}$
$\displaystyle=\bigg{(}-\gamma\Sigma_{v}(t)+\frac{\gamma
k_{B}T(t)}{m}\bigg{)}dt,$ (3b)

where $\Sigma_{x}$ and $\Sigma_{v}$ are the variances of $X_{t}$ and $v_{t}$,
respectively. The averaged expressions also satisfy the first law of
thermodynamics $d\mathcal{E}_{t}=d\bar{}\hskip
1.99997pt\mathcal{Q}_{t}+d\bar{}\hskip 1.99997pt\mathcal{W}_{t}$ where

$\mathcal{E}_{t}=\frac{1}{2}q(t)\Sigma_{x}(t)+\frac{1}{2}m\Sigma_{v}(t)$

is the average total energy of the particle.

The governing equations for the variance are obtained from the Langevin
dynamics (1),

$\displaystyle\dot{\Sigma}_{x}(t)$ $\displaystyle=2\Sigma_{xv}(t),$ (4a)
$\displaystyle\dot{\Sigma}_{xv}(t)$
$\displaystyle=\Sigma_{v}(t)-\frac{q(t)}{m}\Sigma_{x}(t)-\frac{\gamma}{m}{\Sigma}_{xv}(t),$
(4b) $\displaystyle\dot{\Sigma}_{v}(t)$
$\displaystyle=-\frac{2q(t)}{m}{\Sigma}_{xv}(t)-\frac{2\gamma}{m}{\Sigma_{v}}(t)+\frac{2\gamma
k_{B}}{m^{2}}T(t),$ (4c)

where $\Sigma_{xv}(t)$ is the correlation between $X_{t}$ and $v_{t}$.

Although the differential equations for the evolution of the variance are
explicit and in closed-form, the analysis for maximizing efficiency is still
challenging due to the presence of the three coupled ODEs. As a result, it is
common to consider simplifying approximations that hold in different physical
regimes.

Our analysis concerns the underdamped (or low friction) regime that is studied
in [15, 16]. When the temperature $T(t)$ and the harmonic potential intensity
$q(t)$ are constant, the low friction regime holds when the particle in (1)
undergoes many oscillations before dissipation effects become significant.
Mathematically, this amounts to $\sqrt{\frac{q}{m}}\gg\frac{\gamma}{m}$. In
the case where $T(t)$ and $q(t)$ are time varying periodic functions, it is
required in addition that their power spectrum mostly lies below the natural
frequency of oscillations $\sqrt{q/m}$.

Under the low friction assumptions, the equipartition condition

$\frac{1}{2}q(t)\Sigma_{x}(t)=\frac{1}{2}m\Sigma_{v}(t)$ (5)

holds approximately [16, Section 4.2]. As a result, it is possible to express
the dynamics of our system as a function of $\Sigma_{v}(t)$ uniquely. More
precisely, we can use (4a) to express $\Sigma_{xv}(t)$ as follows:

$\Sigma_{xv}(t)=\frac{1}{2}\dot{\Sigma}_{x}=\frac{m}{2q(t)}\dot{\Sigma}_{v}-\frac{m\dot{q}(t)}{2q(t)^{2}}{\Sigma_{v}},$

where the second equality comes from the equipartition condition (5). Now, we
can rewrite $\Sigma_{xv}(t)$ in (4c), to obtain

$\dot{\Sigma}_{v}(t)=\Sigma_{v}(t)\bigg{(}\frac{\dot{q}(t)}{2q(t)}-\frac{\gamma}{m}\bigg{)}+\frac{\gamma
k_{B}}{m^{2}}T(t).$ (6)

For this equation to be meaningful, we require $q(t)$ to be (weakly)
differentiable. Equation (6) is the underlying model for our analysis through
the rest of the paper.

##  II Analysis of maximum power

We consider a stochastic thermodynamic engine driven by a periodic temperature
$T(t)$ with period $t_{f}$. The averaged power output over a cycle is equal to

$\displaystyle\mathcal{P}=-\frac{1}{t_{f}}\int_{0}^{t_{f}}d\bar{}\hskip
1.99997pt\mathcal{W}_{t}.$

Using the first law of thermodynamics $d\mathcal{E}_{t}=d\bar{}\hskip
1.99997pt\mathcal{Q}_{t}+d\bar{}\hskip 1.99997pt\mathcal{W}_{t}$, the power
can be expressed in terms of heat exchange as follows,

$\displaystyle\mathcal{P}$
$\displaystyle=\frac{1}{t_{f}}\int_{0}^{t_{f}}(d\bar{}\hskip
1.99997pt\mathcal{Q}_{t}-d\mathcal{E}_{t})$
$\displaystyle=\frac{\gamma}{t_{f}}\int_{0}^{t_{f}}\bigg{(}\frac{k_{B}T(t)}{m}-\Sigma_{v}(t)\bigg{)}dt+\frac{1}{t_{f}}(\mathcal{E}(t_{f})-\mathcal{E}(0)),$
(7)

where the definition (3b) is used and $\mathcal{E}_{t}\equiv\mathcal{E}(t)$,
for clarity. To account for possible discontinuities in $T(t)$ and $q(t)$, as
it is often the case in the underdamped setting [27, 15, 28], we consider a
cycle from $0^{+}$ to $t_{f}^{+}$. Thus, $q,\Sigma_{v}$ and $\mathcal{E}$,
which are periodic, satisfy

$q(0^{+})=q(t_{f}^{+})\quad\mbox{and}\quad~{}\Sigma_{v}(0^{+})=\Sigma_{v}(t_{f}^{+}),$
(8)

and $\mathcal{E}(0^{+})=\mathcal{E}(t_{f}^{+})$. The latter condition removes
the boundary terms in (7), while conditions (8) impose an integral constraint
on $\Sigma_{v}(t)$. Specifically, if we divide (6) by $\Sigma_{v}(t)$ and
integrate from $0^{+}$ to $t_{f}^{+}$ we obtain

$\log\left(\frac{\Sigma_{v}(t_{f}^{+})}{\Sigma_{v}(0^{+})}\right)=\log\left(\frac{q(t_{f}^{+})}{q(0^{+})}\right)^{\frac{1}{2}}-\frac{\gamma}{m}t_{f}+\frac{\gamma
k_{B}}{m^{2}}\int_{0}^{t_{f}}\frac{T(t)}{\Sigma_{v}(t)}dt.$ (9)

Thus, in view of (8),

$\int_{0}^{t_{f}}\frac{T(t)}{\Sigma_{v}(t)}dt=\frac{mt_{f}}{k_{B}}.$ (10)

This integral constraint is the main ingredient in our analysis for both
maximum power and maximum efficiency.

We first consider maximizing power output, namely,

$\displaystyle\max_{\Sigma_{v}(\cdot)}\frac{\gamma}{t_{f}}\int_{0}^{t_{f}}\bigg{(}\frac{k_{B}T(t)}{m}-\Sigma_{v}(t)\bigg{)}dt,$
(11)

subject to the constraint (10). We observe that due to (10),

$\displaystyle\int_{0}^{t_{f}}\Sigma_{v}(t)dt$
$\displaystyle=\int_{0}^{t_{f}}\left(\Sigma_{v}(t)+c\frac{T(t)}{\Sigma_{v}(t)}\right)dt-c\frac{mt_{f}}{k_{B}}$

for any arbitrary constant $c$. The minimal value for the integrand
$\Sigma_{v}(t)+c\frac{T(t)}{\Sigma_{v}(t)}$, pointwise and with respect to
$\Sigma_{v}(t)$, is attained for $\Sigma_{v}(t)=\sqrt{cT(t)}$. This choice for
$\Sigma_{v}(t)$ also satisfies (10) for
$\sqrt{c}=\frac{k_{B}}{mt_{f}}\int_{0}^{t_{f}}\sqrt{T(t)}dt$. Therefore,

$\Sigma_{v}(t)=\left(\frac{k_{B}}{mt_{f}}\int_{0}^{t_{f}}\sqrt{T(s)}ds\right)\sqrt{T(t)}$
(12)

is the unique maximizer of (11).

The relation between $\Sigma_{v}(t)$ and the protocol $q(t)$ can be obtained
from (6) by direct integration (cf. (9)), which for the optimizing
$\Sigma_{v}(\cdot)$ gives

$q(t)=q_{0}\frac{T(t)}{T(0)}\exp\left(\frac{2\gamma
t_{f}}{m}\left(\frac{t}{t_{f}}-\frac{\int_{0}^{t}\sqrt{T(s)}ds}{\int_{0}^{t_{f}}\sqrt{T(s)}ds}\right)\right).$
(13)

Note that, $q_{0}$ specifies the starting value $q(0)$ of the optimal
protocol. Further, $q(t)$ is continuous at all times when $T(t)$ is, which
agrees with the intuition that jumps in the control are necessitated for
adapting to a sudden jump in temperature [12, 15].

Finally, inserting (12) into the expression for power, we obtain that, under
this protocol, any value for $q_{0}$ gives the same level of power, namely,

$\mathcal{P}^{*}=\frac{\gamma k_{B}}{m}\text{Var}(\sqrt{T})$ (14)

where

$\text{Var}(\sqrt{T})=\frac{1}{t_{f}}\int_{0}^{t_{f}}T(t)dt-\left(\frac{1}{t_{f}}\int_{0}^{t_{f}}\sqrt{T(t)}dt\right)^{2}\;$

quantifies fluctuations in $\color[rgb]{0,0,0}{\sqrt{T(t)}}$.

Thus, expression (14) quantifies the maximum power that can be drawn in terms
of an explicit measure of temperature fluctuations.

It is interesting to consider the shape of temperature fluctuations that
allows maximum power to be drawn under the optimal protocol. In view of the
above, it is easy to see that when $T(t)\in[T_{c},T_{h}]$, the Carnot
temperature profile

$\displaystyle T_{\rm Carnot}(t)=\begin{cases}T_{h},\quad
t\in[0,\frac{t_{f}}{2})\\\ T_{c},\quad
t\in[\frac{t_{f}}{2},t_{f}),\end{cases}$ (15)

allows drawing maximum power, which in this case is

$\frac{\gamma k_{B}}{m}\frac{(\sqrt{T_{h}}-\sqrt{T_{c}})^{2}}{4}.$

This special profile (Carnot) has already been studied in [15].

##  III Analysis of maximum efficiency at fixed power

The classical definition of efficiency is the ratio of work produced over the
amount of heat drawn from the hot bath into the system,

$\eta_{\;{}_{\mathcal{Q}}}=\frac{\color[rgb]{0,0,0}-\mathcal{W}}{\;\mathcal{Q}_{h}}.$
(16)

This definition presumes that the environment that constitutes the hot heat
bath is well defined, and as a consequence that $\mathcal{Q}_{h}$ is well
defined as well. This is clearly the case for a Carnot cycle, $T(t)=T_{\rm
Carnot}(t)$, when the system remains in contact with the bath at $T_{h}$
during the interval $[0,\frac{t_{f}}{2})$.

For this case, at maximum power,

$\mathcal{Q}_{h}=\int_{0}^{\frac{t_{f}}{2}}{\color[rgb]{0,0,0}d\bar{}\hskip
1.99997pt}\mathcal{Q}_{t}=\frac{\gamma
k_{B}}{m}\frac{t_{f}}{4}\sqrt{T_{h}}(\sqrt{T_{h}}-\sqrt{T_{c}}),$

where the second equality is obtained using (3b) and the optimal expression
for $\Sigma_{v}$ in (12). It follows that

$\eta_{\;{}_{\mathcal{Q}}}=1-\sqrt{\frac{T_{c}}{T_{h}}},$

which is precisely the Curzon-Ahlborn efficiency [29, 30]. This result was
obtained in [15].

However, in the setting of an arbitrary periodic temperature profile, it is
not entirely clear what constitutes the hot bath. The same bath serves as hot
and cold, at different times, and therefore $\mathcal{Q}_{h}$ and
$\eta_{\;{}_{\mathcal{Q}}}$ may be given different definitions and be subject
to different interpretations [24, 26].

In the present work, following [26], we adopt

$\eta_{\;{}_{\mathcal{U}}}=\frac{\color[rgb]{0,0,0}-\mathcal{W}}{\;\mathcal{U}},$

as our definition of efficiency, where

$\mathcal{U}=-k_{B}\int^{t_{f}}_{0}S(\rho_{t})\dot{T}(t)dt,$

(also, $\mathcal{U}=k_{B}\int T(t)\dot{S}(\rho_{t})dt$, using integration by
parts) represents the effective uptake of thermal energy, while

$S(\rho_{t})=-\int\int\rho_{t}(x,v)\log(\rho_{t}(x,v))dxdv$

denotes the entropy of the particle-distribution $\rho_{t}$. The value
$\mathcal{U}$ can be thought of as the maximal amount of work that can be
extracted when the dissipation is zero and the transition takes place
quasistatically. In general, over a cycle,

$\mathcal{U}=-\mathcal{W}+\mathcal{W}_{\rm diss},$

where the dissipation

$\mathcal{W}_{\rm diss}=\int
k_{B}T(t)\underbrace{\left(\dot{S}(\rho_{t})dt-\frac{d\bar{}\hskip
1.99997pt\mathcal{Q}}{k_{B}T(t)}\right)}_{\dot{S}_{\rm total}},$

relates to the total entropy production rate $\dot{S}_{\rm total}$ in the
system and the environment. For quasi-static transitions $\mathcal{W}_{\rm
diss}=0$ and $\eta_{\;{}_{\mathcal{U}}}=1$.

In the underdamped low friction regime with quadratic potential, the
distribution $\rho_{t}$ is Gaussian with independent components $(x,v)$ whose
covariance are $\Sigma_{x}(t)$ and $\Sigma_{v}(t)$, respectively. Using the
standard formula for the entropy of Gaussian distributions together with the
equipartition condition (5) expressed in terms of the variances, $\mathcal{U}$
simplifies to

$\mathcal{U}=-\frac{k_{B}}{2}\int^{t_{f}}_{0}\ln\Big{(}(2\pi
e)^{2}\frac{m}{q(t)}\Sigma_{v}^{2}(t)\Big{)}\dot{T}(t)dt.$

In case $T(t)$ is discontinuous, the limits of integration need to be taken as
$0^{+}$ and $t_{f}^{+}$, and similarly in what follows. Integration by parts
gives

$\displaystyle\mathcal{U}=$
$\displaystyle\frac{k_{B}}{2}\int^{t_{f}}_{0}\bigg{(}2\frac{\dot{\Sigma}_{v}(t)}{\Sigma_{v}(t)}-{\frac{\dot{q}(t)}{q(t)}}\bigg{)}{\color[rgb]{0,0,0}T(t)}dt$
$\displaystyle=$
$\displaystyle\frac{k_{B}^{2}\gamma}{m^{2}}\int^{t_{f}}_{0}\frac{T(t)^{2}}{\Sigma_{v}(t)}dt-\frac{k_{B}\gamma}{m}\int^{t_{f}}_{0}T(t)dt$

where we have used (6) to establish the second equality.

We now seek to optimize efficiency at a given power. To this end, it suffices
to minimize the dissipation $\mathcal{U}$ subject to the given power. Thus, we
consider

$\displaystyle\min_{\Sigma_{v}(\cdot)}~{}$
$\displaystyle\int^{t_{f}}_{0}\frac{T(t)^{2}}{\Sigma_{v}(t)}dt,$ (17)

subject to (10) together with

$\frac{\gamma}{t_{f}}\int_{0}^{t_{f}}\bigg{(}\frac{k_{B}T(t)}{m}-\Sigma_{v}(t)\bigg{)}dt=\mathcal{P},$

for any given level of power $\mathcal{P}$. We follow a variational approach.
The Lagrangian for our optimization problem is

$\mathcal{L}=\frac{T(t)^{2}}{\Sigma_{v}(t)}+\lambda\bigg{(}\Sigma_{v}(t)-\frac{k_{B}T(t)}{m}+\frac{\mathcal{P}}{\gamma}\bigg{)}+\mu\bigg{(}\frac{T(t)}{\Sigma_{v}(t)}-\frac{m}{k_{B}}\bigg{)},$

where $\lambda,\mu$ are the Lagrange multipliers. The stationarity condition
for $\Sigma_{v}(t)$ gives

$\Sigma_{v}(t)=\frac{1}{\sqrt{\lambda}}\sqrt{T(t)^{2}+\mu T(t)}.$

Applying the power constraint we obtain

$\sqrt{\lambda}=\frac{\int^{t_{f}}_{0}\sqrt{T(t)^{2}+\mu
T(t)}dt}{\frac{k_{B}}{m}\int^{t_{f}}_{0}T(t)dt-\frac{t_{f}}{\gamma}\mathcal{P}}.$

In a similar manner as before, the optimal protocol can now be obtained by
direct integration of (6), to give

$\displaystyle q(t)=$ $\displaystyle q_{0}\frac{T(t)^{2}+\mu
T(t)}{T(0)^{2}+\mu T(0)}\exp\left(\frac{2\gamma(t-r(t))}{m}\right),$ (18)

where

$r(t)=\frac{\int_{0}^{t_{f}}\sqrt{T(s)^{2}+\mu
T(s)}ds}{\int_{0}^{t_{f}}T(s)ds-\frac{t_{f}m}{\gamma
k_{B}}\mathcal{P}}\int^{t}_{0}\frac{\sqrt{T(s)}}{\sqrt{T(s)+\mu}}ds.$

It now remains to fix $\mu$ by imposing the last constraint, which yields
$r(t_{f})=t_{f}$. This is,

$\displaystyle\int^{t_{f}}_{0}\sqrt{T(t)^{2}+\mu T(t)}$ $\displaystyle
dt\int^{t_{f}}_{0}\frac{\sqrt{T(t)}}{\sqrt{T(t)+\mu}}dt$ $\displaystyle=$
$\displaystyle t_{f}\int^{t_{f}}_{0}T(t)dt-\frac{t_{f}^{2}m}{\gamma
k_{B}}\mathcal{P}.$ (19)

In general, equation (III) cannot be solved explicitly for $\mu$. However, it
is easily seen that when $\mathcal{P}=0$, then $\mu=0$ is a solution. Also, as
$\mathcal{P}\to\mathcal{P^{*}}$, then the solution $\mu\to\infty$. It is
worthwhile to point out that the condition $r(t_{f})=t_{f}$ makes the protocol
(18) continuous for continuous temperature profiles.

Figure 1:  (a) Maximum efficiency at fixed power for the temperature profiles
portrayed in (b). (b) Temperature profiles with colors and patterns matching
the corresponding performance curves in (a).

For notational convenience, we let

$\overline{g(t)}=\frac{1}{t_{f}}\int_{0}^{t_{f}}g(t)dt$

denote the average of a periodic function $g(t)$ over the period. Now, putting
all of the ingredients together, the optimal $\mathcal{U}$ is

$\displaystyle\mathcal{U^{*}}=\frac{t_{f}}{\kappa}\left[\frac{\overline{\sqrt{(T(t)+\mu)T(t)}}}{\overline{T(t)}-\kappa\mathcal{P}}\times\overline{\left(\frac{T(t)^{\frac{3}{2}}}{\sqrt{(T(t)+\mu)}}\right)}-\overline{T(t)}\right]$

where ${\kappa}=m/(\gamma k_{B})$. Consequently, the maximum efficiency at
fixed power is

$\eta_{\;{}_{\mathcal{U}}}^{*}(\mathcal{P})=\frac{\kappa\mathcal{P}}{\left[\frac{\overline{\sqrt{(T(t)+\mu)T(t)}}}{\overline{T(t)}-\kappa\mathcal{P}}\times\overline{\left(\frac{T(t)^{\frac{3}{2}}}{\sqrt{(T(t)+\mu)}}\right)}-\overline{T(t)}\right]}.$
(20)

Figure 1 displays the trade-off between power and efficiency for three
different temperature profiles. Note that for vanishing power we obtain the
maximum possible efficiency, i.e. $\eta_{\;{}_{\mathcal{U}}}^{*}=1$, that
corresponds to vanishing dissipation. On the other hand, as
$\kappa\mathcal{P}\to\kappa\mathcal{P}^{*}=\text{Var}(\sqrt{T})$ (from (14)),
the limiting value from (20) gives that

$\eta_{\;{}_{\mathcal{U}}}^{*}\to\frac{\overline{\sqrt{T(t)}}\text{Var}(\sqrt{T})}{\overline{T(t)^{\frac{3}{2}}}-\overline{T(t)}\times\overline{\sqrt{T(t)}}}.$
(21)

This is precisely the efficiency at maximum power obtained by using the
optimal protocol (13). The limit in (21) can also be expressed in terms of the
third moment of the square root of the temperature

$\mu_{3}(\sqrt{T})=\frac{1}{\;t_{f}}\int_{0}^{t_{f}}\Big{(}\sqrt{T(t)}-\overline{\sqrt{T(t)}}\Big{)}^{3}dt,$

specifically,

$\eta_{\;{}_{\mathcal{U}}}^{*}\to$
$\frac{1}{2+\frac{\mu_{3}(\sqrt{T})}{\text{Var}(\sqrt{T})\overline{\sqrt{T(t)}}}}.$
(22)

Thus, for any temperature profile such that $\mu_{3}(\sqrt{T})=0$, the limit
is $\frac{1}{2}$, which in particular is true for the Carnot-like temperature
profile $T(t)=T_{\rm Carnot}(t)$. We also note that the third moment can be
negative, leading to efficiency that exceeds $\frac{1}{2}$, albeit at the cost
of decreasing maximum power.

##  IV Example: Sinusoidal temperature profile

In the following we compare our results to the ones obtained in the linear
response regime, where it is assumed that $T_{h}-T_{c}\ll T_{h}+T_{c}$ [24,
31]. To this end, we consider the sinusoidal temperature profile

$T(t)=\frac{T_{h}+T_{c}}{2}+\frac{T_{h}-T_{c}}{2}\cos(\omega t),$ (23)

with $\omega\ll\sqrt{\frac{q}{m}}$ so as to ensure that the low friction
assumption holds.

We first use our analytical result to recover the expressions obtained in the
linear response regime. Specifically, assuming $\Delta
T=\frac{T_{h}-T_{c}}{2}$ is much smaller than $\bar{T}=\frac{T_{h}+T_{c}}{2}$,
we have that

$\displaystyle q^{*}(t)$ $\displaystyle=q_{0}+q_{0}\frac{\Delta
T}{\bar{T}}\cos(\omega t)+\mathcal{O}(\Delta T^{2}),$
$\displaystyle\cal{P}^{*}$ $\displaystyle=\frac{\gamma k_{B}}{8m}\frac{\Delta
T^{2}}{\bar{T}}+\mathcal{O}(\Delta T^{3}),$
$\displaystyle\eta_{\;{}_{\mathcal{U}}}^{*}(\mathcal{P^{*}})$
$\displaystyle=\frac{1}{2}+\mathcal{O}(\Delta T^{2}),$

where the first two coincide with the expressions for the optimal protocol and
power that are given in [24, 31], in the limit where $\frac{q_{0}}{m}$ is
large enough to satisfy the requirements of the low friction regime. In fact,
in the linear response regime
$\eta_{\;{}_{\mathcal{U}}}^{*}(\mathcal{P^{*}})=\frac{1}{2}+\mathcal{O}(\Delta
T^{2})$ holds for any temperature profile whose time dependent component is an
odd function of $t$.

Figure 2: Comparison between the optimal protocols $q(t)$ in the low friction
and linear response regimes, for the sinusoidal temperature profile shown with
a red-dotted curve.

Moreover, in order to illustrate the utility of our results, we compare the
optimality of the protocol obtained in this paper, to the one obtained in the
linear response regime. Specifically, we solve numerically the equations for
the variance (4) for two protocols:

  * i)

the optimal protocol obtained with low friction assumption according to (13),

  * ii)

the optimal protocol obtained in linear response regime in [24] and [31, Eqs
29, 30a].

These two protocols, along with the temperature profile, are shown in Figure
2. Then, we use the numerical solution to the equations for the variance to
compute power and efficiency. The result, as the temperature difference
$T_{h}-T_{c}$ is varied, is shown in Figure 3.

In these numerical examples, we made sure
$\frac{\gamma}{m}\ll\sqrt{\frac{q}{m}}$ so that the low friction assumption
remains valid. As it is expected, the optimal protocol under the linear
response regime assumption performs well when the temperature ratio is small,
giving a numerical result that agrees with equations (13) and (20). However,
as the temperature ratio increases, higher order terms become significant and
the protocol is no longer optimal.

Figure 3: Power and efficiency for (1) with a temperature profile (23) as
functions of temperature ratios for two protocols: i) the optimal protocol
obtained under the low friction assumption (13) (solid line), and ii) the
optimal protocol obtained in the linear response regime [31, Eqs 29, 30a]
(dashed line). For comparison, the analytical expressions for maximum power
(14) and efficiency (21) in the low friction regime are also drawn (dotted
curve).

We finally highlight how the numerical value of power and efficiency, using
the optimal protocol (13), deviates from the analytical expressions for
maximum power (14) and efficiency at maximum power (21) with increasing values
of the friction coefficient $\gamma$. This is shown in Figure 4. Specifically,
for small values of $\frac{\gamma}{\sqrt{mq}}$, the analytical and numerical
results for power and efficiency coincide. For larger values of
$\frac{\gamma}{\sqrt{mq}}$ the approximation error becomes significant due to
the fact that the condition $\frac{\gamma}{\sqrt{mq}}\ll 1$, and hence, the
approximate equipartition relation (5) , no longer hold.

##  V Concluding remarks

The present work quantifies maximum power, as well as the efficiency at any
level of power, that can be achieved by a thermodynamic engine in the
(underdamped) low friction regime, modeled by a Langevin equation and in
contact with a heat bath with periodic temperature profile.

It turns out that the maximum power and efficiency are expressed in terms of
finite moments of the square root of the temperature profile. Moreover, when
the temperature changes continuously, the optimal control protocol is
expressed explicitly in terms of the temperature and certain finite moments,
and is continuous in time as well.

Of particular interest is the efficiency at maximum power and, specifically,
the relation between useful work and dissipation. Interestingly, when the
third moment of the square root of the temperature profile is zero (for
instance, due to a suitable time symmetry), at maximum power in the low
friction regime, the losses in dissipation equal the amount of work that can
be extracted by the engine, yielding $\eta_{\;{}_{\mathcal{U}}}=\frac{1}{2}$.
It is of interest to explore possible connections of this result to the
universal bound on the efficiency at maximum power being
$\frac{1}{2}\eta_{C}$, where $\eta_{C}$ is the Carnot efficiency [32, 14].

##  VI Acknowledgements

OMM, RF and AT contributed equally and AT oversaw the technical development.
The research was partially supported by NSF under grants 1807664, 1839441,
1901599, 1942523, and AFOSR under FA9550-20-1-0029.

Figure 4:  Numerical values for power and efficiency as a function of the
friction coefficient for the system governed by equation (4), driven by the
optimal protocol (13) obtained in the low friction regime (solid line). For
comparison, we have included the analytical expressions for maximum power (14)
and the efficiency at maximum power (21) obtained under the low friction
assumption (dotted line).

## References

  * Carnot [1986] S. Carnot,  _Reflexions on the motive power of fire: a critical edition with the surviving scientific manuscripts_ (Manchester University Press, 1986).
  * Adkins [1983] C. J. Adkins,  _Equilibrium thermodynamics_ (Cambridge University Press, 1983).
  * Sekimoto [2010] K. Sekimoto,  _Stochastic energetics_ , Vol. 799 (Springer, 2010).
  * Seifert [2012] U. Seifert, Stochastic thermodynamics, fluctuation theorems and molecular machines, Reports on progress in physics 75, 126001 (2012).
  * Chen  _et al._ [2020] Y. Chen, T. Georgiou, and A. Tannenbaum, Stochastic control and non-equilibrium thermodynamics: fundamental limits, IEEE Transactions on Automatic Control 65, 252 (2020).
  * Jarzynski [2011] C. Jarzynski, Equalities and inequalities: Irreversibility and the second law of thermodynamics at the nanoscale, Annual Review of Condensed Matter Physics 2 (2011).
  * Jarzynski [1996] C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78 (1996).
  * Gallavotti and Cohen [1995] G. Gallavotti and E. G. D. Cohen, Dynamical ensembles in nonequilibrium statistical mechanics, Physical Review Letters 74 (1995).
  * Evans and Searles [1994] D. J. Evans and D. J. Searles, Equilibrium microstates which generate second law violating steady states, Physical Review E 50 (1994).
  * Crooks [1999] G. E. Crooks, Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences, Physical Review E 60 (1999).
  * Hatano and Sasa [2001] T. Hatano and S. Sasa, Steady-state thermodynamics of Langevin systems, Physical Review Letters 86 (2001).
  * Schmiedl and Seifert [2007] T. Schmiedl and U. Seifert, Efficiency at maximum power: An analytically solvable model for stochastic heat engines, Europhysics Letters 81 (2007).
  * Fu  _et al._ [2021] R. Fu, A. Taghvaei, Y. Chen, and T. T. Georgiou, Maximal power output of a stochastic thermodynamic engine, Automatica 123, 109366 (2021).
  * Van den Broeck [2005] C. Van den Broeck, Thermodynamic efficiency at maximum power, Physical Review Letters 95 (2005).
  * Dechant  _et al._ [2017] A. Dechant, N. Kiesel, and E. Lutz, Underdamped stochastic heat engine at maximum efficiency, Europhysics Letters 119 (2017).
  * Zöller [2017] N. Zöller,  _Optimization of Stochastic Heat Engines in the Underdamped Limit_ (Springer, 2017).
  * Blickle and Bechinger [2011] V. Blickle and C. Bechinger, Realization of a micrometre-sized stochastic heat engine, Nature Physics 8 (2011).
  * Abah  _et al._ [2012] O. Abah, J. Roßnagel, G. Jacob, S. Deffner, F. Schmidt-Kaler, K. Singer, and E. Lutz, Single-ion heat engine at maximum power, Physical review letters 109, 203006 (2012).
  * Martínez  _et al._ [2015] I. A. Martínez, E. Roldán, L. Dinis, D. Petrov, and R. A. Rica, Adiabatic processes realized with a trapped brownian particle, Physical Review Letters 114 (2015).
  * Cilek [2009] V. Cilek,  _Earth System: History and Natural Variability-Volume III_ , Vol. 3 (EOLSS Publications, 2009).
  * Yanagida  _et al._ [2007] T. Yanagida, M. Ueda, T. Murata, S. Esaki, and Y. Ishii, Brownian motion, fluctuation and life, Biosystems 88, 228 (2007).
  * Schrödinger [1944] E. Schrödinger, The physical aspect of the living cell, What is Life , 70 (1944).
  * Ornes [2017] S. Ornes, Core concept: how nonequilibrium thermodynamics speaks to the mystery of life, Proceedings of the National Academy of Sciences 114, 423 (2017).
  * Bauer  _et al._ [2016] M. Bauer, K. Brandner, and U. Seifert, Optimal performance of periodically driven, stochastic heat engines under limited control, Physical Review E 93, 042112 (2016).
  * Bauer  _et al._ [2015] M. Bauer, K. Brandner, and U. Seifert, Thermodynamics of micro- and nano-systems driven by periodic temperature variations, Physical Review X 5 (2015).
  * Brandner and Saito [2020] K. Brandner and K. Saito, Thermodynamic geometry of microscopic heat engines, Physical Review Letters 124 (2020).
  * Dechant  _et al._ [2015] A. Dechant, N. Kiesel, and E. Lutz, All-optical nanomechanical heat engine, Physical Review Letters 114 (2015).
  * Gomez-Marin  _et al._ [2008] A. Gomez-Marin, T. Schmiedl, and U. Seifert, Optimal protocols for minimal work processes in underdamped stochastic thermodynamics, Journal of Chemical Physics 129 (2008).
  * Novikov [1958] I. I. Novikov, The efficiency of atomic power stations (a review), Journal of Nuclear Energy 7 (1958).
  * Curzon and Ahlborn [1975] F. L. Curzon and B. Ahlborn, Efficiency of a Carnot engine at maximum power output, Americal Journal of Physics 43 (1975).
  * Fu  _et al._ [2020] R. Fu, O. Movilla Miangolarra, A. Taghvaei, Y. Chen, and T. T. Georgiou, Harvesting energy from a periodic heat bath, in  _2020 59th IEEE Conference on Decision and Control (CDC)_ (2020) pp. 3034–3039.
  * Esposito  _et al._ [2009] M. Esposito, K. Lindenberg, and C. Van den Broeck, Universality of efficiency at maximum power, Physical review letters 102, 130602 (2009).


Mason R. Swofford

Department of Computer Science

Stanford University

<EMAIL_ADDRESS>
https://github.com/mswoff/Image-Completion-on-CIFAR-10

Image Completion on CIFAR-10

###### Abstract

This project performed image completion on CIFAR-10, a dataset of 60,000 32x32
RGB images, using three different neural network architectures: fully
convolutional networks, convolutional networks with fully connected layers,
and encoder-decoder convolutional networks. The highest performing model was a
deep fully convolutional network, which was able to achieve a mean squared
error of .015 when comparing the original image pixel values with the
predicted pixel values. As well, this network was able to output in-painted
images which appeared real to the human eye.

##  1 Introduction

Image completion, or image in-painting, is a fascinating and quickly growing
area of research in machine learning. Given an image with some pixels removed,
the image completion task is to correctly fill in those missing pixels.
Therefore, the input to an image completion algorithm is an image with some
pixels removed, and the algorithm outputs its estimate of the values of those
missing pixels. This has various use cases. For instance, if an image is
partially corrupted, the corrupted pixels can be removed, and an image
completion algorithm can be used to fill in the missing pixels so that the
image looks like the original, uncorrupted version. As well, object removal
can be accomplished with image completion. The pixels corresponding to the
object can be removed, and then an image completion algorithm can fill in the
missing pixels values in accordance with the rest of the image.

This can be a challenging problem, because in some cases a large portion of
the image is removed, and the algorithm must try and recover the missing
information using only the remaining parts of the image. In fact, this is not
a problem most humans perform well at, as it takes artistic abilities to make
an image look natural. My problem will face further challenges, as I chose to
use very small images which therefore contain less information in them. My
approach was to remove a square from the center of each image using a mask and
then apply a variety of different convolutional networks to the masked images,
to see which could best predict the original images. I drew inspiration for my
network architectures from some of the networks presented in the next section.

##  2 Related work

Current state of the art in image completion largely relies on GANs
(Generative Adversarial Networks), which are able to learn a distribution for
images, and then try to fill in the missing pixel values according to this
learned distribution. These networks are able to perform incredibly well even
when a large portion of the original image is removed. However, GANs, and in
particular DCGANs (Deep Convolutional Generative Adversarial Networks),
require huge amounts of data, are typically only trained and used on one class
of images (e.g. faces, scenery, or cats), and are employed on larger images
which contain more pixels and therefore more information [2].

Prior to DCGANs appearance, most image completion algorithms did not rely on
machine learning, and instead focused on mathematical methods for in-painting.
One interesting algorithm performed in-painting by scouring a large corpus of
images to find semantically similar ones, and then in-painted the corrupted
image by taking fragments from these similar corpus images [5]. Other
techniques tried to find similar patches within the same image that could fill
in the missing pixels while satisfying certain constraints [4].

In researching prior work on image completion, I noticed parallels between
image completion and image segmentation, where the input is an image and the
output is a map of which pixels correspond to different objects in the image.
Both of these tasks require high dimensional outputs, which is different from
the most common computer vision task, classification. As well, if an algorithm
can identify which pixels correspond to an object, then I theorized if some of
those pixels are removed, it can more easily in-paint them by using
information from the object. These ideas lead me to consider networks used in
image segmentation.

The Fully Convolutional Network (FCN) used by Long and Shelhamer was one of
the first networks to perform exceptionally well on image segmentation. Unlike
previous work, their network only used convolutional layers and could
therefore be used on inputs of arbitrary sizes [7]. The simplicity of their
approach along with the high performance motivated me to try this model as my
initial approach.

SegNet built on FCN by employing an encoder-decoder architecture where the
image was first transformed into a low dimensional encoding, and then that low
dimensional encoding was up-sampled to produce an output of the same size as
the input [3]. The idea of using encoders and decoders has yielded astounding
results in fields outside of image segmentation as well, most notably in
machine translation.

##  3 Dataset and Features

Figure 1: Sample Training Example

As mentioned, I wanted to see how image completion algorithms would fare on
smaller datasets with both fewer training examples and smaller image sizes. As
well, I wanted the dataset to include multiple object classes, to see how well
an algorithm could perform under tougher conditions. I settled on the CIFAR-10
dataset because it met all of these conditions. CIFAR-10 is a subset of the
Tiny Images dataset which contains 80 million 32x32 RGB images. CIFAR-10
includes just 10 object classes, with 6,000 images per class for a total of
60,000 32x32 RGB images. These classes vary widely, from airplanes to frogs to
trucks [6]. The data was split into 50,000 train images, 5,000 development
images, and 5,000 test images. Since CIFAR-10 was intended for object
classification, some pre-processing was required to transform the dataset.

Instead of feeding each image to the network directly, a mask was applied that
zeroed out the center 8x8 pixels. This 8x8 region was now what the network was
trying to predict. Therefore, the true label for each image was the inverse
mask applied to the image, resulting in an 8x8 RGB label. This label was
flattened into a 192 (8x8x3) dimensional vector, so that it could be more
easily fed into a loss function. The pixel values were also scaled from the
usual [0,255] range down to the [0,1] range. An example of this process is
shown in Figure 1.

##  4  Methods

In total, 6 different networks were trained and compared. While there were 6
total networks, they fall into just 3 different categories: fully
convolutional network, convolutional network with fully connected layers, and
encoder-decoder convolutional network. Regardless of the algorithm used, each
network was given a masked-out image and had to predict the 192-dimensional
output. Mean squared error was used as the loss function for all of the
algorithms. It computes the mean squared difference between each element in
the predicted 192-vector and the label 192-vector. This means that it compares
how close each channel value (red, green, and blue) is in the prediction
versus the ground truth. By optimizing mean squared error, the algorithms
learn to match the prediction to the true value for each pixel. The formula
for mean squared error is as follows:

$MSE=\frac{1}{192}\sum_{i=1}^{192}(y_{i}-\hat{y}_{i})^{2}$

As well, all the networks were trained using mini-batch gradient descent with
the Adam optimization algorithm with variable learning rates which were hand
tuned. All layers (except max pooling and flattening layers) use rectified
linear activations. Since the final per channel output should be in the range
[0,1], sigmoid activations were used at first. However, they proved to work
poorly and instead the rectified linear activation was used even for the final
layer, and it was clipped to be in the range [0,1].

3 fully convolutional networks were trained: a shallow, deep, and super deep
network. By only including convolutional and max pooling layers and no
deconvolutions or fully connected layers, these networks did not have access
to the same amount of global information as the other networks and were
theorized to over-fit the training data less. The deeper networks not only
have more layers, but they also have more filters per layer. This allows them
to learn more features from the input image (e.g. texture patterns, lines,
etc.). The kernels used along with the input size at each layer are as
follows, where the part in brackets is the kernel applied and number of
filters used is in parentheses. For example, 32x32x3 ->[5x5 Convolution
(10)]-> 28x28x10 means that 10 5x5 filters were applied to the 32x32x3 input
and resulted in a 28x28x10 output for that layer. Note: all convolutions are
“valid” convolutions unless otherwise noted, and all max pooling layers are
2x2 filters with strides of 2.

Shallow Network: 32x32x3 ->[5x5 Convolution (10)]-> 28x28x10 ->[5x5
Convolution (20)]-> 24x24x20 ->[Max Pooling]-> 12x12x20 ->[5x5 Convolution
(20)]-> 8x8x20 ->[1x1 Convolution (3)]-> 8x8x3 ->[Flatten]-> 192x1

Deep Network: 32x32x3 ->[5x5 Convolution (20)]-> 28x28x20 ->[5x5 Convolution
(40)]-> 24x24x40 ->[Max Pooling]-> 12x12x40 ->[3x3 Convolution (60)]->
10x10x60 [3x3 Convolution (80)]-> 8x8x80 ->[1x1 Convolution (3)]-> 8x8x3
->[Flatten]-> 192x1

Super Deep Network: 32x32x3 ->[7x7 “Same” Convolution (40)]-> 32x32x40 ->[7x7
“Same” Convolution (40)]-> 32x32x40 ->[7x7 Convolution (40)]-> 26x26x40 ->[7x7
Convolution (60)]-> 20x20x60 ->[5x5 Convolution (60)]-> 16x16x60 ->[Max
Pooling]-> 8x8x60 ->[5x5 Convolution (60)]-> 4x4x60 ->[4x4 Convolution
(192)]-> 1x192

A convolutional network with fully connected layers was also trained. This
model was based on the models typically used in object classification, where
first convolutions are applied, and then the output is flattened and fed into
a series of fully connected layers and finally a softmax layer. However, my
network differed from this in that I did not have a softmax layer, and instead
of having just one output, I had 192. Following the previous conventions, the
network is as follows.

Fully Connected Network: 32x32x3 ->[5x5 Convolution (10)]-> 28x28x10 ->[5x5
Convolution (20)]-> 24x24x20 ->[Max Pooling]-> 12x12x20 ->[5x5 “Same”
Convolution (30)]-> 12x12x30 ->[3x3 “Same” Convolution (40)]-> 12x12x40
->[Flatten]-> 5,760x1 ->[Fully Connected Layer]-> 768x1 ->[Fully Connected
Layer]-> 192x1

The last models tried were 2 encoder-decoder networks. These models first
created a lower dimensional representation of the missing pixels, and then use
deconvolutions to up-sample the representation and create the final output. By
first creating a lower dimensional encoding of the missing pixels, encoder-
decoder networks tend to over-fit less, as they can only store so much
information in the lower dimension. The deconvolutions used were all 3x3
kernels with strides of 2. Following the previous conventions, the networks
are as follows.

Encoder Network: 32x32x3 ->[5x5 Convolution (10)]-> 28x28x10 ->[5x5
Convolution (20)]-> 24x24x20 ->[Max Pooling]-> 12x12x20 ->[5x5 Convolution
(30)]-> 8x8x30 ->[5x5 Convolution (128)-> 4x4x128 ->[3x3 Deconvolution (3)]->
6x6x3 ->[3x3 Deconvolution (3)]-> 8x8x3 ->[Flatten]-> 192x1

Deep Encoder Network: 32x32x3 ->[5x5 Convolution (30)]-> 28x28x30 ->[5x5
Convolution (40)]-> 24x24x40 ->[Max Pooling]-> 12x12x40 ->[5x5 Convolution
(40)]-> 8x8x40 ->->[5x5 Convolution (40)]-> 4x4x40 ->[Flatten]-> 1x640
->[Fully Connected Layer]-> 1x768 ->[3x3 Deconvolution (40)]-> 2x2x40 ->[3x3
Deconvolution (40)]-> 4x4x40 ->[3x3 Deconvolution (3)]-> 8x8x3 ->[Flatten]->
192x1

##  5 Experiments/Results/Discussion

Various hyper-parameters were experimented with over the course of this
project, including the different model architectures, activation functions,
and learning rates. I began with the idea to try at least one of each type of
architecture: fully convolutional network, convolutional network with fully
connected layer, and encoder-decoder network.

For the learning rate, I experimented with different values ranging from .005
to .00001. However, this was all done by looking at the development error and
hand tuning the learning rate for each model. As well, I decreased the
learning rate as the models trained.

The first model tested was the shallow fully convolutional network (Shallow
Network). With this network, I tested the sigmoid verses rectified linear
activation for the final layer, and quickly saw that the sigmoid had poor mean
squared error compared to the rectified linear activation. So, for all
subsequent models, I only used the rectified linear activation (clipped to be
in the range [0,1]); although it is possible the sigmoid activation would have
worked better for some of these networks, it would have doubled the number of
models I had to test. Using the rectified linear activation, the algorithm
achieved a train and development loss that were approximately equal, .0153 and
.0154, respectively. Figure 2a shows the epoch vs loss graph for this model
(the first 5 epochs are removed for all graphs). While I did not know what
Baye’s error would be for this task, after analyzing the images output by the
network, they were not as good as I knew was possible. This lead me to believe
that the algorithm had a bias problem and lead me to my second model.

(a) Shallow Network

(b) Super Deep Network

(c) Deep Network

Figure 2: Epoch vs Loss Graphs

Super Deep Network had orders of magnitude more parameters due to having more
layers and more filters per layer. This allowed the model to learn more
complex features than Shallow Network. However, the network took up to a
minute to run per epoch, compared to around 5 seconds for Shallow Network. As
well, it became apparent that there was not enough data in the training set to
feed this deeper network, as it was able to over-fit the training set with a
train loss of .0161 versus a development loss of .0253 after 320 epochs.
Figure 2b shows the epoch vs loss graph for this model. While the train loss
would likely have continued to decline if I trained for more epochs, it was
clear the model was overfitting the data and the development loss was
increasing every epoch. I believe that more data likely would have reduced a
lot of the variance, and I considered adding the CIFAR-100 dataset (60,000
additional tiny images) to the training set or performing data augmentation on
CIFAR-10. However, since part of my goal was to see how these algorithms
performed on limited data, I decided against it and instead tried a new
network, which was somewhere between Super Deep Network and Shallow Network in
complexity.

The next model, Deep Network, had 5 convolutional layers unlike Shallow
Network’s 4 and Super Deep Network’s 7. As well, it had more filters than
Shallow Network, which I believed to be the biggest reason Shallow Network
could not fit the data. Deep Network was able to perform much better than
either previous network, achieving train and development losses of 0.013620
and 0.014624, respectively, after 870 epochs. Figure 2c shows the epoch vs
loss graph for Deep Network.

Next, I built a model similar to Shallow Network but with 2 fully connected
layers at the end. This model, Fully Connected Network, was able to fit the
training set exceptionally well. Figure 3a shows the epoch vs loss graph for
this model. However, it was unable to achieve development loss as low as Deep
Network. The network clearly had a huge variance problem and lead me to
abandon exploring more fully connected networks. I believed the issue was that
predicting a high dimensional output using fully connected layers was too ripe
for over-fitting.

The final class of models, the encoder-decoder networks, both performed
reasonably well. Encoder Network achieved minimum development loss after 440
epochs, with a train and development loss of 0.015688 and 0.016934,
respectively. Figure 3b shows the epoch vs loss graph for this model. After
training for another 200 epochs, the algorithms development loss increased
slightly and its training loss decreased slightly; however, it appeared to
have almost leveled off, signifying that the model likely was under-fitting
the data.

The final model, Deep Encoder Network, however over-fit the data. It achieved
a training loss as low as .010 and development loss of .018 after 440 epochs.
Figure 3c shows the epoch vs loss graph for this model. This network utilized
a fully connected layer to help make the encoding, which could have
contributed to it overfitting the training data.

(a) Fully Connected Network

(b) Encoder Network

(c) Deep Encoder Network

Figure 3: Epoch vs Loss Graphs

In order to better understand the mistakes the algorithms were making, I
analyzed some of the training and development images output by the networks
(Note: the algorithms only output the center square and I recombined this with
the input image). In Figure 4, I compared the output of Deep Network versus
the output of Fully Connected Network on two sample development set images.
While I used means squared error as my loss function, my end goal was an
algorithm that could fill in the missing pixels to create realistic looking
images. Therefore, here I was looking to see whether a human could tell that
the algorithms’ outputs were not real images. In particular, I was trying to
see if it was obvious that the 8x8 center region had been altered. Looking at
4b, this image appeared realistic enough to be an unmodified image. However,
in 4c, the center 8x8 region clearly had been altered. In comparing 4b and 4c
with the original image, 4a, it was easier to see that both images had been
modified. However, the Deep Network’s output was still pretty good.

One common area where all the algorithms failed was images with sharp edges.
Looking at 4e and 4f, both algorithms blurred the edge between the truck’s cab
and trailer, while the original image 4d had a sharp line here. It is possible
that adding more filters to the initial layers of the networks would allow the
algorithms to better identify edges; however, this would likely have lead to
the over-fitting problems that plagued the deeper networks. This also
showcases how mean squared error can be a bad metric at times. If the
algorithms had output a straight line of red pixels on one side and white
pixels on the other, but shifted the line to the left or white of where it was
originally, the algorithm would have incurred a large MSE, but the image would
have looked realistic. Instead, by blurring the line somewhat, it ensures that
the pixel values are not too far off what they should be, and thus incurs a
lower MSE on average. I looked into alternative loss functions to mitigate
this issue, but did not find one which I felt worked as well as MSE.

(a) Original Image

(b) Deep Network

(c) Fully Connected Network

(d) Original Image

(e) Deep Network

(f) Fully Connected Network

Figure 4: Sample Network Outputs

##  6 Conclusion/Future Work

The goal of this project was to determine whether machine learning image
completion algorithms could be trained on small images and small datasets and
perform well. Figure 5 shows the MSE for the different algorithms on the test
set, where the models are tested at their lowest development loss. From this
chart and the development errors listed previously, Deep Network was the best
performing model, with a test MSE of .015.

Figure 5: Test Losses

After analyzing some of the test images output by Deep Network, many of them
looked realistic. However, images with certain defining features were almost
never predicted correctly (as shown in Figure 4). This problem plagued all the
networks and is one of the key benefits of using GANs. For instance, if an
image of a human face with the eyes removed was fed to Deep Network, it would
likely fill in the missing pixels with normal skin, as it performs a type of
weighted average of the pixels surrounding the removed values. However, a GAN
may have learned the distribution of human faces and recognize that it needs
eyes in those locations to look realistic. This issue is in part due to my
networks being trained on multiple classes of images, while GANs are often
only trained on one class.

While the Deep Network performed the best on this dataset, if I increased the
amount of training data, the deeper, more complex models would have likely
outperformed it. This highlights one of the key issues in deep learning: you
often need more data if you want to train a deeper model.

Given more time, I would like to experiment with this theory and try my models
on a larger dataset with more images and possibly larger images. ImageNet’s
64x64 RGB dataset would be perfect for this task and I think the deeper models
would perform very well. Using larger datasets, I would likely be able to have
more filters without risking over-fitting my training set. These additional
filters could hopefully learn more complex features such as sharp edges, so
that the algorithms would not make the same kinds of mistakes as they did on
CIFAR-10.

As well, I would like to experiment more with encoder-decoder networks. I was
only able to experiment with two different sized encodings, but I would like
to see how larger and small dimensional encoding affect performance. I
theorize that smaller encodings would reduce over-fitting and lead to
interesting results. As well, it would be interesting to build a model that
included encodings and had a skip connection with the input image, so that the
model would have an easier time maintaining important features like sharp
lines.

## References

[1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S.
Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G.
Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg,
D. Mane, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B.
Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F.
Viegas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng.
2015. TensorFlow: Large-scale machine learning on heterogeneous systems.
www.tensorflow.org.

[2] Amos, B. 2018. Image Completion with Deep Learning in TensorFlow.
bamos.github.io/2016/08/09/deep-completion/.

[3] V. Badrinarayanan, A. Kendall and R. Cipolla. 2017. SegNet: A Deep
Convolutional Encoder-Decoder Architecture for Image Segmentation. IEEE
Transactions on Pattern Analysis and Machine Intelligence, vol. 39, no. 12,
pp. 2481-2495, Dec. 1 2017. doi: 10.1109/TPAMI.2016.2644615.

[4] C. Barnes, E. Shechtman, A. Finkelstein, and D. Goldman. 2009. Patchmatch:
a randomized correspondence algorithm for structural image editing. TOG, 2009.

[5] Hays, J., and Efros, A. A. 2007. Scene completion using millions of
photographs. ACM Trans. Graph (SIGGRAPH 2007) 26, 3.

[6] A. Krizhevsky. 2010. Convolutional Deep Belief Networks on CIFAR-10.
www.cs.toronto.edu/ kriz/cifar.html.

[7] Long, J., Shelhamer, E., and Darrell, T. 2014. Fully convolutional
networks for semantic segmentation. CoRR, abs/1411.4038.