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

Wavelet transform methods are important tools in image processing due to their capabilities for
multiresolution analysis and image decomposition~\cite{Mallat99,Burrus98}.
Wavelet-based image processing
find application
in computer graphics~\cite{Schroder96}, including radiosity, global
illumination~\cite{Stollnitz96}, and real volume data~\cite{Muraki93},
and volume rendering~\cite{Roerdink99};
being
routinely
considered as an approach
for texture image decomposition,
image coding,
subband coding,
fast image segmentation~\cite{Aujol2006,Antonini92,Vetterli90,Lim90,Kim2003},
and
3D signal processing~\cite{Taubman94}.
In particular,
wavelet transform coding~\cite{Ohm2004,Bottreau2001,Usevitch2001,DeVore95} has emerged as a
practical and fully-established tool~\cite{Mallat99,Aujol2006}
which
benefits image compressing methods,
such as
the {JPEG} 2000 standard~\cite{Skodras2001,Usevitch2001,Saha2000,Hilton94},
and
multimedia schemes on Internet,
such as
data streaming over {IP}~\cite{Radha2001,Baganne2003}
and
image querying~\cite{Jacobs95}.
Scalable coding for image, audio, and video
largely
adopts wavelet-based methods
as demonstrated in the MPEG-4 codec~\cite{Radha2001,Ohm2004,Bottreau2001,Shaar2003,Sodagar2002}.

Multiwavelets
have been explored
for the assessment of order/disorder
in reconstructed biomedical images~\cite{Zemni2019}.
Further connections between information-theoretical entropy
and
wavelet analysis
were explored in the context of geoscience~\cite{nicolis20152d}.
Recent advances in the field
of image decomposition
include
the proposition
of special filters for
spherical harmonics modeling~\cite{Jallouli2019b}
capable
of
multi-level decomposition
suitable for 3D images
with the introduction
of the concept of spherical harmonics entropy~\cite{Jallouli2019a}.

In the same vein of exploring connections between different research fields,
Ashmead reported~\cite{Ash2012}
a link between
quantum mechanics and wavelets.
Despite the potentially deep mathematical meaning of such link,
it has not been significantly explored in
the
context of wavelet decomposition and image analysis.
Indeed,
the wave nature of light can be deduced from the phenomenon of interference, the photoelectric effect, however, it seems to suggest a corpuscular nature of light.
Theoretical physicists struggled to include observations like the photoelectric effect and the wave-particle duality into their formulations~\cite{Sully2012}.
Erwin Schr\"odinger
employed advanced mechanics to address
such phenomena
and developed an equation that relates the space-time in quantum mechanics. %
Because wavelets are localized in both time and frequency they offer significant advantages for the
analysis of problems in quantum mechanics.

In this paper,
we aim at proposing
an alternative wavelet decomposition scheme.
For such,
we explore
the above discussed link between
wavelets and quantum mechanics
and
shed some light on some of these relations.
Rather than seeking at wavelet features on particles or waves,
we adapted some concepts of quantum mechanics to a novel wavelet decomposition for still images.

The paper is organized as follows.
Section~\ref{sec:motivation}
details the main ideas behind the proposed wavelet system
which
is based on
an interpretation
from particle physics and quantum states~\cite{Eisberg2007,Beiser94}.
We describe the construction of the proposed orbital wavelets
for the image decomposition.
It is shown that the introduced derivation
is naturally
suitable for generating two-dimensional wavelets.
In Section~\ref{sec:implementing}
we submit standard imagery to the proposed wavelet decomposition
scheme
and compare
the results with
standard wavelet decomposition..
Section~\ref{sec:concluding}
concludes the paper.

\section{Orbital Wavelet Decomposition}
\label{sec:motivation}

\subsection{Particle Systems}

The wave functions describing electronic orbitals
can be combined generating
atomic orbitals.
Let us consider
a
two-particle non-interaction systems
with
particles
located at position vectors
$\mathbf{r}_1$
and
$\mathbf{r}_2$,
respectively,
where
each particle
is equipped
with wave functions
$\psi_\alpha(\mathbf{r}_1)$
and
$\psi_\beta(\mathbf{r}_2)$
at states
$\alpha$ and $\beta$,
respectively.
The wave function
that characterizes the orbital interaction of the two particles
is furnished by a combination
$\psi_\alpha(\mathbf{r}_1)$
and
$\psi_\beta(\mathbf{r}_2)$
in two different configurations:
symmetric ($S$)
and
anti-symmetric ($A$)
combinations~\cite{Duck98}.
Such combinations are given by:
\begin{align}
\psi_{S}(\mathbf{r}_1,\mathbf{r}_2)
&
=
\frac{1}{\sqrt{2}}
\Big[
\psi_\alpha(\mathbf{r}_1)
\cdot
\psi_\beta(\mathbf{r}_2)
+
\psi_\alpha(\mathbf{r}_2)
\cdot
\psi_\beta(\mathbf{r}_1)
\Big]
,
\label{eq:4}
\\
\psi_{A}(\mathbf{r}_1,\mathbf{r}_2)
&
=
\frac{1}{\sqrt{2}}
\Big[
\psi_\alpha(\mathbf{r}_1)
\psi_\beta(\mathbf{r}_2)
-
\psi_\alpha(\mathbf{r}_2)
\psi_\beta(\mathbf{r}_1)
\Big]
\label{eq:5}
,
\end{align}
respectively.
The anti-symmetric case
can be conveniently
written
in matrix format
as
$
\psi_{A}(\mathbf{r}_1,\mathbf{r}_2)
=
\frac{1}{\sqrt2}
\operatorname{det}
\left[
\begin{smallmatrix}
\psi_{\alpha}(\mathbf{r}_1) & \psi_{\alpha}(\mathbf{r}_2)\\
\psi_{\beta}(\mathbf{r}_1) & \psi_{\beta}(\mathbf{r}_2)
\end{smallmatrix}
\right]
.
$

A
comparable concept
in the scope of wavelets, also characterized by wave functions,
would
be the combination of different spatial wavelets~\cite{DeO}.

\subsection{Orbital Wavelets}

\subsubsection{Symmetric Case: Standard Wavelets}

Usual wavelet image analysis combines one-dimensional (1D) wavelets to generate
a two-dimensional (2D) wavelet~\cite{Antonini92,Vetterli90}.
This can be done by considering a scaling function $\varphi$ and
a wavelet function $\psi$,
a version for each dimension, abscissa and ordinate~\cite{Burrus98}.
Thus,
we have a 2D scale function
$\varphi_\text{LL}(x,y) = \varphi(x)\cdot\varphi(y)$
and three
2D wavelet functions:
\begin{align}
\label{equation-standard-2d-wavelets}
\begin{split}
\psi_\text{LH}(x,y) &= \varphi(x)\cdot\psi(y)
,
\\
\psi_\text{HL}(x,y) &= \psi(x)\cdot\varphi(y)
,
\\
\psi_\text{HH}(x,y) &= \psi(x)\cdot\psi(y)
.
\end{split}
\end{align}

The wavelets
$\psi_\text{LH}(x,y)$
and
$\psi_\text{HL}(x,y)$
naturally exhibit
reflection symmetry
with respect
to the plane
$x=y$,
i.e.
$
\psi_\text{LH}(x,y)
=
\psi_\text{LH}(y,x)
$
and
$
\psi_\text{HL}(x,y)
=
\psi_\text{HL}(y,x)
$.
Such reflection symmetry
is analogous
to
the symmetric wave function
described in~\eqref{eq:4}.

\subsubsection{Anti-symmetric Case: Orbital Wavelets}

Following this analogy,
considering a single orthogonal wavelet system,
we are compelled
to pursue the definition of wavelets
that
could be regarded as the counterparts
of the
anti-symmetric wave function described in~\eqref{eq:5}.
Thus
the combination of $\varphi(\cdot)$ and $\psi(\cdot)$
should
be arranged to provide anti-symmetry,
i.e.
the
$\hat{\psi}_\text{LH}(x,y)$
and
$\hat{\psi}_\text{HL}(x,y)$
wavelets
should be such that:
\begin{align}
\label{eq:2}
\hat{\psi}_\text{LH}(x,y)
&
=
-
\hat{\psi}_\text{LH}(x,y)
\\
\hat{\psi}_\text{HL}(x,y)
&
=
-
\hat{\psi}_\text{HL}(x,y)
\end{align}
A solution to the above requirement
is
to define an orbital-inspired combination
(cf.~\eqref{eq:5})
of $\varphi$ and $\psi$
according to the following:
\begin{align}
\label{equation-orbital}
\hat{\psi}_\text{LH}(x,y)
\triangleq
\frac{1}{\sqrt{2}}
\Big[
\varphi(x)\cdot\psi(y)-\psi(x)\cdot\varphi(y)
\Big]
.
\end{align}
In order
to illustrate the effect
of this definition
we consider
the case of the Meyer orthogonal wavelet~\cite{Meyer90}.
Figure~\ref{fig:meyer} shows
the surface plots for
the standard  of the discussed functions.

\begin{figure}
\centering

\subfigure[$\psi_\text{LH}(x,y)$]
{\includegraphics[width=0.48\columnwidth]{phipsiB_W.eps}}
\subfigure[$\hat{\psi}_\text{LH}(x,y)$]
{\includegraphics[width=0.48\columnwidth]{phi_psiB_W.eps}}

\vspace*{8pt}

\caption{2D-Meyer decomposing functions of the secondary diagonal:
(a)~$\psi_\text{LH}(x,y)$;
(b)~$\hat{\psi}_\text{LH}(x,y)$.}

\label{fig:meyer}
\end{figure}

Analogously,
we define
the
wavelets
related to the LL and HH
decompositions
according
to the following
expressions:
\begin{align}
\label{equation-orbital-wavelet-LL}
\hat{\varphi}_\text{LL}(x,y)
&
\triangleq
\frac{1}{\sqrt{2}}
\Big[
\varphi^{\ast}(x)
\cdot
\varphi(y)
+
\varphi^{\ast}(y)
\cdot
\varphi(x)
\Big]
,
\\
\label{equation-orbital-wavelet-HH}
\hat{\psi}_\text{HH}(x,y)
&
\triangleq
\frac{1}{\sqrt{2}}
\Big[
\psi^{\ast}(x)
\cdot
\psi(y)
-
\psi^{\ast}(y)
\cdot
\psi(x)
\Big]
.
\end{align}
The above definitions allows
the analysis of images using continuous
complex wavelets~\cite{Kingsbury99}.
For real-valued
wavelets,
the above expressions
collapse to the usual forms
$\varphi(x) \cdot \varphi(y)$
or
$\psi(x)\cdot\psi(y)$
present in standard wavelet analysis.
Thus in the real case,
although
the proposed wavelet
$\hat{\varphi}_\text{LL}(x,y)$
in \eqref{equation-orbital-wavelet-LL}
is well-defined,
the
wavelet
$\hat{\psi}_\text{HH}(x,y)$
would collapse to zero.

An approach to address such degeneracy
is to consider
daughter wavelets at different
scales.
Therefore
let us consider the 1{D} orthogonal~\cite{Lawton91,
Maab96} wavelet mother $\psi(x)$ equipped
equipped
with her daughter wavelets
$\{\psi_{a,b}(x)\}_{a\neq0,b\in R}$.
The formalism shown in~\eqref{equation-orbital-wavelet-HH}
can be extended
by considering
the inclusion
of
two wavelets $\psi_{a_1,b}(\cdot)$ and $\psi_{a_2,b}(\cdot)$
resulting in the following definition.

