\section{Introduction}

The advanced network of ground-based gravitational wave (GW)
interferometers is designed to detect and study GW emission from
events such as the mergers of binary systems composed of neutron stars
and/or black holes to distances of hundreds of Mpc (see
\citet{ligo2013} and references therein).  In mergers containing at
least one neutron star, counterpart electromagnetic radiation is expected,
potentially ranging from a short-duration gamma-ray burst through
optical/near-IR emission from the radioactive decay of r-process
nuclei to radio emission from ejecta interacting with the circumbinary
medium (e.g.,
\citealt{li98,nakar2011,metzger2012,barnes2013,2013ApJ...775..113T,2014ApJ...780...31T,2014ApJS..211....7A,berger2014,cowperthwaite2015}).
The detection of an electromagnetic counterpart will provide critical
insight into the physics of the event,
helping to determine the distance scale, energy scale, and the
progenitor environment, as well as insight into the behavior of matter
post-merger (e.g., the production of jets and outflows).

With this motivation, we recently began an observational program
using the wide-field Dark Energy Camera (DECam; \citealt{flaugher2015}) on the
Blanco 4-m telescope at Cerro Tololo Inter-American Observatory to
search for optical counterparts to GW triggers from the new advanced
GW detectors (LIGO, \citet{ligo2009}; Virgo, \cite{virgo2009}).
This program was awarded three target of opportunity nights to observe LIGO-triggered events during the 2015B semester; observations were coordinated with and managed by the Dark Energy Survey (DES).
Our program is optimized for detection of kilonovae, the hypothesized
optical counterparts
of mergers involving neutron stars, which would appear as
red transients with expected
decay timescale of about a week  (for an overview of our program see \citet{des2016}).

