\section{Introduction}
Parametrically-driven systems and parametric resonance occur in many different
physical systems, ranging from the mechanical domain to the electronic,
microwave, electromechanic, optomechanic, and quantum domains.
In the mechanical domain we have Faraday waves \cite{faraday1831},
inverted pendulum stabilization, stability of boats, balloons, and parachutes \cite{ruby1996}.
A comprehensive review of applications in electronics and microwave cavities
spanning from the early twentieth century up to 1960 can be found in
Ref.~\cite{mumford1960}.
A few relevant recent applications, in
micro and nano systems, include quadrupole ion guides and ion traps \cite{paul90},
linear ion crystals in linear Paul traps designed as
prototype systems for the implementation of quantum computing
\cite{raizen1992ionic, drewsen1998large, kielpinski2002architecture},
magnetic resonance force microscopy \cite{dougherty1996detection},
tapping-mode force microscopy \cite{moreno2006parametric}, axially-loaded
microelectromechanical systems (MEMS) \cite{requa06electromechanical},
torsional MEMS \cite{turner98nature}.
In the quantum domain we could mention  wideband superconducting parametric
amplifiers \cite{eom2012nature} and squeezing
in optomechanical cavities below the zero-point motion \cite{szorkovszky2011}.

Parametric pumping has had many applications in the field of MEMS, which have
been used primarily as accelerometers, for measuring small forces and as ultrasensitive mass detectors
since the mid 80's \cite{binnig86}.
An enhancement to the detection techniques in MEMS was developed by Rugar and
Grütter \cite{rugar91} in the early 90's that uses mechanical parametric amplification (before transduction)
to improve the sensitivity of measurements.
This amplification method works by driving the parametrically-driven resonator
on the verge of parametric unstable zones.
They were
looking for means of reducing noise and increasing precision in a detector
for gravitational waves, when they experimentally found classical
thermomechanical quadrature squeezing, a phenomenon which is reminiscent of quantum squeezed states.
The classical version is characterized
by oscillating levels of the response of the parametric oscillator to noise at
the frequency of the parametric pump, in such a way that the product between
maximum and minimum output noise levels is constant. They observed that when
the pump was turned on, the noise increased in one quadrature, while on the
other it decreased. No theoretical model was proposed by them to explain the effect though.
Subsequently, DiFilippo et al. \cite{difilippo92prl}
and Natarajan et al. \cite{natarajan95prl} proposed
theoretical explanations for this noise squeezing phenomenon, but their models
did not treat noise directly in the equations of motion.

Here, we study  a parametrically-driven
oscillator in the presence of noise with the objective of understanding what
causes thermal squeezing and heating in the stable zone near the transition
line of the first parametric instability.
The one-degree of freedom model studied here may be applied for instance to the
fundamental mode of a doubly-clamped beam resonator that is axially loaded, in which case
the one degree of freedom represents the amount of deflexion of the middle of
the beam from the equilibrium position. The present model can also be applied
to the linear response of ac driven nonlinear oscillators to noise (such as
transversally-loaded beam resonators), see for example Ref.~\cite{almog2007noise}.

One of the objectives of the present investigation is to extend and improve on
recently obtained analytical quantitative estimates of the amount of quadrature noise squeezing
and heating in a parametrically-driven oscillator \cite{batista2011cooling}.
Here we use the Green's function approach, previously developed to solve the
Langevin equation, aligned with averaging techniques up to second order, to
obtain more precise analytical estimates of the thermal fluctuations in the
parametrically-driven oscillator with added noise.
We further show, using an approximate Floquet theory based on first and
second-order averaging approximations, that thermal squeezing and heating are
related to the onset of real-valued Floquet multipliers (FMs) with different
magnitudes. It is shown that one FM grows while one gets closer in parameter
space to the first transition line to instability while the other FM decreases.
As a consequence, one gets two different effective dissipation rates,
while at the same time the input power due to noise remains constant as the pump
amplitude is increased.
We show below that these effects account for most thermal squeezing and heating
observed.
Furthermore, first-order analytical estimates of heating and the amount of
squeezing are also provided.

The contents of this paper are organized as follows. In Sec.~(\ref{theory})
we present our theoretical model, in Sec.~(\ref{numerics}) we present and
discuss our numerical results, and in Sec.~(\ref{conclusions}) we draw our
conclusions.
\section{Theory}
\label{theory}
\noindent
The equation for the parametrically-driven  oscillator (in dimensionless format) is given by
the damped Matthieu's equation
\begin{equation}
    \ddot x+\omega^2_0x=-\gamma\dot x+F_p\cos(2\omega t)\;x,
    \label{parOsc}
\end{equation}
in which $\gamma$ and $F_p\sim O(\varepsilon)$, where $\varepsilon<<1$.
Since we want to apply the averaging method (AM) \cite{verh96, Guck83} to
situations in which we have detuning, it is convenient to rewrite
Eq.~(\ref{parOsc}) in a more appropriate form with the
notation $\Omega=\omega^2_0-\omega^2$, where we also have $\Omega\sim O(\varepsilon)$.
With this substitution we obtain
$\ddot x+\omega^2x=-\Omega x-\gamma\dot x+F_p\cos(2\omega t)\;x$.
We then rewrite this equation in the form
$\dot x = y$, $\dot y = -\omega^2x+f(x,y,t)$,
where $ f(x,y,t)= -\Omega x+F_p\cos(2\omega t)\;x-\gamma y$.
We now set the above equation in slowly-varying form with the transformation to a slowly-varying frame
\begin{equation}
\left(\begin{array}{c}
x \\ y
\end{array}
\right)
=
\left(
\begin{array}{cc}
    \cos\omega t & -\sin \omega t\\
    -\omega\sin\omega t & -\omega\cos\omega t
\end{array}
\right)
\left(\begin{array}{c}
{\mathcal U} \\ {\mathcal V}
\end{array}
\right)
\label{slow_trans}
\end{equation}
   and obtain
\begin{eqnarray}
\left(\begin{array}{c}
\dot {\mathcal U} \\ \dot {\mathcal V}
\end{array}
\right)
&=&
\left(
\begin{array}{cc}
    \cos \omega t & -\frac{1}{\omega}\sin\omega t\\
    -\sin\omega t & -\frac{1}{\omega}\cos\omega t\nonumber
\end{array}
\right)
\left(\begin{array}{c}
0\\  f(x,y,t)
\end{array}
\right)\nonumber\\
&=&
-\frac{1}{\omega}\left(\begin{array}{c} \sin\omega t f(x,y,t)\\\cos\omega t
f(x,y,t)\end{array}\right)\equiv F({\mathcal U}, {\mathcal V}, t)=DF({\mathcal U}, {\mathcal V}, t)\left(\begin{array}{c}
{\mathcal U} \\ {\mathcal V}
\end{array}
\right).
\label{slow_frame}
\end{eqnarray}

The components of the jacobian matrix $DF({\mathcal U}, {\mathcal V}, t)$ of the above flow is
given by
\begin{eqnarray}
DF_{11} &=&
\frac{-1}{2\omega}\left(
\gamma\omega  (1-\cos(2\omega
    t))-\Omega\sin(2\omega t) +\frac{F_p}{2}\sin(4\omega t)\right)\nonumber\\
DF_{12} &=&\frac{-1}{2\omega}\left(
    \Omega(1-\cos(2\omega t)) +\frac{F_p}{2}(1-2\cos(2\omega
    t)+\cos(4\omega t))+\gamma\omega\sin(2\omega t)\right)\nonumber\\
DF_{21} &=& \frac{-1}{2\omega}\left(
    -\Omega(1+\cos(2\omega t))
    +\frac{F_p}{2}[1+2\cos(2\omega t)+\cos(4\omega t)]+\gamma\omega \sin(2\omega
    t)\right)\nonumber\\
DF_{22} &=&\frac{-1}{2\omega}\left(\Omega\sin(2\omega
    t)+\gamma\omega\left(1+\cos(2\omega t)\right)-\frac{F_p}{2}\sin(4\omega t)
    \right)
\end{eqnarray}
After application of the AM to first order (in which, basically,
we filter out oscillating terms at $2\omega$ and $4\omega$ in the above
equation), we obtain
\begin{eqnarray}
    \dot { u} &=& \frac{-1}{2\omega}\left[\gamma\omega  u
    +\left(\Omega +\frac{F_p}{2}\right) v\right],\nonumber\\
    \dot { v} &=& \frac{-1}{2\omega}\left[\left(-\Omega +\frac{F_p}{2}\right) u+\gamma\omega v\right],
    \label{1stOrdAvg}
\end{eqnarray}
where the functions ${\mathcal U}(t)$ and ${\mathcal V}(t)$ are related to their slowly-varying
averages $u(t)$ and $v(t)$, respectively, by the transformation
\begin{equation}
\left(\begin{array}{c}
{\mathcal U} \\ {\mathcal V}
\end{array}
\right)=\left(\begin{array}{c}
 u +{\mathcal W}_1\\  v+{\mathcal W}_2
\end{array}
\right).
\end{equation}
According to the averaging theorem \cite{hol81, bat00}, the vector
${\mathcal W}$ obeys $\partial_t{\mathcal W}(u,v,t)=g (u,v,t)$, where the vector $g(u, v, t)$
corresponds to the explicitly-varying components of the right side of
Eq.~(\ref{slow_frame}).
Namely, we have
{
\begin{eqnarray}
&&\frac{\partial{\mathcal W}(u,v, t)}{\partial
t}=
\left(
\begin{array}{c}
g_1 (u, v, t) \\
g_2 (u,v, t)
\end{array}
\right)=\nonumber\\
&=&\frac{1}{2\omega}\left(
\begin{array}{c}
\left[(\Omega+F_p){v}+\omega\gamma{u}\right]\cos(2\omega
t)
+\left(\Omega{u}-\gamma\omega{v}\right)\sin(2\omega
t)-\frac{F_p}{2}\left[{v}\cos(4\omega
t)+{u}\sin(4\omega t)\right]\\
\left[(\Omega-F_p){u}-\omega\gamma{v}\right]\cos(2\omega
t)
-\left(\Omega{v}+\gamma\omega{u}\right)\sin(2\omega
t)-\frac{F_p}{2}\left[{u}\cos(4\omega t)-{v}\sin(4\omega
t)\right]
\end{array}
\right).\nonumber
\end{eqnarray}
}
Upon integration we find

\begin{flushleft}
${\mathcal W}(u,v, t)=\frac{1}{4\omega^2}\times$
\[
\left[
\begin{array}{ll}
\omega\gamma\sin(2\omega t)-\Omega\cos(2\omega t)+\frac{F_p}{4}\cos(4\omega t)
&
(\Omega+F_p)\sin(2\omega t)+\gamma\omega\cos(2\omega
t)-\frac{F_p}{4}\sin(4\omega t)
\\
(\Omega-F_p)\sin(2\omega t)+\gamma\omega\cos(2\omega
t)-\frac{F_p}{4}\sin(4\omega t) & -\omega\gamma\sin(2\omega t)+\Omega\cos(2\omega t)-\frac{F_p}{4}\cos(4\omega t)
\end{array}
\right]\]
$\times\left[\begin{array}{c} u \\ v\end{array}\right],$
\end{flushleft}
in which the integration constants are set to zero.
The averaging theorem \cite{Guck83} states that
these two sets of functions, namely $(u(t), v(t))$  and $({\mathcal U}(t),\, {\mathcal V}(t))$,
will be close to each other to order $O(\epsilon)$ during a time scale of $O(1/\epsilon)$ if they have initial conditions within an initial distance of $O(\epsilon)$.
So by studying the simpler averaged system, one may obtain very accurate
information about the corresponding more complex non-autonomous original system.
Using the transformations $ u(t)= e^{-\gamma t/2}\tilde u(t)$ and $v(t)=
e^{-\gamma t/2}\tilde v(t)$ in Eqs.~\eqref{1stOrdAvg}, we obtain
\begin{eqnarray}
    \dot{\tilde u} &=& \frac{-1}{2\omega}\left(\Omega+\frac{F_p}{2}\right)\tilde v,\nonumber\\
    \dot{\tilde v} &=&
    \frac{-1}{2\omega}\left(-\Omega+\frac{F_p}{2}\right)\tilde u.
    \label{eq_av1}
\end{eqnarray}
Upon integration of
Eqs.~(\ref{eq_av1}), one finds the solution
\begin{eqnarray}
u(t)&=&e^{-\gamma t/2}\left[u_0 \cosh (\kappa
t)+\frac{\beta -\delta}{\kappa}v_0\sinh(\kappa t)\right]\nonumber,\\
v(t)&=&e^{-\gamma t/2}\left[v_0 \cosh (\kappa
t)+\frac{\beta+\delta}{\kappa}u_0 \sinh(\kappa t)\right],
\label{1st-avg_sol}
\end{eqnarray}
where $\kappa=\sqrt{\beta^2-\delta^2}$, $\beta=-F_p/4\omega$, and
$\delta=\Omega/2\omega$. Hence, we find that the first parametric resonance, i.e. the boundary between
the stable and unstable responses, is given by
\begin{equation}
(\gamma\omega)^2=(F_p/2)^2-\Omega^2.
\label{eq_zones}
\end{equation}
This result is valid for $\omega\approx \omega_0$ even in the presence of
added noise. In Fig.~\ref{zones} we find very good agreement between the
boundary obtained from numerical integration of Eq.~(\ref{parOsc}) and the boundary given by the
averaging technique.

From Eqs. (\ref{slow_trans}) and (\ref{1st-avg_sol}) we obtain the approximate
fundamental matrix $\Phi(t)$ (also known as the time evolution operator)