\begin{definition}
\label{definition-2d-orbital}
The function \mbox{2D}-orbital at the scales $\{a_1,a_2\}$
is defined by:
\begin{align}
\hat{\psi}_\text{HH}(x,y)\triangleq\frac{1}{\sqrt{2}}\operatorname{det}\begin{bmatrix}\psi_{a_1,b}^{*}(x) & \psi_{a_1,b}^{*}(y)\\
\psi_{a_2,b}(x) & \psi_{a_2,b}(y)
\end{bmatrix},
\end{align}

which can be rewritten as:
\begin{align}
\hat{\psi}_\text{HH}(x,y)
=
&
\frac{1}{\sqrt{2|a_1||a_2|}}
\psi^{\ast}\left(\frac{x-b}{a_1}\right)
\cdot
\psi\left(\frac{y-b}{a_2}\right)
\\
&
-
\frac{1}{\sqrt{2|a_1||a_2|}}
\psi\left(\frac{x-b}{a_2}\right)
\cdot
\psi^{\ast}\left(\frac{y-b}{a_1}\right)
.
\end{align}
\end{definition}

The condition $a_1\neq a_2$ eliminates the degeneration
$\hat{\psi}_\text{HH}(x,y)=0$.
This is to some extent in connection to the
Pauli Exclusion Principle~\cite{Duck98},
which
states that with a single-valued many-particle wave function is equivalent to requiring the wave function to be antisymmetric.
An antisymmetric two-particle state is represented as a sum of states in which one particle is in state $\alpha$ and the other in state $\beta$. Besides, the relationship $\hat{\psi}_\text{HH}(y,x)=-\hat{\psi}_\text{HH}(x,y)$ ensures the desired asymmetry.
Here, we use the same wavelet-mother,
but at different scales.
The 2D decomposition stated in
Definition~\ref{definition-2d-orbital} results in
a strict 2D wavelet.

\subsection{Mathematical Properties}

\begin{proposition}
The previously defined 2D-orbital function has oscillatory behavior
satisfying the following properties:
\begin{enumerate}[(i)]
\item
$\intop_{-\infty}^{\infty}\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!x=0$,
\item
$\intop_{-\infty}^{\infty}\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!y=0$,
\item
$\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!x\operatorname{d}\!y=0$.
\end{enumerate}

\end{proposition}

\begin{proof}It follows that

\begin{align}
\intop_{-\infty}^{\infty}\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!x=\frac{1}{\sqrt{2}}\psi_{a_2,b}(y)\cdot\overline{\psi_{a_1,b}^{\ast}(x)}
-\frac{1}{\sqrt{2}}\overline{\psi_{a_2,b}(x)}\cdot\psi_{a_1,b}^{\ast}(y),
\end{align}
where
\begin{align}
\overline{\psi_{a,b}(x)}
\triangleq
\intop_{-\infty}^{\infty}
\psi_{a,b}(x)\operatorname{d}\!x
.
\end{align}
Therefore,
the property~(i) derives from the fact that $\psi_{a,b}(x)$,
$a=\{a_1,a_2\}$ be individual wavelets.
The demonstration for property~(ii)
is similar, considering that
\begin{align}
\intop_{-\infty}^{\infty}\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!y
=
\frac{1}{\sqrt{2}}
\overline{\psi_{a_2,b}(y)}
\cdot
\psi_{a_1,b}^{\ast}(x)
-
\frac{1}{\sqrt{2}}
\psi_{a_2,b}(x)
\cdot
\overline{\psi_{a_1,b}^{\ast}(y)}.
\end{align}
The condition
$\intop_{-\infty}^{\infty}
\intop_{-\infty}^{\infty}
\hat{\psi}_\text{HH}(x,y)\operatorname{d}\!x\operatorname{d}\!y=0$
follows from Fubini's theorem~\cite{DeFig2000},
regardless of the order of integration.
\end{proof}

Hereafter
we assume
an orthogonal wavelet system
with unitary energy.
In other words,
the following conditions hold true:
\begin{enumerate}[(i)]

\item
the inner product
$\langle\psi_{a_1,b},\psi_{a_2,b}\rangle=0$,
i.e., the following integrals cancel out $\forall a_1\neq a_2$:
\begin{align}
\intop_{-\infty}^{\infty}
\psi_{a_1,b}(x)\cdot\psi_{a_2,b}^{*}(x)
\operatorname{d}\!x
=
\intop_{-\infty}^{\infty}
\psi_{a_1,b}^{*}(x)\cdot\psi_{a_2,b}(x)
\operatorname{d}\!x=0
,
\end{align}
and
\item
$\int_{-\infty}^\infty | \psi_{a_1,b}(x) |^2 \operatorname{d}\!x = 1$.

\end{enumerate}
It is also noteworthy that
\begin{align}
\langle\psi_{a_1,b},\psi_{a_2,b}\rangle^{*}
=
\langle\psi_{a_2,b},\psi_{a_1,b}\rangle
.
\end{align}

\begin{proposition}
The 2D-orbital functions have normalized energy.
\end{proposition}
\begin{proof}
Expanding the expression
$|\hat{\psi}_\text{HH}(x,y)|^2=\hat{\psi}_\text{HH}(x,y)\cdot\hat{\psi}_\text{HH}^{\ast}(x,y)$
yields
cross-product terms
in the following form
\begin{align}
\operatorname{cross}(x,y)
\triangleq
-\psi_{a_1,b}(x)
\cdot
\psi_{a_2,b}^{\ast}(y)
\cdot
\psi_{a_1,b}(y)
\cdot
\psi_{a_2,b}^{^{\ast}}(y)
\end{align}
and its complex conjugate
$\operatorname{cross}^\ast(x,y)$.
Performing the integration with respect
to $x$ and $y$
yields:
\begin{align}
\intop_{-\infty}^{\infty}
\operatorname{cross}(x,y)
\operatorname{d}\!x
=
-2
\cdot
\psi_{a_2,b}^{\ast}(y)\psi_{a_1,b}(y)
\cdot
\intop_{-\infty}^{\infty}
\psi_{a_1,b}^{\ast}(x)\psi_{a_2,b}(x)
\operatorname{d}\!x
,
\end{align}
and
\begin{align}
\intop_{-\infty}^{\infty}
\operatorname{cross}(x,y)
\operatorname{d}\!y
=
-2
\cdot
\psi_{a_2,b}^{\ast}(x)\psi_{a_1,b}(x)
\cdot
\intop_{-\infty}^{\infty}
\psi_{a_1,b}^{\ast}(y)\psi_{a_2,b}(y)\operatorname{d}\!x
.
\end{align}
Invoking the orthogonality condition
we obtain
that
all cross terms are void.
The remaining possibly nonnull terms
are:
\begin{align}
|\hat{\psi}_\text{HH}(x,y)|^2
=
\frac{1}{2}|\psi_{a_1,b}(x)|^2
\cdot
|\psi_{a_2,b}(y)|^2
\cdot
|\psi_{a_1,b}(y)|^2
\cdot
|\psi_{a_2,b}(x)|^2
\end{align}
and therefore,
because of the normalized energy condition,
we have:
\begin{align}
\intop_{-\infty}^{\infty}\intop_{-\infty}^{\infty}|\hat{\psi}_\text{HH}(x,y)|^2\operatorname{d}\!x\operatorname{d}\!y=1
.
\end{align}
\end{proof}

It is possible (more
easily) to combine orthogonal 1D-wavelets and use them to build a
new 2D-wavelet.

\begin{proposition}
The 2D-orbital function is a 2D wavelet.\end{proposition}

\begin{proof}
Let $\Psi(\omega)$
and
$\Psi_{a,b}(\omega)$ be the Fourier transforms
of
the wavelet~$\psi(t)$
and
the
daughter wavelets~$\psi_{a,b}(t)$,
respectively.
If
the admissibility condition holds~\cite{Boggess09,Burrus98},
\begin{align}
\intop_{-\infty}^{\infty}
\frac{|\Psi(\omega)|^2}{|\omega|}
\operatorname{d}\!\omega < \infty
\end{align}
then
$\intop_{-\infty}^{\infty}
\frac{|\Psi_{a,b}(\omega)|^2}{|\omega|}
\operatorname{d}\!\omega
<\infty
$,
since
$\Psi_{a,b}(\omega)=\sqrt{|a|}\,\Psi(a\omega)e^{-j \omega b}$~\cite{Mallat99}.
Let us now evaluate the condition for the 2D case. If the Fourier
transform pair
$\hat{\psi}_\text{HH}(x,y)\leftrightarrow\hat{\Psi}_\text{HH}(u,v)$ does exist,
then the 2D spectrum of $\hat{\psi}_\text{HH}~$ can be computed
in terms of the Fourier spectrum
of $\psi$:
\begin{align}
\hat{\Psi}_\text{HH}(u,v)
=
\frac{\sqrt{|a_1a_2|}}{\sqrt{2}}
\Big[
\Psi(a_1u)\Psi^{\ast}(a_2v)-\Psi(a_2u)\Psi^{\ast}(a_1v)
\Big]
.
\end{align}

From the generalized Parseval-Plancherel energy theorem~\cite{DeFig2000,
Burrus98}, the cross-terms vanish due to the orthogonality.
Thus, we have
\begin{align}
|\hat{\Psi}_\text{HH}(u,v)|^2
=
\frac{|a_1a_2|}{2}
|\Psi(a_1u)|^2
\cdot
|\Psi(a_2v)|^2
+
\frac{|a_1a_2|}{2}
|\Psi(a_2u)|^2
\cdot
|\Psi(a_1v)|^2
.
\end{align}
Then,
from the marginal admission conditions of the 1D daughter wavelets,
we
have that
\begin{align}
\int_{-\infty}^{\infty}
\int_{-\infty}^{\infty}
\frac{|\hat{\Psi}_\text{HH}(u,v)|^2}{|u|\cdot|v|}
\operatorname{d}\!u
\operatorname{d}\!v
<
\infty.
\end{align}
\end{proof}

\section{Orbital Wavelet-based  Image Decomposition}
\label{sec:implementing}

The standard 2D wavelet decomposition
effects coefficient matrices representing
vertical,
horizontal,
and diagonal structures
denoted by
sub-images
$\{
L_1L_1,
L_1H_1,
H_1L_1,
H_1H_1
\}$
and
$\{
L_2L_2,
L_2H_2,
H_2L_2,
H_2H_2
\}$
for
the first and second level decompositions~\cite{Strang1997}.
Figure~\ref{table-standard-decomposition}
illustrate the structures.
Noticed that the terms
$L_{1}L_{1}$,
$L_{2}L_{2}$,
$H_{1}H_{1}$,
and
$H_{2}H_{2}$ on the main diagonal correspond to a part of the standard
wavelet decomposition of the image into two levels.
The proposed wavelet analysis
results in a similar structure
with sub-images
$\{
\widehat{L_1L_1},
\widehat{L_1H_1},
\widehat{H_1L_1},
\widehat{H_1H_1}
\}$
and
$\{
\widehat{L_2L_2},
\widehat{L_2H_2},
\widehat{H_2L_2},
\widehat{H_2H_2}
\}$
in a two-level decomposition.
Despite the similarity,
the
sub-images are computed
according to the orbital wavelets.
The resulting scheme is depicted in Figure~\ref{table-orbital-decomposition}.

\begin{figure}
\centering

\subfigure[First level]{\includegraphics{std-wav-level-1.ps}}
\subfigure[Second level]{\includegraphics{std-wav-level-2.ps}}
\caption{Image decomposition scheme according to the standard \mbox{2D} wavelet decomposition.}
\label{table-standard-decomposition}
\end{figure}

\begin{figure}
\centering