On 2015 September 14 at 09:50:45 UT the Advanced LIGO interferometer network detected
a high significance candidate GW event designated GW150914
\citep{ligo2016} and two
days later provided spatial location information in the form of
probability sky maps via a private GCN circular (\#18330; \cite{GCN18330}).
We initiated observations with DECam, a 3 deg$^2$ field-of-view
instrument, on 2015
September 18 in an effort to identify an optical counterpart.  Here
we describe the observations and provide the results of the
three-epoch search.  These DECam observations  are
the deepest search for an optical counterpart to GW event GW150914
\citep{gw-em-overivew}.

\section{DECam Observations of GW150914}
\label{sec:data}

The detection of GW150914 was triggered by the cWB (coherent
WaveBurst; \citet{2008CQGra..25k4029K})
unmodeled burst analysis during real-time data processing.
On 2015 September 16, the LIGO Virgo Collaboration (LVC) provided two
all-sky localization probability maps for the event, generated from
the cWB and LALInferenceBurst (LIB; \citet{veitch2014}) analyses.
The cWB online trigger
analysis makes minimal assumptions about signal morphology by
searching for coherent power across the LIGO network.  The LIB
analysis is a version of the LALInference analysis Bayesian
forward-modeling-based follow up tool that uses a Sine-Gaussian signal
morphology instead of models of compact binary mergers
\citep{veitch2014}; for information on both algorithms see
\citet{essick2015}.  The maps provided initial spatial localization
of 50\% and 90\% confidence regions encompassing about 200 and 750
deg$^2$, respectively.

Our first observations with DECam took place on 2015
September 18 UT.  Overall, we imaged 102 deg$^2$ covering
38\% of the total probability in the initial cWB map;
see Table~\ref{datatable} for a summary of our DECam observations.
As shown in Figure~\ref{datafig}, 18 deg$^2$, were centered on the LMC.
For the remaining 84 deg$^2$
we obtained 3 separate epochs of imaging.
At each epoch we acquired one
   90-sec exposure in $i$ band  and two 90-sec exposures
   in $z$ band.
  The first
epoch spanned 4--5 days post-GW trigger (2015 September 18--19 UT),
the second epoch 7 days post-GW trigger (2015 September 21 UT), and the
third was obtained 24 days post-GW trigger (2015 October 08 UT).

Subsequently, in January 2016, the LVC released a revised sky map of
localization probabilities from a LALInference analysis that used the
assumption that the signal arises from a compact binary coalescence
(CBC). That analysis also showed that the
data are most consistent with models of a binary black hole merger
(BBH). The LALInference-based map is considered the most accurate and
authoritative localization for this event.
Our 102 deg$^2$  cover a total of $11\%$ probability in this new map,
as the localization region has
shifted significantly southward (see Figure~\ref{datafig}) relative
to the initial cWB map.

Our single-epoch exposures achieve median $5\sigma$ point-source
limiting magnitudes of $i = 22.5$ and $z = 21.8$ with an rms variation
among the images of $\pm 0.5$ mag.  This value is a consequence of
night-to-night variations in the observing conditions (see
Table~\ref{datatable}) and of a strong gradient in stellar density and
extinction along the major axis of the region imaged (see
Figure~\ref{datafig}).

\begin{deluxetable*}{lrrrccccc}
\tabletypesize{\footnotesize}
\tablecaption{Summary of Observations
\label{datatable}}
\tablehead{
    \colhead{Program} &
    \colhead{Night}   &
    \colhead{MJD}     &
    \colhead{$\Delta t$\tablenotemark{a}} &
    \colhead{$\langle$PSF(FWHM)$_i$$\rangle$}    &
    \colhead{$\langle$airmass$\rangle$}  &
    \colhead{$\langle$depth$_i\rangle$}   &
    \colhead{$\langle$depth$_z\rangle$}  &
    \colhead{$A_{\mathrm{eff}}$\tablenotemark{b}}
\\
    \colhead{ } &
    \colhead{(UT)} &
    \colhead{ } &
    \colhead{(days)} &
    \colhead{(arcsec)} &
    \colhead{ } &
    \colhead{(mag)} &
    \colhead{(mag)} &
    \colhead{(deg$^2$)}
}
\startdata
Main, 1$^{\mathrm{st}}$ epoch & 2015-09-18 & 57383 &  3.88 & 1.38 & 1.50 & 22.71 & 22.00 & 52.8 \\
                     & 2015-09-19 & 57384 &  4.97 & 1.35 & 1.46 & 22.82 & 22.12 & 14.4 \\
Main, 2$^{\mathrm{nd}}$ epoch & 2015-09-21 & 57286 &  6.86 & 2.17 & 1.51 & 22.18 & 21.48 & 67.2 \\
Main, 3$^{\mathrm{rd}}$ epoch & 2015-10-08 & 57303 & 23.84 & 1.46 & 1.40 & 22.33 & 21.63 & 67.2 \\[2pt]
\tableline \\[-5pt]
LMC, initial           & 2015-09-18 & 57383 &  3.98 & 1.14 & 1.30 & 21.32 & 20.62 & 14.4 \\
LMC, extension    & 2015-09-27 & 57292 & 12.96 & 1.21 & 1.28 & 20.91 & 20.21 & 33.6 \\[-5pt]
\enddata
\tablecomments{Summary of the observations performed in the ``main''
  search program, described in this paper, and the ``LMC'' program,
  described in the companion paper \citet{Annis:2016}.  We observed at
  high airmass because the region of interest was rising at the
  end of the night. The PSF FWHM, and therefore the actual depth
  achieved, are partly affected by these high airmass conditions.
  The reported depth corresponding to 5-$\sigma$ point source detection in the search images.
  Variations in cloud conditions are also responsible for the
  variation in depth.  The effective area
  imaged in the main program corresponds to 28 camera fields. The area
  covered in the LMC program totaled 20 fields.}
\tablenotetext{a}{Time elapsed between the trigger time and the time
  stamp of the first image of the night.}  \tablenotetext{b}{Effective
  area imaged, considering that approximately 20\% of the 3 deg$^2$
  field of view of DECam is lost due to chip gaps (10\%), 3 dead CCDs
  (5\%) and masked edge pixels (5\%).
}
\end{deluxetable*}

\subsection{Observing Strategy}

We chose the location and sequence of DECam observations using an
automated observing strategy algorithm.
The algorithm utilizes the GW localization map, an
estimate of the event distance, and a model of the expected optical
emission (e.g., \citealt{barnes2013}).  This information is folded in
with observational information, including a map of sky brightness
(using the DES sky brightness model;
\citealt{2012ASPC..461..201N}), the atmospheric transmission (using
information on airmass and the interstellar dust extinction from
Planck; \citealt{planck2014}), the expected seeing (from scaling laws
with airmass and wavelength), and the confusion-limit probability
(based on stellar density maps) to produce a full source-detection
probability as a function of sky location.
We used this map to observe the highest probability region that
    included area both inside and outside the DES footprint.

In the case of GW150914 the localization region intersected the
Large Magellanic Cloud (LMC), so we designed a separate set of short
observations to observe the brightest LMC stars.  We
obtained 5-sec $i$ and $z$ band exposures covering 18 deg$^2$ centered
on the LMC on 2015 September 18 and 27.  This shallower data set was used to
search for a potental failed supernova in the LMC; the
results are reported in a separate paper \citep{Annis:2016}.
Figure~\ref{datafig} shows a sky map computed for the end of the first
night of observations, zoomed in to the region of interest and
detailing the fields observed in each of the three epochs in red.

\begin{figure}[h!]
  \includegraphics[width=\linewidth,,trim=10 0 10 0, clip]{fig1.pdf}
  \caption{The color image shows the estimated limiting point-source
    magnitude for a 90-sec $i$ band exposure as a function of sky
    position for our first night of DECam observations just before
    sunrise.  In this area and for this time of night, the variations
    are mostly due to interstellar dust extinction.  The dotted
    contours show the initial (September 2015) {\tt
      skyprobcc\_cWB\_complete} map, while the solid contours are for
    the final (January 2016) {\tt LALInference\_skymap}.  There is an
    island of significant probability in the Northern hemisphere in
    the {\tt skyprobcc\_cWB\_complete} map, not present in the {\tt
      LALInference\_skymap}, so the dotted contours do not show the
    complete 50\% or 90\% areas.  The hexagonal DECam fields observed
    are shown, with red for the main search and orange for the short
    exposure LMC data.  Fields located on the west (left) side of the
    region of interest overlap with the DES area (footprint boundary
    shown in light-gold).  The excluded region (dark grey) is beyond
    the horizon limit that could be observed with DECam at that time.
    The total area inside the camera pointings is about 102
    deg$^2$. We covered about 11\% of the total localization
    probability in the final map, and 38\% of the initial map.
    The projection shown is an equal-area McBryde-Thomas flat-polar quartic projection.}
 \label{datafig}
\end{figure}

\subsection{Image Processing}

Our data analysis relies on subtracting earlier template images from
the science images taken for this program.  In the area that overlaps
the DES footprint (25\% of the total), we used DES images from the
first two seasons of the survey as templates.
In the 75\% of the area outside of the DES footprint, we used publicly
available DECam data from the NOAO Science Archive
({\tt\url{portal-nvo.noao.edu}}), requiring exposures of at least 30
sec in $i$ and $z$ bands.

We processed the DECam search and template images using the DES Data
Management single-epoch image processing software
\citep{2012ApJ...757...83D,DESDM2012,sevilla2011,Gruendl:2016}. Its output images were
 used as input to the difference imaging pipeline, which we
developed from the DES Supernova pipeline \citep{kessler2015}.  The
main adaptation of the pipeline for our purposes was to generalize to
the case of search and template images with arbitrary relative
alignment.
A candidate requires two SExtractor \citep{1996A&AS..117..393B}
detections in
the first epoch in both $i$ and $z$ bands.
To reduce the large number of detected artifacts,
     each detection must satisfy quality requirements
     (Table 3 of \citealt{kessler2015}) and be selected by our automated scanning
     program \citep{2015AJ....150...82G}.
     For each of the 2349 candidate locations, ``forced'' PSF-fitted
     fluxes and uncertainties are obtained at every epoch regardless
     of whether or not there was a detection.

\section{Analysis}\label{sec:analysis}

While a BBH merger is not expected to result in an optical signature,
it is nevertheless of interest to search for a possible optical counterpart.
As our first epoch of observations occurred 4
days after the trigger, our prior on the search is that any candidate
shall be fading slowly enough to be detectable 7 days after the event,
but not 24 days after the event.

Of the 84 deg$^2$ area outside of the LMC, about 20\% is lost due to
camera fill-factor (see Table \ref{datatable} for details) resulting in
an effective area of
67.2 deg$^2$. In addition,
30\% of the area is lost due to
sparse availability of templates outside of the DES footprint.
Another 10\% loss arises from  processing issues.
This results in 40 deg$^2$ which were used in this analysis.

Based on an analysis of a sample of fake point sources injected into
the images in this area, we find that the typical $80\%$ source detection
completeness in the subtracted images is at $i\approx 22.1$ and
$z\approx 21.2$ mag. In the first epoch, where the observing conditions
were better, we achieve that level of completeness at $i\approx 22.7$
and $z\approx 21.8$, comparable to the 5$\sigma$ point source depth for
those images.
The fakes were in all the images we processed,
thus the completeness depth reflects the variation in conditions as well.

\subsection{Sample Selection}

For the selection criteria described below, multiple observations
   per night (primarily in $z$ band) are combined into a single
   weighted-average flux:

\begin{enumerate}

\item Second-epoch signal-to-noise ratio (S/N) above 2 in both $i$ and $z$
(to enable flux change determination with respect to the first epoch);

\item $\ge 3\sigma$ decline in both $i$ and $z$ fluxes from the first epoch to the second (to
  isolate fading sources; $\sigma$ is defined by the quadrature sum of the
  flux errors in the first two epochs);

\item  S/N $\le 3\sigma$  in both $i$ and $z$ third epoch
(at 24 days post-trigger, to reject long-timescale transients such as supernovae).

\end{enumerate}

\subsection{Results}

In Table~\ref{cuts} we show the impact of our selection criteria on
the sample of candidates as a function of the first epoch $i$-band magnitude.
The decaying light curve requirement has the most impact in reducing the sample
size. None of the candidates pass all the selection criteria.
The area analyzed, 40 deg$^2$, covers
3\% of the localization probability in the final LALInference map
(though it covered 12\% in the initial cWB map).

To interpret these results some caveats are required.
Because our selection criteria impose demands on significance
in the second epoch, the actual first epoch search depth  depends
on the decline rate and $i-z$ color of the source model.
In addition, we have not yet accounted for the degraded sensitivity
to candidates located in bright galaxies.

For a particular source model, we can estimate the search
depth. We applied our selection criteria to a sample of
fake sources randomly placed in our search images before processing
with our difference imaging pipeline.
The fakes have a constant decay rate of 0.3 mags/day
and are red, with $(i-z) \approx 1$, as expected from kilonova models.
The magnitude at which we recover
50\% of the fakes, $m_{50\%}$, is about 1 magnitude brighter
than the 5$\sigma$ point source limiting
magnitude reported in Table~\ref{datatable},
i.e., $m_{50\%} - m_{5\sigma} \approx -1$.
Simulations with bluer models show that
for sources  with $(i-z) = 0$ the search depth
is $m_{50\%} - m_{5\sigma}  \approx -1.4$; for
 $(i-z) = -1$, the search depth is
$m_{50\%} - m_{5\sigma}  \approx -2.4$.
We therefore achieve magnitude limit of $i= 21.5,21.1,20.1$ for object
colors $(i-z)=1,0,-1$, respectively.

\begin{deluxetable}{lrrrr}
\tablecolumns{5}
\tablecaption{Number of selected events \label{cuts}}
\tablehead{
    \colhead{mag($i$)} &
    \colhead{raw} & \colhead{cut 1} &
    \colhead{cut 2} & \colhead{cut 3} }
\startdata
18.0--18.5 &   84  &   1  & 0 & 0 \\
18.5--19.0 &   177  &   1  & 0 & 0 \\
19.0--19.5 &   291 &   2  & 0 & 0 \\
19.5--20.0 &   227  &   2  & 1 & 0 \\
20.0--20.5 &   156 &   17  & 2 & 0 \\
20.5--21.0 &   225 &  42  & 3 & 0 \\
21.0--21.5 &   334  & 84  & 2 & 0 \\
21.5--22.0 &   756 & 159  & 1 & 0 \\
22.0--22.5 &  1099 & 183  & 0 & 0 \\
\tableline\\[-5pt]
total      &  2349 & 491 & 9 & 0 \\[-5pt]
\enddata
\end{deluxetable}

\section{Conclusions}
\label{sec:conclusion}

We presented our search for an optical counterpart to the first
gravitational wave event, GW150914, using the wide-field
DECam instrument.  Our observations cover $102$ deg$^2$ corresponding to 11\%
of the total probability map.
The search images used in this analysis reach median 5$\sigma$
point source depth of $i=22.5$ and $z=21.8$
mag.
Our DECam/Blanco observations are the deepest optical follow-up
for this GW event.

Using  selection criteria which
isolate fading transients over the analysis region covering
3\% of the total localization probability,
we find no candidate counterparts.
We are still investigating improved background
rejection criteria using information such as:
matching against a galaxy catalog
to remove transients associated with high-redshift galaxies,
angular separation between $i$ and $z$ exposures to reduce
asteroids, and detailed simulations of supernovae
and source models
    to better optimize selection requirements as well
    as the search strategy for future events.

Although these results are not
surprising given the partial areal coverage and the likely
BBH merger nature of the
event, our search is a crucial first step and
demonstrates the viability of DECam
for deep optical follow-up of GW events.

\acknowledgments

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain,
the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing
Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago,
the Center for Cosmology and Astro-Particle Physics at the Ohio State University,
the Mitchell Institute for Fundamental Physics and Astronomy at Texas A\&M University, Financiadora de Estudos e Projetos,
Funda{\c c}{\~a}o Carlos Chagas Filho de Amparo {\`a} Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cient{\'i}fico e Tecnol{\'o}gico and
the Minist{\'e}rio da Ci{\^e}ncia, Tecnologia e Inova{\c c}{\~a}o, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energ{\'e}ticas,
Medioambientales y Tecnol{\'o}gicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh,
the Eidgen{\"o}ssische Technische Hochschule (ETH) Z{\"u}rich,
Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ci{\`e}ncies de l'Espai (IEEC/CSIC),
the Institut de F{\'i}sica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universit{\"a}t M{\"u}nchen and the associated Excellence Cluster Universe,
the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth,
SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A\&M University.

The DES data management system is supported by the National Science Foundation under Grant Number AST-1138766.
The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2012-39559, ESP2013-48274, FPA2013-47986, and Centro de Excelencia Severo Ochoa SEV-2012-0234.
Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) including ERC grant agreements
 240672, 291329, and 306478.

This research uses services or data provided by the NOAO Science Archive.
NOAO is operated by the Association of Universities for
Research in Astronomy (AURA), Inc. under a cooperative agreement
with the National Science Foundation.


\section{Introduction}
A ubiquitous class of problems found in many thermodynamics textbooks is the one where the student is required to find the equilibrium temperature of a thermodynamics system and to verify that the entropy change in the equilibrating process is consistent with the Second law \cite{Reif1965}. Here, we essentially turn the problem on its head and ask whether one can state and prove useful inequalities inspired by the form of the second law in specific physical situations. Thus, the Second law provides the motivation for studying certain integral inequalities which are expected to hold independently as mathematical results.
Treating the inequalities as purely mathematical results we generalize them to a larger class of functions and to higher number of dimensions. Throughout this paper we will adopt a strategy of generalization illustrated by the figure below:
\begin{figure}[ht]
	\begin{centering}
		\includegraphics[width = \textwidth]{flowchart.eps}
		\par\end{centering}
	\caption{Flow chart for generalizing Theorem 1
		\label{fig:1} }
\end{figure}

The inequalities derived are important due to their utility as inequality generators. The generalized versions also prove to be useful in extending well known results such as the weighted means inequality and in proving results like the Jensen inequality. We then finally show how these results can be interpreted as a proof of the equivalence of the Planck statement of the Second Law and the positivity of specific heats in certain special cases.

\section{MOTIVATION FOR THE REVERSIBLE ONE-DIMENSIONAL INEQUALITY}
Consider a closed system of N identical blocks each with a positive heat capacity $f(T)$ (some function of temperature) at initial temperatures $T_i$ where the index $i$ runs from 1 to N. If now all pairs of blocks are connected via a Carnot engines that reside inside the system then the system will finally reach an equilibrium temperature $T_0$ dictated by the constraint of no net creation of entropy in the system.\\
Symbolically:\\
$T_i>0,i = 1$ to $N$ \& $f(T) > 0 $ $\forall$ $T >0$
\begin{equation}
\sum_{i=1}^N\int_{T_i}^{T_o} \frac{{f(T)}}{T}dT = 0
\end{equation}

The fact that the net-work output of the engines must be non-negative leads to:
\begin{equation}
\sum_{i=1}^N\int_{T_i}^{T_o}f(T)dT \geq 0
\end{equation}
If this result is interpreted mathematically as an inequality concerning the integral for a positive function then questions about the nature of function for which the inequality works can be asked. Also, it is reasonable to demand a proof of the result that does not use any physics whatsoever for a pre-stated class of functions $f(T)$. We must also prove that the final equilibrium temperature lies between the maximum and the minimum temperatures.

\section{THE REVERSIBLE ONE-DIMENSIONAL INEQUALITY:}

\subsection*{Theorem 1}
Suppose ${x_i },i=1$ to $n$ is a sequence of non-decreasing positive real numbers and let  $f(x):[x_1,x_n] \to \mathbb{R}^+$ be a step-wise continuous function.
Define
\begin{equation}
F(y) \equiv \sum_{i=1}^N\int_{y}^{x_i} \frac{{f(x)}}{x}dx ; y \in [x_1,x_n] \nonumber
\end{equation}
Then there exists $x_0$ such that -
\begin{equation}
x_0 \in [x_1,x_n] \& F(x_0) = \sum_{i=1}^N\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx = 0
\end{equation}
And the following holds -
\begin{equation}
\sum_{i=1}^N\int_{x_0}^{x_i} f(x) dx \geq 0
\end{equation}
\subsubsection*{Proof:}
Since, every step-wise continuous function is Riemann- integrable on its domain, f(x) is also Riemann- integrable over $[x_1,x_n]$ \cite{Pugh2002}. Using the fact that indefinite integral of a Riemann-integrable function is continuous \cite{Rudin1976}, one can say that each integral inside the sum in eq(3) and hence F (y) is a well-defined continuous function of $y$. Since, \noindent f(x)/x$>0\ \forall \ x\epsilon [x_1,x_n]\ $

we have:
	\[F\left(x_n\right)\ge 0\ \&\ F\left(x_1\right)\le 0\]

And thus by application of Bolzano’s theorem \cite{Apostol1967} to $F(x)$  on $[x_1,x_n]$ we obtain eq(3). Now, $x_0$  partitions ${x_i}$ such that -
\begin{eqnarray}
	x_p\leq x_0 \leq x_{p+1};1\leq p \leq n-1 \nonumber  \\
	F(x_0) &=& \sum_{i=1}^N\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx \nonumber \\
	F(x_0) &=& \sum_{i=1}^p\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx + \sum_{i=p+1}^n\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx \nonumber \\
	\Rightarrow \sum_{i=1}^p\int_{x_i}^{x_0} \frac{{f(x)}}{x}dx &=& \sum_{i=p+1}^n\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx \\
	\because \sum_{i=p+1}^n\int_{x_0}^{x_i} \frac{{f(x)}}{x}dx &\leq& \sum_{i=p+1}^n\frac{1}{x_0}\int_{x_0}^{x_i}\frac{{f(x)}}{x}dx \nonumber \\
	\because \sum_{i=1}^p\int_{x_i}^{x_0} \frac{{f(x)}}{x}dx &\geq& \sum_{i=1}^p\frac{1}{x_0}\int_{x_i}^{x_0} \frac{{f(x)}}{x}dx \nonumber
\end{eqnarray}
We get from eq(5)
\begin{eqnarray}
	\sum_{i=p+1}^n\frac{1}{x_0}\int_{x_0}^{x_i}f(x)dx &\geq& \sum_{i=1}^p\frac{1}{x_0}\int_{x_i}^{x_0}f(x)dx \nonumber \\
	\Rightarrow	\sum_{i=p+1}^n\int_{x_0}^{x_i}f(x)dx &\geq& \sum_{i=1}^p\int_{x_i}^{x_0}f(x)dx \nonumber \\
	\Rightarrow \sum_{i=1}^n\int_{x_0}^{x_i}f(x)dx &\geq& 0 \nonumber
\end{eqnarray}

\subsubsection*{Comments:}
\begin{enumerate}
	\item The fact that only piece-wise continuity is required for the proof to go through has a nice physical consequence. If the heat capacity of the blocks is a step-function of temperature then blocks undergo phase transformation/s,which is a situation to which this analysis can be applied by the assumption of piece-wise continuity of $f(x)$.  Despite this eq(1) \& eq(2) will continue to hold and so will eq(3) \& eq(4).
	\item Notice that the proof also works if we replace the $x$ in the denominator with any continuous increasing function $g(x)$ which leads to theorem below.
	\item If the function $g(x)$ above is decreasing the sign of the inequality in eq(7) is reversed.
\end{enumerate}

\subsection*{Theorem 2}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals and \\ let  $f\left(x\right):[x_1,x_n]\to {\mathbb{R}}^+$ be a step-wise continuous function.

Let $g\left(x\right):[x_1,x_n]\to {\mathbb{R}}^+$be a continuous function such that:
\[g\left(a\right)\ge g\left(b\right)\forall \ a,b\epsilon \left[x_1,x_n\right]\ \&\ a>b\]

Define
\[F\left(y\right)\equiv \sum^n_{i=1}{\int^{x_i}_y{\frac{f\left(x\right)}{g(x)}dx\ ;\ y\epsilon [x_1,x_n]\ \ }}\]

Then there exists a $\ x_0$ such that,
\begin{equation}
	F\left(x_0\right)=\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{g(x)}dx=0\ \ \ \ \&\ \ }}x_0\epsilon [x_1,x_n]
\end{equation}

And the following holds:
\begin{equation}
	\sum^n_{i=1}{\int^{x_i}_{x_0}{f(x)dx\ge 0}}
\end{equation}

\section{PROOF OF RESTRICTED JENSEN INEQUALITY}

We now present a proof of the following restricted form of Jensen's inequality \cite{Jensen1906} using Theorem 2.
\subsection*{Restricted Jensen's inequality}
Let $h:{\mathbb{R}}^+\to {\mathbb{R}}^+$ be a doubly differentiable function such that
\begin{equation}
	h^{''}\left(x\right),h^{'}\left(x\right)>0\ \forall \ x>0
\end{equation}

And let${\{y}_j$$\}$, j=1 to m be a sequence of positive reals, then the following holds:
\begin{equation}
\frac{\sum^m_{j=1}{h(y_j)}}{m}\ge h\left(\frac{\sum^m_{j=1}{y_j}}{m}\right)
\end{equation}

\subsection*{Proof}
Consider a new sequence of positive reals$\{x_j$$\}$, j= 1 to m which is just a permutation of${\{y}_j$$\}$, in a non-decreasing order. Choose $f\ \&\ g\ $in Theorem 2 such that:
\begin{equation}
f\left(x\right)=g\left(x\right)=h'\left(x\right)\ \forall \ x\epsilon [x_{1,}x_m]
\end{equation}
Now $h^{''}\left(x\right)\ $ exists and is greater than 0 so that:
\begin{equation}
	f^{'}\left(x\right)=g^{'}\left(x\right)=h^{''}\left(x\right)>0\ \forall \ x\epsilon \left[x_{1,}x_m\right]
\end{equation}

Thus, such a choice of $f\ \&\ g$ can always be made consistent with the conditions of Theorem 2. We can say so because $f\ \&\ g$ are continuous since differentiability of a function implies continuity (from (10)) and because $g$ is increasing (from (10)). Apply Theorem 2 to$\{x_j$$\}$ with the given choice of functions.\\
By Theorem 2:
\begin{eqnarray}
	F\left(y\right)\equiv \sum^m_{j=1}{\int^{x_j}_y{\frac{f\left(x\right)}{g(x)}dx=\sum^m_{j=1}{(}\ x_j-y);\ y\epsilon [x_1,x_n]\ \ }} \nonumber \\   \end{eqnarray}
Thus
\begin{eqnarray}
 \ x_0=\ \frac{\sum^m_{j=1}{y_j}}{m} \nonumber
\end{eqnarray}

Now,
\begin{equation}
	\sum^m_{j=1}{\int^{x_j}_{x_0}{f\left(x\right)dx\ge 0}}
\end{equation}

But, from eq(10)
\begin{equation}
\sum^m_{j=1}{\int^{x_j}_{x_0}{f\left(x\right)dx=\ \sum^m_{j=1}{(h\left(x_j\right)-h(x_0)}}})
\end{equation}

Therefore, from eq(12) and eq(13), we get
\begin{equation}
	\frac{\sum^m_{j=1}{h(y_j)}}{m}\ge h\left(\frac{\sum^m_{j=1}{y_j}}{m}\right) \nonumber
\end{equation}

\section{THE ODD-DIMENSIONAL VERSION OF THE REVERSIBLE ONE-DIMENSIONAL INEQUALITY}

\subsection*{Theorem 3}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals. Consider n    k-tuples of reals:
\[(x_i,x_i,\dots ,x_i)\ ;i=1\ to\ n\ \&\ k\ is\ odd\ as\ well\ as\ positive\ \]

And let  $f:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be a continuous function.

Let $g:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be a continuous function that is non-decreasing with respect to each of its arguments.

Define
\[F\left(y\right)\equiv \sum^n_{i=1}{\int^{x_i}_y{\dots .}\int^{x_i}_y{\ \ }\int^{x_i}_y{\frac{f}{g}dz_1\ dz_2\dots dz_k;\ y\epsilon [x_1,x_n]\ \ }}\]

Then there exists $x_0$ such that
\begin{equation}
	F\left(x_0\right)=0\ \&\ x_0\epsilon [x_1,x_n]
\end{equation}

And the following holds:
\begin{equation}
	\sum^n_{i=1}{\int^{x_i}_{x_0}{\dots .}\int^{x_i}_{x_0}{\ \ }\int^{x_i}_{x_0}{fdz_1\ dz_2\dots dz_k\ \ge 0\ \ }}
\end{equation}

\subsubsection*{Proof:}
The proof of Theorem 3 is identical to the proof of Theorem 2 except for the fact that multiple integrals replace the single integral. The oddness of k is exploited while writing the expression corresponding to eq (6).

\subsubsection*{Comments:}
\begin{enumerate}
	\item  The theorem does not hold for even values of k.

	\item  The sign of inequality is reversed if g is decreasing with respect to each of its arguments.
\end{enumerate}
\subsection*{Counter-example in even dimensions}
Take $\ k=n=2$ and$\ f=g=z_1$.

Then the equation for $x_0\ $turns out to be:
\[2x^2_0-2x_0\left(x_1+x_2\right)+\left(x^2_1+x^2_2\right)=0\]

Which has no real roots if$\ x_2\neq x_1$. Thus $x_0$ fails to exist and so the theorem collapses.

\section{GENERALIZATION OF  WEIGHTED MEANS INEQUALITY}
We propose the following generalization of weighted means inequality \cite{Hardy1952}:
\subsection*{Corollary 1}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals. Let k be an odd positive number and b \& a be any two reals such that$\ a\ \ge b\ \&\ a,b\neq 0,-1$ . Then the following holds:
\[\exists \ x_0\ \in {\mathbb{R}}^+\ such\ that\ \sum^n_{i=1}{{(x^b_i-x^b_0)}^k\ =0\ }\]
And
\[\sum^n_{i=1}{{(x^a_i-x^a_0)}^k\ \ge 0\ }\]
The sign of inequality flips if$\ a\ \le b$.
\subsubsection*{Proof:}
Choose $\ f\ \&   g$ in Theorem 3 as follows:
\[f\left(x_1,x_2,\dots .,x_n\right)=x^{a-1}_1x^{a-1}_2\dots .x^{a-1}_n\]
And
\[g\left(x_1,x_2,\dots .,x_n\right)=x^{a-b}_1x^{a-b}_2\dots .x^{a-b}_n\]
Clearly f and g have all the properties required for Theorem 3 to work and thus we obtain Corollary 1.

\noindent If $a\ \le b$ then g turns decreasing in each of its arguments and thus the sign of inequality flips as mentioned in Comment number 2 on the proof of Theorem 3.

\noindent The case where either a or b equals 0 or -1 can be handled similarly as above however because the integrals now involve natural logarithms instead on just polynomials the final form of the result changes from Corollary 1. The case where k is unity corresponds to the Weighted means inequality.

\section{GENERALIZATIONS OF THEOREM 1 TO CASES WITH IRREVERSIBLE HEAT EXCHANGE}

\subsection*{Theorem 4}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals and let  $f\left(x\right):[x_1,x_n]\to {\mathbb{R}}^+$ be a step-wise continuous function and $x_0\ $be positive real

Then the following holds:
\begin{equation}
	\frac{1}{x_0}\sum^n_{i=1}{\int^{x_i}_{x_0}{f(x)dx\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f(x)}{x}}}}}\ dx
\end{equation}

\subsubsection*{Proof:}
We split the theorem into 3 cases.
\begin{itemize}
	\item \textbf{Case 1: $x_0<x_1\ $}
	\begin{eqnarray}
		\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}}dx&\ge& 0 \nonumber \\
		\Longrightarrow x_0\int^{x_i}_{x_0}{\frac{f(x)}{x_0}}dx&\geq& x_0\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}}dx\ \forall \ i=1\ to\ n\ \nonumber \\
		\Longrightarrow \frac{1}{x_0}\sum^n_{i=1}\int^{x_i}_{x_0}f(x)dx &\geq& \sum^n_{i=1}\int^{x_i}_{x_0}\frac{f(x)}{x}dx
	\end{eqnarray}

	\item \textbf{Case 2: $x_n<x_0$ }
	\begin{eqnarray}
	\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}}dx&\le& 0 \nonumber \\
	\Longrightarrow x_0\int^{x_i}_{x_0}{\frac{f(x)}{x_0}}dx&\geq& x_0\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}}dx\ \forall \ i=1\ to\ n\ \nonumber \\
	\Longrightarrow \frac{1}{x_0}\sum^n_{i=1}\int^{x_i}_{x_0}f(x)dx &\geq& \sum^n_{i=1}\int^{x_i}_{x_0}\frac{f(x)}{x}dx
	\end{eqnarray}

	\item \textbf{Case 3:} $x_1<x_0<x_n$ \\
	Now, $x_0\ $partitions $\left\{x_i\right\}$ such that
	\[{x_p\le x}_0\le x_{p+1\ };1\le p\le n-1\]

	\begin{eqnarray}
	\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}dx\ \ \ \ }}
	&=&\sum^p_{i=1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}dx+\sum^n_{i=p+1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}dx\ \ \ }}\ \ \ }} \nonumber \\
	\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}dx\ \ \ \ }}&\le& \ \frac{1}{x_0}\sum^p_{i=1}{\int^{x_i}_{x_0}{f\left(x\right)dx}}+\frac{1}{x_0}\sum^n_{i=p+1}{\int^{x_i}_{x_0}{f\left(x\right)dx}}\ \ \ \ (By\ \left(18\right)\ \&\ \left(19\right)) \nonumber \\
	\Longrightarrow \frac{1}{x_0}\sum^n_{i=1}\int^{x_i}_{x_0}f(x)dx&\ge& \sum^n_{i=1}\int^{x_i}_{x_0}\frac{f(x)}{x}\ dx\
	\end{eqnarray}

\end{itemize}

\subsubsection*{Physical interpretation}
This describes the physical situation where the N identical blocks in the system are connected to one another and allowed to attain thermal equilibrium. It is not necessary that they be connected by Carnot engines. Since this process is irreversible, the entropy of the universe must increase \cite{Planck1903}. This can be proven form eq(17) as follows:

\[\Delta Q=0\Longrightarrow \sum^n_{i=1}{\int^{x_i}_{x_0}{f\left(x\right)dx}}=0\]
Using eq(17),
\[ -\Delta S=\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f\left(x\right)}{x}}}dx\ \le 0\ \ \]
\subsubsection*{Comments:}
\begin{enumerate}
	\item  Notice that the proof also works if we replace the x in the denominator with any continuous increasing function g(x) which leads to theorem below.
	\item  If the function g(x) above is decreasing then sign of the inequality in eq(21) is reversed.
	\end{enumerate}
\subsection*{Theorem 5}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals and let  $f\left(x\right):[x_1,x_n]\to {\mathbb{R}}^+$ be a step-wise continuous function and $x_0\ $be positive real. Let $g\left(x\right):[x_{min},x_{max}]\to {\mathbb{R}}^+$be a continuous function such that:

$g\left(a\right)\ge g\left(b\right)\forall \ a>b\ $

In its domain

Where $x_{min}={\mathrm{min} \left\{x_0,x_n\right\}\ }\ \&\ \ x_{max}={\mathrm{max} \left\{x_0,x_n\right\}\ }$

Then the following holds:
\begin{equation}
\frac{1}{g(x_0)}\sum^n_{i=1}{\int^{x_i}_{x_0}{f(x)dx\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f(x)}{g(x)}}}}}\ dx
\end{equation}

The above theorem can be generalized to odd dimensions as follows:
\subsection*{Theorem 6}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals. Consider n    k-tuples of reals:
\[(x_i,x_i,\dots ,x_i)\ ;i=1\ to\ n\ \&\ k\ is\ odd\ as\ well\ as\ positive\ \]

And let  $f:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be a continuous function.

Let $g:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be a continuous function that is non-decreasing with respect to each of its arguments. Let $x_0\ $be positive real, then following holds:
\begin{equation}
\sum^n_{i=1}{\int^{x_i}_{x_0}{\dots .}\int^{x_i}_{x_0}{\ \ }\int^{x_i}_{x_0}{\frac{f}{g(x_0,x_0,\dots ,x_0)}dz_1\ dz_2\dots dz_k\ \ }}\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\dots .}\int^{x_i}_{x_0}{\ \ }\int^{x_i}_{x_0}{\frac{f}{g}dz_1\ dz_2\dots dz_k\ \ }}
\end{equation}

\subsubsection*{Proof:}
The proof of Theorem 6 is identical to the proof of Theorem 4 except for the fact that multiple integrals replace the single integral. The oddness of k is exploited while writing the expression corresponding to$\int^{x_i}_{x_0}{\frac{f\left(x\right)}{g(x)}}dx\le 0$ in the second case where  $x_n<x_0.$

\subsubsection*{Comments}
\begin{enumerate}
	\item  The theorem does not hold for even values of k.

	\item  The sign of inequality is reversed if g is decreasing with respect to each of its arguments.
\end{enumerate}
\subsection*{Theorem 7}
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals and $x_0\ $be positive real let  $f_i\left(x\right):[\mathrm{min}\mathrm{}(x_0,x_i),\mathrm{max}\mathrm{}(x_0,x_i)]\to {\mathbb{R}}^+$ be n step-wise continuous functions

Then the following holds:
\begin{equation}
\frac{1}{x_0}\sum^n_{i=1}{\int^{x_i}_{x_0}{f_i\left(x\right)dx\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{x}}}}}\ dx
\end{equation}

If $g\left(x\right):[x_{min},x_{max}]\to {\mathbb{R}}^+$be a continuous function such that:

$g\left(a\right)\ge g\left(b\right)\forall \ a>b\ $  In its domain

Where $x_{min}={\mathrm{min} \left\{x_0,x_n\right\}\ }\ \&\ \ x_{max}={\mathrm{max} \left\{x_0,x_n\right\}\ }$

Then the following holds:
\begin{equation}
	\frac{1}{g(x_0)}\sum^n_{i=1}{\int^{x_i}_{x_0}{f_i\left(x\right)dx\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{g(x)}}}}}\ dx
\end{equation}

\subsubsection*{Proof:}
The result follows from the fact that  \\
\[\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{x_0}}dx\ge \int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{x}}dx\ \] \& \[\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{g{(x}_0)}}dx\ge \int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{g(x)}}dx\ \ \ \forall \ i=1\ to\ n\ \]

Which holds irrespective of whether $x_i$ is greater than, less than or equal to$\ x_0$.

\subsubsection*{Physical interpretation}
Physically this corresponds to the situation in which the different blocks labelled by index i and having different heat capacities ($f_i\left(x\right))$ and initial temperatures ($x_i)\ $ are connected to an infinite reservoir at temperature  $x_{0.}$ . The net entropy of universe is expected to increase from the second law of thermodynamics and this can be proven using eq(23).
\begin{eqnarray}
	\Delta S\ for\ reservoir &=&\ \frac{1}{x_0}\sum^n_{i=1}{\int^{x_i}_{x_0}{f_i\left(x\right)dx}} \nonumber \\
\&\ \Delta S\ for\ all\ blocks&=&-\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{x}}}dx\ \nonumber \\
\Longrightarrow \ \Delta S\ for\ universe &=&
\frac{1}{x_0}\sum^n_{i=1}{\int^{x_i}_{x_0}{f_i\left(x\right)dx}}-\sum^n_{i=1}{\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{x}}}dx\ge 0\ (from\ 23) \nonumber
\end{eqnarray}

\subsection*{Theorem 8 }
Suppose $\left\{x_i\right\},\ i=1\ to\ n$ is a sequence of non-decreasing positive reals. Consider n    k-tuples of reals:
\[(x_i,x_i,\dots ,x_i)\ ;i=1\ to\ n\ \&\ k\ is\ odd\ as\ well\ as\ positive\ \]

And let  $f_i:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be n continuous functions.

Let $g:[\mathrm{(}x_1,x_1,\dots .,x_1),\mathrm{(}x_n$$,x_n,\dots .,x_n)]\to {\mathbb{R}}^+$ be a continuous function that is non-decreasing with respect to each of its arguments. Let $x_0\ $be positive real, then following holds:
\[\sum^n_{i=1}{\int^{x_i}_{x_0}{\dots .}\int^{x_i}_{x_0}{\ \ }\int^{x_i}_{x_0}{\frac{\mathrm{\ }f_i}{g(x_0,x_0,\dots ,x_0)}dz_1\ dz_2\dots dz_k\ \ }}\ge \sum^n_{i=1}{\int^{x_i}_{x_0}{\dots .}\int^{x_i}_{x_0}{\ \ }\int^{x_i}_{x_0}{\frac{\mathrm{\ }f_i}{g}dz_1\ dz_2\dots dz_k\ \ }}\]
\subsubsection*{Proof:}
The proof of Theorem 8 is identical to the proof of Theorem 7 except for the fact that multiple integrals replace the single integral. The oddness of k is exploited while writing the expression corresponding to
\begin{equation}
\int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{g{(x}_0)}}dx\ge \int^{x_i}_{x_0}{\frac{f_i\left(x\right)}{g(x)}}dx \nonumber
\end{equation}
for the case where  $x_i<x_0.$
\section{ EQUIVALENCE OF POSITIVITY OF SPECIFIC HEATS AND THE SECOND LAW OF THERMODYNAMCIS:}

 We will now prove the equivalence of the positivity of the specific heats and of the Planck's statement of the second law of thermodynamics in some special cases

 Consider the two statements P \& Q as follows:

P: all statistical systems have positive specific heat.

Q: Every process occurring in nature proceeds in the sense in which the sum of the entropies of all bodies taking part in the process is increased.

\subsubsection*{Proof of Equivalence of P \& Q:}

  Proof that$\ \ P\to Q$:

 The theorems until now imply that in the scenarios considered the statement P implies the statement Q.

Proof that$\ Q\to P$:

 Instead of proving that $Q\to P\ $we will prove that $\sim P\to \sim Q$ by performing a thought experiment as follows:

Since P is not true consider a block with negative heat capacity, let say --C (where C $>$ 0) at initial temperature  $T_1$ . Put it in thermal contact with a block of heat capacity 2C at initial temperature$\ T_2$. Now close the system consisting of two blocks using adiabatic walls. The system will reach an equilibrium temperature\cite{Bell1998} which can be calculated using the first law of thermodynamics:

 Heat lost/gained by block with heat capacity 2C: 2$C(T_{eq}-T_2)$

 Heat lost/gained by block with heat capacity -C: -$C(T_{eq}-T_1)$

Therefore, we get from the first law of thermodynamics:
\begin{equation} \label{GrindEQ__24_}
T_{eq}=2T_2-T_1
\end{equation}
Change in entropy of the block with heat capacity -C:
\begin{equation} \label{GrindEQ__25_}
\Delta S_1=\int^{T_{eq}}_{T_1}{\frac{-C}{T}}dT\ \
\end{equation}
Change in entropy of the block with heat capacity 2C:
\begin{equation} \label{GrindEQ__26_}
\Delta S_2=\int^{T_{eq}}_{T_2}{\frac{2C}{T}}dT
\end{equation}

 Thus, the total change in entropy of the system is:

\[\Delta S=Cln({\left(2T_2-T_1\right)T_1}/{T^2_2)\ }\]

 Which is negative when  $T_2>(\sqrt{2}-1)T_1$ . Therefore, the negation of statement P implies the negation of statement Q, which completes the proof of equivalence of P and Q.

\section{Discussion}
We would like to remark that the Theorems work even if the function f is just Lebesgue integrable and thus the given versions are a restricted form of the inequalities; however the fact that the function might be discontinuous over a set of measure zero hardly makes any difference as regard to its intended use.

We have been unable to figure out any physical interpretation of either the odd dimensional inequalities or of the Lebesgue integrable form of Theorems suggested above. We would like to welcome everyone to tackle this issue of interpretation.

In each of the above results, we have assumed that : `All statistical systems have positive specific heat' \cite{Sch1952} \ \cite{Andersen2014},holds and then proved what was expected on the basis of the Planck's statement of Second Law which states: `Every process occurring in nature proceeds in the sense in which the sum of the entropies of all bodies taking part in the process is increased'.   Thus, the results derived here can be interpreted as a mathematical proof of the equivalence of Second law of thermodynamics and of
positivity of specific heats in various cases. The theorems are also a remarkable way of churning out new inequalities by feeding them with the appropriate functions as has been demonstrated in the paper.

\begin{acknowledgments}

The authors would especially like to thank Dr.Subroto Mukerjee for discussions on sections VII and VIII of the paper.The authors would like to thank Dr.Gautam Bharali and Mr.Vivek Wagh for their valuable suggestions. We would also like to thank Mr.Salil Apte for typesetting the draft and Mr.Samyak Jha for proofreading it.
\end{acknowledgments}


\section{Introduction}
\label{sec:Intro}
\input{Intro.tex}

\section{Convolutional Tables Ensemble}
\label{sec:Method}
\input{Method.tex}

\section{Empirical results}
\label{sec:EmpiricalResults}
\input{Empiric.tex}

\section{Conclusions and further work}
\label{sec:Conclusions}

We introduced improvements to the convolutional tables framework in terms of bit functions used, word calculator structure, calculator optimization and global optimization. We have shown that for highly computational constrained tasks CTE may provide accuracy higher than CNNs. A natural direction for future research is to replace the flat structure of CTEs with a layered approach, in order to try and enjoy the accuracy of CNNs with the speed of CTEs.

{\small
\bibliographystyle{ieee}

\subsection{Implementation and data details}
\label{sub:Setup}

The experiments were conducted on $4$ publicly available datasets: MNIST, CIFAR-10, SVHN and 3-HANDPOSE. The first three are standard recognition benchmarks in gray-scale (MNIST) or RGB (CIFAR-10,SVHN), with $10$ classes each. 3-HANDPOSE are a $4$-class dataset, with $3$ hand poses and a fourth class of 'other', and its images contain depth and IR channels. The image sizes are between $28 \times 28$ (MNIST) and $36 \times 36 $ (3-HANDPOSE). The training set size ranges from $50000$ (CIFAR-10) to $604000$ (SVHN).

CTE training code was written in Matlab, with some routines using code from the packages~\cite{PMT2006,Fan2008,vedaldi15matconvnet}. The test time classier was implemented and optimized in C.
For ferns we implemented algorithm~\ref{alg:CTE_test} with SSE operations. Words are computed for $8$ neighboring pixels together, and voting is done for $8$ classes at once.
For trees
we implemented a program generating efficient code of the bit computation loop for a specific tree, so the tree parameters are part of the code. This obtained an acceleration factor of ~$2X$ over standard C code.
We also thread-parallelized the algorithm over the convolutional tables, with good a speed-up of $3.6\times$ obtained from $4$ cores. However, we report and compare single thread performance to keep the methodology as simple as possible.

CNN models were trained using MatConvNet~\cite{vedaldi15matconvnet}. The implementation is efficient, reported to be comparable to Caffe~\cite{Caffe2014} in~\cite{vedaldi15matconvnet}, with the convolutional and global layers reduced to matrix multiplication done using an SSE-optimized BLAS package. When measuring execution time, we measured net run time of the convolutional, pooling and global layers alone, without Matlab overhead.
Time measurements were made on a Lenovo Thinkpad W530 quad core laptop, with i7-3720QM core running at 2.6Ghz.

\begin{figure*}
\begin{tabular}{cccc}
\includegraphics[scale=0.09]{MNist_1Thread_AccuracySpeed}  &
\includegraphics[scale=0.09]{Cifar10_1Thread_AccuracySpeed}     &
\includegraphics[scale=0.09]{SVHN_1Thread_AccuracySpeed}  &
\includegraphics[scale=0.09]{Durango_1Thread_AccuracySpeed}     \\
MNIST & CIFAR-10 & SVHN & 3-HANDPOSE
\end{tabular}
 \caption{\footnotesize \textbf{Speed-accuracy trade-offs curves.}
X axis states classifier run time in microseconds, in log scale base $10$. Each point states the (speed,acuarcy) result of a single trained classifier. the Lines are the lower envelope of classifiers from the CTE and CNN family. The number of classifiers trained over the data is 172 (MNIST), 283 (CIFAR-10), 89 (SVHN) and 111 (3-HANDPOSE).
}
\label{fig:TradeOffs}
\end{figure*}

\subsection{Comparison and variation}
\label{sub:Ablations}

\textbf{Comparison with DFE:} The Discriminative Ferns Ensemble (DFE) was suggested in~\cite{DFE} for classification of 3-HANDPOSE, and can be seen a baseline for CTE, which enhances it in many aspects. The first two columns in Table~\ref{table:Ablations} present errors of DFE and CTE on the $4$ datasets, using $50$ ferns for MNIST, SVHN, 3-HANDPOSE and $100$ for CIFAR-10. MNIST was trained with softmax distillation loss (see below for details), and the others with SVM loss. The aggregation area ${\{A^m\}}_{m=1}^M$ were chosen to be identical for all tables in a classifier, forming a centered square occupying most of the image.
To enable the comparison, M-class error rates are extracted from DFE (in~\cite{DFE} such errors are not reported, and $3$-class average true positive rates are reported instead). It can be seen that CTE base provides significant improvements of $24-45\%$ error reduction over DFE, with $28\%$ obtained for 3-HANDPOSE, where DFE was originally applied. Note that the CTE base is not the best choice for 3-HANDPOSE. With additional parameter tuning result of $2\%$ can be obtained with $50$ ferns, which is an improvement of $38\%$ over DFE.

\textbf{Ablation experiments:} The accuracy obtained by a CTE is influenced by many small incremental improvements related to structural and algorithmic variations. Columns 3-9 in Table~\ref{table:Ablations} show the contribution of some ingredients by removing them from the baseline CTE. For MNIST, where the effects are small due to the low error, results were averaged over $5$ experiments varying in their random seed, a with seed-induced std of $0.1\%$. It can be seen that these ingredients consistently contribute to accuracy for non-depth data.

\textbf{Trees/Ferns trade-off:} The trade-off between ferns and trees for MNIST and CIFAR-10 is presented in Figure~\ref{fig:MoreTables}(Right). For MNIST, the results were averaged over $5$ experiments, with a seed induced std of $0.028\%$. It can be seen that trees provide better accuracy. However, the speed cost of using trees is significant, due to the inability to efficiently vectorize their implementation.

\textbf{Distillation experiments:} We experimented with knowledge distillation from a CNN to a CTE using the method suggested in~\cite{Hinton2015}. In such experiments, soft labels are taken from our best CNN model, and a CTE is trained to optimize a convex combination of the standard softmax loss and the Kullback-Leibler distance from the CNN-induced probabilities.
We attempted this for MNIST and CIFAR-10 using our best CNN models, providing $0.31\%$ and $8.6\%$ error  respectively as distillation sources. For MNIST, this training methodology proved to be successful. Averaging over $5$ seeds, the accuracy of a $50$-fern CTE optimized for softmax was $0.66\%$(the std was $0.025\%$) without distillation, and $0.45\%(0.029\%)$ with distillation. For comparison, an SVM-optimized CTE with the same parameters obtained $0.61\%(0.04\%)$ error. For CIFAR-10 distillation did not consistently improve the results.

\subsection{Speed-Accuracy trade-off}
\label{sub:AS}

We are interested in the trade-off or Pareto curves~\cite{Climacao1997}, showing the best accuracy obtainable for a specific speed constraint and vice versa. Since the design space for variations of $CNNs$ and $CTE$ algorithms is huge, and the training time of the algorithms is considerable, we needed to sample it wisely to get a good curve approximation. Our sampling technique is based on two stages. In stage 1 we searched for the most accurate classifiers for CTE and CNN with loose speed constraints, so even slow classifiers
were considered. We then used the few top accuracy variants of each architecture as baselines and accelerated them by systematically varying certain design parameters.

Our CNN baseline architectures are variations of DeepCNiN(l,k)~\cite{Graham2014}, with $l=3-4$ convolutional layers and $k=60-100$, implying usage of $i\cdot k$ maps at the $i$-th layer. It was shown in~\cite{Graham2014} that higher $l,k$ values provide better accuracy, but such architectures are much slower than $1$ CPU millisecond and so they are outside our domain of interest. We experimented with dropout~\cite{Srivastava2014}, parametric RELU units~\cite{He2015}, affine image transformations following~\cite{Graham2014}, and HSV image transformations following~\cite{Snoek2015}. Acceleration of the baseline architectures used three main parameters. The first was reducing parameter $k$ controlling the network width. The second was reduction of the number of maps in the output of the NIN layers. This reduces the number of input maps for the next layer, and can dramatically save computation with relatively small loss of accuracy. The third was raising the convolution stride parameter from $1$ to $2$. For CTEs, our exploration space was sketched in Section~\ref{sec:Method}, and it includes both ferns and trees. The best performing configurations were then accelerated using a single parameter: the number of tables in the ensemble.

Trade-off graphs for the $4$ datasets are shown in Figure~\ref{fig:TradeOffs}. Classification speed in microseconds is displayed along the $X$-axis in log scale with base $10$. For all datasets, there is a high speed regime where CTEs provide better accuracy than CNNs. Specifically CTEs are preferable for all datasets when less than $100$ microseconds are available for computation.
Starting from $1000$ microseconds and up CNNs are usually better, with CTEs still providing comparable accuracy for MNIST and 3-HANDPOSE at the $1-10$ milliseconds regime.
Viewed conversely, for a wide range of error rates, if the error rate is obtainable by a CTE, it is obtainable with significant speedups over CNNs. Some examples of this phenomenon are given in Figure~\ref{fig:MoreTables}(Left). Note that while a working point of $0.25$ error for CIFAR-10 may seem high, the majority of the one-versus-one errors of such a classifier are lower than $0.05$, which may be good enough for many purposes.

\subsection{Notation and classifier structure}
\label{sub:Struct}

A convolutional table ensemble is a classifier $f:I \to {\{1,..,C\}}$ where $I\in \mathbb{R}^{ S_x\times S_y \times D}$ and $C$ is the number of classes. The image $I$ may undergoes a preparation stage where additional feature channels ('maps') may be added to it, so it is transformed to $I^e \in \mathbb{R}^{ S_x\times S_y \times D_e}$ with $D_e\ge D$.
After preparation, the ensemble sums the votes of $M$ convolutional tables, each in turn sums the votes of a word calculator over all pixels in a aggregation area. We now explain this process bottom-up.

\textbf{Word calculator:} For an image location $p=(x,y) \in \{1,..,S_x\} \times \{1,..S_y\}$, we denote
its neighborhood by $N(p)$, and by $I_{N(p)}$ the patch centered at location $p$.
A word calculator is a feature extractor applied to such a patch and returning a $K$-bit index, i.e. a function $W:\mathbb{R}^{l\times l \times D_e} \to \{0,1\}^K$, where $l$ is the neighborhood size. In a classifier we have $M$ calculators denoted by $B^m = B(\cdot ~; \Theta^m)$ where $m=1,..,M$ and $\Theta^m$ are the calculator parameters. The calculator computes its output by applying $K$ bit-functions to the patch, each producing a single bit.
Several types of bit functions are discussed in Section~\ref{sub:StructDetails}.

\textbf{ Convolutional table:} Each word calculator is applied to all locations in an integration area, and each word extracted casts votes for the $C$ output classes. Convolutional table $m$ is hence a triplet $(B^m,A^m,W^m)$ where $A^m \subset \{1,..,S_x\} \times \{1,..S_y\}$ is the integration area of word calculator $B^m$, and $W^m \in \mathbb{M}_{C \times 2^K }$ is its weight matrix. Denote by $b_p^m =B(I^e_{N(p)} ; \Theta^m)$ the word calculated at location $p$. We gather a histogram $H^m=(H^m_0, .., H^m_{2^K-1})$ counting word occurrences; i.e.,

\begin{equation}
H^{m}_b=\sum_{p\in A^{m}}\delta(b_{p}^{m}-b)     \label{eq:Table}
\end{equation}
with $\delta$ the discrete delta function. The class support of the convolutional table is the $C$-element vector $W^m (H^m)^t$.

\textbf{ Convolutional tables ensemble:} The ensemble classification is done by accumulating the class support of all convolutional tables into a linear classifier with a bias term. Let $H = [ H^1, .., H^M]$ , and $W=[ W^1,..W^M] \in \mathbb{R}^{C \times M2^K}$. The classifier's decision is given by
\begin{equation}
C^{*} =  \arg\max_c W H^t - T^t    \label{eq:classifier}
\end{equation}
where $T=(T^1,..T^C)$ is a vector of class biases. Algorithm box~\ref{alg:CTE_test} shows the classifier's test time flow. Note that the histograms $H^m$ are not accumulated in practice, and instead each word computed directly votes for all classes.

\begin{algorithm}[t]
\caption{\textbf{\label{alg:training}} Convolutional Tables Ensemble: Training}

{

\textbf{\small{{Input}}}{\small{{: A labeled training set $S=\{I_{i},y_{i}\}_{i=1}^{N}$}}}{\small \par}

{\small{{$~~~~~~~~~~~$ Parameters $M,K,C, {\{A^m\}}_{m=1}^M$ , convex loss $L(S ; W,T)$ }}}{\small \par}

\textbf{\small{{Output}}}{\small{{: A classifier $(\Theta^{m},A^{m},W^{m,c},T^c)_{m=1,c=1}^{M,C}$      }}}{\small \par}

{\small{{Initialization: $g_i^c=1/|\{I_{i}|y_{i}^c=1\}|~~$ if $y_{i}^c=1$,}}}{\small \par}

{\small{{$~~~~~~~~~~~~~~~~~~~~~~$ $g_i^c=-1/|\{I_{i}|y_{i}^c=-1\}|$
if $y_{i}^c=-1$}}}{\small \par}

{\small{{For $m=1,..,M$}}}{\small \par}

{\small{{$~~~~~~$ Table addition: choose  $\Theta^{m}$ to optimize: }}}{\small \par}

{\small{{$~~~~~~$    $\Theta^{m} = \arg\max_{\Theta}  \sum_c \sum_{b\in \{0,1\}^K} | \sum_{\{i,p: b_{i,p}(\Theta)=b\}} g_i^c |$  }}}{\small \par}

{\small{{$~~~~~~$ Update representation: $~~~ \forall i=1,..N, b\in{\{0,1\}}^K $  }}}{\small \par}

{\small{{$~~~~~~~~~~~~$  $H^m_b[I_i]= \sum_{p}\delta(b^m_{p,i}=b)$, $H=[ H, H^m ]$     }}}{\small \par}

{\small{$~~~~~~$   Global optimization:  train $W,T$ by solving  }}{\small \par}

{\small{$~~~~~~~~~~~~~$  $ \arg\min_{W,T} L( {\{H_i,y_i\}}_{i=1}^N ; W,T)$  }}{\small \par}

{\small{$~~~~~~$   If $m<M$ get loss gradients:  $ g_i^c = \frac{dL}{dS_i^c}  $    }}}{\small \par}

\end{algorithm}

\subsection{Training}
\label{sub:Train}

In~\cite{DFE,Ren_2014_CVPR} instances of convolutional tables ensemble were discriminatively optimized for specific tasks and losses (hand pose recognition using SVM in~\cite{DFE}, and face alignment using $l^2$ regression in~\cite{Ren_2014_CVPR}). The main idea behind these methods is to iterate between solving a convex problem for a fixed representation, and augmenting the representation based on gradient signals from the obtained solution. Here we adapt these ideas to linear $M$-classification with an  arbitrary $l_2$-regularized convex loss function, using techniques from~\cite{Mason2000,BarHillel2005}. Assume a labeled training sample with fixed representation ${\{(H_i,y_i)\}}_{i=1}^{N}$ where $H_i \in \mathbb{R}^{m2^K}$, $y_i \in \{1,..,C\}$, and denote the $c$-th row of the weight matrix $W$ by $W_c$. We want to learn a linear classifier of the form $C^{*}=\arg\max_c s^c$ with $s^c=W_c H^t-T^c$ by minimizing a sample loss function of the form
\begin{equation}
L({\{H_i,y_i\}}_{i=1}^N) = \frac{1}{2}||W||^2 + \sum_{i=1}^N l({\{s_i^c\}}_{c=1}^C, y_i)    \label{eq:Loss}
\end{equation}

with $l({\{s^c\}}_{c=1}^C, y)$ a convex function of $s^c$.
$L$ is strictly convex with a single global minimum, hence solvable using known techniques. Once the problem has been solved for the fixed representation $H_i$,
 we want to extend the representation by incorporating a new table, effectively adding $2^K$ new features. In order to choose the new features wisely, we consider how the loss changes if a new feature $f^{+}$ is added to the representation with small class weights, regarded as a small perturbation of the existing model.

Denote by $f^{+}_i$ the value of a new feature candidate for example $i$.After incorporating the new feature, example $i$'s representation changes from $H_i$ to $H_i^+ = [ H_i, f^{+}_i]$ and weights vectors $W_c$ are augmented to $[W_c, w^+_c]$ with $w^+_c\in \mathbb{R}$. Class scores $s_i^c$ are updated to $s_i^{c,+} = W^+_c (H_i^{+})^t -t^c = s_i^c + w^+_c f^{+}$. Finally, the loss is changed to $L^+ = L({\{H_i^+,y_i\}}_{i=1}^N) = \frac{1}{2}||W||^2 + \frac{1}{2}\sum_{c=1}^C {w^+_c}^2 + \sum_{i=1}^N l({\{s_i^{c,+}\}}_{c=1}^C, y_i)$.
Denote $W^{+}=[w^+_1,..,w^+_C]$ the new weights vector. We assume that the new feature is added with small weights; i.e., $w^+_c\le \epsilon$ for all $c$. $L^+$ can be Taylor approximated around $W^{+}=0$, with the gradient $\frac{dL^+}{dW^{+}} |_{W^{+}=0}$:

\begin{equation}
\begin{aligned}
	\frac{dL^+}{dW^{+}} \bigg|_{W^{+}=0}  & = w^{+} + \sum_{i=1}^N \frac{dl({\{s_i^{c,+}\}},y_i)}{ds_i^{c,+}} f^{+}_i 	\bigg|_{W^{+}=0}   \\
	 &	=			\sum_{i=1}^N \frac{dl({\{s_i^{c}\}},y_i)}{ds_i^{c}} f^{+}_i  				\label{eq:LossGrad}
\end{aligned}
\end{equation}

Using the gradient in a Taylor approximation of $L^+$ gives
\begin{equation}
\begin{aligned}
	L^+  &= L + W^+ 	(\frac{dL^+}{dW^{+}})^t + O(||W^+||^2) \\			\label{eq:LossApprox}
	    &= L + \sum_{c=1}^C w^+_c \sum_{i=1}^N  \frac{dl({\{s_i^{c}\}},y_i)}{ds_i^{c}} f^{+}_i  + O(||W^+||^2)
\end{aligned}
\end{equation}

Denote $g_i^c = \frac{dl({\{s_i^{c}\}},y_i)}{ds_i^{c}}$. For loss minimization we want to minimize $\sum_{c=1}^C w^+_c \sum_{i=1}^N  g_i^c f^{+}_i$ over $W^+$ and $f^+$. For fixed $f^+$ minimizing over $W^+$ is simple. Denoting $R^c(f^+)=\sum_{i=1}^n g_i^c f_i^+$, we have to minimize $\sum_{c=1}^C w^+_c R^c(f^+)$ under the constraint $w^+_c\le \epsilon,~ \forall c$. We can minimize each term in the sum independently to get $w^{+,c}_{opt} = -\epsilon sign(R^c)$, and the value of the minimum is $-\epsilon \sum_{c=1}^C |R^c(f^+)|$. Hence, for a single feature addition, we need to maximize the score $\sum_{c=1}^C |R^c(f^+)|$.

To return to our scenario, we add $2^K$ features at once, generated by a new word calculator $B$.
The derivation above can be done for each of them independently, so for the addition of the features ${\{H^+_b\}}_{b\in \{0,1\}^K}$ we get
\begin{equation}
\begin{aligned}
	L^+  & \approx L - \sum_{b \in {\{0,1\}}^K}  -\epsilon \sum_{c=1}^C |R^c(H^+_b) | \\ \label{eq:gradientAccum}
			 & =  L -  \epsilon \sum_{c=1}^C \sum_{b \in {\{0,1\}}^K}  \boldsymbol{|} \sum_{i=1}^N g_i^c \sum_{p \in A^+} \delta(b_{i,p}^+ = b)  \boldsymbol{|} \\
			 & = L - \epsilon \sum_{c=1}^C \sum_{b \in {\{0,1\}}^K}  \boldsymbol{|} \sum_{\{ i,p: b_{i,p}^+=b\}} g_i^c  \boldsymbol{|} ~\accentset{\Delta}{=} ~L -\epsilon R(B)
\end{aligned}
\end{equation}

where we used Equation~\ref{eq:Table} for $H^+_b$ and denoted $b_{i,p}^+=B(I_{i,N(p)})$. The resulting training algorithm, summarized in algorithm box~\ref{alg:training}, iterates between global classifier optimization and greedy optimization of the next convolutional table by maximizing $R(B; \Theta)$.

\subsection{Structural variants}
\label{sub:StructDetails}

The word calculator concept described in Section~\ref{sub:Struct} is very general. Here we describe the bit functions and word calculator types we have explored.

\textbf{Bit functions and input preparation:} Word calculators compute an index descriptor of a patch $P\in \mathbb{R}^{l\times l \times D_e}$ by applying $K$ bit functions, each producing a single bit. Each such function is composed of a simple comparison operation, with a few parameters stating its exact operation.  Specifically we use the following bit function forms:

\begin{itemize}
\item {\bf One pixel:} $F(P) = \sigma(P(x,y,d)-t)$
\item {\bf Two pixels:} $F(P) = \sigma(P(x_1,y_1,d) - P(x_2,y_2,d)-t)$
\item {\bf Get Bit $l$:} $F(P) = ( P(0,0,d) << l ) \& 1$
\item {\bf Integral channel bit:} $F(p) = \sigma(P(x_1,y_1,d) - P(x_1,y_2,d) - P(x_2,y_1,d) + P(x_2,y_2,d)-t)$
\end{itemize}

where $\sigma$ is the Heaviside step function. The first two bit function types can be applied to any input channel $d\in \{1,..,D_e\}$, while the latter two are meaningful only for specific channels. The channels we consider are as follows:

\begin{itemize}
\item {\bf Original image channels:} Gray scale and color channels, or depth and IR (multiplied by a depth-based mask) for depth images.
\item {\bf Gradient-based channels:} Two kinds of gradient maps are computed from the original channels following~\cite{Dollar2009}. A normalized gradient channel includes the norm of the gradient for each pixel location. In oriented gradient channels the gradient energy of a pixel is softly quantized into $N^O$ orientation maps.
\item {\bf Integral channels:} Integral images~\cite{Viola01} of channels from the previous two forms, again following~\cite{Dollar2009}. Applying integral channel bits to these channels allows fast calculation of channel area sums.
\item {\bf Spatial channels:} Two channels stating the horizontal and vertical location of a pixel in the image. These channels state the quantized location, using $N^H=\lfloor \log_2 S_x \rfloor$ and $ N^V=\lfloor \log_2 S_y \rfloor$ bits respectively.
\end{itemize}

After preparation, the channels are optionally smoothed by a convolution with a triangle filter. Spatial channels enable the incorporation of patches' position in the word computed. They are used with a 'Get Bit $l$' bit function type, with $l$ referring to the higher bits. This effectively puts a spatial grid over the image, thus turning the global summation pooling into local summation using a pyramid-like structure~\cite{Lazebnik06}. For example using two bit functions, checking for the $N^H$-th horizontal bit and the $N^V$-th vertical bit, effectively puts a $2\times2$ grid over the image where words are summed independently and get different weights for each quarter. Similarly using $4$ spatial bits one gets a $4\times4$ pyramid, etc. We found that enforcing a different number of spatial bits in each convolutional table improves feature diversity and consequently the accuracy.

\textbf{Word calculator structure:} The main design decision in this respect is the choice between ferns and trees.
Ferns include only $K$ bit functions, so the number of parameters is relatively small and over-fitting during local optimization is less likely. Trees are a much larger hypothesis family
with up to $2^K-1$ bit functions in a full tree.
Thus they are likely to
enable higher accuracy, but also be more prone to overfit. We explored this trade-off using a 'long tree' structure enabling a gradual interplay between the fern and full tree extremes.

In a long tree the $K$ bits to compute are divided into $N^S$ stages, with $K_s$ bits computed at stage $s = 1,..,N^S$, so $\sum_{s=1}^{N^S} K_s =K$. A tree
of depth $N^S$ is built, where a node in stage $s$ contains $K_s$ bit functions computing a $K_s$-bits word.
A node in stage $s=1,..,N^S-1$ has $q_s$ children, and it has a child-directing table of size $2^{K_s}$, with entries containing child indices in $1,..,q_s$.
Computation starts at stage $1$ at a root node, and after computation of the $K_s$ bits in a node the produced word is used as an index to the child-directing table, whose output is the index of the child node to descend to. The tree structure is determined by the vectors $(K_1,..,K_{N^S})$ and $(q_1,..,q_{N^S-1})$ of stage size and stage split factors respectively.

When speed is considered, the most important point is that ferns can be efficiently implemented using vector operations (like SSE), constructing the word in several locations at the same time. The efficiency arises because computing the same bit function for several contiguous patches involves access to contiguous pixels, which can be done without expensive gather operations. Conversely, for trees different bit functions are applied at contiguous patches
so the accessed pixels are not contiguous in memory. As will be seen in Section~\ref{sec:EmpiricalResults}, trees can be more accurate, but ferns provide considerably better accuracy-speed trade-off.

\begin{figure*}
\begin{tabular}{| l |c |c| c| c| c |c| c| c| c | }
\hline
Dataset    &    DFE  & CTE base & \textbackslash Opt TH     & \textbackslash Ftr Norm      & \textbackslash WC opt & \textbackslash Chnls & \textbackslash Smooth & \textbackslash Spatial & \textbackslash Sp. Enforce \\
\hline
MNIST      &  0.77   & 0.45     & 0.48       & 0.58          & 0.48          & 0.7      & 0.66   &  0.51   & 0.48                \\
CIFAR-10   &  31.3   & 20.3     & 21.3       & 21.9          & 21.8          & 22.0     & 22.2   &  21.5   & 21.0               \\
SVHN		   &  11.9   & 6.5		  & 7.1			   & 7.1           & 10.5				   & 11.6		  & 7.2	   & 13.2	  & 7.6				          \\
3-HANDPOSE &   3.2   & 2.3      & 2.1		  	 & 2.5           & 3.5					 & 4.4		  & 2.7 	 & 2.2	  & 2.2  	        \\
\hline
\end{tabular}
 \caption{\footnotesize \textbf{Comparison and ablation:} Columns one and two present errors of a DFE~\cite{DFE} and a baseline CTE. The next columns show results obtained when the CTE baseline is {\it \bf ablated} in a single aspect. \textbf{ Opt TH:} bit function thresholds are randomized instead of optimally chosen. \textbf{ Ftr Norm:} histogram features are not multiplicatively normalized before training. \textbf{ WC opt} bit functions in the word calculator are randomly chosen, and not optimized. \textbf{ Chnls:} appearance channels beyond original ones are not used. \textbf{ Smooth:} input channels are not smoothed during preprocessing. \textbf{ Spatial:} Spatial bits are not used. \textbf{ Sp. Enforce:} spatial bits are not enforced (When enforced, a random number in $\{1,..,5\}$ of spatial bits is enforced for each table).
}
\label{table:Ablations}
\end{figure*}

\subsection{Training variants}
\label{sub:TrainDetails}

As stated in algorithm~\ref{alg:training}, training iterates between gradient based word calculator optimization and global optimization of table weights. We now describe the methods we explored for these two components.

\textbf{Word calculator optimization:} We consider several mechanisms for the optimization of $R(B;\Theta)$, including forward bit function selection, optimal threshold finding, and iterative bit function replacement/refinement.

In forward selection, we optimize $R(B)$ by adding one bit after the other. For fern growing there are $K$ such stages. At stage $l=1,..K$, ${\{F^j\}}_{j=1}^{N_c}$ candidate bit functions are generated, with their type and parameters drawn from a prior distribution. For each $j$, we augment the current word calculator $B$ to $B^+=[B,F^j]$ and choose the one with the highest score. However, we found that simple greedy computation of $R(B^+)$ at each stage is not the best way to optimize $R(B)$, and an auxiliary score which additively normalizes the newly-introduced features does a better job. Denote the patch features of a word calculator $B$ by $\delta^b(P)=\delta(B(P)=b)$, by $\delta^b_{i,p}$ the value of $\delta^b$ for pixel $p$ in image $i$ and by $R^c(f(P)) \accentset{\Delta}{=} \sum_{i,p} g^c_i f_{i,p}$ the score $R$ induced by a patch feature $f$. The addition of a new bit effectively replaces the feature $\delta^b$ for $b\in\{0,2^l-1\}$ with $2$ new features $\delta^{(b,0)}$ and $\delta^{(b,1)}$. If the gradients in cell $b$ are not balanced; i.e., $\sum_{i,p} g^c_i \delta^b_{i,p}= C_0 \ne 0$, as is often the case, a feature $\delta^{(b,0)}$ may get a good $R^c(\delta^{(b,0)})$ score of $C_0$ even if the new bit function is constant, or otherwise uninformative. To handle this, we score a normalized version of the new features, with an average value of 0, which more effectively measures the added information in the new features. The following lemma shows that this is a valid, as well as computationally effective strategy:

\begin{lemma}\label{lemma:L1Delta} Let $\bar{\delta}^{(b,u)} = \delta^{(b,u)} - \frac{\#\delta^{(b,u)}}{\#\delta^b}\delta^{b}$ for $u=0,1$ and $\#\delta^a = \sum_{i,p} \delta^a_{i,p}$. The following properties hold
\begin{enumerate}
\item \label{lem:equivalence} Using $\bar{\delta}^{(b,1)}$, $\delta^b$ in a classifier is equivalent to using $\delta^{(b,0)}$,$\delta^{(b,1)}$; i.e, for any weight choice $w_0,w_1$ there are $w_b,w_\Delta$ such that $w_0 \delta^{(b,0)} + w_1 \delta^{(b,1)} = w_b\delta^b + w_\Delta \bar{\delta}^{(b,1)}$
\item \label{lem:half}  $R^c(\bar{\delta}^{(b,0)}) = R^c(\bar{\delta}^{(b,1)})$

\item \label{lem:dual}  $R^c(\bar{\delta}^{(b,0)})= \sum_{i,p}(g^c_i - E[g^c_i|b])\delta^{(b,0)}$ with $E[g^c_i|b]  \accentset{\Delta}{=} \frac{\sum_{i,p}g^c_i\delta^b}{\#\delta^b}$
\end{enumerate}
\end{lemma}

The proofs are rather simple and appear in appendix~\ref{App:LemmaProof}.
Property~\ref{lem:equivalence} shows that we may score $\delta^b,\bar{\delta}^{(b,1)}$ features instead of $\delta^{(b,u)}$ features. Since only $\bar{\delta}^{(b,1)}$ is affected by the new candidate bit, we can score only those terms when selecting among candidates.
Property~\ref{lem:dual} shows that we can normalize the gradient instead of the feature candidates, which is cheaper (as there are $N_c$ candidates but only a single gradient vector). In summary, we optimize the next bit selection by maximizing
\begin{equation}
\label{eq:L1Delta}
R_\Delta ([B,F^j]) \accentset{\Delta}{=} \sum_{c=1}^C \sum_{ b\in \{0,1\}^{l}}  \boldsymbol{|} \sum_{\{i,p: b_{i,p}=(b,1)\}} (g_i^c - E[g^c_i|b]) \boldsymbol{|}
\end{equation}

over the choice of $F^j$. The calculation requires a single histogram aggregation sweep over all patches $(i,p)$.

Most of the bit functions obtain their bit by comparing an underlying patch measurement to a threshold $t$. For such functions, the optimal threshold parameter can be found with a small additional computational cost. This is done by sorting the underlying values of $F^j$ and computing the sum over $i,p$ in Equation~\ref{eq:L1Delta} by running with the sorted order. This way, a running statistic of the $R_{\Delta}$ score can be maintained, computing the score for all possible thresholds and keeping the best.

For a long tree a similar algorithm is employed, with ferns internal to nodes optimized as full ferns, but tree splits requiring special treatment. Assume we are splitting a node in stage $s$, so the current word calculator has already computes a $L_s=\sum_{i=1}^s K_i$-bit word, among which $K_s$ were computed in the current node. We now choose the first bit functions of all the children, as well as the redirection table, to optimize $R$. Since different prefixes of the current calculator $B$ are augmented by different bit functions we need to decompose the $R$ score. Denote by $a$ the index set $\{L_s-K_s+1,..,L_s\}$ of bits computed by the current node, and by b(a) the limitation of a binary word $b$ to indices $a$.
For a $K_s$-bit word $z\in{\{0,1\}}^{K_s}$, we define the component of $R$ contributed by words with $b(a)=z$ by

\begin{equation}
R_{b(a)=z}(B) \accentset{\Delta}{=} \sum_{c=1}^C \sum_{ \substack{b\in \{0,1\}^{L_s} \\ b(a)=z} } \boldsymbol{|} \sum_{\{i,p: b_{i,p}=b\}} g_i^c \boldsymbol{|}
\end{equation}

For the tree split we draw a large set $\mathbb{F}$ of candidate bits, and choose the first bits of the $q_s$ children by optimizing

\begin{equation}
\max_{ \substack{G \subset \mathbb{F} \\ |G|=q_s}}  \sum_{z \in {\{0,1\}}^{K_s}} \max_{F_{z}\in G} R_{b(a)=z}([B,F_{z}])
\end{equation}

with $G$ the set of chosen bits for the children and entry $z$ in the redirection table set to the index of the child containing $F_z$. For this optimization we compute the score matrix $S \in \mathbb{M}_{2^{K_s} \times |\mathbb{F}|}$ with $S(i,j)=R_{b(a)=i}([B,F_j])$. Given a choice of $G$, amounting to a choice of column subset in $S$, the optimization over $F_z$ is trivial and the score is easy to compute. We optimize over $G$ by exhaustively trying all choices of $G$ for $|G|=2$, and greedily adding columns to $G$ until it contains $q_s$ members.

In addition to forward bit selection, we implemented iterative bit replacement and refinement stages. The rationale for this is the observation that while the last bit functions in a fern are chosen to complement the previous ones, the bits chosen at the beginning are not optimized to be complementary and may be suboptimal in a long word calculator. The bit replacement algorithm operates after forward bit selection. It runs over all the bit functions several times and attempts to replace each function with several randomly drawn candidates. A replacement step is accepted if it improves the $R_{\Delta}(B)$ score. In a similar manner, a bit refinement algorithm attempts to replace a bit function by small perturbations of its parameters, thus effectively implementing a local search. For trees, bit replacement/refinement is done only for bits inside a node, and once a split is made the node parameters are fixed.

\begin{figure*}
\begin{tabular}{cc}
\begin{tabular}{| l | c | c | c |}
\hline
Dataset and Error & CTE ($\mu S$) & CNN ($\mu S$) & speedup   \\
\hline
MNIST 0.01		  		 &  4.8  			 &  63.9  		 & 13.1$\times$      \\
CIFAR-10 0.25   		 & 168.8 			 &  882  		   & 5.2$\times$     \\
SVHN  0.15 		  		 & 18.4  			 &  88  		   & 4.7$\times$     \\
3-HANDPOSE 0.035		 & 6.3				 &  1250       & 199.3$\times$		 \\
\hline
\end{tabular}  &
\begin{tabular}{| l | c | c | c |}
\hline
Dataset &  Tree form   								 & Error(\%)  & speed ($\mu S$) \\
\hline
MNIST    &  Fern				 					     & 0.45    		& 106           \\
				 &  Depth 3, 4-way splits		   & 0.43  	  	& 398           \\
				 &  Depth 6, 2-way splits 	   & 0.39  	  	& 498           \\
\hline
CIFAR    &  Fern  										 & 21.0  			& 800         \\
	~~-10	 &  Depth 3, 4-way splits			 & 19.7   		& 3111        \\
				 &	Depth 4, 3-way splits      & 19.3       &	3544			  \\
\hline
\end{tabular}
\end{tabular}
 \caption{\footnotesize \textbf{Left:CTE-CNN speed differences:} Some examples of required accuracy points where a CTE can meet the accuracy while being considerably faster than a CNN.
\textbf{Right:Ferns/Trees trade-off:} Accuracy and speed for several fern/tree configurations on MNIST and CIFAR-10. The fern classifiers are the baseline classifiers whose results are reported in Figure~\ref{table:Ablations}. For trees the bits are split evenly between the $N^S$ stages (for MNIST, which uses $11$-bit tables, the last stage gets one less bit), and all the split parameters are equal.
}
\label{fig:MoreTables}
\end{figure*}

\textbf{Global optimization:}
We considered two global loss functions in our classification experiments: an SVM-based loss, and a softmaxloss as typically used in neural networks optimization. In the SVM loss, we take the sum of $C$ SVM programs, each minimizing a one-versus-all error. Let $y_{i,c}=2*\delta(y_{i},c)-1$ be binary class labels. The loss is

\begin{equation}
L_{SVM} = \frac{1}{2}||W||^2 + \Lambda \sum_{c=1}^C \sum_{i=1}^N max(1-y_{i,c}s^c_i, 0)
\end{equation}

The loss aims for class separation in $C$ independent classifiers. Its advantage lies in the availability of fast and scalable methods for solving large and sparse SVM programs~\cite{Shalev07,Hsieh2008}. The loss gradients are $g^c_i = -y_{i,c}$ if example $i$ is a support vector, and $0$ otherwise. In~\cite{BarHillel10} a first order approximation for $\min_W L_{SVM}$ is derived for new feature addition, in which the example gradients are $-\alpha_i y_{i,c}$ with $\alpha_i$ the dual SVM variables at the optimum. The two expressions are similar and we did not find noticeable difference between them empirically.
The softmax loss is
\begin{equation}
L_{LR} = \frac{1}{2}||W||^2 - \Lambda \sum_{i=1}^N \log \frac{exp(s_i^{y_i})}{\sum_c exp(s^c_i)}
\end{equation}
This loss provides a direct minimization of the $M$-class error.
The gradients are $g^c_i=exp(s_i^{y_i})/ \sum_c exp(s^c_i) - \delta(y_{i},c)$.
Conveniently, it can be extended to a distillation loss~\cite{Hinton2015}, which enables guidance of the classifier using an internal representation of a well-trained CNN classifier.

Features in a word histogram have significant variance,
as some words appear in large quantities in a single image. Without normalization such words may be arbitrarily preferred due to their lower regularization cost- they can be used with lower weights. Denote the column of a feature across all examples by $Col^m_b= (H^m_{b,1},..H^m_{b,N})$. We found that normalizing each features column by the expected count of active examples $L_1(Col^m_b)/ L_0(Col^m_b)$  improved accuracy and convergence speed in many cases.