\begin{eqnarray}
\left(\begin{array}{c}
x(t) \\ y(t)
\end{array}
\right)
&=&\Phi(t)\left(\begin{array}{c}
x_0 \\ y_0
\end{array}
\right)\nonumber\\
&\approx&
e^{-\gamma t/2}\times\left(
\begin{array}{cc}
    \cos\omega t & -\sin \omega t\\
    -\omega\sin\omega t & -\omega\cos\omega t
\end{array}
\right)
\left(\begin{array}{cc}
\cosh(\kappa\,t)  & -\frac{\beta-\delta}{\omega\kappa}\sinh(\kappa\,t)\\
\frac{\beta+\delta}{\kappa}\sinh(\kappa\,t)&-\frac{1}{\omega}\cosh(\kappa\,t) \\
\end{array}
\right)\left(\begin{array}{c}
x_0 \\ y_0
\end{array}
\right)\nonumber,
\end{eqnarray}
where $\Phi(0)=I$, in which $I$ is the identity matrix.
From Floquet theory \cite{verh96} we know that $\Phi(t)=P(t)e^{Bt}$, where
$P(t)$ is a periodic matrix with period $T=\pi/\omega$. We also know that $P(0)=I$.
The eigenvalues of $e^{BT}$ are known as the Floquet multipliers.
We rewrite the fundamental matrix in the following form
\begin{equation}
\Phi(t)
\approx
\left(
\begin{array}{cc}
    \cos\omega t & \frac{1}{\omega}\sin \omega t\\
    -\omega\sin\omega t & \cos\omega t
\end{array}
\right)
\left(\begin{array}{cc}
e^{-\gamma t/2}\cosh(\kappa\,t)  &
-\frac{\beta-\delta}{\omega\kappa}e^{-\gamma t/2}\sinh(\kappa\,t)\\
-\frac{\omega(\beta+\delta)}{\kappa}e^{-\gamma
t/2}\sinh(\kappa\,t)&e^{-\gamma t/2}\cosh(\kappa\,t) \\
\end{array}
\right).
\end{equation}
Hence, we notice that from the approximate solution of the fundamental matrix,
via first-order averaging, we can find the approximate Floquet multipliers.
They are given by
\begin{equation}
\lambda_\pm=-e^{-(\frac{\gamma}{2}\pm\kappa)T}.
\end{equation}

Further improvements can be made by going to second order averaging.
According to Ref.~\cite{bat00},
the second-order corrections are given by the time average of

\begin{eqnarray}
&&\overline{DF({\mathcal U}, {\mathcal V}, t){\mathcal W}({\mathcal U}, {\mathcal V})}=\nonumber\\
&=&\frac{-1}{8\omega^3}\left(
\begin{array}{cc}
-\gamma\omega F_p& -\gamma^2\omega^2-\Omega(\Omega+ F_p)-F_p^2/8\\
\gamma^2\omega^2+\Omega(\Omega-F_p) +F_p^2/8& \gamma \omega F_p
\end{array}
\right)
\left(
\begin{array}{c}
u \\ v
\end{array}
\right).
\nonumber
\end{eqnarray}
Hence, the second-order approximation to Eq.~\eqref{slow_frame} becomes
\begin{eqnarray}
    \dot { u} &=& \frac{-1}{2\omega}\left[\gamma(\omega - \frac{F_p}{4\omega}) u
    +\left(\Omega+\frac{F_p}{2}-\frac{\gamma^2}{4}-\frac{\Omega(\Omega+F_p)}{4\omega^2}-\frac{F_p^2}{32\omega^2}\right)
    v\right],\nonumber\\
    \dot { v} &=& \frac{-1}{2\omega}\left[\left(-\Omega
    +\frac{F_p}{2}+\frac{\gamma^2}{4}+\frac{\Omega(\Omega-F_p)}{4\omega^2}+\frac{F_p^2}{32\omega^2}\right)
    u+\gamma(\omega +\frac{F_p}{4\omega})v\right].
    \label{eq_av2}
\end{eqnarray}
The solution is given by
\begin{eqnarray}
u(t) &=& e^{-\gamma t/2}\left[\left(
\cosh(\xi t)+a\frac{\sinh(\xi t)}{\xi}\right)u_0+(b-c)\frac{\sinh(\xi t)}{\xi}
v_0\right],\nonumber\\
v(t) &=& e^{-\gamma t/2}\left[(b+c)\frac{\sinh(\xi t)}{\xi} u_0+\left(
\cosh(\xi t)-a\frac{\sinh(\xi t)}{\xi}\right)v_0\right],
\end{eqnarray}
with $\xi=\sqrt{a^2+b^2-c^2}$, $a=-\frac{\gamma\beta}{2\omega}$,
$b=\beta-\delta\beta/\omega$, and
$c=\delta-\gamma^2/(8\omega)-\delta^2/(2\omega)-\beta^2/(4\omega).$
Thus, We can write the fundamental matrix in second-order approximation as
\begin{equation}
\Phi(t)
\approx
e^{-\gamma t/2}\left(
\begin{array}{cc}
    \cos\omega t & \frac{1}{\omega}\sin \omega t\\
    -\omega\sin\omega t & \cos\omega t
\end{array}
\right)
\left(\begin{array}{cc}
\cosh(\xi\,t) +a\frac{\sinh(\xi\,t)}{\xi} &
-\frac{b-c}{\omega\xi}\sinh(\xi\,t)\\
-\frac{\omega(b+c)}{\xi}\sinh(\xi\,t)&\cosh(\xi\,t)-a\frac{\sinh(\xi\,t)}{\xi}
\end{array}
\right).
\nonumber
\end{equation}
After some simple algebraic operations, we find the Floquet multipliers
(eigenvalues of $\Phi(T)$) to be given by
\begin{equation}
\lambda_\pm=-e^{-(\gamma/2\pm\xi )T}.
\label{2ndOrderFM}
\end{equation}
Hence, the transition line to instability in second-order averaging is
given by
\begin{equation}
\gamma=2\xi.
\label{2ndOrderTransLine}
\end{equation}