\subfigure[First level]{\includegraphics{orb-wav-level-1.ps}}
\subfigure[Second level]{\includegraphics{orb-wav-level-2.ps}}
\caption{Image decomposition scheme according to the 2D orbital wavelet decomposition.
A two-level wavelet-orbital decomposition.}
\label{table-orbital-decomposition}
\end{figure}

Considered the symlet wavelet of order~4 (\texttt{symlet4})~\cite{Daub92},
we applied
the proposed decomposition
to the standard image \texttt{woman}~\cite{Misiti2002};
the resulting sub-images are shown
in Figure~\ref{fig:woman}.
For a qualitative comparison,
we included the sub-images obtained from the usual wavelet decomposition.
The computation was performed
in the Matlab environement~\cite{Misiti2002}.
It is worth nothing that the subtraction of the images
resulting from
$\psi_\text{LH}(x,y)$
and
$\psi_\text{HL}(x,y)$
results
in the image obtained by the wave function
in~\eqref{equation-orbital}.

\begin{figure}
\centering

\subfigure[standard wavelet decomposition]
{
\begin{tabular}{cc}
\includegraphics[width=0.2\columnwidth]{womanSym4_B_W1.eps}
&
\includegraphics[width=0.2\columnwidth]{womanSym4_B_W3.eps}
\\
\includegraphics[width=0.2\columnwidth]{womanSym4_B_W2.eps}
&
\includegraphics[width=0.2\columnwidth]{womanSym4_B_W4.eps}
\end{tabular}
}
\subfigure[Proposed decomposition]
{
\begin{tabular}{cc}
\includegraphics[width=0.2\columnwidth]{womanoSym4_B_W5.eps}
&
\includegraphics[width=0.2\columnwidth]{womanoSym4_B_W6.eps}
\\
\includegraphics[width=0.2\columnwidth]{womanoSym4_B_W7-negative.eps}
&
\includegraphics[width=0.2\columnwidth]{womanoSym4_B_W8.eps}
\end{tabular}
}

\vspace*{8pt}

\caption{First-level decomposition woman image using \texttt{symlet4}
wavelet as defined in Matlab
according to
(a)~the standard wavelet decomposition
and
(b)~the proposed decomposition.}
\label{fig:woman}
\end{figure}

\section{Concluding Remarks and Future Work}
\label{sec:concluding}

This paper offers an alternative
approach for image
decomposition engendered by asymmetric orthogonal wavelets
whose definition
is inspired by the wave function theory
from particle physics.
Despite the focus being essentially on still image analysis,
the proposed
approach allows a fully scalable multimedia decomposition.
It remains
to be investigated the potential of this methodology
for
image compressing~\cite{DeVore95},
in 3D processing and scalable coding for multimedia
schemes~\cite{Ohm2004}.
Applications in other scenarios such as
wavelet-based watermarking~\cite{JHU2002}
or
steganography~\cite{Carrion2008}
also deserve an investigation.
To the best of our knowledge,
this is the first work to link
the exchanging formalism
from particle wave functions to wavelets analysis.

\section*{Acknowledgments}
The second and third authors
acknowledge
the partial support
from the Brazilian agencies
CAPES and CNPq,
respectively.

{\small
\singlespacing
\bibliographystyle{siam}


\section{Introduction}

A typical machine learning paradigm for discriminative tasks is to learn the distribution of an output, $y$ given a complete set of features, $x \in \mathbb{R}^d$: $p(y \mid x)$. Although this paradigm is successful in a multitude of domains, it is incongruous with the expectations of many real-world intelligent systems in two key ways: first, it assumes that a complete set of features has been observed; second, as a consequence, it also assumes that no additional information (features) of an instance may be obtained at evaluation time. These assumptions often do not hold; human agents routinely reason over instances with incomplete data and decide when and what additional information to obtain. Take for instance a doctor making a diagnosis on a patient. The doctor usually has not observed all the possible measurements (such as blood samples, x-rays, etc.) for the patient. He/she is also not forced to make a diagnosis based on the observed measurements; instead, he/she may dynamically decide to take more measurements to help determine the diagnosis. Of course, the next measurement to make (feature to observe), if any, will depend on the values of the already observed features; thus, the doctor may determine a different set of features to observe from patient to patient (instance to instance) depending on the values of the features that were observed. For example, a low value from a blood test may lead a doctor to ask for a biopsy, whereas a high value may lead to an MRI. Hence, not each patient will have the same subset of features selected (as would be the case with typical feature selection).
In order to more closely match the needs of many real-world applications, we propose a dynamic feature acquisition model that not only makes predictions with incomplete/missing features, but also determines what next feature would be the most valuable to obtain for a particular instance.

Our proposed model is deployed in an $\mathbb{R}^d$ feature space, where we have a subset (perhaps empty) of observed features $o \subseteq \{1,\ldots,d\}$ with values $x_o \in \mathbb{R}^{|o|}$ and we can acquire new features from the unobserved subset $u = \{1,\ldots,d\} \setminus o$. To simplify the problem, we utilize a sequential acquisition strategy, where only one feature is acquired at each step. There are two desired properties for acquiring a new feature: the feature should be informative for the target variable $y$; the feature contains non-redundant information which cannot be inferred from $x_o$. We leverage the conditional mutual information, $I(x_i;y \mid x_o)$, to quantify the dependency between the target variable $y$ and each candidate feature $i \in u$. The feature with maximum mutual information will be the next feature to acquire. Before we actually perform the acquisition, we do not know the exact value of $x_i$. After acquiring its value, we add the newly acquired feature $i$ to the observed subset and proceed to the next acquisition step if necessary.

Estimating mutual information in general is non-trivial \cite{kraskov2004estimating}. We need to capture the joint distribution, $p(x_i, y \mid x_o)$, and the marginal distributions $p(x_i \mid x_o)$ and $p(y \mid x_o)$. Our paradigm introduces another level of complexity. Since we require our model to be able to reason regardless of the set of previously observed features, both $x_o$ and $x_u$ could be arbitrary subsets. Not only does this mean the dimensionality of $x_o$ and $x_u$ can be arbitrary, but the type of dependencies we need to capture also scales exponentially. A naive method where different models $p_o(\cdot \mid x_o)$ are built for each different conditioning $x_o$ will quickly fail since it requires an exponential number of models with respect to the dimensionality $d$. In this work, we leverage a recently proposed generative model, ACFlow \cite{li2019flow}, to learn the arbitrary conditional distributions with one unified model.

In summary, our model sequentially selects features to acquire their actual values and extend the observed feature subset accordingly until it reaches the stopping criterion (e.g., acquisition budget). Then, the model will make a final prediction using the current acquired feature subset.
Our model is different from conventional feature selection, where a \emph{fixed} subset of features is selected for all instances; in contrast,  our model selects a specialized subset for each instance \emph{dynamically}.

Since the observed subset is updated as the acquisition process runs, we need to recalculate $I(x_i ; y \mid x_o)$ at each acquisition step. For an $\mathbb{R}^d$ feature space, each acquisition step has $O(d)$ time complexity. We further propose to leverage the Bayesian network structure over $x$ and $y$ to reduce the search space. We also propose a modified Grow-Shrink \cite{margaritis2000bayesian} algorithm to learn the Bayesian network structure for non-Gaussian data.

In addition to general real-valued data, we also apply our proposed method on time series, where the acquired features must follow the chronological order. Inspired by Thompson sampling \cite{thompson1933likelihood}, we integrate our knowledge about the order with a prior distribution and perform the acquisition based on the posterior. We also consider time-critical applications, where one may want to make a prediction as early as possible and acquire features consecutively but stop acquiring when prediction reaches a specified confidence threshold.

Our contributions are as follows: 1) We extend flow models, specifically ACFlow, to perform the dynamic feature acquisition task, which has not been explored previously. 2) We demonstrate the advantage of leveraging Bayesian network to reduce the searching space for feature acquisition. 3) We propose a modified Grow-Shrink algorithm to learn the Bayesian network structure from data, which leverages a deep generative model for conditional independence tests. 4) We propose a time series feature acquisition strategy to integrate our chronological prior knowledge. 5) We achieve state-of-the-art performance among information based dynamic acquisition methods on both synthetic and real-world datasets.

\section{Method}
In this section, we first formally describe the dynamic feature acquisition problem. Then, we introduce our method for acquiring new features. We describe our methods in detail for various settings. We also briefly review the ACFlow model \cite{li2019flow}, which we utilize to learn arbitrary conditional distributions.

\subsection{Problem Formulation}
Consider a real-valued feature vector, $x \in \mathbb{R}^d$, and a target variable, $y$. We do not have access to all the entries of $x$. Instead, we only have an observed feature subset $o$ with values $x_o \in \mathbb{R}^{|o|}$, where $o \subseteq \{1,\ldots,d\}$ (may be empty), and we can acquire more feature values $x_u$, where $u = \{1,\ldots,d\} \setminus o$ to further improve our prediction. To simplify the problem, we acquire one feature $i$ at each acquisition step, where $i \in u$. After, the acquired feature $i$ is observed, added to the observed set, and one proceeds to the next acquisition step if necessary. The goal of dynamic feature acquisition is to acquire as few features as possible while predicting $y$ as accurately as possible. In next section, we describe our method to determine which feature, $i$, to acquire at each acquisition step.

\subsection{Dynamic Feature Acquisition (DFA)}
In this section, we assume access to the arbitrary conditional distributions and leverage them for dynamic feature acquisition. Those conditional distributions can be estimated from a pretrained ACFlow model. We defer the training process to the next section.

The goal of feature acquisition is to increase the accuracy of our prediction, therefore, we would like the acquired feature to contain as much information about target variable as possible. We use the conditional mutual information to measure the amount of information. Specifically, we estimate the mutual information between each $x_i$ and $y$ conditioned on current observed set $x_o$, i.e., $I(x_i ; y \mid x_o)$. At each acquisition step, the feature with the highest mutual information will be our choice. Note that before we actually acquire the feature, we do not know the exact value of $x_i$, and we never observe the target variable $y$. Therefore, we need the arbitrary conditional distributions to infer these quantities. After the acquisition, the newly acquired feature, as well as its value, are added to the observed subset and the model proceeds to acquire the next feature. See Algorithm \ref{alg:dfa} for pseudo code of our acquisition approach.

\begin{algorithm}[t]
\SetAlgoLined
\KwIn{Training dataset $(\mathbf{X_t},Y_t)$ with all features and labels available; Test dataset $\mathbf{X_e}$ with no features available at the beginning; Total feature dimension $d$;}
1. Train ACFlow on training dataset by optimizing \eqref{eq:obj} \\
2. Dynamic feature acquisition and prediction: \\
\ForEach{$x \in \mathbf{X_e}$}{
    $o \leftarrow \emptyset; u \leftarrow \{1,\ldots,d\}; x_o \leftarrow \emptyset$ \tcp*[f]{initialize the observed set as empty}\\
    \Repeat{Stopping criterion reached}{
    estimate $I(x_i ; y \mid x_o)$ for $\forall i \in u$ \tcp*[f]{compute reward for each candidate}\\
    $i^* \leftarrow \argmax_i I(x_i ; y \mid x_o)$ \tcp*[f]{select the one with the highest reward}\\
    acquire the feature value $x_{i^*}$ \tcp*[f]{acquire the value of selected feature}\\
    $o \leftarrow o \cup i^*; u \leftarrow u \setminus i^*; x_o \leftarrow x_o \cup x_{i^*}$ \tcp*[f]{update observed and unobserved set}\\
    predict $y = \argmax p(y \mid x_o)$ \tcp*[f]{predict using current acquired subset}\\
    \tcp*[f]{stop if prediction is good enough or out of acquisition budget}
    }
}
\caption{Dynamic Feature Acquisition with ACFlow}
\label{alg:dfa}
\end{algorithm}

The conditional mutual information can be factorized as follows:
\begin{align}
    I(x_i ; y \mid x_o) &= I(x_i ; y, x_o) - I(x_i ; x_o) \label{cmi:0}\\
    &= H(x_i \mid x_o) - \mathbb{E}_{y \sim p(y \mid x_o)} H(x_i \mid y, x_o) \label{cmi:1}\\
    &= H(y \mid x_o) - \mathbb{E}_{x_i \sim p(x_i \mid x_o)} H(y \mid x_i, x_o) \label{cmi:2},
\end{align}
where $H(\cdot)$ represents the entropy.
From \eqref{cmi:0}, we see our acquisition metric prefers features that contain more information about $y$ without redundant information already in $x_o$.
According to \eqref{cmi:1}, when $x_o$ is observed, we prefer $x_i$ with higher entropy, that is, if we are already certain about the value of $x_i$, we do not need to acquire it anymore. Since $H(y \mid x_o)$ is fixed at each acquisition step, from \eqref{cmi:2}, we prefer feature that can reduce the most of our uncertainty about the target $y$.
In the following sections, we will describe how to estimate the conditional mutual information in different settings.

\subsubsection{Classification and Regression}
The real-valued target variable $y$ can be concatenated with $x$ as inputs to the ACFlow model so that $p(y \mid x_o)$ can be predicted by the same ACFlow as well.
The conditional mutual information, by definition, is
\begin{equation}\label{cmi:reg}
    I(x_i ; y \mid x_o) = \mathbb{E}_{p(x_i, y \mid x_o)} \log \frac{p(x_i, y \mid x_o)}{p(x_i \mid x_o)p(y \mid x_o)} = \mathbb{E}_{p(x_i, y \mid x_o)} \log \frac{p(y \mid x_i,x_o)}{p(y \mid x_o)}.
\end{equation}
We then perform a Monte Carlo estimation by sampling from $p(x_i, y \mid x_o)$. Note that $p(y \mid x_i, x_o)$ is evaluated on sampled $x_i$ rather than the exact value, since we have not acquired its value yet.