\subsection{Green's function method}
\noindent
The equation for the Green's function of the parametrically-driven oscillator is
given by
\begin{equation}
    \left[\frac{\partial^2}{\partial t^2} +\omega_0^2+\gamma \frac{\partial}{\partial t}-F_p\cos(2\omega t)\right]\;G(t,t')=\delta(t-t').
    \label{green_eq}
\end{equation}
Since we are interested in the stable zones of the parametric oscillator,
for $t<t'$ $G(t,t')=0$ and by integrating the above equation near $t=t'$,
we obtain the initial conditions when $t=t'+0^+$, $G(t,t')=0$ and $\frac{\partial}{\partial t}G(t,t')=1.0$.
\subsubsection{1st-order averaging}
Although Eq. \eqref{green_eq} may be solved exactly by using Floquet
theory \cite{wies85}, one obtains very complex
 solutions.
Instead, we find fairly simple analytical approximations to the Green's
functions and, subsequently, to the statistical averages of fluctuations using
the averaging method.
From Eq.~\eqref{slow_trans} we obtain the approximate Green's function is
$G(t,t')= \cos(\omega t){\mathcal U}(t)-\sin(\omega t){\mathcal V}(t)= \cos(\omega t)[u(t)+{\mathcal W}_1(u(t), v(t), t)]-\sin(\omega t)[(v(t)+{\mathcal W}_2(u(t), v(t),
t)]$, the functions $u(t)$ and $v(t)$ are given by the solution of the system of
coupled differential equations \eqref{1st-avg_sol}, where the time $t$ is
replaced by $t-t'$ and the initial conditions set at $t=t'$ are given by
$u(t')=-\sin(\omega t')/\omega$ and $v(t')=-\cos(\omega t')/\omega$.
For simplicity we set ${\mathcal W}_1(t')={\mathcal W}_2(t')=0$.
We find the approximate Green's function to be
\begin{eqnarray}
G(t,t')&\approx& -\frac{e^{-\gamma (t-t')/2}}{\omega}\left[\cos(\omega
t)\left(\cosh(\kappa\,(t-t'))\sin(\omega t')+\frac{\beta-\delta}{\kappa}\sinh(\kappa\,(t-t'))\cos(\omega t')\right)\right.\nonumber\\ &&\left.-\sin(\omega
t)\left(\frac{\beta+\delta}{\kappa}\sinh(\kappa\,(t-t'))\sin(\omega
t')+\cosh(\kappa\,(t-t')) \cos(\omega t')\right)\right]\nonumber\\
&=&
\frac{e^{-\gamma (t-t')/2}}{\omega} \left\{
\cosh(\kappa\,(t-t'))\sin(\omega (t-t'))
+\frac{\delta}{\kappa}\sinh(\kappa\,(t-t'))\cos(\omega (t-t'))\right.\nonumber\\
&&\left.-\frac{\beta}{\kappa}\sinh(\kappa\,(t-t'))\cos(\omega
(t+t')) \right\},
\label{green_avg}
\end{eqnarray}
for $t>t'$ and $G(t,t')=0$ for $t<t'$. In the stable zone of the
parametrically-driven oscillator, when $|\beta|>|\delta|$ , we can rewrite the
Green's function replacing the initial conditions and using simplifying
trigonometrical identities. The change of variables $t'=t-\tau$ leads to
\begin{eqnarray}
G(t,t-\tau)&\approx& \frac{e^{-\gamma \tau/2}}{\omega} \left\{
\cosh(\kappa\,\tau)\sin(\omega \tau)
+\frac{\delta}{\kappa}\sinh(\kappa\,\tau)\cos(\omega \tau)\right.\nonumber\\
&&\left.-\frac{\beta}{\kappa}\sinh(\kappa\,\tau)\cos(\omega (2t-\tau)) \right\}.
\label{green_avg1}
\end{eqnarray}
\subsubsection{2nd-order averaging}
\label{2nd-orderGF}
The Green's function obtained by second-order averaging is given by
\begin{eqnarray}
&&G(t,t')\approx -\frac{e^{-\gamma
(t-t')/2}}{\omega}\times\nonumber\\
&&\left\{\cos(\omega t) \left[\left(
\cosh(\xi (t-t'))+a\frac{\sinh(\xi
(t-t'))}{\xi}\right)\sin(\omega t')+(b-c)\frac{\sinh(\xi (t-t'))}{\xi}
\cos(\omega t')\right]\right.\nonumber\\
&&\left.-\sin(\omega t)
\left[(b+c)\frac{\sinh(\xi (t-t'))}{\xi} \sin(\omega t')+\left( \cosh(\xi (t-t'))-a\frac{\sinh(\xi
(t-t'))}{\xi}\right)\cos(\omega t')\right]\right\}\nonumber,
\end{eqnarray}
for $t>t'$ and $G(t,t')=0$ for $t<t'$.
This can be rewritten in a shorter format as
\begin{eqnarray}
&&G(t,t-\tau)\approx \frac{e^{-\gamma\tau/2}}{\omega}\left\{
\cosh(\xi \tau)\sin(\omega\tau) +c\frac{\sinh(\xi \tau)}{\xi}
\cos(\omega\tau)\right.\nonumber\\
&&\left.-\frac{\sinh(\xi\tau)}{\xi}\left[a\sin(\omega
(2t-\tau))+b\cos(\omega (2 t-\tau)\right]\right\},
\label{green_avg2}
\end{eqnarray}
for $\tau>0$ and $G(t,t-\tau)=0$ for $\tau<0$.
We notice that this second-order approximate Green's function can be put in the
same format as the one of the first-order approximation given by
Eq.~\eqref{green_avg1}.
\begin{eqnarray}
&&G(t,t-\tau)\approx \frac{e^{-\gamma\tau/2}}{\omega}\left\{
\cosh(\xi \tau)\sin(\omega\tau) +c\frac{\sinh(\xi \tau)}{\xi}
\cos(\omega\tau)\right.\nonumber\\
&&\left.-\beta'\;\frac{\sinh(\xi\tau)}{\xi}\cos(\omega[2(t-t_0)-\tau])\right\},
\label{green_avg2n}
\end{eqnarray}
where $\beta'= \beta\sqrt{(1-\delta/\omega)^2+\gamma^2/(4\omega^2)}$,
$a=\beta'\sin(2\omega t_0)$, and $b=\beta'\cos(2\omega t_0)$.
In Eq.~\eqref{green_avg1} we have to replace $\kappa$ by $\xi$, $\delta$ by
$c$, $\beta$ by $\beta'$ and $t$ by $t-t_0$ to obtain the second-order
Green's function given by Eq.~\eqref{green_avg2n}.
Note that the Green's functions given here represent the
first two steps of a Green's function
renormalization procedure based on the averaging method.
\subsection{Thermal fluctuations}
We will now investigate the effect of noise on the parametric oscillator
\cite{batista2011cooling}.
We start by adding noise to Eq.~(\ref{parOsc}) and obtain
\begin{equation}
    \ddot x=-\omega^2_0 x-\gamma\dot x+F_p\cos(2\omega t)\;x+R(t),
    \label{DuffingNoise}
\end{equation}
where $R(t)$ is a random function that
satisfies the statistical averages
$\langle R(t)\rangle=0$
and $\langle R(t)R(t')\rangle= 2T\gamma\delta(t-t')$, according to the fluctuation-dissipation theorem \cite{kubo66}.
The temperature of the heat bath in which the oscillator (or resonator) is
embedded i $T$. Once we integrate these equations of motion we can show how
classical mechanical noise squeezing and heating occur.
We shall now review the analytical method developed in Refs.
\cite{batista2011cooling, batista2011signal} to study the parametrically-driven
oscillator with noise as given by Eq.~(\ref{DuffingNoise}).

Using the Green's function we obtain the solution $x(t)$ of
Eq.~(\ref{DuffingNoise}) in the presence of noise $R(t)$
\begin{align}
x(t)&= x_h(t)+\int_{-\infty}^{\infty}dt'\;G(t,t')R(t'),\\
\dot x(t)&=v(t)= v_h(t)+\int_{-\infty}^{\infty}dt'\;\frac{\partial}{\partial t}G(t,t')R(t'),
\end{align}
where $x_h(t)$ is the homogeneous solution, which in the stable zone decays
exponentially with time; since we assume the pump has been turned on for a long
time, $x_h(t)=0$.
By statistically averaging the  fluctuations of $x(t)$
we obtain
\begin{eqnarray}
\langle x^2(t)\rangle &=& \iint_{-\infty}^{\infty}dt'\;dt''G(t,t')G(t,t'')\langle R(t')R(t'')\rangle=2T\gamma\int_{0}^{\infty}d\tau\;G(t,t-\tau)^2,
\label{x2t_avg}\\
\langle v^2(t)\rangle
&=&2T\gamma\int_{-\infty}^{t}dt'\;\left[\frac{\partial}{\partial
t}G(t,t')\right]^2,
\label{v2_avg}
\end{eqnarray}
where $\tau=t-t'$.

 We notice that by varying the pump amplitude $F_p$ and the detuning
$\Omega$, we can create a continuous family of classical thermo-mechanical
squeezed states, generalizing the experimental results of Rugar and Gr\"utter
\cite{rugar91}.
An estimate of the time average of the statistically averaged thermal
fluctuations, when $|\beta|>|\delta|$, is given by
\begin{eqnarray}
 \overline{\langle x^2(t)\rangle}&=&
 \frac{2T\gamma}{\omega^2}\int_0^\infty e^{-\gamma\tau}
 \left\{\left[\cosh(\kappa\tau)\sin(\omega\tau)+\frac{\delta}{\kappa}\sinh(\kappa\tau)\cos( \omega\tau)\right]^2
 +\frac{\beta^2}{2\kappa^2}\sinh^2(\kappa\tau)\right\}\,d\tau\nonumber\\
&=&\frac{2T\gamma}{\omega^2} \left[I_1+I_2+I_3+I_4\right],
\label{x2_avg_dc2}
\end{eqnarray}
where the integrals are given by
\begin{eqnarray}
I_1&=&\frac{\beta^2}{2\kappa^2}\int_0^\infty e^{-\gamma\tau}\sinh^2(\kappa\tau)
d\tau=
\frac{\beta^2}{\gamma(\gamma^2-4\kappa^2)},\nonumber\\
I_2&=&\int_0^\infty e^{-\gamma\tau}\cosh^2(\kappa\tau) \sin^2(\omega\tau) d\tau\nonumber\\
&=&\frac{1}{2}\left\{\frac{1}{2\gamma}-\frac{\gamma}{2(\gamma^2+4\omega^2)}
+\frac{\gamma}{2(\gamma^2-4\kappa^2)}-\frac{1}{4}\mbox{Re}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}\right]\right\},\nonumber\\
I_3&=&\frac{\delta^2}{\kappa^2}\int_0^\infty e^{-\gamma\tau}\sinh^2(\kappa\tau)
\cos^2(\omega\tau) d\tau\nonumber\\ &=&\frac{1}{\kappa^2}\left\{\frac{\kappa^2}{\gamma(\gamma^2-4\kappa^2)}+\frac{1}{8}\mbox{Re}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}\right]-\frac{\gamma}{4(\gamma^2+4\omega^2)}\right\},\nonumber\\
\label{I_3}
I_4&=&\frac{\delta}{2\kappa}\int_0^\infty
e^{-\gamma\tau}\sinh(2\kappa\tau)\sin(2\omega\tau)
d\tau=\frac{1}{4}\mbox{Im}\left[
\frac{1}{\kappa(\gamma-2\kappa-2i\omega)}-\frac{1}{\kappa(\gamma+2\kappa-2i\omega)}\right].\nonumber
\end{eqnarray}

An estimate of the statistically averaged thermal fluctuations, when $|\beta|>|\delta|$, is given by

\begin{eqnarray}
\langle x^2(t)\rangle&\approx&\overline{\langle x^2(t)\rangle}+A_{2\omega}\cos(2\omega t)+B_{2\omega}\sin(2\omega t)+A_{4\omega}\cos(4\omega t)
+B_{4\omega}\sin(4\omega t),
\label{x2_avg2}
\end{eqnarray}
where
\begin{eqnarray}
A_{2\omega}&=& -\frac{4\beta T\gamma}{\omega^2}(K_1+K_2),\nonumber\\
B_{2\omega}&=& -\frac{4\beta T\gamma}{\omega^2}(K_3+K_4),\nonumber
\end{eqnarray}
with
\begin{eqnarray}
K_1&=&
\frac{1}{8}\mbox{Im}\left[ \frac{1}{\kappa(\gamma-2\kappa-2i\omega)}-\frac{1}{\kappa(\gamma+2\kappa-2i\omega)}\right],\nonumber\\
K_2&=&\frac{\delta}{\kappa^2}\left\{\frac{\kappa^2}{\gamma(\gamma^2-4\kappa^2)}-\frac{\gamma}{4(\gamma^2+4\omega^2)}+\frac{1}{8}\mbox{Re}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}\right]\right\},\nonumber\\
K_3&=&\frac{1}{8\kappa}\left[\frac{4\kappa}{(\gamma^2-4\kappa^2)}+\mbox{Re}\left( \frac{1}{\gamma-2\kappa-2i\omega}-\frac{1}{\gamma+2\kappa-2i\omega}\right)
\right],\nonumber\\
K_4&=&\frac{\delta}{8\kappa^2}\mbox{Im}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}-\frac{2}{\gamma-2i\omega}\right]\nonumber\\
&=&\frac{\delta\omega}{4\kappa^2}\left[\frac{1}{(\gamma+2\kappa)^2+4\omega^2}+\frac{1}{(\gamma-2\kappa)^2+4\omega^2}-\frac{2}{\gamma^2+4\omega^2}\right].\nonumber
\end{eqnarray}
The remaining coefficients of Eq.~(\ref{x2_avg2}) are given by
\begin{eqnarray}
A_{4\omega} &=&\frac{\beta^2 T\gamma}{4\omega^2\kappa^2}\mbox{Re}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}-\frac{2}{\gamma-2i\omega}\right],\nonumber\\
B_{4\omega} &=&\frac{\beta^2 T\gamma}{4\omega\kappa^2}\mbox{Im}\left[\frac{1}{\gamma-2\kappa-2i\omega}+\frac{1}{\gamma+2\kappa-2i\omega}-\frac{2}{\gamma-2i\omega}\right].\nonumber
\end{eqnarray}
Notice that when one gets close to the zone of instability, we obtain a far
simpler expression for the average fluctuations. It is given approximately by

\begin{equation}
\langle x^2(t)\rangle \approx \frac{2T}{\omega^2}
\left[ \beta^2+\frac{\gamma^2}{4}+\delta^2
\right]\frac{1}{\gamma^2-4\kappa^2}-\frac{4\beta
T}{\omega^2(\gamma^2-4\kappa^2)}\left[\delta\cos(2\omega
t)+\frac{\gamma}{2}\sin(2\omega t)\right].
\label{x2_approx}
\end{equation}
It is easy to verify that the minimum of the position fluctuation is given by
\[
\sigma_{\mbox{min}}^2\equiv\langle x^2(t)\rangle_{\mbox{min}} \approx
\frac{T}{2\omega^2}\frac{\sqrt{\delta^2+\gamma^2/4}-|\beta|}{\sqrt{\delta^2+\gamma^2/4}+|\beta|}
\]
and the maximum is given by
\[
\sigma_{\mbox{max}}^2\equiv{\langle x^2(t)\rangle}_{\mbox{max}} \approx
\frac{T}{2\omega^2}\frac{\sqrt{\delta^2+\gamma^2/4}+|\beta|}{\sqrt{\delta^2+\gamma^2/4}-|\beta|},
\]
such that
\begin{equation}
{\sigma}_{\mbox{min}} {\sigma}_{\mbox{max}}\approx
\frac{T}{2\omega^2}.
\label{squeezing}
\end{equation}
This is a verification that the classical phenomenon of thermal squeezing indeed
occurs near the transition line to instability.
These approximate expressions for squeezing are not valid with detuning
($\delta \neq 0$) when $\beta=0$, although they correctly predict no squeezing
in such situation and the equipartition theorem is correct to $O(\epsilon)$ if
$\Omega=O(\epsilon)$.

From Eqs. \eqref{green_avg2n} and \eqref{x2t_avg}  we obtain the dc
contribution to fluctuations in the second-order approximation
\begin{eqnarray}
\overline{\langle x^2(t)\rangle} &\approx&
\frac{2\gamma T}{\omega^2}\int_{0}^{\infty}d\tau\;e^{-\gamma\tau}\left\{
\left[\cosh(\xi \tau)\sin(\omega\tau) +c\frac{\sinh(\xi \tau)}{\xi}
\cos(\omega\tau)\right]^2+\beta'^2\frac{\sinh^2(\xi\tau)}{2\xi^2}\right\}\nonumber\\
&=&\frac{2T\gamma}{\omega^2} \left[\tilde I_1+\tilde I_2+\tilde
I_3+\tilde I_4\right],
\label{x2_avg2_approx}
\end{eqnarray}
where the integrals are given by
\begin{eqnarray}
\tilde I_1&=& \frac{\beta'^2}{\gamma(\gamma^2-4\xi^2)},\nonumber\\
\tilde I_2&=&\frac{1}{2}\left\{\frac{1}{2\gamma}-\frac{\gamma}{2(\gamma^2+4\omega^2)}
+\frac{\gamma}{2(\gamma^2-4\xi^2)}-\frac{1}{4}\mbox{Re}\left[\frac{1}{\gamma-2\xi-2i\omega}+\frac{1}{\gamma+2\xi-2i\omega}\right]\right\},\nonumber\\
\tilde I_3&=&\frac{1}{\xi^2}\left\{\frac{\xi^2}{\gamma(\gamma^2-4\xi^2)}+\frac{1}{8}\mbox{Re}\left[\frac{1}{\gamma-2\xi-2i\omega}+\frac{1}{\gamma+2\xi-2i\omega}\right]-\frac{\gamma}{4(\gamma^2+4\omega^2)}\right\},\nonumber\\
\label{I2_3}
\tilde I_4&=&\frac{1}{4}\mbox{Im}\left[
\frac{1}{\xi(\gamma-2\xi-2i\omega)}-\frac{1}{\xi(\gamma+2\xi-2i\omega)}\right].\nonumber
\end{eqnarray}
All the other first-order approximation results given at
Eqs.~(\ref{x2_avg2}-\ref{squeezing}) are the same in second-order
approximation, except for the parameter replacements given at
the end of subsubsection (\ref{2nd-orderGF}).
It is noteworthy to mention that the squeezing condition given at first-order
approximation in Eq.~\eqref{squeezing} is still valid in the second-order
approximation.
This implies that a renormalization procedure based on the averaging method
is possible for the Green's function of the parametric oscillator with
parameters set near the onset of the first instability zone.
\subsection{Energy balance}
It is important to verify that in the stable zone of
the parametric oscillator, on average, the power input from the external noise
and from the internal pump is balanced by the dissipated power.

 The instantaneous input power due to the additive noise is given by
\begin{equation}
P_{noise}(t)=\dot x(t)R(t)=
v_h(t)R(t)+\int_{-\infty}^{\infty}dt'\;\frac{\partial}{\partial
t}G(t,t')R(t')R(t).
\end{equation}
Hence, we obtain the statistically averaged  noise input power
\begin{equation}
\langle P_{noise}(t)\rangle=
\int_{-\infty}^{t}dt'\;\frac{\partial}{\partial
t}G(t,t')\langle R(t')R(t)\rangle=\gamma T \left[\frac{\partial}{\partial
t}G(t,t')\right]_{t'=t^-}=\gamma T.
\label{noise}
\end{equation}
The statistically-averaged dissipated power is given by
\begin{eqnarray}
\langle P_{diss}(t)\rangle &=& -\gamma\langle \dot x(t)^2\rangle=
-\gamma \int_{-\infty}^{t}dt'\int_{-\infty}^{t}dt''\;
\frac{\partial}{\partial t}G(t,t')\frac{\partial}{\partial
t}G(t,t'')\langle R(t')R(t'')\rangle\nonumber\\
&=&-2\gamma^2 T
\int_{-\infty}^{t}dt'\;\left[\frac{\partial}{\partial
t}G(t,t')\right]^2.
\label{dissPower}
\end{eqnarray}
The statistically-averaged pump power is given by
\begin{eqnarray}
\langle P_{pump}(t)\rangle &=& F_p\cos(2\omega t)\langle x(t)\dot
x(t)\rangle\nonumber\\
&=& F_p\cos(2\omega
t)\int_{-\infty}^{t}dt'\;\int_{-\infty}^{t}dt''\;G(t,t')\frac{\partial}{\partial
t}G(t,t'') \langle R(t')R(t'')\rangle\nonumber\\
&=&2\gamma T F_p\cos(2\omega
t)\int_{-\infty}^{t}dt'\;G(t,t')\frac{\partial}{\partial t}G(t,t')\nonumber\\
&=&
\gamma T F_p\cos(2\omega
t)\frac{\partial}{\partial t}\int_{-\infty}^{t}dt'\;G(t,t')^2\nonumber\\
&=&
\frac{F_p}{2}\cos(2\omega
t)\frac{\partial}{\partial t}\langle x(t)^2\rangle.
\end{eqnarray}
From Eq.~(\ref{x2_avg2}) we obtain approximately the time-averaged
statistically-averaged pump power.
 It is given by
\begin{equation}
\overline{ \langle P_{pump}(t)\rangle}\approx\frac{F_p\omega}{2}B_{2\omega}.
\end{equation}
Near the threshold to instability in first-order approximation it becomes
\[
\overline{ \langle P_{pump}(t)\rangle}
\approx \frac{F_p^2
T\gamma}{4\omega^2(\gamma^2-4\kappa^2)},
\]
while on second-order approximation it is
\[
\overline{ \langle P_{pump}(t)\rangle}\approx
\frac{F_p^2
T\gamma}{4\omega^2(\gamma^2-4\xi^2)}\sqrt{(1-\delta/\omega)^2+\gamma^2/(4\omega^2)},
\]

In the stable zone of the parametric oscillator, when the stationary point is reached,
one gets the energy balance on average, that is
\begin{equation}
\overline{\langle P_{noise}(t)\rangle}+\overline{\langle
P_{diss}(t)\rangle}+\overline{\langle P_{pump}(t)\rangle}=0.
\label{energyBalance}
\end{equation}
One can use the above expression to obtain the time average of the velocity
fluctuations
\[
\overline{\langle \dot x(t)^2\rangle}= \left(T+\overline{\langle
P_{pump}(t)\rangle}/\gamma\right)=\left(T+\frac{F_p\omega}{2\gamma}B_{2\omega}\right).
\]
By comparison of the above equation with Eq.~\eqref{x2_avg2_approx} one notices
that the equipartion of energy breaks down once pumping is on.
\section{Results and Discussion}
\label{numerics}
In Figs.~(\ref{fmw0.9}, \ref{fmw1}, \ref{fmw1.1}) we show the dependance of the
magnitude of the Floquet multipliers (FMs) on the pump amplitude $F_p$.
When the FMs given by Eq.~(\ref{2ndOrderFM}) become real, they branch off in two
magnitudes.
When this occurs one of the FMs ($\lambda_-$) becomes larger as $F_p$ is
increased, eventually becoming larger than one, when the system given by
Eq.~(\ref{parOsc}) becomes unstable, while the other FM ($\lambda_+$)
becomes smaller. From Eq.~(\ref{2ndOrderFM}) this implies into two effective dissipation rates,
along the stable manifolds of the Poincaré first-return map (with period
$\pi/\omega$, i.e. $\Phi(\pi/\omega)$).
This phenomenon causes both heating and squeezing of the parametric oscillator
with added noise, since the average power input due to thermal noise remains constant
as is shown in Eq.~(\ref{noise}) and one effective dissipation rate is
decreased.
It also causes quadrature thermal squeezing since in one direction,
the $\lambda_-$-stable manifold of $\Phi (\pi/\omega)$ in the phase space of $u$
and $v$, there is less effective dissipation and, consequently, more
fluctuations, while in another direction, the $\lambda_+$-stable manifold of
$\Phi (\pi/\omega)$ in the phase space of $u$ and $v$, there is more effective
dissipation and , consequently, less fluctuations.
This inbalance in the effective dissipation rates, we claim, is the main cause
of thermal squeezing in the parametric oscillator with added noise.
The dependance on pump amplitude of the effective dissipations can be seen on
Fig.~\ref{e-diss} (second-order averaging result) and on
Fig.~\ref{e-dissN} (Floquet theory numerical result).

In Fig.~\ref{greenfs} we show several squared Green's functions with initial
conditions spread out evenly in time during one period of the pump
($\pi/\omega$). They are vertically spaced only for clarity, since all of their assymptotes are
zero.
The Green's functions are shown squared because that is the way they contribute
to the thermal fluctuations in Eq.~(\ref{x2t_avg}).
One notices that the second-order approximation Green's functions yield a much
better approximation to the numerical Green's functions than the first-order
approximation.
This is specially evidenced the closer one gets to first transition line to
instability, a consequence of the fact that the second-order analytical
expression for the transition line given in Eq.~(\ref{2ndOrderTransLine}) is
a better approximation than the first-order expression in Eq.~(\ref{eq_zones}).

In Fig.~\ref{x2_t} we show a comparison between time series of the
fluctuation $\langle x^2(t)\rangle$ given by Eq. (\ref{x2t_avg}) in which the Green's functions are
given either by the second-order approximation expression from
Eq.~(\ref{green_avg2}) or by the numerical Floquet theory Green's functions.
We use several different Green's functions with negative detuning
($\omega=0.95$), in resonance ($\omega=1.0$) and positive detuning
($\omega=1.05$) all with the same pump amplitude $F_p=0.185$.
We observe that the numerical and approximate Green's functions are very
similar and that the squeezing amplitude and heating (proportional to the time
average of $\langle x^2(t)\rangle$) are very dependant on detuning from
resonance. From Fig. \ref{zones} one sees that the off resonance results are
slightly below the threshold of the strong heating and squeezing zone.
In Fig.~\ref{x2_floquet} we show again several time series of the fluctuation
$\langle x^2(t)\rangle$, but this time the pump amplitudes are chosen such that
the parameters are inside the heating zone (as given in Fig.~\ref{zones}) and
close to the transition line to instability. One sees then  considerably higher dynamical temperatures and squeezing
amplitudes for the detuned fluctuations than in the previous figure.

In Fig.~\ref{heating} we show a logarithmic plot of the dc component of
the mean-square displacement over the heat bath temperature.
Most of the heating occurs inside the heating zone in which the Floquet
multipliers are real, each one with a different amplitude, one that increases
and the other that decreases as the pump amplitude is increased.

In Fig.~\ref{tempcont} we show level sets in decibels of the dc
component of the mean-square displacement  over the heat bath
temperature in second-order approximation. One sees that most of the heating
occurs inside the heating zone in which the Floquet multipliers are real.
In Fig.~\ref{sqzAmp} we show level sets of the squeezing amplitude of
the mean-square displacement over the heat bath temperature in second-order
approximation. This attests that most of the squeezing also occurs inside the
heating zone as claimed before.

\section{Conclusions}
\label{conclusions}
Here, we studied  a parametrically-driven
oscillator with added noise with the objective of understanding what
causes heating and thermal squeezing in the stable zone near the
transition line of the first parametric instability.
We improved on our previous work \cite{batista2011cooling}
and obtained a more accurate expression for the Green's functions
and, consequently, obtained more precise estimates of the amount of
heating and quadrature noise squeezing  in the parametrically-driven oscillator.
Furthermore, we used an approximate Floquet theory,
based on first and second-order averaging approximations, to explain why
heating and thermal squeezing occur in the parametrically-driven oscillator
with added noise investigated here.
These phenomena are related to the onset of real-valued Floquet multipliers
(FMs) with different magnitudes. It was shown that as one FM grows while one gets closer
 (in parameter space) to the first transition line to instability the other FM
decreases. As a consequence, one gets two different effective dissipation rates,
while at the same time the input power due to noise remains constant as the pump
amplitude is increased.
We showed that these effects account for most thermal squeezing and
heating observed.
We showed as well that the second-order Green's function of the parametric
oscillator has the same form as the first-order Green's function, which implies
that with a simple parameter change given in the text all first-order results
for the amount of heating and thermal squeezing also apply, with increased
accuracy, to second-order approximation.
This indicates that a Green's function renormalization is very well feasible
if one goes to higher orders of approximation in the averaging method.

The one-degree of freedom model studied here may be applied for instance to the
dynamics of the fundamental mode of an axially loaded doubly-clamped beam
resonator. The present model could also be applied
to the linear response of ac driven nonlinear oscillators,
such as transversally-loaded beam resonators, to noise.
We note further that this model can be applied as well
to the investigation of the dynamics of ions in quadrupole RF ion guides or
traps \cite{paul90} in the presence of thermal noise.
Most importantly, the
notion that we can have two different effective dissipation rates, one for each stable Floquet multiplier,
could be used to create control schemes similar to those of the OGY method
\cite{ott1990controlling} to reduce the effects of noise even further in
parametric amplifiers and in RF ion traps.


\section{Introduction}

It has been a long standing problem how to properly define quantum mechanics
on a surface in three dimensional Euclidean space. On one hand, Dirac stressed
that in his \emph{Principles} on the canonical quantization assumption that
"is found in practice successful only when applied with the dynamic
coordinates and momenta referring to a Cartesian system of axes and not to
more general curvilinear coordinates." \cite{dirac1} On the other hand, there
is in textbooks a routine recipe proposed by DeWitt by hypothesizing the
quantum kinetic energy operator to be proportional to Laplace-Beltrami
operator $\Delta_{LB}$ on the surface, \cite{dewitt
\begin{equation}
T=-\frac{\hbar^{2}}{2m}\Delta_{LB}.
\end{equation}
Is there a way to start from Dirac to reach DeWitt? Certainly, Dirac put
forward a theory for systems of second-class constraints \cite{dirac2}\ which
really encompasses the DeWitt's hypothesis as a special case, \cite{1990} but
contains much more than what was expected.

When we use the tensor covariant and contravariant components and the Einstein
summation convention, the so-called standard parametrization $\mathbf{r
(q^{1},q^{2})$ of the 2D surface is given by,
\begin{equation}
\mathbf{r}(q^{1},q^{2})\mathbf{=(}\text{ }x(q^{1},q^{2}),y(q^{1
,q^{2}),z(q^{1},q^{2})\text{ }\mathbf{).} \label{standform
\end{equation}
In differential geometry, ($q^{1},q^{2}$) is generally denoted by $q^{\mu}$
and $q^{\nu}$ with lowercase greek letters $\mu,\nu$ taking values $1,2$, and
$\mathbf{r}^{\mu}=g^{\mu\nu}\mathbf{r}_{\nu}=g^{\mu\nu}{\partial_{\nu
}\mathbf{r}$ $=g^{\mu\nu}\partial\mathbf{r/}q^{\nu}$ with $g_{\mu\nu
}={\partial_{\mu}}\mathbf{r}\cdot{\partial_{\nu}}\mathbf{r}$ being the metric
tensor. At this point $\mathbf{r}$, $\mathbf{n=(}n_{x},n_{y},n_{z}\mathbf{)}$
is the normal and $M\mathbf{n}$ symbolizes the mean curvature vector field, a
geometric invariant. \cite{liu07} In physics, this two dimensional (2D)
surface can more realistically be considered as a 3D shell whose thickness is
negligible in comparison with the dimension of the whole system. Then, there
are two ways to performing the calculus on the surface: Explicitly, when the
2D curved surface is conceived as a limiting case of a curved shell of equal
thickness $d$, where the limit $d\rightarrow0$ is then taken, great
discrepancies present as firstly taking limit $d\rightarrow0$ then defining
the derivatives on the surface, and as firstly defining derivatives in bulk
then letting $d\rightarrow0$. The second order is named as the \emph{confining
}$\emph{procedure}$ for studying motion on 2D surface embedded in 3D.
\cite{jk,dacosta,FC,BP} This kind of exploration was initialized in 1971,
\cite{jk} fundamentally finished in 1981, \cite{dacosta} and with correct
inclusion of electromagnetic field in 2008 \cite{FC} etc. Remarkably, as the
\emph{confining procedure} is applied to the momentum operator $\mathbf{p
=-i\hbar\nabla$, we find that the resultant momentum on the surface is with
$M$ denoting the mean curvature, \cite{liuPRA,liu11}
\begin{equation}
\mathbf{p}=-i\hbar(\mathbf{r}^{\mu}\partial_{\mu}+M\mathbf{n}), \label{GM
\end{equation}
which was originally discovered in 2007 \cite{liu07} by an entirely
independent development on the quantization of the momentum on 2D surface
embedded in 3D flat space. This momentum corresponds to the so-called standard
parametrization $\mathbf{r}(q^{1},q^{2})$ of the 2D surface (\ref{standform})
in mathematics therefore should be preferable over other forms of momentum
such as the generalized momenta ($p_{q^{1}},p_{q^{2}}$) canonically conjugated
to parameters ($q^{1},q^{2}$). Paralleling to the confining procedure-induced
\emph{geometric potential } $V_{gp}\equiv-\hbar^{2}/2m(M^{2}-K)$ with $K$
being the gaussian curvature, \cite{Szameit,Schultheiss} we call (\ref{GM})
\emph{geometric momentum}. \cite{liu11} This scheme of building up quantum
mechanics on the surface echoes the historic comments of Dirac on the
canonical quantization in his \emph{Principles}. \cite{dirac1}

In 2010, with help of the femtosecond laser writing technology, the optical
analogue of the quantum \emph{geometric potential} is experimentally realized
and its experimental effects on optical wave packets constrained on curved
surfaces are demonstrated. \cite{Szameit} In 2012, the \emph{geometric
potential} effects on the electronic properties of materials such as
Tomononaga-Luttinger liquids are directly confirmed with an observation of the
\textit{in situ} high-resolution ultraviolet photoemission spectra of a
one-dimensional metallic $C_{60}$ polymer with an uneven periodic
peanut-shaped structure. \cite{epl} These two experimental verifications may
have influences on further developments of physics and mathematics for the 2D
curved surfaces, for the \emph{geometric momentum }and the\emph{ geometric
potential }are, upon two constant factors, the gradient and Laplacian
operator, respectively, as pointed out in Refs. \cite{liuPRA,liu11}.

The principal purpose of this study is to explicitly show that, with use of
the Monge parametrization of the 2D surface, the geometric momentum is
compatible with Dirac's theory for systems of second-class constraints all around.

\section{Geometric momentum with Cartesian variables ($x,y$)}

By Monge parametrization, we mean that a 2D surface given by the form
$z=f(x,y)$ where ($x,y,z$) are Cartesian variables. For a sphere of radius $r$
in $R^{3}$, we have the so-called standard form
\begin{equation}
\mathbf{r}(x,y)=(x,y,\sqrt{r^{2}-x^{2}-y^{2}}).
\end{equation}
The covariant derivatives $\mathbf{r}_{\mu}$ and contravariant derivatives
$\mathbf{r}^{\mu}$ can be easily computed and the results are respectively
\begin{equation}
\left(
\begin{array}
[c]{c
\mathbf{r}_{_{x}}\\
\mathbf{r}_{_{y}
\end{array}
\right)  =\left(
\begin{array}
[c]{lll
1, & 0, & -x/\sqrt{r^{2}-x^{2}-y^{2}}\\
0, & 1, & -y/\sqrt{r^{2}-x^{2}-y^{2}
\end{array}
\right)  ,
\end{equation
\begin{equation}
\left(
\begin{array}
[c]{c
\mathbf{r}^{x}\\
\mathbf{r}^{y
\end{array}
\right)  \equiv\left(
\begin{array}
[c]{c
g^{x\mu}\mathbf{r}_{\mu}\\
g^{y\mu}\mathbf{r}_{\mu
\end{array}
\right)  =\frac{1}{r^{2}}\left(
\begin{array}
[c]{lll
r^{2}-x^{2}, & -xy, & -x\sqrt{r^{2}-x^{2}-y^{2}}\\
-xy, & r^{2}-y^{2}, & -y\sqrt{r^{2}-x^{2}-y^{2}
\end{array}
\right)  .
\end{equation}
The normal $\mathbf{n}$ and the mean curvature $M$ are given by respectively
\begin{equation}
\mathbf{n=}\frac{1}{r}(x,y,\sqrt{r^{2}-x^{2}-y^{2}}),\text{ }M=-\frac{1}{r}.
\end{equation}
Then, the geometric momentum operators $p_{i}$ $(i=x,y,z)$ are
\begin{align}
p_{x}  &  =-i\hbar\frac{1}{r^{2}}(\left(  r^{2}-x^{2}\right)  \frac{\partial
}{\partial x}-xy\frac{\partial}{\partial y}-x),\label{px}\\
p_{y}  &  =-i\hbar\frac{1}{r^{2}}(-xy\frac{\partial}{\partial x}+\left(
r^{2}-y^{2}\right)  \frac{\partial}{\partial y}-y),\label{py}\\
p_{z}  &  =i\hbar\frac{\sqrt{r^{2}-x^{2}-y^{2}}}{r^{2}}(x\frac{\partial
}{\partial x}+y\frac{\partial}{\partial y}+1). \label{pz
\end{align}
In flat space, we take commutator $[x_{i},p_{j}]=i\hbar\delta_{ij}$ for
granted. But we can easily verify that one the sphere the correct results turn
out to be $\left[  x_{i},p_{j}\right]  =i\hbar\left(  \delta_{ij}-x_{i
x_{j}/r^{2}\right)  $ with use of (\ref{px})- (\ref{pz}). In next section, we
will show that we need the Dirac's theory for systems of second-class
constraints which accounts for this fact.

\section{Dirac's theory for systems of second-class constraints}

On the sphere in the Monge paramerization, the primary Hamiltonian $H_{p}$ is,
\cite{dirac2
\begin{equation}
H_{p}=\frac{p_{i}^{2}}{2m}+\lambda\left(  \sqrt{r^{2}-x^{2}-y^{2}}-z\right)
+up_{\lambda}, \label{ham
\end{equation}
where $\lambda$ is the Lagrangian multiplier enforcing the constrained of
motion on the surface, and $u$ is also a Lagrangian multiplier guaranteeing
that this Hamiltonian is defined on the symplectic manifold, and $p_{i}$
$\left(  i=x,y,z\right)  $ and $p_{\lambda}$ are respectively the canonical
momenta conjugate to variables $x_{i}$\ and $\lambda$. The Poisson bracket is
defined by.
\begin{equation}
\left\{  f,H_{p}\right\}  \equiv\frac{\partial f}{\partial x_{i}
\frac{\partial H_{p}}{\partial p_{i}}+\frac{\partial f}{\partial\lambda
\frac{\partial H_{p}}{\partial p_{\lambda}}-(\frac{\partial f}{\partial p_{i
}\frac{\partial H_{p}}{\partial x_{i}}+\frac{\partial f}{\partial p_{\lambda
}\frac{\partial H_{p}}{\partial\lambda}). \label{poisson
\end{equation}
The equations of motion for ($x,y,z,\lambda$) are given by
\begin{equation}
p_{i}=m\dot{x},\text{ }p_{\lambda}=0. \label{plamb
\end{equation}
The primary constraint is then
\begin{equation}
\varphi_{1}=p_{\lambda}\approx0,
\end{equation}
hereafter symbol "$\approx$" implies a weak equality. After all calculations
are finished, weak equality takes back the strong one. The secondary
constraints (not confusing with second-class constraints) are then determined
by
\begin{equation}
\left\{  \varphi_{i},H_{p}\right\}  \approx0.
\end{equation}
And the complete secondary constraints are
\begin{align}
\varphi_{2}  &  =\sqrt{r^{2}-x^{2}-y^{2}}-z\approx0,\label{fir}\\
\varphi_{3}  &  =\frac{xp_{x}+yp_{y}}{m\sqrt{r^{2}-x^{2}-y^{2}}}+\frac{p_{z
}{m}\approx0,\label{sec}\\
\varphi_{4}  &  =\frac{\left(  r^{2}-x^{2}-y^{2}\right)  (p_{x}^{2}+p_{y
^{2}+p_{z}^{2})-r^{2}\sqrt{r^{2}-x^{2}-y^{2}}m\lambda}{m^{2}\left(
r^{2}-x^{2}-y^{2}\right)  ^{3/2}}\approx0,\label{thi}\\
\varphi_{5}  &  =\frac{p_{z}(p_{x}^{2}+p_{y}^{2}+p_{z}^{2})+r^{2}m^{2}u
{m^{3}\left(  r^{2}-x^{2}-y^{2}\right)  }\approx0. \label{for
\end{align}
Eqs. (\ref{thi}) and (\ref{for}) determine the Lagrangian multipliers
$\lambda$ and $u$ respectively. With introduction of the Dirac bracket instead
of the Poisson one for the canonical variables $A$ and $B,$\cite{dirac2
\begin{equation}
\left\{  A,B\right\}  _{D}=\left\{  A,B\right\}  -\left\{  A,\varphi_{\alpha
}\right\}  C_{\alpha\beta}^{-1}\left\{  \varphi_{\beta},B\right\}  ,
\end{equation}
which the matrix elements $C_{\alpha\beta}\left(  \alpha,\beta=1,2,3,4\right)
$ is defined by
\begin{equation}
C_{\alpha\beta}=\left\{  \varphi_{\alpha},\varphi_{\beta}\right\}  ,
\end{equation}
the positions $x_{i}$ and the momenta $p_{i}$ satisfy following Dirac bracket,
\cite{1968,1983,1985,1990,1992,1996-1,1996-2,1998,2000,2009}
\begin{align}
\left\{  x_{i},x_{j}\right\}  _{D}  &  =0,\left\{  x_{i},p_{j}\right\}
_{D}=\delta_{ij}-\frac{x_{i}x_{j}}{r^{2}},\nonumber\\
\left\{  p_{i},p_{j}\right\}  _{D}  &  =-\frac{1}{r^{2}}(x_{i}p_{j}-x_{j
p_{i}), \label{xp
\end{align}
and other Dirac brackets between $x_{i}$ and $p_{j}$ vanish. The equation of
motion is in genera
\begin{equation}
\dot{f}=\left\{  f,H_{p}\right\}  _{D}, \label{geneqm
\end{equation}
from which we have for $x_{i}$ and $p_{i}$ respectively
\begin{align}
\dot{x}_{i}  &  =\left\{  x_{i},H_{p}\right\}  _{D}=\frac{p_{i}}{m
,\label{xh}\\
\dot{p}_{i}  &  =\left\{  p_{i},H_{p}\right\}  _{D}=-\frac{x_{i}p_{i}^{2
}{mr^{2}}. \label{ph
\end{align}
Note that in these calculations (\ref{geneqm}), (\ref{xh}) and (\ref{ph})
where the constraints are of second-class, we need to deal with $H$ instead of
$H_{p}$ for we have
\begin{equation}
\dot{f}=\left\{  f,H_{p}\right\}  _{D}=\left\{  f,H\right\}  _{D}.
\end{equation}

In quantum mechanics, the Hamiltonian is, \cite{jk,dacosta}
\begin{align}
H  &  =-\frac{\hbar^{2}}{2m}\nabla^{2}=-\frac{\hbar^{2}}{2m}\frac{1}{\sqrt{g
}\partial_{\mu}\sqrt{g}g^{\mu\nu}\partial_{\upsilon}+V_{gp}\nonumber\\
&  =-\frac{\hbar^{2}}{2m}(\frac{r^{2}-x^{2}}{r^{2}}\frac{\partial^{2
}{\partial x^{2}}-\frac{2xy}{r^{2}}\frac{\partial}{\partial x}\frac{\partial
}{\partial y}-\frac{2x}{r^{2}}\frac{\partial}{\partial x}-\frac{2y}{r^{2
}\frac{\partial}{\partial y}+\frac{r^{2}-y^{2}}{r^{2}}\frac{\partial^{2
}{\partial y^{2}})+V_{gp}\text{ } \label{h
\end{align}
where the factor $g\equiv\det(g_{\mu\upsilon})$ is the determinant of the
matrix $g_{\mu\upsilon}$, and the geometric potential $V_{gp}=0$ for
$M^{2}-K=0$ on sphere. The quantum commutator $\left[  A,B\right]  $ of two
variables $A$ and $B$ is attainable by direct correspondence of the Dirac
bracket as $\left[  A,B\right]  /i\hbar\rightarrow\left\{  A,B\right\}  _{D}$,
and the fundamental commutators are
\begin{align}
\left[  x_{i},x_{j}\right]   &  =0,\label{dxx}\\
\left[  x_{i},p_{j}\right]   &  =i\hbar\left(  \delta_{ij}-\frac{x_{i}x_{j
}{r^{2}}\right)  ,\label{dxp}\\
\left[  p_{i},p_{j}\right]   &  =-\frac{i\hbar}{r^{2}}(x_{i}p_{j}-x_{j}p_{i}),
\label{dpp
\end{align}
There is no operator ordering problem in the right-hand side of Eq.
(\ref{dpp}) because the commutator must satisfy the Jacobian identity. It is
easily to verify that operators $p_{i}$ (\ref{px})-(\ref{pz}) satisfy
relations (\ref{dxx})-(\ref{dpp}). The second category of the fundamental
commutators is given by quantization of (\ref{xh}) and (\ref{ph}
\begin{align}
\left[  x_{i},H\right]   &  =i\hbar\frac{p_{i}}{m},\label{qxh}\\
\left[  p_{i},H\right]   &  =-i\hbar\frac{x_{i}H+Hx_{i}}{mr^{2}}. \label{pH
\end{align}
Strikingly, the geometric momentum (\ref{px})-(\ref{pz}) satisfies all
commutators (\ref{dxp})-(\ref{pH}) above, not only (\ref{dxp})-(\ref{dpp}). As
pointed out in Ref. \cite{liuPRA}, the usual canonical momentum $p_{\theta}$
violates the fundamental commutator (\ref{pH}).

\section{The self-adjointness of the Geometric Momentum}

By a self-adjoint operator, we mean that all eigenvalues of it are real, and
eigenfunctions corresponding to distinct eigenvalues are mutually orthogonal
and they form a complete set. But direct demonstration of the self-adjointness
of the geometric momentum (\ref{px})-(\ref{pz}) is relatively difficult. With
variable transform $(x,y,z)\rightarrow(r,\theta,\varphi)$ with $\theta
\in(0,\pi),\varphi\in(0,2\pi)$ as,
\begin{equation}
\mathbf{r}\left(  \theta,\varphi\right)  =\left(  r\sin\theta\cos\varphi
,r\sin\theta\sin\varphi,r\cos\theta\right)
\end{equation}
being made, the task becomes easy. The geometric momentum operators in terms
of ($\theta,\varphi$) take following forms, \cite{liu07,liu03,liu10
\begin{align}
p_{x}  &  =-\frac{i\hbar}{r}\left(  \cos\theta\cos\varphi\frac{\partial
}{\partial\theta}-\frac{\sin\varphi}{\sin\theta}\frac{\partial}{\partial
\varphi}-\sin\theta\cos\varphi\right)  ,\\
p_{y}  &  =-\frac{i\hbar}{r}\left(  \cos\theta\sin\varphi\frac{\partial
}{\partial\theta}-\frac{\cos\varphi}{\sin\theta}\frac{\partial}{\partial
\varphi}-\sin\theta\sin\varphi\right)  ,\\
p_{z}  &  =-\frac{i\hbar}{r}\left(  -\sin\theta\frac{\partial}{\partial\theta
}-\cos\theta\right)  .
\end{align}
Their complete solutions to eigenvalue equations $p_{j}(\theta,\varphi
)\psi_{p_{j}}(\theta,\varphi)=p_{j}\psi_{p_{j}}(\theta,\varphi)$ can be easily
determined. To note that the eigenvalues $p_{j}$ on the right hand side of
these equations differ from the operators $p_{j}(\theta,\varphi)$ on the left
hand side. Explicitly, the solutions are, \cite{liu12}
\begin{align}
\psi_{p_{x}}\left(  \theta,\varphi\right)   &  =f_{x}\left(  \frac
{1+\sin\theta\cos\varphi}{1-\sin\theta\cos\varphi}\right)  ^{irp_{x}/2\hbar
}\frac{\sqrt{\cos\theta\sin\theta\sin\varphi}}{1-\sin^{2}\theta\cos^{2
\varphi},\label{f1}\\
\psi_{p_{y}}\left(  \theta,\varphi\right)   &  =f_{y}\left(  \frac
{1+\sin\theta\sin\varphi}{1-\sin\theta\sin\varphi}\right)  ^{irp_{y}/2\hbar
}\frac{\sqrt{\cos\theta\sin\theta\cos\varphi}}{1-\sin^{2}\theta\sin^{2
\varphi},\label{f2}\\
\psi_{p_{x}}\left(  \theta,\varphi\right)   &  =C_{z}\left(  \cot\frac{\theta
}{2}\right)  ^{irp_{z}/\hbar}\frac{1}{\sin\theta}, \label{f3
\end{align}
where $f_{x}$ and $f_{y\text{ }}$ are two arbitrary functions of the same
variable $\tan\theta\sin\varphi$. In terms of variables $(x,y,z)$, we have
from (\ref{f1})-(\ref{f3})
\begin{align}
\psi_{p_{x}}\left(  x,y\right)   &  =f_{x}\left(  \frac{r+x}{r-x}\right)
^{irp_{x}/2\hbar}\frac{r}{r^{2}-x^{2}}\sqrt{y\sqrt{r^{2}-x^{2}-y^{2}}},\\
\psi_{p_{y}}\left(  x,y\right)   &  =f_{y}\left(  \frac{r+y}{r-y}\right)
^{irp_{y}/2\hbar}\frac{r}{r^{2}-y^{2}}\sqrt{x\sqrt{r^{2}-x^{2}-y^{2}}},\\
\psi_{p_{z}}\left(  x,y\right)   &  =C_{z}\left(  \frac{r+\sqrt{r^{2
-x^{2}-y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)  ^{irp_{z}/\hbar}\frac{r
{\sqrt{x^{2}+y^{2}}},
\end{align}
where $f_{x}$ and $f_{y\text{ }}$ are two arbitrary functions of variable
$y/\sqrt{r^{2}-x^{2}-y^{2}}$. One can then verify that $p_{j}(x,y)\psi_{p_{j
}(x,y)=p_{j}\psi_{p_{j}}(x,y)$ are satisfied with geometric momentum of form
(\ref{px})-(\ref{pz}).

\section{Remarks and summary}

Two dimensional surface can be considered as three dimensional shell whose
thickness is negligible in comparison with the dimension of the whole system.
We can study the quantum mechanics on surface by first formulating it in the
bulk, and then taking the limit of vanishing thickness, the gradient operator
and the Laplace operator originally defined in flat space converges to the
geometric ones. The presence of the geometric momentum and geometric potential
well reflects the Dirac's penetrating insight into the canonical quantization.

On the two dimensional sphere embedded\ in three dimensional flat space, the
geometric momentum in the Monge paramerization is extensively explored in this
paper. The apparent commutator $[x_{i},p_{j}]=i\hbar\delta_{ij}$ does not hold
true any more, and we must resort to the Dirac's theory on second-class
constrained motion. The correct results turn out to be $\left[  x_{i
,p_{j}\right]  =i\hbar\left(  \delta_{ij}-x_{i}x_{j}/r^{2}\right)  $. This
geometric momentum is geometric invariant under parameters transformation, and self-adjoint.

\begin{acknowledgments}
This work is financially supported by National Natural Science Foundation of
China under Grant No. 11175063.
\end{acknowledgments}


\section{Introduction} \label{sec:intro}

Boosting, and in particular AdaBoost \cite{schapire2012,freund1996,freund1997}, is an effective method of aggregating classifiers.
AdaBoost produces a reliable binary classifier and often avoids overfitting.
Nevertheless, it can be sensitive to ``noisy'' data and may severely underperform as a result.
In this paper, we embed binary classification in a Bayesian model and show how it interfaces with the boosting paradigm.
With this model, we can address the vulnerability to noise in a principled way.

Real-world data will almost always include noise, even with binary labels.
In the U.S. Presidential election of 2000, the country was kept in suspense for more than a month while votes were recounted in the state of Florida.
During the recount it had become apparent that the use of the ``butterfly ballot'' had confused voters \cite{toobin2002}.
Probabilistically, we can model a confused voter as one who casts a vote that is \textit{independent} of his/her actual intention.
These votes---borne out of confusion---are considered ``noisy'' and any attempt to learn a voter-to-vote connection, e.g., via boosting, becomes difficult.
However, this does not preclude the extraction of important noise information.
If we can \textit{detect and quantify} the noise properties of a given dataset, then it should be reflected in our expectations of constructing a good classifier.

In addressing label noise, we have chosen to interpret aggregating classifiers in a fully-Bayesian model.
Once in place, this model lets us incorporate additional latent variables to account for noise.
In our context, noise means that the true label is ignored and randomly reassigned, i.e., it may be inverted.
A learning algorithm such as AdaBoost is sensitive to this type of label perturbation because it focuses on the examples that pose a greater difficulty of classification.
Using this augmented model, we construct an algorithm that performs approximate inference of the posterior distribution associated with the latent variables.
Although the intent is inference, the algorithm is able to produce a binary classifier accompanied by noise statistics that reflect the quality of the learned classifier.
We also show that the algorithm---in its simplest form---reduces to a smoothed version of AdaBoost.

In developing a Bayesian model for aggregating binary classifiers, we begin with the logistic regression model proposed by \cite{friedman2000}.
Given a set of base classifiers, the latent variables of the model are the weights placed on these base classifiers.
We then introduce variables to account for label perturbations.
Finally, we use variational inference to estimate the posterior distributions.

Our ideas lead to a new boosting-like algorithm called \textit{VIBoost}---boosting stemming from \underline{v}ariational \underline{i}nference.
AdaBoost employs a greedy search for incorporating new base classifiers.
Similarly, in VIBoost each main-loop iteration introduces a new base classifier, which induces a new model.
With this new model, variational inference is applied using previous values for a warm start.
In the process, noise statistics are cultivated.
Our experiments reveal that VIBoost performs on par with AdaBoost and supplies meaningful characterizations of the label perturbations.

Much has been done to cast boosting in a statistical setting.
Friedman~et~al.~\cite{friedman2000} leveraged the logistic regression model and then used a functional gradient to derive a boosting update.
Collins~et~al.~\cite{collins2002} used information geometry to derive AdaBoost and algorithms emerged with exponential and logistic loss objectives.
Lebanon~\&~Lafferty~\cite{lebanon2002} solidified the relationship between AdaBoost and maximum likelihood via duality.
These ideas led to a Bayesian perspective of boosting and provided a way to incorporate prior knowledge \cite{schapire2002}.

There have been many approaches for handling noise.
For example, Servedio addressed label noise in a PAC learning framework and developed SmoothBoost \cite{servedio2003}.
Through a statistical formulation, Krause~\&~Singer~\cite{krause2004} addressed noise in the context of symmetric, random label inversions, and devised algorithms to alleviate the resulting adverse effects.
In one of these algorithms they used expectation maximization to construct a classifier while simultaneously updating a noise parameter.
Building on this work, we use variational inference and address label noise in the process.

The paper is organized as follows:
the initial groundwork for the Bayesian model is given in \S\ref{sec:core}.
In \S\ref{sec:logistic} we introduce two probability distributions that will play a role in the model.
The proposed model is presented in \S\ref{sec:model} and variational inference is applied in \S\ref{sec:vi}.
We discuss the connection to AdaBoost in \S\ref{sec:adaboost}.
We give experimental results in \S\ref{sec:exp} and we conclude in \S\ref{sec:con}.

\section{The Core Model} \label{sec:core}

In the binary classification problem we are given a set of $N$ labeled examples $\{(\mathbf{x}_n,y_n)\}_{n=1}^N$.
Each example is an element of some space $\mathcal{X}$ and the labels are elements of $\{-1,+1\}$.
In addition to the labeled examples, we also have a set of $M$ base classifiers $\mathcal{F} = \{f_1,\ldots,f_M\}$.
Each element of $\mathcal{F}$ is a function that maps $\mathcal{X}$ to $\{-1,+1\}$.
Additionally, we assume that (a) $h \in \mathcal{F} \Rightarrow -h \notin \mathcal{F}$ and (b) $h_1, h_2 \in \mathcal{F} \Rightarrow \exists \ i \in \{1,\ldots,N\}$ such that $h_1(\mathbf{x}_i) \neq h_2(\mathbf{x}_i)$.
These assumptions ensure a finite number of classifiers and prevent identifiability problems.

For a fixed $\mathbf{x} \in \mathcal{X}$, suppose the logarithm of the ``$+1$''-to-``$-1$'' label odds is given by $F(\mathbf{x}) = \log \frac{ p(y=+1 \mid \mathbf{x}) }{ p(y=-1 \mid \mathbf{x}) }$.
Thus, we can form the conditional probability mass function for the labels as ${ p(y \mid \mathbf{x}, F) = \frac{ 1 }{ 1 + \exp(-y F(\mathbf{x})) } }$.

A label sampled in this way shall be called a \textit{true} label.
Logistic regression models the spatially-variant log-odds-ratio as a weighted sum over all base classifiers, i.e.,  $F(\mathbf{x}) = \sum_{m=1}^M c_m f_m(\mathbf{x})$.

Consider a model defined by the following generative process:
\begin{samepage}
\begin{ENUM}
  \item Draw $c_m \stackrel{\text{iid}}{\sim} \mathcal{P}_C$ \ ($m=1,\ldots,M$) \ .
  \item Construct $F = \sum_{m=1}^M c_m f_m$ \ .
  \item Draw $\mathbf{x}_{1:N}$ independently according to some distribution over $\mathcal{X}$ \ .
  \item Draw $y_n \in \{-1,+1\}$ independently according to ${  p(y_n \mid \mathbf{x}_n, F) = \frac{1}{1 + \exp[-y_n F(\mathbf{x}_n) ]} }$ \ .

\end{ENUM}
\end{samepage}
The graphical model is shown in Figure~\ref{fig:core}.

The latent variables are the base classifier weights $c_{1:M}$.
Using the labeled examples, we seek the posterior distribution over the weights.
A similar approach was posed by Minka \cite{minka2001} with the Bayes Point Machine \cite{herbrich2001}.
In contrast with our work, the author considered (i) a linear classifier without the notion of base classifiers, (ii) expectation propagation as opposed to variational inference, and (iii) a Gaussian prior for the weights.

The posterior distribution over the weights reflects a compromise of the observed data (\,$\{\mathbf{x}_n,y_n\}_{n=1}^N$\,) with our prior beliefs (\,$\mathcal{P}_C$\,).
It also has the potential of yielding a classifier via the $M$-dimensional mean or mode, for example.
Combining prior beliefs with observed data is made easier through conjugacy, which is how we propose a distribution for $\mathcal{P}_C$.
This is the subject of the next section.

\begin{figure}[t!]
\begin{minipage}[b]{0.48\linewidth}
\centering
                \includegraphics[width=0.50\columnwidth]{coregraphicalmodel8.pdf}
                \caption{The core graphical model for the boosting problem. Each label depends on the example and log-odds-ratio function. The only latent variables are the base classifier weights.}
                \label{fig:core}
\end{minipage}
\hspace{0.4cm}
\begin{minipage}[b]{0.48\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{logoneplusz.pdf}
                \caption{For a given versatile logistic with unit multiplicities, the negative logarithm of one of the product terms is shown above ($\beta\!=\!1$ , $\gamma\!=\!2$). Each curve of this form is tightly lower-bounded by a piecewise-linear function with one knot at $z\!=\!\gamma$. The slopes are $0$ and $\beta$.}
                \label{fig:log1pez}
\end{minipage}
\end{figure}

\section{The Versatile Logistic \& Binary Logistic Distributions} \label{sec:logistic}

We use two conjugate distributions to specify the model described in Figure~\ref{fig:core}.
The first distribution is used as a prior for the weights $c_{1:M}$, and the second is associated with label generation.
For vectors ${\boldsymbol{\beta}},{\boldsymbol{\gamma}} \in {\mathbb{R}}^K$ and ${\boldsymbol{\mu}} \in {\mathbb{R}}^K_+$, we define the density over the reals
\begin{align}
  p(z) &\propto \textstyle{\prod_{k=1}^K} \left( \frac{1}{1+\exp[\beta_k (z - \gamma_k)]} \right)^{\mu_k}
  \label{eq:vlog}
\end{align}
to be the \textit{Versatile Logistic Distribution}---written $\text{v-Log}({\boldsymbol{\beta}},{\boldsymbol{\gamma}},{\boldsymbol{\mu}})$--- with \textit{slope vector} ${\boldsymbol{\beta}}$, \textit{knot vector} ${\boldsymbol{\gamma}}$, and \textit{multiplicity vector} ${\boldsymbol{\mu}}$.
Figure~\ref{fig:log1pez} provides the motivation behind this nomenclature.
Define $\bar{\mathbf{u}} \triangleq \left[ \begin{smallmatrix} +1 \\ -1 \end{smallmatrix} \right] \in {\mathbb{R}}^2$.
A familiar density is $\text{v-Log}(\bar{\mathbf{u}},{\mathbf{0}},{\mathbf{1}})$, which is a logistic distribution.

The density described in \eqref{eq:vlog} is valid if and only if there exists both a positive and negative slope with corresponding positive multiplicity.
Consequently, we must have $K \geq 2$.
Additionally, this distribution is unimodal, so it is reasonable to estimate its mean with an approximate mode.
We prove these facts in \S\ref{supp:sec:vlog}.
The product represented in \eqref{eq:vlog} relates to a Product of Experts \cite{hinton2002}; however, each factor by itself does not correspond to a valid density.

We now define a probability mass function for the binary random variable $Y$ taking values in $\{-1,+1\}$.
For scalars $z$, $\beta$, and $\gamma$ we define
\begin{align}
  p(y) &= \frac{1}{1 + \exp[-y \beta (z - \gamma) ]}
  \label{eq:binarylogistic}
\end{align}
to be the corresponding \textit{Binary Logistic Distribution}, written $\text{b-Log}(z,\beta,\gamma)$.
In comparing \eqref{eq:binarylogistic} to label generation in our model, we see that $\beta$ and $\gamma$ encode base classifier information.
The versatile logistic and binary logistic are conjugate in the following way:
if $z \sim \text{v-Log}({\boldsymbol{\beta}},{\boldsymbol{\gamma}},{\boldsymbol{\mu}})$ and $y_n | z \sim \text{b-Log}(z,\theta_n,\phi_n)$---drawn independently for $n=1,\ldots,N$---then the posterior of $z$ given $y_{1:N}$ is also a versatile logistic with parameters
\begin{align}
  {\boldsymbol{\beta}}' &= [ \beta_1,\ldots,\beta_K,-y_1\theta_1,\ldots,-y_N\theta_N]^T \in {\mathbb{R}}^{K+N} \\
  {\boldsymbol{\gamma}}' &= [ \gamma_1,\ldots,\gamma_K,\phi_1,\ldots,\phi_N]^T \in {\mathbb{R}}^{K+N} \\
  {\boldsymbol{\mu}}' &= [ \mu_1,\ldots,\mu_K,1,\ldots,1]^T \in {\mathbb{R}}^{K+N} \ .
\end{align}
In the binary classification problem, the b-Log--v-Log conjugacy relationship helps with posterior inference.
By construction, the posterior distribution of the weights is a versatile logistic.

\begin{figure}[t!]
\begin{minipage}[b]{0.48\linewidth}
\centering
                \includegraphics[width=0.95\columnwidth]{graphicalmodel5.pdf}
                \caption{The graphical model for Bayesian boosting.}
                \label{fig:gm}
\end{minipage}
\hspace{0.4cm}
\begin{minipage}[b]{0.48\linewidth}
\centering
                \includegraphics[width=0.95\columnwidth]{graphicalmodel8.pdf}
                \caption{The dynamic graphical model for Bayesian boosting.}
                \label{fig:gmu}
\end{minipage}
\end{figure}

\section{Incorporating Noise} \label{sec:model}

We will now build upon the model presented in \S\ref{sec:core}.
Suppose we fix an instance $\mathbf{x} \in \mathcal{X}$ and repeatedly generate labels from $p(y \mid \mathbf{x} , F)$.
If the labels are \textit{true}, then the empirical ratio of plus-to-minus labels will converge to $\exp[F(\mathbf{x})]$, i.e., the odds ratio.

On the other hand, according to our model, if the labels are \textit{noisy}, the empirical odds ratio converges to some fixed value, which is independent of $\mathbf{x}$.
Let $e^\xi$ be this noise-related odds ratio.
Equivalently, $\xi$ is an instance-independent, static log-odds-ratio, which we refer to as the \textit{noise grade}.
For example, noise grades of $-\infty$, $0$, and $+\infty$, translate to random label assignments of $+1$ with probability $0$, $1/2$, and $1$, respectively.

Let $w$ take on values in $\{0,1\}$ and encode whether a label is true or noisy.
We can merge the two label types into the following conditional label probability:
\begin{align}

  \!\!\!
  p(y \mid w, \mathbf{x}, \xi, F) &= \frac{1}{1+\exp[-y(w F(\mathbf{x}) + (1{-}w) \xi)]}
  = \frac{w}{1+\exp[-y F(\mathbf{x})]} + \frac{1{-}w}{1 + \exp[-y \xi]} .
  \label{eq:merge}
\end{align}
The role of $w$ selects the label type: true or noisy.
Treating $w$ as a latent variable, we embellish the model of \S\ref{sec:core}:
\begin{samepage}
\begin{ENUM}
  \item Draw $c_m \stackrel{\text{iid}}{\sim} \text{v-Log}(\bar{\mathbf{u}},{\mathbf{0}},\mu_0 {\mathbf{1}})$\ ($m=1,\ldots,M$) \ .
  \item Construct $F = \sum_{m=1}^M c_m f_m$ \ .
  \item Draw $\mathbf{x}_{1:N}$ independently according to some distribution over $\mathcal{X}$ \ .
  \item Draw $\xi \sim \text{v-Log}(\bar{\mathbf{u}},{\mathbf{0}},\mu_0' {\mathbf{1}})$ \ .
  \item Draw $\theta \sim \text{Beta}(\zeta_1,\zeta_2)$ \ \ \ (${\boldsymbol{\zeta}} \in {\mathbb{R}}^2_+$) \ .
  \item Draw $w_n | \theta \stackrel{\text{iid}}{\sim} \text{Bernoulli}(\theta)$ \ ($n=1,\ldots,N$) \ .
  \item Draw $y_n \in \{-1,+1\}$ independently according to
	\begin{align}
	  p(y_n \mid w_n, \mathbf{x}_n, \xi, F) &=  \frac{1}{1+\exp[-y_n(w_{n}F(\mathbf{x}_n) + (1{-}w_{n})\xi)]}  \ .

	\end{align}
\end{ENUM}
\end{samepage}
There is now a prior assigned to $c_m$, and the new steps (4-6) model noise.
The graphical model is depicted in Figure~\ref{fig:gm}.
Although not immediately apparent, this model subsumes label inversion as a form of noise (\S\ref{supp:sec:inversion}).

From the classification standpoint, the primary latent variables of the above generative process are still $c_{1:M}$, or the \textit{weights}.
The latent $w_{1:N}$, or \textit{type selectors}, are responsible for the type of label generated.
They are drawn independently from $\theta$, the \textit{type prior}.
We can reason that $\theta/(1-\theta)$ represents a signal-to-noise ratio (SNR).
This stems from the expected value of $N\theta$ true labels and $N(1{-}\theta)$ noisy labels (see \S\ref{supp:sec:inversion} for full details).
Alternatively, we can use the prior of $\theta$ for the SNR estimate, yielding $\mathbb{E}\{\Theta\}/\mathbb{E}\{1{-}\Theta\} = \zeta_1/\zeta_2$.

\section{Variational Inference} \label{sec:vi}

With a Bayesian model in place, our focus turns to the posterior distribution of the latent variables.
This allows us to construct a classifier by estimating the mean or mode of the posterior weights ($c_m$).
We accomplish this with stagewise variational inference.\footnote{The common approach in this situation is to use a Gibbs sampler. The Gibbs sampler for the graphical model of Figure~\ref{fig:gm} is given in \S\ref{supp:sec:gibbs}. We used Adaptive Rejection Sampling \cite{gilks1992} to sample a v-Log. In practice, this approach was too time consuming, which is why we turned to inference.}

\paragraph{Full Variational Inference.}

Before motivating our stagewise approach, we first review variational inference.
In our model we have the observed variables $\langle \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} \rangle$ and the latent variables $\langle c_{1:M} , \xi , w_{1:N} , \theta \rangle$.
We are interested in the posterior $p(c_{1:M} , \xi , w_{1:N} , \theta \,|\, \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} )$, which is proportional to the joint $p(c_{1:M} , \xi , w_{1:N} , \theta , \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} )$.
In variational inference, we introduce a distribution $q(c_{1:M} , \xi , w_{1:N} , \theta)$ to bound the log of the marginal probability of the observations \cite{wainwright2008},
\begin{align}
  \!\!
   \log  p(  \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} \! )

	  &\!\geq\! \textstyle{\int} q(c_{1:M} , \xi , w_{1:N} , \theta) \log p(c_{1:M} , \xi , w_{1:N} , \theta , \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} \! ) \text{d} c_{1:M} \text{d} \xi \text{d} w_{1:N} \text{d} \theta \notag \\
  & \qquad \quad - \textstyle{\int} q(c_{1:M} , \xi , w_{1:N} , \theta) \log q(c_{1:M} , \xi , w_{1:N} , \theta) \text{d} c_{1:M} \text{d} \xi \text{d} w_{1:N} \text{d} \theta \label{eq:elbo}
				  \ .
\end{align}
The right-hand side of \eqref{eq:elbo} is referred to as the \underline{e}vidence \underline{l}ower \underline{bo}und (ELBO).
Using the KL-divergence, we can also write
\begin{align}

\log  p(  \mathcal{F} , \mathbf{x}_{1:N} , y_{1:N} \! ) = \textsc{kl}(\ q(c_{1:M} , \xi , w_{1:N} , \theta) \parallel \text{posterior}\ ) + \textsc{elbo} \ .
\end{align}
The KL divergence provides a measure of closeness between the auxiliary distribution and the posterior.
We maximize the ELBO with respect to the parameters of $q$, thereby minimizing the KL divergence to the posterior.
We use \textit{mean-field variational inference}, i.e., we assume a factorized $q$:
\begin{align}
  q(c_{1:M} , \xi , w_{1:N} , \theta) = \textstyle{\prod_{m=1}^M} q(c_m) \cdot q(\xi) \cdot \textstyle{\prod_{n=1}^N} q(w_n) \cdot q(\theta) \ .
  \label{eq:factor}
\end{align}
Each component of the factorized variational distribution has a form and variational parameters.
For example, a reasonable form of $q(c_m)$ is a versatile logistic with variational parameters given by some slope, knot and weight vectors.
Typically, we optimize the parameters with coordinate ascent, updating each in turn, holding the others fixed.
In our model, this yields the following updates \cite{bishop2006}:
\begin{align}
  \log q^*(c_m) &\leftarrow \mathbb{E}_q [ \log p( C_{1:m-1}, c_m, C_{m+1:M}, \Xi, W_{1:N}, \Theta, \mathcal{F}, \mathbf{x}_{1:N}, y_{1:N} ) ] +\text{const} \label{eq:qc} \\
  \log q^*(\xi) &\leftarrow \mathbb{E}_q [ \log p( C_{1:M}, \xi, W_{1:N}, \Theta, \mathcal{F}, \mathbf{x}_{1:N}, y_{1:N} ) ] +\text{const} \\
  \log q^*(w_n) &\leftarrow \mathbb{E}_q [ \log p( C_{1:M}, \Xi, W_{1:n-1}, w_n, W_{n+1:N}, \Theta, \mathcal{F}, \mathbf{x}_{1:N}, y_{1:N} ) ] +\text{const} \\
  \log q^*(\theta) &\leftarrow \mathbb{E}_q [ \log p( C_{1:M}, \Xi, W_{1:N}, \theta, \mathcal{F}, \mathbf{x}_{1:N}, y_{1:N} ) ] +\text{const}  \ .
\end{align}
Each term on the right is a leave-one-out expectation over the latent variables, resulting in a function of the corresponding left-out latent variable.
Running the variational inference algorithm repeatedly cycles through these updates.

This algorithm is not convenient.
The chosen form of the approximate posterior weight distributions is a versatile logistic.
From conjugacy, the number of parameters required to specify each distribution is linear in the number of examples ($N$).
Additionally, we hope to use a large number of base classifiers, even for small datasets.
Thus, for our classification problem, cycling through all auxiliary weight distributions is impractical because integrating over the weights is too much of a computational burden.

\begin{algorithm}[tb]
   \caption{VIBoost}
   \label{alg:viboost}
\begin{algorithmic}
   \STATE {\bfseries Input:} $\{(\mathbf{x}_n,y_n)\}_{n=1}^N$, $\mathcal{F}$, $\mu_0\in{\mathbb{R}}_+$, $\mu_0'\in{\mathbb{R}}_+$, ${\boldsymbol{\zeta}}\in{\mathbb{R}}_+^2$

   \STATE Initialize $H : \mathcal{X} \to \{-1,+1\}$ to the zero function
   \STATE Initialize ${\boldsymbol{\eta}} \in {\mathbb{R}}_+^2$, ${\boldsymbol{\omega}} \in {\mathbb{R}}_+^2$, and ${\boldsymbol{\phi}} \in [0,1]^N$
   \STATE Define ${\boldsymbol{\beta}}(h) \triangleq [+1,-1,-y_1 h(\mathbf{x}_1),\ldots,-y_N h(\mathbf{x}_N)]^T$
   \STATE Define ${\boldsymbol{\gamma}}(H,h) \triangleq [0,0,-H(\mathbf{x}_1) h(\mathbf{x}_1),\ldots,-H(\mathbf{x}_N) h(\mathbf{x}_N)]^T$
   \STATE Define ${\boldsymbol{\mu}}({\boldsymbol{\phi}}) \triangleq [\mu_0,\mu_0,\phi_{1},\ldots,\phi_{N}]^T$
   \FOR{$t=1$ {\bfseries to} $T$}
   \vskip 0.4ex

   \STATE $h_t \leftarrow \argmax_{h \in \mathcal{F}} \ \vert \text{Mode}[\text{v-Log}({\boldsymbol{\beta}}(h),{\boldsymbol{\gamma}}(H,h),{\boldsymbol{\mu}}({\boldsymbol{\phi}}))] \vert$
   \vskip 0.4ex
   \WHILE{\textsc{ELBO} increases significantly}
   \vskip 0.4ex
   \STATE $\alpha_t \leftarrow \text{Mode}[\text{v-Log}({\boldsymbol{\beta}}(h_t),{\boldsymbol{\gamma}}(H,h_t),{\boldsymbol{\mu}}({\boldsymbol{\phi}}))]$
   \vskip 0.4ex
   \STATE $\omega_1 \leftarrow \mu_0' + \textstyle{\sum_{n=1}^N} (1-\phi_{n}) \mathds{1}} %\text{I\!\!\!\!1}\{y_n = -1\}$
   \vskip 0.4ex
   \STATE $\omega_2 \leftarrow \mu_0' + \textstyle{\sum_{n=1}^N} (1-\phi_{n}) \mathds{1}} %\text{I\!\!\!\!1}\{y_n = +1\}$
   \vskip 0.4ex
   \STATE $\kappa_n \stackrel{1:N}{\leftarrow} \displaystyle{\frac{ \exp\left[ \uppsi(\eta_1) - \uppsi(\eta_2) + \uppsi(\omega_0) - \uppsi(\omega_2)\mathds{1}} %\text{I\!\!\!\!1}{\{y_n=+1\}} - \uppsi(\omega_1)\mathds{1}} %\text{I\!\!\!\!1}{\{y_n=-1\}} \right] }{ 1 + \exp[-y_n(H(\mathbf{x}_n) + \alpha_t h_t(\mathbf{x}_n))]}}$
   \vskip 0.4ex
   \STATE $\phi_{n} \stackrel{1:N}{\leftarrow} \kappa_n / (1 + \kappa_n)$
   \vskip 0.4ex
   \STATE $\eta_1 \leftarrow \zeta_1 + \textstyle{\sum_{n=1}^N} \phi_{n}$
   \vskip 0.4ex
   \STATE $\eta_2 \leftarrow \zeta_2 + \textstyle{\sum_{n=1}^N} (1-\phi_{n})$
   \vskip 0.4ex
   \ENDWHILE
   \STATE $H \leftarrow H + \alpha_t h_t$
   \ENDFOR
   \STATE {\bfseries Output:} classifier\ \ $\sign\{H(\cdot)\}$
\end{algorithmic}
\end{algorithm}

\paragraph{Stagewise Variational Inference.}

To address these issues, we propose a dynamic model over the current static one: with a current estimate of $F = \sum_m c_m f_m$, we introduce a single base classifier and then run variational inference on the latent $\langle c,\xi, w_{1:N}, \theta \rangle$.
The regression counterpart would be Forward Stagewise Regression, a greedy algorithm which finds a sparse subset of covariates and is structurally similar to AdaBoost \cite{friedman2009}.

In each main loop iteration, let $H(\cdot)$ be the current estimate of the true log-odds-ratio $F(\cdot)$ and suppose we have a ``promising'' candidate $h \in \mathcal{F}$ that we wish to merge with $H$.
This promising classifier is found greedily (details below) and, once found, becomes a fixed variable in the model.
Now, rather than $M$ latent weights, we have a single latent weight $c$ corresponding to $h$.
Every update of $H$ induces a new model to which we apply variational inference.
This new, time varying graphical model is featured in Figure~\ref{fig:gmu}.

Let $\mathcal{D}$ denote the evidence, i.e., the observed variables $\mathbf{x}_{1:N}$, $y_{1:N}$, $H$, and $h$.
At each stage we assume the following distributions:
\begin{align}
  p(c \mid \mathcal{D}) &\approx q( c \mid {\boldsymbol{\beta}},{\boldsymbol{\gamma}},{\boldsymbol{\mu}} )
                        \sim \text{v-Log}({\boldsymbol{\beta}},{\boldsymbol{\gamma}},{\boldsymbol{\mu}})
						&
  p(\xi \mid \mathcal{D}) &\approx q( \xi \mid {\boldsymbol{\omega}} )
                        \sim \text{v-Log}(\bar{\mathbf{u}},{\mathbf{0}},{\boldsymbol{\omega}}) \\
  p(w_n \mid \mathcal{D}) &\approx q( w_n \mid \phi_n )
                        \sim \text{Bernoulli}(\phi_n)
						&
  p(\theta \mid \mathcal{D}) &\approx q( \theta \mid {\boldsymbol{\eta}} )
                        \sim \text{Beta}({\boldsymbol{\eta}}) \ .
\end{align}

The variational updates and the ELBO are derived in \S\ref{supp:sec:vi} and \S\ref{supp:sec:elbo}, respectively.
The general approach is to isolate the terms of the log-likelihood that feature the variable of interest---all other terms will cancel after normalization and are extraneous.
We then take expectations and attempt to identify a distribution.

The resulting algorithm, \textit{VIBoost}, is presented in Algorithm~\ref{alg:viboost} (\ $\uppsi(\cdot)$ is the digamma function ).
Possible modifications include (a) fixing the number of variational inference iterations so that ELBO calculations are avoided, and (b) setting the $\alpha_t$ once and skipping its update in the variational inference loop.

We greedily select the next base classifier by finding the v-Log posterior with maximal mode.
A large mode suggests that the corresponding weight possesses discriminative classification strength.
We opted for the mode rather than the mean; we now justify this choice.

A v-Log distribution with more than two slope/knot/multiplicity terms has the advantage of being a one-dimensional density, but is cumbersome when evaluating statistics of interest.
Computing the normalization constant, mode and mean require iterative techniques, which can bog down any algorithm.
However, if we replace $\mu_k \log(1 + e^{\beta_k(z-\gamma_k)})$, a summand of the log-density, with the single-tail approximation $\mu_k e^{\beta_k \tau (z - \gamma_k)}$  ($\tau > 0$) we arrive at the modal estimate of
\begin{align}
  \textstyle{
  \alpha =
  \frac{1}{2\tau\beta} \log \left( \frac{
           \sum_{k:\beta_k<0} \mu_k e^{\tau\beta\gamma_k}
     }{
           \sum_{k:\beta_k>0} \mu_k e^{-\tau\beta\gamma_k}
		   } \right) } \ ,
		\label{eq:appxmode}
\end{align}
where $\beta = \vert \beta_k \vert$ is constant (\S\ref{supp:sec:cupdate} and \cite{friedman2000}).
For a unimodal distribution, this closed-form expression provides an efficient way of estimating expectations.
Thus, in avoiding numerical integration, Algorithm~\ref{alg:viboost} is performing approximate variational inference.

\section{Relation to AdaBoost} \label{sec:adaboost}

We now compare our algorithm to AdaBoost.
Consider the simpler model of \S\ref{sec:core}, a true-label dataset with prior assignments (Figure~\ref{fig:core}).
This leaves the greedy step of finding the maximal mode in Algorithm~\ref{alg:viboost} and updating $H$ without the variational inference.
We now investigate the approximation supplied by \eqref{eq:appxmode} with $\mu_0 = \tau=1$.
Let $Z=\sum_{n=1}^N e^{-y_n H(\mathbf{x}_n)}$ and $d_n=e^{-y_n H(\mathbf{x}_n)} /Z$ so that $\sum_{n=1}^N d_n=1$.
The approximate mode $\alpha$ is
\begin{align}
  \textstyle{
  \frac{1}{2}
  \log \left( \frac{ 1/Z + (1 - \varepsilon) }{ 1/Z + \varepsilon } \right) }
  \ ,  \label{eq:1oz}
\end{align}
where $\varepsilon =  \sum_{n} d_n \mathds{1}} %\text{I\!\!\!\!1}\{h(\mathbf{x}_n) \neq y_n\} $ is a weighted error ascribed to the new classifier (\S\ref{supp:sec:modeappx}).
When compared to AdaBoost, the update is identical when the $1/Z$ term is not present.
Effectively, the $1/Z$ term results in a shrinkage of the assigned weights (see Figure~\ref{fig:1Z}).

The variable $Z$ is equal to the current exponential loss.
If $Z$ is small then $1/Z$ is large which leads to a dampened weight assignment (and vice versa).
Assuming the exponential loss decreases with more iterations, the algorithm acts like AdaBoost early on and then becomes more conservative with each iteration.
Quinlan \cite{quinlan1996} incorporated similar smoothing for AdaBoost and described it as ``necessarily ad-hoc''.
In the proposed model, this smoothing results from the prior assignment.
We also note that AdaBoost selects the base classifier that minimizes $\varepsilon$.
From \eqref{eq:1oz}, this coincides with the largest approximation-based mode

The slopes ${\boldsymbol{\beta}}(h)$ as defined in Algorithm~\ref{alg:viboost} contain individual [mis]matches of the base classifier with the labels, whereas the knots contain individual, weighted [mis]matches of the base classifier with the current log-odds-ratio estimate.
The prior effectively augments the data by inserting two phantom examples.
Each example lies in the zero level set of $H$ as indicated by a knot of $0$ ($H(\mathbf{x})h(\mathbf{x}) = 0 \Rightarrow H(\mathbf{x}) = 0$).
The slopes of $\pm 1$ presume that the base classifier succeeds in correctly labeling one of the pseudo-examples, while failing with the other.

Finally, leveraging $d_{1:N}$ we can rewrite Algorithm~\ref{alg:viboost} to use the weighted error $\varepsilon$ rather than a mode search.
Using these errors for ranking the base classifiers---as done in AdaBoost---decreases computation time significantly when searching for a new candidate base classifier.
The greedy search in VIBoost would then closely match AdaBoost's in computation time, thereby leading to an efficient algorithm with a similar runtime to AdaBoost.

\begin{figure}[t]
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{oneoverz.pdf}
                \caption{\small The effects of the versatile logistic prior on classifier weights ($\mu_0\!=\!w_{n}\!=\!\tau\!=\!1$)}
                \label{fig:1Z}
\end{minipage}
\hspace{0.03\linewidth}
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{spamerr.pdf}
                \caption{\small Spam dataset, classification error \\ \ }
                \label{fig:spamerr}
\end{minipage}
\hspace{0.03\linewidth}
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{techerr.pdf}
                \caption{\small State dataset, classification error \\ \ }
                \label{fig:texterr}
\end{minipage}
\end{figure}

\section{Experiments}  \label{sec:exp}

We studied VIBoost on real and synthetic data.
We found that VIBoost works as well as AdaBoost for binary classification.
More importantly, we show that the variables accounting for label noise are a meaningful diagnostic of misfit.
For all experiments, our VIBoost initialization was $\mu_0\!=\!\mu'_0\!=\!\phi_{n}\!=\!\zeta_j\!=\!\eta_j\!=\!\tau\!=\!1$.
Setting $\tau\!=\!1$ provides the closest means of comparison with AdaBoost.
For all experiments we used decision stumps as our base classifiers.
Using the variational parameters of the type prior, we use $\eta_1/\eta_2$ for the SNR.
For the noise grade we use $\log(\omega_2/\omega_1)$, the mode associated with the approximate posterior.
All results presented are average values calculated over $40$ runs.

As VIBoost outputs a classifier, we investigate classifier quality on two real-world, text datasets.
The first dataset is the $57$-feature \textit{spam} dataset \cite{ucimlr}.
With 6,401 examples, each run trained on a random $10\%$ and tested on the remaining $90\%$.
The second dataset is a \textit{state} dataset \cite{gabrilovich2004} comprising 145 documents with 22,648 features (bag of words).
Instead of the word count, however, we used a present/absent binary value.
Each document relates to Illinois or Michigan.
We trained on $30\%$ and tested on the remaining $70\%$ (random splits).
Error results are featured in Figures \ref{fig:spamerr} and \ref{fig:texterr} and reveal that VIBoost and AdaBoost performed similarly.

In addition to a classifier, VIBoost also provides noise statistics.
Using a synthetic dataset, we now look at the algorithm's estimate of the posterior SNR ($\eta_1/\eta_2$) and the posterior noise grade (\,$\log(\omega_2/\omega_1)$\,) after $50$ iterations.
We simulated 100 examples on the real line with $\mathcal{X} = \{-99,-97,-95,\ldots,+99\}$.
Following the generative process of \S\ref{sec:model}, we constructed the \textit{step} dataset with $F(x)=+\infty$ for $x$ positive and $-\infty$ for $x$ negative ($+1$ label for $x$ positive and $-1$ label for $x$ negative).
With a noise grade of $\log 3 \approx 1.1$, we varied the type prior, $\theta$, of the generative process.
Figures \ref{fig:stepsnr} and \ref{fig:stepxi} respectively show the SNR and noise grade with varying $\theta$.
In a pure-noise situation ($\theta = 0$) the SNR is at its lowest and the noise grade is best estimated.
Conversely, in the absence of noise ($\theta = 1$) the SNR is at its greatest, rendering the noise grade estimate irrelevant.

The last dataset we considered in this paper was also simulated.
Inspired by \cite{long2010}, we constructed a 1,200-example, 31-feature \textit{Long-Servedio} dataset---a dataset that provably ``breaks'' AdaBoost and many other algorithms with a convex loss minimization.
The details can be found in \cite[\S12.3]{schapire2012} and Matlab code is included in \S\ref{supp:sec:longserv}.
Each run comprised $200$ training examples and $1000$ testing examples (random splits).
Following \cite{schapire2012}, the noise level was set to $0.20$.
In our context, this translates to a type prior of $0$ (always reassigning a random label) and a noise grade of $\approx -1.4$.
Not surprisingly, AdaBoost and VIBoost do not succeed in finding a decent classifier for this set.
However, the SNR values produced by VIBoost indicate that poor classification should be expected (Figure~\ref{fig:lssnr}).
As a result, the algorithm shifts its focus from classification to noise quantification.

As variational inference navigates through a vast set of auxiliary distributions, our only verifiable means of efficacy is provided by the ELBO (\S\ref{supp:sec:elbo}).
Empirically, we have noticed that ELBO increases are larger in the beginning main-loop iterations.
As the algorithm progresses, changes in the ELBO are quite small and sometimes negative (and small).
The small changes are expected because the composite classifier's accuracy is improving.
We hypothesize that the negative changes stem from the modal approximation used in a $w_n$-update expectation.
Alternatively, the ELBO requires a v-Log normalization constant, which we compute numerically and may be inexact.

\begin{figure}[t]
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{stepsnr1.pdf}
                \caption{\small Step dataset, Signal-to-noise ratio}
                \label{fig:stepsnr}
\end{minipage}
\hspace{0.03\linewidth}
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{stepxi.pdf}
                \caption{\small Step dataset, Noise Grade}
                \label{fig:stepxi}
\end{minipage}
\hspace{0.03\linewidth}
\begin{minipage}[b]{0.30\linewidth}
\centering
                \includegraphics[width=1.0\columnwidth]{lslogsnr.pdf}
                \caption{\small Long-Servedio dataset, Signal-to-noise ratio}
                \label{fig:lssnr}
\end{minipage}
\end{figure}

\section{Conclusion}  \label{sec:con}

We have developed a new boosting-like algorithm.
VIBoost attempts to fit a posterior distribution by applying variational inference to a dynamic model.
We began with a model centered around the binary classification problem and augmented it hierarchically to account for noise.

We did not set out to improve AdaBoost.
In addition to forming a binary classifier, the Bayesian model facilitated a label noise extension and we were able to extract information beyond classification.
We have observed experimentally that a good classifier is accompanied by a large SNR.
The SNR may explain why a poor classifier is returned by the learning algorithm.
We demonstrated this by analyzing the Long-Servedio dataset.

This model and accompanying algorithm are fertile ground for future work.
This paper did not address multi-class problems or regression.
We can also extend our model by forming connections between instances, base classifiers, and classifier weights (currently, these three features are conditionally independent given the labels).
We can also form dependencies between instances and label types, modeling varying levels of noise throughout the instance space.

\nocite{jordan1995}
\nocite{chen1995}
{ \small