For a discrete target variable $y$, we employ Bayes's rule:
\begin{equation}\label{eq:bayes_rule}
    P(y \mid x_i,x_o) = \frac{p(x_i,x_o \mid y)P(y)}{\sum_{y'} p(x_i,x_o \mid y')P(y')} = \softmax_y(\log p(x_i,x_o \mid y') + \log P(y')).
\end{equation}
$P(y \mid x_o)$ can be similarly computed. The benefit of conditioning on $y$ is that we do not need an external classifier for prediction; we can use the same ACFlow model for both feature acquisition and prediction. To estimate conditional mutual information for classification problems, we can further simplify \eqref{cmi:reg} as
\begin{equation}\label{cmi:cls}
    I(x_i ; y \mid x_o) = \mathbb{E}_{p(x_i \mid x_o)P(y \mid x_i, x_o)} \log \frac{P(y \mid x_i,x_o)}{P(y \mid x_o)}
    = \mathbb{E}_{p(x_i \mid x_o)} \text{KL}[P(y \mid x_i, x_o) \| P(y \mid x_o)],
\end{equation}
where the KL divergence between two discrete distributions can be analytically computed. Note $x_i$ is sampled from $p(x_i \mid x_o)$ as before. We again use Monte Carlo estimation for the expectation.

\subsubsection{Pruning the Searching Space}
\begin{wrapfigure}{r}{0.22\textwidth}
\vspace{-15pt}
\begin{center}
    \includegraphics[width=0.21\textwidth]{bn.pdf}
\end{center}
\caption{An example of Bayesian Network.}
\label{fig:bn}
\end{wrapfigure}

As shown in Algorithm \ref{alg:dfa}, we need to estimate the conditional mutual information, $I(x_i ; y \mid x_o)$, for all possible candidates at each acquisition step. This induces a $O(d)$ time complexity at each acquisition step in $\mathbb{R}^d$ feature space. We propose to leverage the Bayesian Network to prune the searching space. Since our ultimate goal is to predict $y$, if observing $x_o$ makes part of $x_u$ independent of $y$, then we do not need to acquire those features. Note that acquiring new features might also introduce dependencies. For instance, in Fig.~\ref{fig:bn}, if we observe $x_2$, then $x_4$ is independent of $y$, therefore we do not need to acquire $x_4$ anymore; if we observe $x_1$, then $x_3$ becomes dependent, so we will need to reconsider $x_3$ as a candidate.

This is different from using the Markov Blanket of $y$ for feature selection. For instance, in Fig.~\ref{fig:bn}, the Markov Blanket will constrain the searching space to $\{x_1, x_2, x_3\}$, and at each step a subset of the Markov Blanket is considered as candidates. However, in our method, the candidates are dynamically chosen; it can add or delete some features based on the current acquired features, which means our method potentially selects over fewer candidates.

\subsubsection{Bayesian Network Structure Learning}
As shown in the last section, the Bayesian network can be leveraged to reduce the searching space. However, the groundtruth BN is generally not known with real-world data. Therefore, we would like to learn the Bayesian network structure from data. Conventional methods typically consider either discrete or Gaussian distributed data. Here, we propose to utilize the dependencies captured by an ACFlow model to infer the Bayesian network in order to deal with non-Gaussian data.

\begin{wrapfigure}{r}{0.55\textwidth}
    \centering
    \vspace{-15pt}
    \subfigure[]{\includegraphics[width=0.11\textwidth]{bnlearn.pdf}\label{fig:bnlearn:a}}
    \hspace{3mm}%
    \subfigure[]{\includegraphics[width=0.11\textwidth]{bnlearn_step1.pdf}\label{fig:bnlearn:b}}
    \hspace{3mm}%
    \subfigure[]{\includegraphics[width=0.11\textwidth]{bnlearn_step2.pdf}\label{fig:bnlearn:c}}
    \hspace{3mm}%
    \subfigure[]{\includegraphics[width=0.11\textwidth]{bnlearn_step3.pdf}\label{fig:bnlearn:d}}
    \caption{

    Structure learning.
    (a) The groundtruth Bayesian network. (b) The moral graph where each node is linked to its Markov Blanket with undirected edges. (c) PDAG returned by GS algorithm. Note the edge between A and B is undirected. (d) Resolve the orientation for remaining undirected edges. Here, the orientation is randomly chosen, since both directions do not introduce v-structures.}
    \label{fig:bnlearn}
    \vspace{-7pt}
\end{wrapfigure}

Our method is based on the Grow-Shrink (GS) algorithm \cite{margaritis2000bayesian}. Grow-Shrink is a constraint based structure learning algorithm. It starts by identifying the Markov Blanket for each node (feature). After getting the Markov Blanket, we can build a moral graph \cite{lauritzen1988local}, where parents and children are linked with undirected edges, spouses are linked together as well (Fig.~\ref{fig:bnlearn:b}). Then conditional independence tests are performed to detect and delete spouse links. Edges involved in v-structures can further be oriented by performing conditional independence tests. The GS algorithm returns a partial directed acyclic graph (PDAG) where some edges cannot be oriented and remain undirected (Fig.~\ref{fig:bnlearn:c}).  We refer readers to \cite{margaritis2000bayesian,pellet2008using} for details and pseudo code of the GS algorithm. The final step is to orient the remaining undirected edges. Following \cite{chickering2002learning}, we orient those edges in a way that they will not introduce v-structures (Fig~\ref{fig:bnlearn:d}).

We first propose a way of identifying Markov Blanket for both features $x$ and target $y$ using ACFlow. For notation succinctness, we concatenate $y$ into $x$ and denote them as $\mathrm{x}$. Specifically, we demonstrate the following theorem (see Appendix~\ref{sec:proof} for proof):
\begin{theorem}\label{thm:bn}
Let $\mathrm{x} \in \mathcal{X}^d$ (either discrete or continuous or mixed), $MB(\cdot)$ represents the Markov Blanket of the input, $I(\mathrm{x}_i ; \mathrm{x}_j \mid \mathrm{x}_{\widebar{\{i,j\}}}) = 0 \iff \mathrm{x}_i \notin MB(\mathrm{x}_j) \land \mathrm{x}_j \notin MB(\mathrm{x}_i)$, where  $\mathrm{x}_{\widebar{\{i,j\}}}$ represents all other variables except $\mathrm{x}_i$ and $\mathrm{x}_j$.
\end{theorem}
The Markov Blanket for a specific feature $i$ can be identified by estimating the conditional mutual information $I(\mathrm{x}_i;\mathrm{x}_j \mid \mathrm{x}_{\widebar{\{i,j\}}})$ in a similar fashion as \eqref{cmi:reg} and \eqref{cmi:cls} for each $j \neq i$. If $I(\mathrm{x}_i;\mathrm{x}_j \mid \mathrm{x}_{\widebar{\{i,j\}}})$ is larger than a specified threshold $\epsilon$, we take $\mathrm{x}_i$ and $\mathrm{x}_j$ as conditionally dependent given $\mathrm{x}_{\widebar{\{i,j\}}}$. We hence add $\mathrm{x}_i$ and $\mathrm{x}_j$ to the Markov Blanket of $\mathrm{x}_j$ and $\mathrm{x}_i$ respectively.

Typical applications of GS algorithm assume that features are Gaussian distributed, therefore $\chi$-squared tests are used to test conditional independence. To alleviate the Gaussian assumption, we leverage the ACFlow model for these tests. We again threshold the conditional mutual information between corresponding variables to test the conditional independence. For instance, in order to test if $\mathrm{x}_i$ and $\mathrm{x}_j$ are independent conditioned on $\mathrm{x}_c$, we estimate $I(\mathrm{x}_i;\mathrm{x}_j \mid \mathrm{x}_c)$. If $I(\mathrm{x}_i;\mathrm{x}_j \mid \mathrm{x}_c) \leq \epsilon$, we regard them as conditionally independent and vice versa.

\vspace{-5pt}
\subsubsection{Time Series}\label{sec:ts}
\vspace{-5pt}
In this section, we apply our method on time series data. For example, consider a scenario where sensors are deployed in the field with very limited power. We would like the sensors to decide when to put themselves online to collect data. The goal is to make as few acquisitions as possible while still making an accurate prediction. In contrast to ordinary vector data, the acquired features must follow a chronological order, i.e., the newly acquired feature $i$ must occur after all elements of $o$ (since we may not go back in time to turn on sensors). In this case, it is detrimental to acquire a feature that occurs very late in an early acquisition step. For example, for a time series with total $T$ time steps, if we acquire feature $x_{T-1}$ as the first feature, then we will lose the opportunity to observe features ahead of it.

Inspired by Thompson sampling \cite{thompson1933likelihood,russo2017tutorial}, we employ a prior distribution to encode our chronological constraint. Specifically, we set the prior as a Dirichlet distribution that is biased towards the selection of earlier time steps:
$
    \pi(\rho) = \text{Dir}\left[\alpha(T- (\max(o)+1)),\ldots, \alpha(T-(T-1)) \right](\rho),
$
where $\alpha$ is a hyperparameter, $T$ is the total time steps, $\max(o)$ represents the latest time step already acquired, and $\rho$ is a distribution for acquisition over the remaining future time steps.
However, we still desire that the acquired features are informative for target $y$. Hence, we update the prior to a posterior using time steps $V$ that are drawn according to how informative they are:
\begin{equation}
    p(V_n=t) \propto \exp(I(x_t ; y \mid x_o)),\ t \in \{\max(o)+1,\ldots, T-1 \},\ n \in \{1,\ldots,N\},
\end{equation}
where $N$ is the total number of samples.
Due to conjugacy, the posterior is also a Dirichlet distribution
\begin{equation}
    p(\rho \mid V) = \text{Dir}\left[\alpha(T- (\max(o)+1)) + \sum_{n=1}^{N} \mathbb{I}\{V_n=\max(o)+1\}, \ldots  \right](\rho).
\end{equation}
Samples from posterior represent the probabilities of choosing each candidate, which now prefer both earlier time steps and informative features. We draw a sample from posterior and select the most likely time step at each acquisition step.

In some cases, the number of sensor readings is not as important as making a prediction as early as possible. That is, sensors may be collecting data at each time step towards a prediction, but one would like to make a final prediction as soon as a confidence threshold is hit. This setting is especially applicable in time-critical applications such as an autonomous vehicle predicting if a nearby pedestrian will cross the street. In this setting, we propose to use the prediction probability, $\max_y p(y \mid x_o)$, as the stopping criterion. When the probability reaches the specified threshold, we stop acquiring more features (time steps) and predict the target variables. Note that deep models usually underestimate uncertainty \cite{gal2016uncertainty}. Therefore, the actual accuracy may be lower than the specified threshold. We perform uncertainty calibration following \cite{kumar2019verified}. Note we calibrate uncertainty for each time step separately by utilizing a held-out validation set.

\vspace{-8pt}
\subsection{Arbitrary Conditional Flow}
\vspace{-5pt}
Above we assumed access to the arbitrary conditional distributions via a pretrained ACFlow model \cite{li2019flow}. Here, we give a brief introduction to the ACFlow model and our modified training procedure in different settings.
ACFlow extends flow models to learn the arbitrary conditional distributions $p(x_u \mid x_o)$.
Since both $x_u$ and $x_o$ could be an arbitrary subset of $x$, one needs to model an exponential number of conditionals with respect to the dimensionality of $x$. \cite{li2019flow} utilizes a masking mechanism to train ACFlow in a multi-task fashion; i.e., \emph{all} conditionals are captured by a single model. The multi-task training mechanism allows the model to generalize to unseen combinations of $x_u$ and $x_o$.

The original training process of ACFlow selects two non-overlapping subsets of features at random as $x_u$ and $x_o$ for each training example. Then, the training objective is to maximize the arbitrary conditional log likelihoods $\log p(x_u \mid x_o)$. In our case, we train ACFlow to optimize
\begin{equation}\label{eq:obj}
    \mathcal{L} = \log p(x_i, y \mid x_o) = \log p(x_i \mid x_o) + \log p(y \mid x_i, x_o),
\end{equation}
where $x_o$ is an arbitrary subset, while $x_i, i \in u$ is an arbitrary element in unobserved subset. Note that the actual value of $x_i$ is used during training.

For real-valued target variable, we concatenate $x$ and $y$ so that $p(x_i \mid x_o)$ and $p(y \mid x_i, x_o)$ can both be evaluated by ACFlow.
For discrete target variable, we train a conditional version of ACFlow conditioned on $y$. Then $P(y \mid x_i, x_o)$ can be evaluated by Bayes rule as in \eqref{eq:bayes_rule}.  $p(x_i \mid x_o)$ can be computed by marginalizing out $y$:
$
    p(x_i \mid x_o) = \frac{\sum_{y'} p(x_i,x_o \mid y') P(y')}{\sum_{y'} p(x_o \mid y')P(y')},
$
where $p(x_i, x_o \mid y')$ and $p(x_o \mid y')$ are conditional likelihoods from the conditional ACFlow.

\vspace{-5pt}
\section{Related Works}
\vspace{-5pt}
\paragraph{Active Learning}
Traditional active learning \cite{settles2009active} performs instance level acquisition, where new instances are selected to acquire their labels to train a better model. Information based acquisition functions play an important role in this area \cite{mackay1992information,houlsby2011bayesian}. In contrast, DFA does not consider a labeling oracle at training time. Instead, we perform a feature level acquisition, where new features are acquired for the same instance at evaluation time. In addition, we require the model to make prediction based on partial observations, while active learning typically assumes all features are observed for selected instances.

\vspace{-7pt}
\paragraph{Feature Selection} Conventional methods for feature selection eliminate redundant features, which can help reduce computation complexity and improve generalization \cite{chandrashekar2014survey}. The same selected subset is then applied to each instance. In contrast, in DFA
a personalized subset of features is selected for each instance depending on the specific values that observed features held. Thus, DFA allows for an instance specific strategy to select informative features. It is worth noting that DFA may be applied after an initial feature selection preprocessing step to reduce the search space.

\vspace{-7pt}
\paragraph{Dynamic Feature Acquisition}
Acquiring features dynamically is of great use in many real-world decision-making scenarios. Several methods have been proposed to actively acquire new features in cost-sensitive setting \cite{chai2004test,ling2004decision,sheng2006feature,nan2014fast}. In \cite{zubek2002pruning,ruckstiess2011sequential}, the dynamic feature acquisition problem is cast as a Markov Decision Process. \cite{shim2018joint} further propose a deep reinforcement learning based method. As in active learning, using an information-based acquisition metric is a promising direction, but it imposes a requirement to learn the arbitrary conditional distributions. A recent work EDDI \cite{ma2018eddi} utilizes a VAE based generative model to learn arbitrary marginal distributions $p(x_o)$. Since the VAE assumes the likelihood $p(x_i \mid x_o)$ for each feature $x_i$ is conditionally independent, they hence evaluate $p(x_i \mid x_o)$ by $p(x_i \mid z)$, where $z$ is the latent code. In our method, we directly model the arbitrary conditional distributions $p(x_i \mid x_o)$ with a flow based generative model.

\vspace{-7pt}
\paragraph{Bayesian Network Structure Learning} Learning the Bayesian network structure has been studied for decades. There are two primary classes of algorithms: constraint-based algorithms, where conditional independence tests are utilized to infer the dependency structure, and score-based algorithms, which maximize goodness-of-fit scores \cite{scutari2018learns,scanagatta2019survey}. Hybrid methods that combine both approaches have also been considered \cite{scutari2018learns,scanagatta2019survey}. In this work, we extend the constraint-based Grow-Shrink algorithm \cite{margaritis2000bayesian} with deep generative models to exploit the benefits of better distribution modeling.

\vspace{-5pt}
\section{Experiments}
\vspace{-5pt}
In this section, we evaluate our model in different settings by comparing to EDDI \cite{ma2018eddi}, a state-of-the-art information based feature acquisition method. We evaluate two acquisition strategies: the dynamic feature acquisition (DFA), where specialized acquisition orders are applied for different instances, and the static feature acquisition (SFA), where the same acquisition order is applied for each instance. The static order is similarly decided by the conditional mutual information metric, but averaged over all test instances. For estimating the conditional mutual information, we draw 10 samples from the corresponding distributions. We also compare to a decision tree based approach, where multiple trees with different depth are fitted. The depth controls the number of acquired features.

\vspace{-5pt}
\paragraph{Classification and Regression}
\begin{figure}
    \vspace{-0.1in}
    \centering
    \subfigure[example of acquired features]{\includegraphics[width=0.52\textwidth]{mnist_exp.png}\label{fig:mnist:a}}
    \subfigure[test accuracy]{\includegraphics[width=0.4\textwidth]{mnist_cls.png}\label{fig:mnist:b}}
    \caption{Classification on MNIST. (a) two examples of the acquisition process. The blue mask indicates unobserved features. We also plot the prediction probability at each acquisition step. (b) prediction accuracy after each acquisition step averaged over the test set.}
    \label{fig:mnist}
    \vspace{-0.2in}
\end{figure}

\begin{figure}
    \begin{minipage}{0.48\textwidth}
    \centering
    \subfigure[synthetic]{\includegraphics[width=0.49\textwidth]{syn_cls.png}}
    \subfigure[gas]{\includegraphics[width=0.49\textwidth]{gas_cls.png}}
    \caption{Test Accuracy. Higher is better.}
    \label{fig:uci_cls}
    \end{minipage}
    \vspace{-0.1in}
    \hfill
    \begin{minipage}{0.5\textwidth}
    \centering
    \subfigure[housing]{\includegraphics[width=0.48\textwidth]{house_reg.png}}
    \subfigure[whitewine]{\includegraphics[width=0.49\textwidth]{wine_reg.png}}
    \caption{Test RMSE. Lower is better.}
    \label{fig:uci_reg}
    \end{minipage}
    \vspace{-0.1in}
\end{figure}

In this section, we demonstrate that our model can achieve accurate classification and regression by acquiring only a small subset of features.
We first conduct experiments on MNIST dataset \cite{lecun1998gradient}. To reduce the number of total features, we downsample the images to $16 \times 16$. Data preprocessing and architecture details are listed in Appendix~\ref{sec:cls_reg}. The results as well as several examples are shown in Fig.~\ref{fig:mnist}. From Fig.~\ref{fig:mnist:a}, we can see our model acquires different features for different instances, which demonstrates that our model is capable of identifying specialized feature subsets for different instances. After acquiring only a small subset of features, the model is already certain about the target variable. It is also worth noting that our model is capable of exploring the spatial correlation between features. For instance, the acquired features present a checkerboard-like pattern, which is consistent with our impression that nearby pixels usually contain redundant information. Figure.~\ref{fig:mnist:b} compares the test accuracy after each acquisition step. Our model consistently outperforms EDDI at each acquisition step. We can see the DFA policy is always better than the SFA policy for both EDDI and ACFlow, which verifies the benefits of acquiring specialized subsets. We also notice that ACFlow with a static policy can already outperform EDDI with a dynamic policy, which we conjecture is due to the superior conditional distribution modeling ability of ACFlow. First, it improves the prediction performance since $p(y \mid x_o)$ is better captured by ACFlow. Second, it improves the conditional mutual information estimation, which gives better acquisition policy.

We then evaluate our model on real-valued vector data from both synthetic dataset and UCI datasets \cite{ucilichman}. Please refer to Appendix~\ref{sec:cls_reg} for experimental details. We employ a validation set to select the best architecture and hyperparameters for EDDI. Figures~\ref{fig:uci_cls} and~\ref{fig:uci_reg} present the results for classification and regression. We report the test accuracy and root mean squared error (RMSE) after each acquisition step for classification and regression tasks, respectively. We see our model achieves superior performance compared to EDDI. Additionally, the dynamic acquisition policy tends to give better results compared to a static alternative. Although the decision tree based method utilizes multiple models, our method exceeds it using only one unified model.

\begin{figure}
    \vspace{-0.1in}
    \begin{minipage}{0.48\textwidth}
    \centering
    \subfigure[asia]{\includegraphics[width=0.49\textwidth]{asia_bn.png}}
    \subfigure[sachs]{\includegraphics[width=0.49\textwidth]{sachs_bn.png}}
    \caption{On synthetic datasets, our learned BN performs equally well compared to groundtruth.}
    \label{fig:bn_syn}
    \end{minipage}
    \hfill
    \vspace{-0.1in}
    \begin{minipage}{0.48\textwidth}
    \centering
    \subfigure[housing]{\includegraphics[width=0.48\textwidth]{house_bn.png}}
    \subfigure[whitewine]{\includegraphics[width=0.49\textwidth]{wine_bn.png}}
    \caption{On UCI datasets, our learned BN performs better than the one learned from GS.}
    \label{fig:bn_reg}
    \end{minipage}
\end{figure}

\begin{figure}
    \vspace{-0.1in}
    \begin{minipage}{0.48\textwidth}
    \centering
    \subfigure[digits]{\includegraphics[width=0.49\textwidth]{digits_ts.png}}
    \subfigure[pedestrian]{\includegraphics[width=0.49\textwidth]{ped_ts.png}}
    \caption{Time series classification accuracy.}
    \label{fig:ts_cls}
    \end{minipage}
    \hfill
    \vspace{-0.1in}
    \begin{minipage}{0.48\textwidth}
    \centering
    \subfigure[digits]{\includegraphics[width=0.48\textwidth]{digits_uty.png}}
    \subfigure[pedestrian]{\includegraphics[width=0.49\textwidth]{ped_uty.png}}
    \caption{Accuracy and calibrated probability.}
    \label{fig:ts_uty}
    \end{minipage}
\end{figure}

\vspace{-5pt}
\paragraph{Bayesian Network}
In this section, we utilize the Bayesian network structure to accelerate the acquisition process. We first conduct experiments on two synthetic datasets from Bayesian network repository \cite{bnrepository}, where we know the underlying Bayesian network structure. Data are sampled from conditional Gaussian distributions. Please refer to Appendix~\ref{sec:bn} for data generation process. The results are presented in Fig.~\ref{fig:bn_syn} by comparing our learned BN with the groundtruth BN. We also show the results without using BN for comparison. We notice that utilizing BN does not affect the performance and our learned BN performs equally well compared to the groundtruth cases. Compared to ACFlow DFA without BN, on average our method reduces the candidates set size by 19.6\% and 13.5\% respectively for \textsl{asia} and \textsl{sachs} datasets. We then evaluate the learned BN on UCI datasets for regression in Fig.~\ref{fig:bn_reg}. Compared to the BN learned by the original GS algorithm, our learned BN tends to achieve lower RMSE with the same acquisition budget. Our learned BN also reduces the candidates set size by 28.2\% and 3.6\% for \textsl{housing} and \textsl{whitewine} respectively.

\vspace{-5pt}
\paragraph{Time Series}
In this section, we test our method on time series data. We compare to EDDI by applying the same acquisition policy described in Sec.~\ref{sec:ts}. Figure~\ref{fig:ts_cls} shows the results on two datasets from UEA \& UCR time series classification repository \cite{bagnall16bakeoff}. We can see our model outperforms EDDI by a large margin. Our model reaches the plateau using fewer features, indicating that our model is more efficient for acquisition. We then test the consecutive acquisition case where we use the prediction probability as the stopping criterion. In Fig.~\ref{fig:ts_uty}, we plot the accuracy at each time step (as opposed to acquisition step) and the corresponding prediction probability. We calibrate the probability for each time step separately based on a held-out validation set. We see that after calibration, the estimated probability matches the true accuracy. Thus, one can use the probability as the stopping criterion to make a prediction as soon as possible on a per-instance basis.
We see that by setting a threshold at 90\%, on average, we would only need to acquire 5.01 and 10.69 time steps for \textsl{digits} and \textsl{pedestrian} datasets respectively, and obtain the requested 90\% accuracy.

\vspace{-5pt}
\section{Conclusion}
\vspace{-5pt}
In this work, we present a framework for dynamic feature acquisition which acquires new features sequentially to improve prediction. We utilize an information theoretic metric, conditional mutual information, to determine which feature to acquire next. The conditional mutual information is estimated by leveraging a flow based generative model, ACFlow, to capture the arbitrary conditional distributions. We conduct experiments in different settings and demonstrate superior performance over strong baselines. To overcome the potential high time complexity for high-dimensional data, we propose to learn the Bayesian network from data using the same ACFlow model. The learned Bayesian network can help reduce the searching space and thus reduce the time complexity. In future work, we will explore this framework in other settings, such as spatial-temporal feature acquisition.

\section*{Broader Impact}
Although human agents routinely reason on instances with incomplete and muddied data (and weigh the cost of obtaining further features), much of machine learning (ML) is devoted to the unrealistic, sterile environment where all features are observed, and further information on an instance is obviated. Thus, the current paradigm in ML is woefully unprepared for the future of automation, which will be highly interactive and able to obtain further information when making predictions. To enable ML for a future of automation that makes decisions in realistic, uncertain environments, we propose to develop dynamic feature acquisition (DFA) methods that can intelligently make assessments on instances in the face of incomplete and costly features.

The impact of DFA would be broadly felt in a myriad of domains across the sciences and business. For instance, DFA could inform autonomous vehicles when it is possible to accurately make an early prediction based on the limited data collected. In these time-critical applications, DFA has the potential to save lives. DFA could also make a direct impact on applications where the computer is actively interacting with a user to gather information and make an assessment. Examples include surveys and costumer service chat-bots, where DFA would enable such systems to determine the most informative next question and to quickly and accurately make predictions based on the limited answers collected. Related applications of interest are AI education systems that personalize curricula and materials to each student. DFA would enable an education system to quickly assess the student’s current knowledge and abilities with a limited number of questions, which would allow the system to expedite the curriculum and focus on concepts not well understood by the student. DFA would also have a deep impact for applications that collect data out in the field. For instance, in environmental applications agents must physically collect samples from various locations to assess if a region has been contaminated. DFA would not only potentially reduce the number of locations sampled, but may also be used to suggest the most informative locations to query with autonomous agents. DFA would also have impact in medical and health applications. For example, as aforementioned, DFA could sharply reduce the number of invasive biopsies collected on a patient when diagnosing in health-related applications.

Since DFA strives to reduce the number of acquired features and thus the cost of acquisition, a failure of the system would only decrease the efficacy. At worst, it will just degenerate to the case where we acquire all features.
However, like any other model, DFA might unintentionally learn spurious correlations and biases within the dataset. Thus, we encourage practitioners to carefully design the training set or utilize other debiasing techniques.

\bibliographystyle{unsrt}


\section{Introduction}
The incessant evolution of hydrogen bond network makes water as one of the most interesting liquid.\cite{Bagchi2011,Ball2015} In spite of its smaller size, water molecules are capable of behaving hydrogen bond donor and acceptor simultaneously. As a result a gigantic hydrogen bond network is abundant in bulk water. The ultrafast evolution of this massive hydrogen bond network leads to many unique features of water.\cite{Nandi2000,Bagchi2005,Bellissent-Funel2016} The presence of external solutes perturbs the length scale and the time scale of evolution of the hydrogen bond network that makes the systems even more fascinating.\cite{Ikeguchi1998,Kropman2001,Raschke2005,Bakulin2011,Grdadolnik2017}
\par
The hydrophobic effect is the manifestation of the interaction of non-polar moieties with water. The hydrophobic effect is typically inferred at two levels. First, the interaction between a non-polar molecule and the surrounding water molecules is called \textit{hydrophobic hydration}. Second is known as \textit{hydrophobic interaction} or \textit{pair hydrophobicity} which describes the interaction between two non-polar molecules in water as a function of their separation distance. The hydrophobic effect shows pivotal role in biological processes such as protein folding, formation of cell membrane, formation of vesicles and lipid bilayer, assembly of proteins into functional complexes, and many more.\cite{Tanford1973,Srinivas2002,Srinivas2003,Mason2004,Chandler2005,Berne2009,Hazra2014} The significant role in Biology and Chemistry makes the study of hydrophobic effect a subject of major interest till date. In spite of several earlier conclusive experimental and theoretical studies, the microscopic origin and the length scale of hydrophobic interaction are still not understood completely and yet remain an active area of research.
\par
The presence of many body interactions could give rise to different structural arrangement surrounding the hydrophobes.
\cite{Bakulin2011,Grdadolnik2017,Chandler2005,Wiggins1997,Lum1999,Huang2000} It is considered that the presence of hydrophobic solutes modifies the shape of the solvation shell into tiny \textit{icebergs}.
\cite{Frank1945,Kauzmann1959} This type of solvent arrangement significantly alters the hydrogen bond strength in the close vicinity of the hydrophobic moieties.\cite{Raschke2005,Bakulin2011,Grdadolnik2017,Koh2000,Scatena2001,Davis2012,Montagna2012} It has also been argued that for sufficiently weak solute-water attraction, a large smooth hydrophobic surface might be enclosed by a microscopically thin film of water vapor.\cite{Stillinger1973} The recent finding also shows that the hydrophobic solutes reinforce the water hydrogen bonds of solvation shell water molecules.\cite{Grdadolnik2017,Koh2000}
As the vibrational frequencies of water are highly sensitive to the local microscopic solvation configuration, infrared (IR) spectroscopy is a tempting method for exploration of these systems.\cite{Fecko2003,Biswas2012} IR spectroscopy of O-H stretching mode is usually employed as the most trustworthy and sensitive method for estimating relative strengths of hydrogen bonds.\cite{Davis2012,Fecko2003,Hecht1992,Hecht1993,Sharp2001,Auer2007,Laage2009,Stiopkin2011,Biswas2013,Biswas2016}
\par
In spite of its extensive use, meaningful interpretation of IR spectroscopic data face the formidable challenge of establishing the connection among the experimental spectra to molecular structures in the bulk phase. Experimental spectra represent a superposition of different transient solvation structures, hence illustrate the macroscopic response. Besides, the existence of strong anharmonic couplings leads to further delocalization through different vibrational modes mixing. On the contrary, the microscopic resolution of computer simulation aided spectroscopy modeling empowers us to investigate these systems at the molecular level.\cite{Auer2007,Biswas2013,Biswas2016,Corcelli2004,Corcelli2005,Roberts2009,Biswas2017,Samanta2018} Therefore, the theoretical spectroscopy offers a window to decipher the experimental data in rather quantitative fashion.\par
In spite of several earlier computational studies, the structure spectrum correlation of aqueous hydrophobic system needs further quantification using appropriate spectroscopic modeling. As in bulk water, vibrational modes of O-H oscillators are significantly coupled and give rise to very broad response, which makes it challenging to differentiate contribution from diverse transient configurations. To reduce the complexity of the problem, isotope dilution strategies are often used.\cite{Auer2007,Biswas2016,Corcelli2005,Samanta2018} The isotope dilute system eventually allows us to study the system in lower dimensionality by isolating a local oscillator. Thus, a small percentage of H$_2$O in a fully deuterated solution can provide isolated O-H oscillators, which can function as a local probe to different solvation structures.\par
In a number of earlier works, mixed quantum-classical (MQC) models have been extensively utilized to explain the O-H stretching vibrations in various isotope dilute aqueous systems.\cite{Auer2007,Biswas2016,Corcelli2004,Corcelli2005,Roberts2009,Samanta2018} Usually, these models identify a collective coordinate from classical molecular dynamics (MD) trajectories to depict the influence of the solvent on the quantum mechanical spectroscopic coordinate. The electronic structure calculation coupled with MD simulation is the fundamental basis of the MQC approach. The MQC models are widely used to construct the required spectroscopic maps, which can essentially be used to generate the trajectories of time-dependent transition frequencies and dipole moments.\par
In this work, we investigate the effect of hydrophobicity in water structure by employing computer simulation aided spectroscopy modeling which employs the microscopic environmental sensitivity of vibrational frequencies. The rest of this article is structured as follows: In section II, we will elaborate the spectroscopic modeling and simulation details, data analysis and major findings are discussed in section III and conclusions are given in section IV.
\section{methods}
\subsection{Classical molecular dynamics}
We perform molecular dynamics simulation using GROMACS version 2019.1.\cite{Abraham2015} The system consists of a cubic box of 255 SPC/E\cite{Berendsen1987} water molecules and 1 OPLS-AA methane molecule.\cite{Jorgensen1996} For bulk water, we take 256 water molecules. Periodic boundary conditions are applied in all three directions.\cite{Frenkel2002} We perform energy minimization of the systems using steepest decent method. Thereafter, the systems are equilibrated in \textit{NPT} ensemble for 1ns. Finally, the data acquisitions are done in \textit{NVT} ensemble over 5 ns long trajectory. We use Berendsen thermostat\cite{Berendsen1984} with a relaxation time of 0.1 ps and Parinello-Rahman barostat\cite{Parrinello1981,Nose1983} with a relaxation time of 1.0 ps for maintaining constant temperature at 300 K and pressure at 1 bar, respectively. We use a 9 \text{\AA} cutoff radius for neighbor searching and non-bonded interactions and all the bonds are kept fixed using LINCS.\cite{Hess1997} The long-ranged electrostatic interactions are calculated using Particle Mesh Ewald (PME) with FFT grid spacing of 1.6 \AA.\cite{Darden1993}

\subsection{Mixed quantum mechanical calculation}
We use mixed quantum classical approach to build the single oscillator O-H stretch spectroscopic modeling of methane in water. This technique is based on semi-empirical cluster-based mapping method, which has been widely used in simulating the O-H vibrational spectroscopy of isotope diluted aqueous systems.\cite{Biswas2016,Corcelli2004,Corcelli2005,Roberts2009} We utilize electronic structure calculations on instantaneous small frozen clusters extracted from the classical MD trajectory to construct the spectroscopic map. In this, we correlate the quantum mechanical O-H transition frequency and transition dipole moment onto a collective coordinate which can easily be calculated from the classical trajectory. Subsequently, by using this classical analogue of the collective variable, we generate the transition frequency and transition dipole trajectories, which can further be utilized to calculate linear and nonlinear IR spectra from respective time-domain response functions.\par
From the classical trajectory, we extract small clusters by identifying a central H atom that belongs to the water molecule closest to the methane molecule and include any molecules having its oxygen within a 7.0 \text{\AA} radius. We use such 165 methane-water clusters. For bulk water response, we choose 130 small clusters extracted in the similar fashion from bulk water trajectory. These clusters contain average of $\sim44$ water molecules, which is adequate to approximately resemble the bulk like environment around the central O-H bond, and slight fluctuation in the number of water molecules in a cluster does not affect our findings. Furthermore, these clusters are chosen in such a way so that we can sample all relevant configurations. A more detail description concerning the cluster selection is provided in the supplementary information.\par
We get the O-H stretch response for each selected cluster, from the one-dimensional adiabatic potential energy surface (PES). We construct the quantum mechanical PES, by stretching the central O-H bond from $r_{\mathrm{OH}} = 0.7$ to 1.6 \text{\AA} with a grid spacing of 0.08 \text{\AA}, while keeping the remaining degrees of freedom frozen. For getting the single point energy, we perform DFT calculations employing B3LYP hybrid functional\cite{Becke1993,Lee1988,Vosko1980,Stephens2002} and 6-311++G(d,p) basis set in Gaussian 16 package.\cite{Frisch2016} The choice of B3LYP functional is inspired by the earlier work of Skinner and co-workers, which demonstrated that potential energy obtained using B3LYP functional along the O-H stretch coordinate is in promising agreement with the coupled cluster CCSD(T) prediction.\cite{Gruenbaum2013} Although the classical simulations are executed with isotopically pure systems, we consider that all hydrogens but the proton of the selected O-H bond are deuterons. We obtain the eigenstates, $|n \rangle $, and eigenvalues, $E_n$, of PES $U(r_{\mathrm{OH}})$ by solving the one dimensional Schrodinger equation using the discrete variable representation (DVR)\cite{Colbert1992,Groenenboom1993} with a grid spacing of 0.01 \text{\AA} and reduced mass of the O-H vibration of the HOD molecule (0.954426 amu).\cite{Roberts2009} Then from the energy eigenvalues, we obtain the transition frequencies $\omega_{nm} = (E_n-E_m)/ \hbar$ between different vibrational states.\par
We calculate the transition dipole moments in a manner similar to Corcelli and Skinner.\cite{Corcelli2005} The transition dipole matrix elements for the transition between vibrational states n and m is expressed as:
\begin{equation}
   \vec{\mu}_{nm} = \langle n | \vec{\mu} |m \rangle
\end{equation}
After expanding the dipole operator about the minimum of the O-H stretching potential ($r_{eq}$ to the first order term in $r_\mathrm{{OH}}$, one can approximate $\vec{\mu}_{nm}$ as
\begin{equation}
    \vec{\mu}_{nm} = \langle n | \vec{\mu}_0 + r_{\mathrm{OH}}\left(\frac{d\vec{\mu}}{dr_{\mathrm{OH}}}\right)_{r_{\mathrm{OH}}=r_{eq}} |m \rangle \approx \mu^\prime r_{nm} \hat{\textbf{u}}
\end{equation}
where $r_{nm} = \langle n | r_{\mathrm{OH}} |m \rangle$ are the matrix elements obtained from the eigenstates of the DVR calculation, $\vec{\mu}_{nm}$ can be divided into a magnitude $\mu$ and a direction $\hat{\textbf{u}}$, which we consider as lying along the O-H bond axis and $\mu^\prime$ is the dipole moment derivative. We acquire $\mu^\prime$ for each chosen cluster by calculating $\vec{\mu}_{nm} \cdot \hat{\textbf{u}}$ at five $r_{\mathrm{OH}}$  displacements separated by 0.01 \text{\AA} about $r_{eq}$, and then numerically differentiate with respect to $r_{\mathrm{OH}}$ . Finally, we construct correlation maps for $\omega_{nm}$ and $\mu_{nm}$ against a collective coordinate.\par
Although it was reported earlier that the projected electric field on the O-H stretch could serve as a potential collective variable, which can describe the spectroscopy of isotope diluted aqueous system.\cite{Auer2007,Corcelli2005} However, we find it is not the case for the water methane system (see supplementary material). Therefore, we use an alternative collective variable as used in the case of aqueous hydroxide and aqueous proton systems.\cite{Roberts2009,Biswas2016} We use solvation energy as the collective variable defined as the difference in the potential energy at two reference points along $r_{\mathrm{OH}}$ .
\begin{equation}
    \Delta E_{\mathrm{DFT}}=U(r_f)-U(r_i)
\end{equation}
As the potential energy is highly sensitive to the local environment, of the solvation energy accurately captures the alteration in the PES as a function of the hydrogen bond strength. After rigorous tweaking, we find that the $r_i = 1.0 $ \text{\AA} and $ r_f = 1.4 $ \text{\AA} combination provides a better correlation with the spectroscopic variables.

\section{Results and Discussion}
\subsection{Spectroscopy maps}
We have mentioned earlier that hydrophobic solutes alter the water hydrogen bond network and thus change the local electronic environments of O-H oscillators. The intense local environmental sensitivity makes vibrational spectroscopy of O-H stretching mode $(\nu_{OH})$ a most reliable and sensitive approach for investigating relative strengths of H-bonds. The alteration of H-bond strength will result in shifting of observed frequency for $(\nu_{OH})$ mode. Thus, enhancements of the H-bond strength will eventually be reflected in frequency redshift of the $(\nu_{OH})$ mode when compared the spectra of bulk water with those of the water molecules perturbed by hydrophobic solute. Furthermore, structural order can also be investigated by examining the spectral line width. As more structural ordering will have less variations in the local structures and will eventually results into a narrower line width.\par
\begin{figure}
	\includegraphics[width=1.0\linewidth]{figure/fig_1.eps}
	\caption{\label{fig1}(a) Correlation between fundamental O-H transition frequency and solvation coordinate obtained from DFT calculation  using clusters taken from the classical MD. Black color filled circles represents the bulk water data, whereas all other colored filled symbols represent the data for methane-water cluster. The cyan solid line represents the fourth order polynomial fit. (b) Definition of the angle between methane hydrogen, water hydrogen and water oxygen.}
\end{figure}
The general observation is that the observed one-dimensional O-H stretch potential energy surfaces becomes less anharmonic in presence of methane, which results in O-H stretch frequency in the high frequency range of the bulk-water response. The correlation of the fundamental transition frequency and solvation energy obtained from the quantum mechanical calculation is presented in Figure \ref{fig1}. The O-H stretch frequencies are monotonic and highly correlated (with correlation coefficient 0.9964) through a nonlinear relationship with the collective solvation coordinate within the 1000 cm-1 range spanned by the different cluster configurations. Although, we select all the instantaneous configurations without any bias, we have further investigated the effect of orientation of the tagged O-H bond with the methane molecule. Note that, this orientational dependencies are investigated only for the O-H oscillators, which are directly facing to the methane molecule. This was achieved by sampling the configurations having angle $\angle \mathrm {O_WH_WH_{Me}}$ (as defined in Figure \ref{fig1}(b)) within a particular range. We find there exist minimal or almost no effect of the orientation of the O-H oscillator with respect to the methane moiety on the spectroscopy map. We use a $4^{th}$ order polynomial fitting function to get the empirical relations of the transition frequencies with solvation coordinate (given in Table \ref{table1}). We furthermore explore the spectroscopic maps for other transition frequencies and these are shown in the supplementary material (Figure \ref{figs3}).\par
\begin{figure}[b]
	\includegraphics[width=0.8\linewidth]{figure/fig_2.eps}
	\caption{\label{fig2}Correlation between fundamental O-H stretching frequency ($\omega_{10}$) and transition dipole matrix elements. Cyan line represents the fourth order polynomial fit.}
\end{figure}
\begin{table*}[]
\centering
\caption{Empirical relation between the transition frequency and the solvation energy obtained by fitting the quantum mechanical data as shown in Figure \ref{fig1}. The fitting function is $\omega_{10}=a_0+a_1 \Delta E_{\mathrm{DFT}}+a_2\Delta E_{\mathrm{DFT}}^2+a_3\Delta E_{\mathrm{DFT}}^3+a_4\Delta E_{\mathrm{DFT}}^4$. The parameters for other transition frequencies are given in supplementary information (Table S2).}
\begin{tabular}{llllll}
\hline
$a_0$    & $a_1$   & $a_2$        & $a_3$          & $a_4$         & Correlation Coefficient \\
\hline
\hline
1090.37 & 0.2636 & $4.7833\times 10^{-7}$ & $-9.0792\times10^{-10}$ & $3.0082\times10^{-14}$ & 0.9964\\
\hline
\end{tabular}
\label{table1}
\end{table*}
Subsequently, we construct the empirical maps for transition dipole moments. We follow the strategy as explained earlier. It was shown earlier that empirical relation can be established between the dipole moment derivative $\mu^\prime$ and the electric field projected along $r_{\mathrm{OH}}$  in the $\mathrm{HOD/D_2O}$ system or the solvation coordinate for the aqueous hydronium and hydroxide systems. In the present work, we find the best correlations between $\mu^\prime$ to $\Delta E_{\mathrm{DFT}}$ and $\mu_{nm}$ to $\omega_{nm}$. The correlation between the fundamental transition dipole moment $\mu_{10}$ and the fundamental transition frequency $\omega_{10}$ is shown in Figure \ref{fig2}. The fundamental transition dipole moments are monotonic and highly correlated (with correlation coefficient 0.9999) through a nonlinear relationship with the fundamental transition frequency (Figure \ref{fig2}). We use a $4^{th}$ order polynomial fitting function to get the empirical relation between fundamental transition dipole moment and fundamental transition frequency (fitting parameters are shown in Table \ref{table2}).\par
\begin{table*}
\centering
\caption{Empirical relationships of the fundamental O-H stretching frequency with transition dipole matrix element using $\mu_{nm}=a_0+a_1 \omega_{nm}+a_2\omega_{nm}^2+a_3\omega_{nm}^3+a_4\omega_{nm}^4$.}
\begin{tabular}{llllll}
\hline
$a_0$    & $a_1$   & $a_2$        & $a_3$          & $a_4$         & Correlation Coefficient \\
\hline
\hline
21.2154	& -0.0258	& $1.1895\times10^{-5}$ &	$-2.4285\times10^{-9}$	& $1.8459\times10^{-13}$	& 0.9999\\
\hline
\end{tabular}
\label{table2}
\end{table*}
Since the quantum mechanical calculation of the solvation energy at each time step for a tagged O-H bond is computationally expensive; this necessitates an appropriate quantity that will be adequate for estimating the transition frequencies from the classical MD trajectory. To overcome this, we construct a correlation map between the quantum mechanical solvation coordinate i.e. $\Delta E_{\mathrm{DFT}}$ and that calculated from the classical simulation force field $\Delta E_{\mathrm{MD}}$ (Figure \ref{fig3}). The classical solvation coordinate $\Delta E_{\mathrm{MD}}$ is also computed in similar fashion by using $\Delta E_{\mathrm{MD}} = U_{\mathrm{MD}}(r_f)-U_{\mathrm{MD}}(r_i)$. Although the classical forcefield is not designed to compute the correct full PES for O-H stretch, our calculation nevertheless depicts that these two variables are linearly correlated with a correlation coefficient of 0.83 (Figure \ref{fig3}). The fitting parameters are represented in Table \ref{table3}.
\begin{figure}
\includegraphics[width=0.8\linewidth]{figure/fig_3.eps}
\caption{\label{fig3}Correlation between solvation energy obtained from DFT calculation and MD calculations. Black color filled circles represents the bulk water data, whereas rests represent the data for methane-water cluster. The cyan solid line represents the linear fit.}
\end{figure}

\begin{table}[]
\centering
\caption{Linear empirical relation between solvation energy obtained from DFT calculation and MD calculations. We use $\Delta E_{\mathrm{MD}}=a_0+a_1 \Delta E_{\mathrm{DFT}}$ fitting function and the parameters and correlation coefficient are shown below.}
\begin{tabular}{lll}
\hline
$a_0$   & $a_1$    & Correlation Coefficient \\
\hline
\hline
11044.4	& 0.1326 &	0.83\\
\hline
\end{tabular}
\label{table3}
\end{table}

\subsection{Structure-spectrum correlations}
 To investigate the microscopic origin of different observed frequencies, we explore the structure-spectrum relationships. This was achieved by analyzing the static frequency distributions for distinct local environments sampled by a tagged O-H bond. We examine the effect of hydrophobicity on the transition frequency of O-H stretch by investigating the water molecule nearest to the methane moiety. Furthermore, we select the O-H bond of that selected water, which is pointing towards the methane. In Figure \ref{fig4}, we represent the static histogram distribution of fundamental transition frequency for bulk water and methane-water systems. The methane-water response is further decomposed into three sub-ensembles: $1^{st}$ solvation shell, $2^{nd}$ solvation shell and $3^{rd}$ solvation shell. The criteria for solvation shell decompositions are shown in supplementary information. It is evident from the Figure \ref{fig4} that the frequency distributions of closest water that are facing towards methane moiety and $1^{st}$ solvation shell water are significantly different than that of bulk water. The response of water molecules, which are directly facing the methane moiety, shows a two-ensemble blueshifted distribution (Figure \ref{fig4}(a)). Although the observed blueshift is apparent from the fact that the water molecules that are directly facing the methane moiety are lacking in number of hydrogen bonds; the two-ensemble picture needs further exploration. In case of $1^{st}$ solvation shell water, the distribution is redshifted with respect to the bulk water distribution. We indeed find that the distribution maximum shifted from bulk water response i.e. $\sim 3460 cm^{-1}$ to $\sim 3445 cm^{-1}$ in case of $1^{st}$ shell water (Figure \ref{fig4}(c)). This observation is similar to the earlier findings that the water in the first solvation shell of methane has enhanced structural ordering.\cite{Raschke2005,Grdadolnik2017,Koh2000} However, the water in the $2^{nd}$ and $3^{rd}$ solvation shell of methane almost behave like bulk water (Figure \ref{fig4}(a) and 5(b)).
 \par
\begin{figure*}[t!]
\includegraphics[width=.9\linewidth]{figure/fig_4.png}
\caption{\label{fig4}(a) Comparison of stretching frequency  of bulk water (black), closest water to methane (cyan), first solvation shell (red), second solvation shell (magenta) and third solvation shell (green). (b) Comparison of stretching frequency  of bulk water, second and third solvation shell water of methane-water system. (c) Comparison of stretching frequency  of bulk water and first solvation shell water of methane-water system.}
\end{figure*}

To understand the observed blueshift in the O-H stretch frequency of the water molecule nearest to the methane moiety, we perform bulk water simulation at T=373 K. In Figure \ref{fig5}(a) we represent the O-H frequency distributions of bulk water, water molecule closest to methane and water at 373 K. We indeed observe that the frequency distribution in the case of nearest water molecule is much more similar to that of water at T=373 K (Figure \ref{fig5}(a)). This suggests that the nearest water molecule behaves like low-density water, in fact it is almost similar to the dangling water in case of air-water interface.\cite{Du1994} This observation supports the suggestion made by Stillinger as well.\cite{Stillinger1973} It is clear that the O-H bond pointing towards methane moiety is experiencing extremely low hydrogen bond environment which is originating because of the alteration of hydrogen bond network by the hydrophobic methane moiety. The presence of methane moiety makes the neighboring environment similar to that of surface water.\cite{Stillinger1973}
\par
\begin{figure}
\includegraphics[width=.8\linewidth]{figure/fig_5.png}
\caption{\label{fig5} a) Comparison of stretching frequency $(\omega _{10})$ bulk water (black), closest water to methane (cyan) and water at 373 K (red). b) Stretching frequency $( \omega _{10} )$ of closest water to methane bimodal distribution divided based on angle.}
\end{figure}
However, this does not explain the two-ensemble nature of the distribution. To inspect that, we correlate the observed frequency with the $\angle \mathrm{O_WH_WH_{Me}}$ (Figure \ref{fig5}(b)). We find that the high frequency ensemble with larger population is coming from the configurations in which $\angle \mathrm{O_WH_WH_{Me}} > 120 \degree $ and the low frequency ensemble is originating from the configurations having $\angle \mathrm{O_WH_WH_{Me}} < 120 \degree $ (Figure \ref{fig5}(b)). Thus, the blueshift is increasing when the configurations are having more linear arrangement of the O-H bond of the nearest water molecule with the methane C—H bond. With increase in linearity, the propensity of interacting the O-H bond with the neighboring polar water molecules reduces, which results in higher blueshift in the observed transition frequency.\par
We have already discussed that the $1^{st}$ solvation shell water molecules (except the water molecules which directly faces the methane moiety) of methane show redshifted O-H stretch frequency distribution (Figure \ref{fig4}). To explore the observed redshift in the O-H stretch frequency of the $1^{st}$ shell water, we perform bulk water simulation at T=290 K, 280 K and 270 K. In Figure \ref{fig6}, we represent the O-H stretch frequency distribution in case of first shell water, bulk water at T=300 K, 290 K, 280 K and 270 K. It is evident from the figure that the first shell water behaves like more structured low temperature water.
\par
\begin{figure}
\includegraphics[width=.8\linewidth]{figure/fig_6.eps}
\caption{\label{fig6}Comparison of frequency distributions of first solvation shell water and bulk water at T=300K, 290K, 280K and 270K.}
\end{figure}

\begin{table}[]
\centering
\caption{Average hydrogen bond numbers per water molecule in bulk water, $1^{st}$ shell, $2^{nd}$ shell and $3^{rd}$ shell water of methane molecule.}
\begin{tabular}{ll}
\hline
Contributing water ensemble	& Average $n_{HB}$\\
\hline
\hline
Bulk water	& 3.55\\
1st shell	& 3.51\\
2nd shell	& 3.54\\
3rd shell	& 3.54\\
\hline
\end{tabular}
\label{table4}
\end{table}

In order to understand the above trends in the observed O-H stretch frequency, we further investigate the hydrogen bond distribution using standard geometrical criteria.\cite{Luzar1993} A water molecules is considered to form hydrogen bond with another water molecule if the inter-oxygen distance is less than 3.5 \text{\AA}, the hydrogen acceptor distance is less than 2.6 \text{\AA}, and H$_d$-O$_d$-O$_a$ angle is less than 30\degree, where the subscript ``$d$'' and ``$a$'' symbolize donor and acceptor respectively. The average hydrogen bond numbers in the different solvation shell water molecules are presented in Table \ref{table4}. We show the distribution of hydrogen bond number $(n_{HB})$ in each solvation shell along with the bulk water in Figure \ref{fig7}. It is evident from the figure is that the water molecules in the first solvation shell are having a much wider distribution spanning from $(n_{HB}) \sim 3.0$ to $(n_{HB}) \sim 4.0$. Although the average hydrogen bond number does not correlate with the observed frequency data, however the distribution makes a clear revelation of the microscopic picture. Thus, more structured water molecules are indeed present in the first solvation shell which eventually is reflected in the redshift of the observed O-H frequency. The extended low $(n_{HB})$ tail in the distribution of first shell water also supports the blueshift of the transition frequency in case of some of the first shell water molecules. On further inspection, we find that the later water ensemble mainly consist the water molecules, which are directly facing the methane moiety. In case of second and third shell water molecules, the hydrogen bond distributions are much narrower than that of the first shell water. These distributions are in fact slightly different than that of the bulk water. The observed frequency distributions in these two solvation shells also reflect that water in these two solvation shells behave almost similar like bulk water.
\begin{figure}
\includegraphics[width=0.8\linewidth]{figure/fig_7.eps}
\caption{\label{fig7}Comparison of hydrogen bond $(n_{HB})$ distribution in bulk water and different solvation shells of methane. All symbols represent the data and solid line represents the respective Gaussian fit.}
\end{figure}

\section{CONCLUSIONS}
In this work, we develop a semiempirical quantum-classical spectroscopy map to investigate the effects of hydrophobic solute in the water solvation structure. We find a strong correlation with the single oscillator level O-H stretch frequency with a collective variable solvation energy. Our model predicts that the water hydrogen bond network gets significantly modified in presence of small hydrophobic molecule such as methane. A bimodal frequency distribution with a blue shifted population of transition frequency illustrates presence of more like low-density water molecules and dangling water molecules in the hydration shell of methane. We also find that there exist a strong correlation between the transition frequency of nearest water O-H bond pointing towards methane moiety and the orientation of that O-H bond with that of  C-H  bond of methane. Temperature dependent simulation depicts that the water molecules facing to the methane molecule behave similar to the high temperature water. The solvation shell-wise decomposition of the O-H stretch frequencies further established that there exists a significant increase in ordering of the first solvation water except the water molecules, which are directly facing the methane molecule. This is manifested in redshift of the observed transition frequencies. Temperature dependent simulation also depicts that the water molecules in the first solvation shell except the water molecules, which are facing to the methane molecule behave more like cold water.

\section{SUPPLEMENTARY MATERIAL}
See supplementary information in Appendix \ref{sinfo} for complete information regarding cluster selection, correlation between electric field and transition frequency, spectroscopic maps and empirical relationships for other vibrational transitions and construction of different solvation shells of methane hydration.

\section{Acknowledgements}
RB acknowledges IIT Tirupati for support through new faculty seed grant and also the computational support.