\section{Introduction}
\label{sec:intro}
Contemporary cosmology consists of a few hypothesis and paradigms. The crucial one is the cosmological principle claiming that
the Universe is approximately homogeneous and isotropic on large scales, and assuming the matter and geometry of spacetime are
represented by smooth function of a position.
On the other hand, our Universe is fulfilled with the cosmic microwave background (CMB) radiation, as observed today, almost homogenous and isotropic.
The tiny fluctuations of CMB radiation are of order $10^{-5}$ \cite{smoot}.
Those deviations of CMB radiation appeared after the recombination epoch of our Universe when the photons decoupled from the electrons
making the Universe transparent. It constitutes radiation of free photons with a perfect blackbody spectrum \cite{ruth} which has been very well verified
observationally \cite{sta18}.
Measurements of the anisotropies in CMB radiation provide us information about the cosmological model describing the contents and evolution of the Universe.
The data delivered by them led to a standard cosmological model with only six parameters. The key CMB observables like temperature, polarization
anisotropies play the main role in the recent understanding of the initial conditions and structure formations inside the Universe.
The fluctuations on the large scale seen in the CMB spectrum allow us to determine the deviations from the
Friedmann-Robertson-Lemaitre-Walker (FRLW) geometry.

The question arises (the so-called horizon problem) why we observe the same background temperature in regions not being in causal contact before the CMB photons could freely
move through the spacetime. The evolution of the Universe provided by the $\Lambda$CDM model does not give enough time for creating so similar
patches with the same temperature. Another ones are the flatness problem and monopole problem (unwanted relics) - the latter is
predicted by models of particle physics giving that at very high temperatures there appear some artifacts which are not visible in our Universe presently.

The inflationary cosmology, being the standard paradigm, proposed as a solution to the horizon, flatness and monopole problems \cite{starob, guth},
delivers a very generous mechanism for generating primordial perturbations, which result in large-scale structures and temperatures
anisotropies of CMB radiation. The observed departures from cosmological homogeneity and isotropy observed in CMB radiation
are so tiny that one may conclude that they should be dominated by a Gaussian probability distribution.

The inflationary phase, characterized by the very rapid Universe expansion, is assumed to be dominated by a scalar field (the most popular is given by the so-called inflaton field \cite{LL}), which quantum fluctuations provide
density fluctuations reproducing the almost scale invariant spectrum of classical density perturbations \cite{komatsu}.

The cosmic of no-hair theorem \cite{gib77} claims that the inflationary Universe asymptotically tends to the de-Sitter spacetime till the end of the inflationary phase.
However, it was recently claimed that the effects of matter and spacetime inhomogeneities to inflation should be taken into account \cite{cho}.

 In the simplest form one deals with the definition of
power spectrum in Fourier space of the primordial density perturbation $\delta(\mathbf{k})$ in the following form\footnote{Later on we will assume
that the density perturbations $\delta(\mathbf{k})=\kappa\chi$.}
\begin{equation}
 \langle\delta(\mathbf{k})\delta(\mathbf{q})\rangle=P(k)(2\pi)^3\delta(\mathbf{k+q}),
\label{dd}
\end{equation}
or
\begin{equation}
 \langle\delta(\mathbf{x})\delta(\mathbf{y})\rangle=\int d^3k e^{i\mathbf{k\cdot(x-y)}} P(k),
\end{equation}
where we set $k=|\mathbf{k}|$ and $P$ is the power spectrum.

The form of the relation (\ref{dd}) is connected with rotational and translational invariance of the two-point correlation of the Fourier transform taken
from the primordial density $\delta( \mathbf{k})$. Namely,
the dependence of $P$ on the magnitude of the wave-vector $k$ informs
about the rotational invariance, while the delta function presence of $\mathbf{k+q}$ is a consequence of the translational
invariance. Now on, when we suppose that there could exist some objects during the inflationary era which violated the invariance, they should leave imprints on the
anisotropy of the CMB radiation caused by such a violation. Therefore, it was suggested \cite{ack} how to compute the power spectrum of fluctuations
in the case of a small violation of statistical isotropy which is characterized by a preferred direction in space along a unit vector $\mathbf{n}$. Thus, the primordial
power spectrum which takes into account the effects of the rotational invariance violation is now expressed (a parity $\mathbf{k\rightarrow-k}$ is assumed)
\begin{equation}
 P'(\mathbf{k})=P(k)\Big(1+g(\mathbf{\hat{k}\cdot n})^2\Big),
\end{equation}
where $\mathbf{\hat{k}}$ is the unit vector along the direction of $\mathbf{k}$ and $g$ is a $k$-independent constant. Further, the similar examination was done for a small
violation of translational invariance \cite{carrol}. Similarly, they assume that an object breaking the invariance during inflation disappears after the end of the
inflationary phase. They obtained expressions for the translational invariance broken by a special point, line or plane with the focus on the first one. The general form of
the density perturbation's two-point correlation breaking statistical translational invariance by a presence of a special point $\mathbf{x}_0$ preserving rotational
invariance about that point is provided by
\begin{equation}
 \langle\delta(\mathbf{x})\delta(\mathbf{y})\rangle=\int d^3k\int d^3q e^{i\mathbf{k\cdot(x-x_0)}} e^{i\mathbf{q\cdot(y-x_0)}} P_t(k,q,\mathbf{k\cdot q}),
\end{equation}
where $P_t$ is symmetric under interchange of $k$ and $q$. Consistent with the data, these violations are small resulting that $P_t$ is strongly peaked
about $\mathbf{k=-q}$. A special line, which one identifies via a point $\mathbf{x}_0$ and a unit tangent vector $\mathbf{n}$
(for simplicity one locates it along the $z$ axis) preserves the rotation invariance about the distinguished axis, while the translational invariance along this preferred
direction is left broken
\begin{equation}
  \langle\delta(\mathbf{x})\delta(\mathbf{y})\rangle=\int dk_z \int d k_\perp \int d^2q_\perp e^{ik_z(x_z-y_z)}e^{i\mathbf{k_\perp(\cdot x_\perp-x_{0\perp})}}
  e^{i\mathbf{q_\perp\cdot (y_\perp-z_{0\perp}}}P_t(k_\perp,q_\perp,k_z,\mathbf{k_\perp\cdot q_\perp}),
\end{equation}
where $P_t$ is symmetric under interchange of $k_\perp$ and $q_\perp$. The similar expression can be obtained for the plane \cite{carrol}. Those relations enable us
to compute directly the correlations between spherical harmonic coefficients of the CMB temperature field, which can be used for the comparison the result of the
CMB observations.

There are many candidates as the sources of anomalies, delivering small violation of rotational and translational invariance in the primordial
density fluctuations, which can be spotted in the CMB radiation pattern.
The impact of the preferred direction during the inflationary idea on the CMB anisotropy were studied in \cite{gum, arm, pereira, gum2, yoko, kan, wat}.
Quantum fluctuations of a free scalar field in the inflationary epoch were already investigated in a different context in the literature.
The point defect (Schwarzschild-de
Sitter black hole) were considered in \cite{cho}, while an infinitely long straight cosmic string was examined in \cite{string}.

The main objective of the work is to shed some light on the problem of a small violation of translational and rotational invariances during the early phase of
our Universe, which can be in principle imprinted on the deviation from the statistical isotropy of the cosmic microwave background radiation.

The paper is organized as follows: in the section \ref{sec2} we derive the two-point correlation function of the massless scalar
field in the expanding universe with a wormhole. The numerical results and our conclusions are given in the section \ref{sec3}.

\section{The two-point correlation function of a massless scalar field}\label{sec2}
The primordial quantum fluctuations occurring during the inflation phase are statistically homogeneous, being the consequence of translation
invariance. However, the presence of compact objects, like cosmological defects, black holes or wormholes can potentially
lead to the violation of the aforementioned pattern, which produces the imprint on CMB radiation, leading to its anisotropy.

Our main aim in this section will be to find the two-point correlation function of scalar massless field which mimics the dilaton field during the inflation epoch.
 We shall be interested in the spacetime of a wormhole in the expanding universe whose existence one supposes in the inflationary era of the Universe
evolution. In our calculations we use the so-called $\langle in - in \rangle$ formalism \cite{Weinberg}, in which the Hamiltonian of the system, in the
interactive picture, will be found in a perturbative way. The unperturbed spacetime will authorize the de Sitter manifold.

The line element describing the object in question is provided by \cite{mir06}
\begin{equation}\label{metric}
 ds^2=-dt^2+a^2(t)\left(1+\frac{b^m(r)}{r^m}\right)^\frac{4}{m}(dr^2+r^2d\theta+dz^2),
\end{equation}
where by $a(t)=e^{Ht}$ an usual scale factor is denoted, while
$H$ stands for the Hubble parameter.
The function
$b(r)$ is the so-called shape function, as it determines the spatial shape of the wormhole.
The $r$-coordinate decreases from infinity to the value $b_0$, being the size of the wormhole throat.
Then it increases from it back to the infinity (to another universe or
to some very distant place from ourselves).
Here, we are interested in the simplest wormhole, that is, when $b(r)=b=\text{const}$ and $m=1$. For
convenience we denote $\alpha=\left(1+\frac{b}{r}\right)$.

The other main ingredient in our consideration is massless scalar field which will mimic the inflaton field. Its action is given by the relation
\begin{equation}
S =- \frac{1}{2}\int d^4x\sqrt{-g} \nabla_\mu\chi \nabla^\mu\chi.
\end{equation}
In the spacetime of interests the Lagrangian density of the dilaton field implies
\begin{equation}
 \mathcal{L}=\frac{ra^3\alpha^6}{2}\left(\dot{\chi}^2-\frac{(\partial_r\chi)^2}{\alpha^4 a^2}-\frac{(\partial_\theta\chi)^2}{\alpha^{4}a^{2}r^{2}}-
 \frac{(\partial_z\chi)^2}{a^{2}\alpha^{4}}
 \right).
\end{equation}
On the other hand, the Legendre transformation enables one to find the Hamiltonian of the system
\begin{equation}
 \mathcal{H}=\int d^3x~
 \frac{ra^2\alpha^6}{2}\left(\dot{\chi}^2+\frac{(\partial_r\chi)^2}{\alpha^4 a^2}+\frac{(\partial_\theta\chi)^2}{\alpha^{4}a^{2}r^{2}}+
 \frac{(\partial_z\chi)^2}{a^{2}\alpha^{4}}
 \right).
\end{equation}
In order to implement the $\langle in - in \rangle$ formalism we expand in a Taylor series the Hamiltonian, near
$b \rightarrow 0$. The result can be written as
\begin{equation}
 \mathcal{H}=\int rdrd\theta dz~
 \frac{a^3r^2\dot{\chi}^2+a(\partial_\theta\chi)^2+ar^2(\partial_r\chi)^2+ar^2(\partial_z\chi)^2}{2r^2}+\mathcal{H}_\text{int},
\end{equation}
where the first term is an unperturbed Hamiltonian in the de Sitter spacetime, while
\begin{equation}\label{hami}
  \mathcal{H}_\text{int}=b\int d^3x \left( 3a^3\dot{\chi}^2+a(\partial_r\chi)^2+\frac{a}{r^2}(\partial_\theta\chi)^2+a(\partial_z\chi)^2
 \right).
\end{equation}
describes the interaction-picture Hamiltonian $\mathcal{H}_\text{int}(t)$, up to the first order in $b$ since in that case we have that
$[\mathcal{H}_0,\mathcal{H}_\text{int}]=0$.

Upon the quantization the inflaton field $\chi$ becomes a quantum operator
\begin{equation}
 \chi_I(\mathbf{x},\tau)=\int\frac{d^3k}{(2\pi)^3}e^{ik_zz+ik_\perp r\cos(\theta-\theta_k)}\Big( \chi_{\mathbf{k}}(\tau)\beta(\mathbf{k}) +\chi^*_{\mathbf{k}}(\tau)\beta^\dag(\mathbf{-k}) \Big)
\end{equation}
where $\tau=-\frac{1}{{H}}e^{-{H} t}$ is a conformal time.
The conformal time varies from $- \infty$ to $0$, while $t$-coordinate changes from $-\infty$ to $\infty$.
$\beta(\mathbf{k})$ annihilates the vacuum state satisfying the standard commutation
relations $[\beta(\mathbf{k}),\beta^\dag(\mathbf{-q})]=(2\pi)^3\delta(\mathbf{k}-\mathbf{q})$. Therefore, we rewrite the Hamiltonian (\ref{hami}) as
\begin{align}\label{int} \nonumber
 \mathcal{H}_\text{int}(\tau')&=b\int d^3x \frac{d^3k}{(2\pi)^3}\frac{d^3q}{(2\pi)^3}e^{i(\mathbf{k}+\mathbf{q})\mathbf{x}} \\ \nonumber
&{} \left\{ 3a^3(\tau')\Big( \dot{\chi}_{\mathbf{k}}(\tau)\beta(\mathbf{k}) +\dot{\chi}^*_{\mathbf{k}}(\tau)\beta^\dag(\mathbf{-k}) \Big)
\Big( \dot{\chi}_{\mathbf{q}}(\tau)\beta(\mathbf{q}) +\dot{\chi}^*_{\mathbf{q}}(\tau)\beta^\dag(\mathbf{-q}) \Big)\right.\nonumber\\
  &\left.+\left(-a(\tau')k_\perp q_\perp\cos(\theta-\theta_k)\cos(\theta-\theta_q)-a(\tau')k_zq_z-(k_xy-k_yx)(q_xy-q_yx)\frac{a(\tau')}{x^2+y^2}\right)\right.\nonumber\\
 &\left.\Big( \chi_{\mathbf{k}}(\tau)\beta(\mathbf{k}) +\chi^*_{\mathbf{k}}(\tau)\beta^\dag(\mathbf{-k}) \Big)
 \Big( \chi_{\mathbf{q}}(\tau)\beta(\mathbf{q}) +\chi^*_{\mathbf{q}}(\tau)\beta^\dag(\mathbf{-q}) \Big)
 \right\}.
\end{align}
In order to compute the two-point correlation function $\langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle$ of the inflaton, we use the formula
\begin{equation}\label{calka}
\langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle\simeq\langle\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)\rangle+
i\int^t_{-\infty}dt'e^{-\epsilon'|t'|}\langle[\mathcal{H}_\text{int}(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)] \rangle,
\end{equation}
where $\epsilon'$ is an infinitesimal parameter that cuts off the early time part of the integration. The relation (\ref{calka}) was originally introduced in
\cite{Weinberg}, the so-called $\langle in - in \rangle$ formalism applied to the calculation of higher order Gaussian and non-Gaussian corrections in cosmology.

To commence with, we calculate
the commutator appearing in the relation (\ref{calka}). Namely, one obtains
\begin{equation}
 [\mathcal{H}_\text{int}(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)]=[\mathcal{H}_1(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)]+
 [\mathcal{H}_2(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)],
\end{equation}
where we have denoted by $\mathcal{H}_1(t')$ the part including $\dot{\chi}$, while by $\mathcal{H}_2(t')$ the rest of (\ref{int}). It implies
\begin{align}\label{wyn1}
 &[\mathcal{H}_1(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)]\nonumber\\&=
 \frac{3}{2}bH\int\frac{d^3x'}{\tau'}\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{i\mathbf{x}'(\mathbf{k}+\mathbf{q})
 -i(\mathbf{kx}+\mathbf{qy})}\left((\tau+\frac{i}{q})(\tau+\frac{i}{k})e^{-i(\mathbf{q+k})(\tau'-\tau)}-\text{h.c.}
 \right),
\end{align}
and for the $\mathcal{H}_2(t')$-part, one deals with the following relation:
\begin{align} \label{above}
 &[\mathcal{H}_2(t'),\chi_I(\mathbf{x},t)\chi_I(\mathbf{y},t)] \\&=
 \frac{bH^3}{2}\int\frac{d^3x'}{\tau'}\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}\frac{e^{i\mathbf{x}'(\mathbf{k}+\mathbf{q})
 -i(\mathbf{kx}+\mathbf{qy})}}{kq}\Big[ \frac{(k_xy-k_yx)(q_xy-q_yx)}{x^2+y^2}+k_zq_z\nonumber\\
 &+k_\perp q_\perp \cos(\theta-\theta_k)\cos(\theta-\theta_q)
 \Big]\left((\tau'-\frac{i}{k})(\tau'-\frac{i}{q})
 (\tau+\frac{i}{k})(\tau+\frac{i}{q})e^{-i(\mathbf{q+k})(\tau'-\tau)}-\text{h.c.}
 \right), \nonumber
\end{align}
where one decomposes the vectors along $z$-axis and the two-dimensional subspace, perpendicular to it.

Let us notice that the first term in the square bracket in the equation (\ref{above}) has been already calculated in \cite{string}. The rest of
the formula in question, after
the change of the time coordinate $dt'=-\frac{d\tau'}{H\tau'}$ on the conformal time, and integrating over $\tau$, yields
\begin{align}\label{wyn2}
 \langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle_{2_\text{part}}
 &=-bH^2\int d^3x'\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{i\mathbf{x}'(\mathbf{k}+\mathbf{q})
 -i(\mathbf{kx}+\mathbf{qy})}\Big[k_zq_z
 +\frac{(\mathbf{k}_\perp\mathbf{x}'_\perp)(\mathbf{q}_\perp\mathbf{x}'_\perp)}{\mathbf{x}'^{2}_\perp}
 \Big]\nonumber\\
 &\times\left(\frac{kq+q^2+k^2+\tau k^2q^2}{k^3q^3(k+q)}
 \right).
\end{align}
On the other hand,
the result of the computations of the first term is given in \cite{string}, by the following:
\begin{align}
 \langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle_{2_\text{part2}}
 &=-bH^2\int d^3x'\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{i\mathbf{k}(\mathbf{x}'-\mathbf{x})
 +i\mathbf{q}(\mathbf{x}'-\mathbf{y})}\nonumber\\
 &\times\left(\frac{kq+q^2+k^2+\tau k^2q^2}{k^3q^3(k+q)}
 \right)\left(\frac{(y'k_x-x'k_y)(y'q_x-x'q_y)}{x'^2+y'^2}\right).
\end{align}

Similarly,  for  the relation (\ref{wyn1}), one obtains
\begin{align}\label{wyn1a}
  \langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle_{1}=2bH
  \int d^3x'\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{i\mathbf{x}'(\mathbf{k}+\mathbf{q})-i(\mathbf{kx}+\mathbf{qy})}
  \frac{kq\tau^2-1}{kq(k+q)}.
\end{align}
Taking into account the result presented in \cite{string} and ours given by the equations  (\ref{wyn2}) and (\ref{wyn1a}), finally we arrive at the formula
\begin{align}
 \Delta\langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle
 &=b\int d^3x'\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{i\mathbf{x}'(\mathbf{k}+\mathbf{q})-i(\mathbf{kx}+\mathbf{qy})}
 \left(3H\frac{kq\tau^2-1}{kq(k+q)}\right.\nonumber\\
 &\left.-H^2\mathbf{k} \cdot\mathbf{q}\Big( \frac{kq+q^2+k^2+k^2q^2\tau^2}{k^3q^3(k+q)} \Big)
 \right).
\end{align}
After integrating over $x'$, $ \Delta\langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle$ is provided by
\begin{align}\label{chi}
  \Delta\langle\chi(\mathbf{x},t)\chi(\mathbf{y},t)\rangle
 &=3Hb\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{-i(\mathbf{kx}+\mathbf{qy})}(2\pi)^\frac{3}{2}\delta(\mathbf{k+q})
 \frac{kq\tau^2-1}{kq(k+q)}\\
 &-bH^2\int \frac{d^3k}{(2\pi)^3}\int\frac{d^3q}{(2\pi)^3}e^{-i(\mathbf{kx}+\mathbf{qy})}(2\pi)^\frac{3}{2}
 \Big[ \frac{kq+q^2+k^2+k^2q^2\tau^2}{k^3q^3(k+q)} \Big]\mathbf{k\cdot q}. \nonumber
\end{align}
We notice that the first term is translational invariant because of the presence of the Dirac delta function while rotational invariance is broken since the power
spectrum depends on both vector fields magnitudes $k$ and $q$. The second term violates the translational invariance by the presence of
the special point $\mathbf{x}_0$ (which is taken at the origin) and preserves the rotational
invariance about this point \cite{carrol}.
Thus, if we approximate the density perturbations as $\delta=\kappa\chi$, the part of the power spectrum violating rotational and translational invariance,
which might be visible now ($\tau=0$) is written as
\begin{equation}\label{result}
 \Delta P(k,q,\mathbf{k\cdot q})= -3(2\pi)^\frac{3}{2}b\kappa^2H\left(\frac{\delta(\mathbf{k+q})}{kq(k+q)} +\frac{1}{3}
 \Big[ \frac{kq+q^2+k^2}{k^3q^3(k+q)} \Big]\mathbf{k\cdot q}
 \right).
\end{equation}
The relation (\ref{result}) constitutes the main result of our considerations. It will be discussed in what follows.

\section{Discussion and conclusions}\label{sec3}
We have computed the influence on the energy density perturbations of our Universe which is caused by the simple wormhole (\ref{metric}), with
the shape parameter
$b(r)=b=\text{constant}$ and
$m=1$,that could exist during the inflation epoch and have passed our horizon volume. Calculating of the two-point correlation function of the free massless scalar field (\ref{calka})
revealed that the effect of the wormhole on the perturbations of the energy density envisages in violating translational and rotational invariance.
It leads to the conclusion that the obtained result (\ref{result})
may be potentially  found in the anisotropy data of CMB.

Let us analyze the behavior of the fluctuation of the power spectrum. The formula (\ref{result}) in the general case depends on the mode wave
vectors $(\mathbf{k,q})$ which for convenience we express in Hubble parameter units. Let us notice that if the inflation duration is assumed to last about $60$
e-folds, the wavelength of the modes with $k/H=1$ is about the size of the Universe.
The dependence on the conformal time is suppressed in (\ref{result}), because of $\tau=0$ in the current evolutionary phase, but $\tau$-dependence can be easily obtained from
(\ref{chi}).
In all the plots presented here we took $\kappa = 10^{-3}$ and $b$ as unity, moreover the Dirac delta function will be approximated by the Lorentzian $\delta ( k ) \simeq \zeta / [ \pi ( k^{2} + \zeta^{2} ) ]$ with
broadening $\zeta = 0.0025$. The influence of the configuration of $\mathbf{(k,q)}$ modes on the function is as expected, that is,
we confirm for $\tau = 0$ that $\Delta P$ is large for small values of $k$ and $q$ due to the term with Dirac delta function
and it decreases rapidly with growth of the arguments. Therefore, the size of the initial peak is strongly dependent on a
spatial configuration of the vectors which is present in the scalar product.
To investigate it, without loss of generality, we set $\mathbf{k}$ to be oriented along the $x$-axes direction in the Fourier space, then we
describe $\mathbf{q}$ vector with spherical coordinate angles $(\theta, \phi)$. By rotating $\mathbf{q}$ we plot a map of values dependent
on the aforementioned angles (figure \ref{fig.modecenter}).

\begin{figure}[h!]
	\centering
	\includegraphics[width=0.80\textwidth]{map.pdf}
	\caption{
	\label{fig.modecenter}
	Map of $\Delta P(k,q,\mathbf{k\cdot q})$ showing the value of the function near initial peak ($k = q = 0.1$) with respect to different
	orientations of $\mathbf{q}$ described in terms of $(\theta, \phi)$ angles. The lowest value of power spectrum perturbation occurs
	for $\theta = 90^{\circ}, \phi = 180^{\circ}$ which corresponds to anti-parallel $\mathbf{k} = -\mathbf{q}$ setting.}
\end{figure}

We find out that the biggest perturbations occur for $\mathbf{k} = -\mathbf{q}$ setting which confirms the statement
from \cite{carrol}. Therefore, we focus on this configuration and we plot the values of $\Delta P$ with increasing value of $k$.

In figure  \ref{fig.k_t=0} we plot the power spectrum function for $\tau = 0$ versus the different values of the momenta $\bf k$. It can be seen that for the
increasingly values of $\bf k$, the function $\Delta P$ rapidly tends to its zero value. The only modes responsible for small values of the momenta,
influence the anisotropy spectrum observed nowadays.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.80\textwidth]{k_t=0.pdf}
	\caption{
	\label{fig.k_t=0}
	$\Delta P(k,q,\mathbf{k\cdot q})$ with $\tau = 0$ for different values of $k$. Function rapidly vanishes with
	increasing argument. Only the modes with small values of $k$ can have an influence on the anisotropy spectrum
	observed nowadays ($\tau = 0$).}
\end{figure}

Finally, in figure \ref{fig.dP_t},  we plot the time dependence of power spectrum perturbation for a few different e-folds. Consequently
we draw a conclusion that the longer inflation lasts, the weaker effects of the studied wormhole we observe. It can be envisaged that
$\Delta P$-function rapidly decreases to its final value after about three e-folds.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.80\textwidth]{P_t.pdf}
	\caption{
	\label{fig.dP_t}
	Power spectrum perturbation versus e-folding number connected with the conformal time $\tau = -\frac{e^{-Ht}}{H}$. Initial fluctuations decrease to their final values
	after approximately three e-folds.}
\end{figure}

One concludes that the observed CMB anisotropies deliver a new research window onto physics and phenomena that may occur during the rapid expansion phase of the
early Universe. The resulting theoretical claims can be compared with the observations, giving us the new opportunity to test
our models of  the early history of the Universe.

\acknowledgments
MR and AW were partially supported by the grant of the National Science Center \\
$DEC-2014/15/B/ST2/00089$. AW is also partially supported by FAPES (Brazil).


\section{Full Algorithm}
\label{sec:alg}

\begin{wrapfigure}{t}{0.55\linewidth}
\centering
\vspace{-0.6cm}
\input{diagram2.tex}
\caption{Overview of CoNAS. Train the one-shot model for the proxy measurements in stage 1. Then CoNAS approximates the one-shot model $f$ via compressive sensing. Stage 2 repeats until the minimizer of $g$ has enough information to construct the architecture.}
\end{wrapfigure}

We propose Compressive Recovering Neural Architecture Search (CoNAS), a new neural architecture search algorithm which approximates the one-shot architecture with Boolean functions via compressive sensing. We explain CoNAS according to three dimensions: search space, search strategy, and performance estimation strategy, suggested by NAS survey~\cite{elsken2018neural} to categorize the NAS methods. We define our search space used in our NAS algorithm in Section 3.1, succeeding the methodology of search space defined in \cite{li2019random}, \cite{liu2018darts}, and \cite{zoph2018learning}, but more considerable possible outcomes. Then we describe a fast proxy measuring technique by training the one-shot model through a weight-sharing by \cite{bender2018understanding} and \cite{li2019random} along with our new search space. Finally, we present our algorithm based on the compressive measuring scheme to recover the approximating function of the one-shot model with a sparse polynomial basis.

\subsection{Search Space}
\label{subsec:search space}

Following the recent NAS methodology as \cite{li2019random}, \cite{liu2018darts}, \cite{pham2018efficient}, \cite{real2018regularized}, and \cite{zoph2018learning}, our algorithm aims to find the best performing cell, a standard unit consists of the architecture, with less searching time compare to preexisting NAS algorithms. Zoph~\textit{et al.}~\cite{zoph2018learning} introduces a method of constructing the architecture by stacking and repeating the cell inspired by preexisting hand manually designed models such as Inception~\cite{szegedy2015going} and ResNet~\cite{he2016deep}. Previous NAS literature (\cite{li2019random}, \cite{liu2018darts}, \cite{pham2018efficient}, \cite{real2018regularized}, \cite{zoph2018learning}) has the following rules for cell construction:

\begin{itemize}
    \item Each cell has two input nodes and one output node. All intermediate nodes can only be connected by predecessor nodes including input nodes.
    \item Intermediate nodes are wired to two predecessor nodes in CNN and one predecessor node in RNN.
\end{itemize}

We propose an expanded search space using the Boolean architecture encoder denoted by $\alpha$. Starting from the continuous relaxation idea from DARTs~\cite{liu2018darts}, we define a direct acyclic graph (DAG) where all predecessor nodes fully connect every intermediate node with all possible operations. The architecture encoder $\alpha$, a binary string with a length of total number of edges in the DAG, contains information of the sub-graph from fully connected architecture. Each binary digit corresponds to the edge in the graph, and $1$ and $-1$ indicate its edge activated and inactivated respectively.

Figure~\ref{fig:diagram} describes how the architecture encoder $\alpha$ samples the sub-architecture of the fully-connected model in case of convolutional neural network. Our algorithm aims to find the best encoder $\alpha^*$ which is "close enough" to the global optimum returning the best validation accuracy by constructing the final model with $\alpha^*$ encoded sub-graph.

The new search space allows exploring a cell with more diverse connectivity patterns by being capable of switching on and off each edge independently. Moreover, the number of possible configuration from our search space exceeds the previous methods (\cite{li2019random}, \cite{liu2018darts}, \cite{pham2018efficient}, \cite{real2018regularized}, \cite{zoph2018learning}) with constrained wiring rules. For example, let the DAG for a convolutional neural network has $N=7$ nodes with five operations. Assuming we activate $\frac{1}{4}$ fractions from the total number of edges equivalent to $E=(2+3+4+5) \cdot 5 \cdot 2 = 140$, the possible combination of generating both \emph{normal cell} and \emph{reduce cell} is equal to $\binom{140}{35} \approx 1.2 \times 10^{33}$ where DARTs with five operations have $((5+1)^{14})^2 \approx 6.1 \times 10^{21}$

\subsection{Search Strategy}
\label{subsec:search strategy}

We propose a compressive measuring strategy to approximate the proxy model with Boolean function with Fourier expansion. Let function $f: \{-1, 1\}^n \rightarrow \mathcal{R}$ maps the sub-graph encoded binary string $\alpha$ to its performance (e.g. loss, accuracy) with given validation data. A similar approach to Cho~\textit{et al.}, Hazen~\textit{et al.}, and Stobbe~\textit{et al.} (\cite{stobbe2012learning}, \cite{hazan2017hyperparameter}, \cite{cho2019reducing}), we approximate the function $g \approx f$ with a unique sparse Fourier expansion via compressive sensing with randomly sampled measurements. Then we solve the $\argmin_{x}g(x)$ which is a highly non-convex optimization problem requiring sparsity $s$ be a small number to calculate all distinct points in a subcube. If the solution of the $\argmin_{x}g(x)$ does not return enough edges to construct the cell, we recover the approximate function $g$ repeating the sparse recovery with \emph{restriction} (Definition~\ref{def:restriction}).

We restrict the approximate function $g$ by reducing the domain with the solution found in the previous recovery. Then we randomly sample the sub-architectures while fixing the bit values in the previous sparse recovery solution. The multiple stages of sparse recovery enable to find the approximate function $g$ with extra basis contrasting to the single stage recovery.

\subsection{Performance Estimation Strategy}
\label{subsec:PES}

The conventional method of gathering the measurements to recover the sparse expression of the model involves the massive computational resource problem since each measurement requires training and evaluation independently. For the sake of better understanding, assume we want to train an architecture constructed from a randomly sampled cell by concatenating eight layers with 16 initial channels for a convolutional neural network (which is the same proxy model setup with DARTs). Suppose we measure the validation loss from the 50 epochs trained model where each epoch takes five minutes of training. The total sampling time in a unit of day to collect 4000 measurements (which is an equivalent number to our experiment), is equal to $5 \cdot 50 \cdot 4000 \cdot \frac{1}{60} \cdot \frac{1}{24} \approx 694$ days.

Our work is the extension of results from Bender et al.~\cite{bender2018understanding} and Li et al.~\cite{li2019random} suggesting the new direction of NAS algorithm not utilizing a controller but constructing a one-shot proxy model with weight-sharing. The experiment from Bender et al.~\cite{bender2018understanding} confirmed the strong correlation between the one-shot architecture and the stand-alone model build on its sub-architecture from the one-shot model. While Bender et al.~\cite{bender2018understanding} stabilizes the one-shot model training with a vraint of ghost batch normalization~\cite{hoffer2017train} and an extra hyperparameter called path dropout rate, Li et al.~\cite{li2019random} suggest the simpler training algorithm with regular backpropagation, only updating the weights of $\alpha_{w}$ correspond to the activated edges from randomly sampled sub-graph for each minibatch. We use the one-shot weight-shared model with the same methodology suggested by Li et al. to accelerate the gathering the measurements to recover the sparse approximation of $f$.

\begin{figure}[t]
\centering
    \input{diagram.tex}
    \caption{Diagram of how the architecture encoder $\alpha$ for the convolutional neural network samples the sub-architecture for $N=5$ nodes (two intermediate nodes). Each components in $\alpha$ maps to the edges one-to-one in all \emph{Choice} block in a \emph{Cell}. If the digit corresponds to $1$, the edge activates while $-1$ turns off its edge. Since CNN search space finds both \emph{normal cell} and \emph{reduce cell}, the length of $\alpha$ is equivalent to $(2+3) \cdot 5 \cdot 2 = 50$.}
    \label{fig:diagram}
\end{figure}

\subsection{Full Algorithm}

CoNAS has two main stages as shown in Algorithm~\ref{alg: CoNAS}. We first train the weight-shared model using the algorithm from Li.et.al~\cite{li2019random}, but applying different search space defined in Section~\ref{subsec:search space}. We randomly sample the sub-graph encoded string $\alpha \in \{-1, 1\}^n$ under $\textrm{Bernoulli}(p)$ for each digit in $\alpha$ and update the weight via the rule described in section~\ref{subsec:PES}.

Our algorithm has two main stages as shown in Algorithm \ref{alg: CoNAS}: the one-shot architecture training and sparse polynomial approximation via compressive sensing. We first train the one-shot architecture with the algorithm proposed by Li~\textit{et al.}. While RSWS followed the exact search space proposed by DARTs, we randomly sample the sub-architecture encoder $\alpha \in \{-1, 1\}^n$ under $Binomial(n, p)$ where $n$ is the total number of edges in the architecture and $p$ is the probability of being $1$ for each binary digits. In the second stage, we collect $m$ measurements of randomly sampled sub-architecture performance denoting $\mathbf{y} = (f(\alpha_1), f(\alpha_2), \ldots, f(\alpha_m))^T$. Then we construct the \emph{sampling matrix} $\mathbf{A} \in \{-1, 1\}^{m \times \sum_{l=0}^{d}\binom{n}{l}}$ with entries
$$
\mathbf{A}_{l, k} = \chi_k(\alpha_l), \:\:\:\:\:\:\: l \in [m], k \in [n], |k| \leq d
$$
where $d$ is the degree of Fourier series. The total number of basis, $\sum_{l=0}^{d}\binom{n}{l}$, is not too big by selecting $d=2$ (we assume that a Boolean function's Fourier spectrum is mostly concentrated at small degree). Then we solve the following convex optimization problem known as Lasso~\cite{tibshirani1996regression}, represented by the equation:

\begin{align}
\label{eq:Lasso}
    \mathbf{x}^* = \argmin_{\mathbf{x} \in \sum_{l=0}^{d}\binom{n}{l}} \|\mathbf{y} - \mathbf{A} \mathbf{x}\|_2^2 + \lambda\|\mathbf{x}\|_1
\end{align}

We define an approximate function $g \approx f$ with Fourier coefficients with top $s$ absolutely largest coefficients from $\mathbf{x}^*$. Since $g$ is the function with sparse polynomial expression, we compute the minimizer $\mathbf{\alpha}^* = \argmin_{\mathbf{\alpha}} g(\alpha)$ calculating all possible points in the subcube (the minimizer computation is feasible with small $s$ value). The multiple stages of sparse recovery of the function $f$ restricting to previous optimal $alpha^*$ enable to approximate the minimizer with extra monomial basis. After finding the best optimal $\alpha^*$ from CoNAS, we obtain the cell to construct the final architecture by only activating the edges corresponding to $\forall_{i \in [n]}\alpha_i^* = 1$.

\begin{algorithm}[t]
\caption{\textsc{CoNAS}}
\label{alg: CoNAS}
\begin{algorithmic}[1]
\State\textbf{Inputs:} Measurements $m$, stage $t$, sparsity $s$, lasso parameter $\lambda$
\While{not converged}
    \State Randomly sample a sub-architecture encoded binary vector $\alpha$
    \State Update weights $w_{\alpha}$ by descending $\nabla_{w_{\alpha}} \mathcal{L}_{train}(w_{\alpha})$
\EndWhile
\For{$i \in \{0, \ldots, t-1\}$}
\State Collect $\mathbf{y} = (f(\alpha_1), f(\alpha_2), \ldots, f(\alpha_m))^\top$.
\State Solve $$\mathbf{x}^* = \argmin_{\mathbf{x}} \|y - \mathbf{A} \mathbf{x}\|_2^2 + \lambda\|\mathbf{x}\|_1$$
\State Let $x_1, x_2, \ldots, x_s$ be the $s$ absolutely largest coefficients of $\mathbf{x}^*$ and $\chi_1, \chi_2, \ldots, \chi_s$ are the Boolean basis corresponding to the coefficient respectively. Let $g$ be a Fourier series with $s$ monomials:
$$g(\alpha) = \sum_{i=1}^{s} x_i \chi_i(\alpha)$$
\State Compute a minimizer $\mathbf{z} = \argmin_{\alpha} g(\alpha)$ and let $\overline{J}$ the set of indices where values in $\alpha^*$ are not zero.
\State $f = f_{J|\mathbf{z}}$
\EndFor
\State Construct the cell by activating the edge where $z_i = 1$ where $i \in [n]$.
\end{algorithmic}
\end{algorithm}
\section{Appendix}

\subsection{Detailed Background of One-Shot Neural Architecture Search}
\label{appendix:background one-shot}

In this section, we explain the one-shot neural architecture search in more details regarding of three directions introduced by \cite{elsken2018neural}.

\paragraph{Search Space.}
First, we start with the search space, defining the principles of constructing neural architecture. As discussed before, simplifying the search space using human prior knowledge can help search algorithms to find an optimal candidate faster. However, it may limit the algorithm to find a novel architecture beyond our knowledge due to human bias.

\paragraph{Search Strategy}
Next, we have a search strategy which is considered as a methodology to find the best neural architecture. In the last two years, different methodologies for search space such as reinforcement learning, evolutionary algorithms, SMBO, Bayesian optimization, bilevel optimization, and randomness are introduced in NAS literature. For instance, one-shot architecture uses random search as a search strategy.

\paragraph{Performance Estimation Strategy.}
The last direction aims to reduce the computational cost efficiently by estimating the validation performance. The conventional performance estimating methods require the massive computational resources as each model need to be trained and evaluated separately. In particular, assume that our goal is to construct a CNN by using a proxy architecture through randomly selected cell samples by concatenating eight layers with 16 initial channels (This is the same proxy model setup in DARTs). Then, we measure the validation loss from 50 epochs trained model where each epoch takes five minutes of training. As a result, the total sampling time to collect 4,000 randomly sampled of the architecture performance (which is the same number of measurements as our experiment using CoNAS) equals to $5 \cdot 50 \cdot 4000 \cdot \frac{1}{60} \cdot \frac{1}{24} \approx 694$ days.

\subsection{Theoretical Support for CoNAS}
\label{sec:theoretical support}

A Boolean function $f$ defined over a binary $n$-bit input may expresses using $2^n$ polynomial basis functions (which is intractable). As discussed in~\cite{hazan2017hyperparameter}, instead, we can focus on all the polynomial basis with degree at most $d$, that is, we have $|\mathcal{P}_d| \equiv \mathcal{O}(n^d)$ basis (parity) functions). On the other hand, the system of linear equations $\mathbf{y} = \mathbf{Au}$ with the graph-sampling matrix $\mathbf{A} \in \{-1, 1\}^{m \times O(n^d)}$, measurements $\mathbf{y} \in \mathbb{R}^m$, and Fourier coefficient $\mathbf{u} \in \mathbb{R}^{O(n^d)}$ is an ill-posed problem when $m \ll O(n^d)$ for large $n$. While the underdetermined system of linear equations is an ill-posed, it has been shown that if the graph-sampling matrix satisfies \emph{Restricted Isometry Property (RIP)}, the sparse coefficients, $\mathbf{u}$ can be recovered.

We recall the definition of RIP:
\begin{definition}
A matrix $\mathbf{A}\in\mathbb{R}^{m\times \mathcal{O}(n^d)}$ satisfies the restricted isometry property of order $s$ with some constant $\delta$ if for every $s$-sparse vector $\mathbf{u}\in\mathbb{R}^{\mathcal{O}(n^d)}$ (i.e., only $s$ entries are non-zero) the following holds:
$$(1-\delta) \|\mathbf{u}\|_2^2\leq\|\mathbf{Au}\|_2^2\leq(1+\delta) \|\mathbf{u}\|_2^2.$$
\end{definition}

Now, we assume that the Fourier support of $f$ is $s$-sparse. In addition, in our case, the columns of graph-sampling matrix $\mathbf{A}$ (the atoms of dictionary $\mathbf{A}$) include the parity function. Hence, the columns of matrix $A$ form a bounded orthonormal system (BOS), i.e., $\sup_{\bm{\alpha} \in \{-1, 1\}^n} |\chi_S(\bm{\alpha})| \leq 1 \:\:\: \textrm{for all }S \in [n]$.

There has been significant research during the last decade in proving upper bounds on the number of rows of bounded orthonormal dictionaries (matrix $\mathbf{A}$) for which $\mathbf{A}$ is guaranteed to satisfy the restricted isometry property with high probability. One of the first BOS results was established by~\cite{candes2006near}, where the authors proved an upper bound scales as $\mathcal{O}(sd^6\log^6 n)$ for a subsampled Fourier matrix. While this result is seminal, it is only optimal up to some \textit{polylog} factors. In fact, the authors in chapter $12$ of \cite{foucart2017mathematical} have shown a necessary condition (lower bound) on the number of rows of BOS which scales as $\mathcal{O}(sd\log n)$. In an attempt to achieve to this lower bound, the result in~\cite{candes2006near} was further improved by~\cite{rudelson2008sparse} to $\mathcal{O}(sd\log^2s\log(sd\log n)\log n)$. Motivated by this result, \cite{cheraghchi2013restricted} has even reduced the gap further by proving an upper bound on the number of rows as $\mathcal{O}(sd\log^3 s\log n)$.  The best known available upper bound on the number of rows appears to be $\mathcal{O}(sd^2\log s\log^2n)$; however with worse dependency on the constant $\delta$, i.e., $\delta^{-4}$ (please see~\cite{bourgain2014improved}). To the best of our knowledge, the best known result with mild dependency on $\delta$ ($\delta^{-2}$) is due to~\cite{haviv2017restricted}, which we can invoke for our setup. In our case, the graph-sampling matrix $\mathbf{A}$ in our proposed CoNAS algorithm satisfies BOS for $K=1$ (equation~\ref{eq:sampling matrix}).

\begin{theorem}
\label{thm:recovery}
Let graph-sampling matrix $A\in\{-1,1\}^{m\times \mathcal{O}(n^d)}$ be constructed by taking $m$ rows (random sampling points) uniformly and independently from the rows of a square matrix $\mathbf{M}\in\{-1,1\}^{\mathcal{O}(n^d)\times \mathcal{O}(n^d)}$. Then normalized matrix $\mathbf{A}$ (multiplied by $\sqrt{\frac{\mathcal{O}(n^d)}{m}}$) with $m = \mathcal{O}(\log^2(\frac{1}{\delta})\delta^{-2}s\log^2(\frac{s}{\delta})d\log n)$ with probability at least $1-2^{-\Omega(d\log n\log(\frac{s}{\delta}))}$ satisfies the restricted isometry property of order $s$ with constant $\delta$; as a result, every $s$-sparse vector $\mathbf{u}\in \mathbb{R}^{\mathcal{O}(n^d)}$ is recovered from the samples
$$
\mathbf{y} = \mathbf{Au} = \big(\sum_{j=1}^{|\mathcal{O}(n^d)|} u_j\mathbf{A}_{i,j}\big)_{i=1}^{m}
$$
by LASSO (equation~\ref{eq:Lasso}). It is worthwhile to mention that we can use any other sparse recovery method such as IHT~\cite{blumensath2009iterative}, or OMP~\cite{tropp2007signal} in our algorithm.

\end{theorem}

Therefore, Theorem~\ref{thm:recovery} leads to successfully guarantee approximation of the function using the Boolean Fourier basis given the sufficient number of measurements.

\subsection{Effect of Multiple Stages Sampling}

Now we talk briefly about the advantage of using multiple stages in CoNAS. As mentioned in the algorithm pseudocode, after each stage of recovery, the algorithm uses multiple stages of recovery process to find the final architecture. This is similar to the idea of de-biasing step in \emph{Hard-Thresholding (HT)} algorithm~\cite{foucart2017mathematical} introduced in the sparse recovery literature where first the support is estimated and then within the estimated support the sparse coefficients are calculated through a least-square computation. However in our setup, we have to solve a non-linear optimization problem to get better estimation of the sparse Fourier coefficients. For instance, in Table~\ref{table:stage sampling}, we have illustrated an experiment which shows that the sparse recovery process with multiple stages can discover the important edges; as a result to find a more accurate  architecture. In this experiment, we have sampled $1000$ measurements and computed the mean and the standard deviation of validation loss/perplexity for each stage. As we can see, the larger number of stages, the smaller loss/perplexity.

\begin{table}[!ht]
  \caption{Changes in Measurements for Each Stage.}
  \label{table:stage sampling}
  \centering
  \vspace{0.3cm}
  \begin{tabular}{c c c c c}
    \toprule
    \multicolumn{1}{l}{Architecture Type}     & Stage 1    & Stage 2  & Stage 3 & Stage 4 \\
    \midrule
    \multicolumn{1}{l}{CNN (Valid loss)} & $8.21 \pm 26.76$ & $3.01\pm 8.48$ & $2.48 \pm 3.55$ & $2.37 \pm 1.13$ \\
    \multicolumn{1}{l}{RNN (Valid perplexity)} & $340.65 \pm 239.84$ & $207.30\pm 56.74$ & - & - \\
    \bottomrule
  \end{tabular}
\end{table}

\newpage
\subsection{Architecture Found from Our Experiment}

\subsubsection{Compressive sensing-based Neural Architecture Search (CoNAS) for RNN}
\label{appendix:RNN CoNAS}
\begin{figure}[!ht]
    \centering
    \includegraphics[width=1.0\linewidth]{figures/CoNAS_recurrent.pdf}
    \vspace{5mm}
    \caption{Recurrent Cell found from CoNAS}
    \label{fig:rnn cell}
\end{figure}

\subsubsection{Random Search with Weight-Sharing (RSWS) for CNN}

\begin{figure}[!ht]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=0.8\linewidth]{figures/RSWS_normal.pdf}} \\
    \subfloat[Reduce Cell]{\includegraphics[width=0.8\linewidth]{figures/RSWS_reduce.pdf}}
    \vspace{5mm}
    \caption{Convolutional cell found from RSWS with search space in Section~\ref{subsec:search space CoNAS}}
    \label{fig:RSWS cnn cell}
\end{figure}

\newpage
\subsubsection{Differentiable Neural Architecture Search (DARTs) for CNN}

\begin{figure}[!ht]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=0.9\linewidth]{figures/DARTS_normal.pdf}} \\
    \subfloat[Reduce Cell]{\includegraphics[width=0.9\linewidth]{figures/DARTS_reduce.pdf}}
    \vspace{5mm}
    \caption{Convolutional Cell found from DARTs with the original setting in \cite{liu2018darts}.}
\end{figure}

\newpage
\subsubsection{Randomly Wired CNN Architectures}
\label{appendix:random wired CNN}
\begin{figure}[!ht]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=1.0\linewidth]{figures/random1_normal.pdf}} \\
    \subfloat[Reduce Cell]{\includegraphics[width=1.0\linewidth]{figures/random1_reduce.pdf}}
    \vspace{5mm}
    \caption{Randomly wired cell corresponding to the first row result to Table~\ref{table:randomly wired model}.}
\end{figure}

\newpage
\begin{figure}[!ht]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=1.0\linewidth]{figures/random2_normal.pdf}} \\
    \subfloat[Reduce Cell]{\includegraphics[width=1.0\linewidth]{figures/random2_reduce.pdf}}
    \vspace{5mm}
    \caption{Randomly wired cell corresponding to the second row result to Table~\ref{table:randomly wired model}.}
\end{figure}

\newpage
\begin{figure}[!ht]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=1.0\linewidth]{figures/random3_normal.pdf}} \\
    \subfloat[Reduce Cell]{\includegraphics[width=1.0\linewidth]{figures/random3_reduce.pdf}}
    \vspace{5mm}
    \caption{Randomly wired cell corresponding to the third row result to Table~\ref{table:randomly wired model}.}
\end{figure}

\pagebreak
\newpage
\subsubsection{Random Search with Weight-Sharing (RSWS) for RNN}
\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.7\linewidth]{figures/RSWS_recurrent.pdf}
    \vspace{5mm}
    \caption{Recurrent cell found from RSWS allocating equal amount of search time to CoNAS}
\end{figure}

\section{Background}

\paragraph{One-shot NAS}
\label{sec:background of one-shot}
We briefly describe one-shot neural architecture search techniques \cite{bender2018understanding,li2019random}; a lengthier description is available in Appendix~\ref{appendix:background one-shot}. Following the recent NAS survey paper~\cite{elsken2018neural}, one-shot NAS approaches have three main components: a search space, a search strategy, and a performance estimation strategy.

\paragraph{Search Space.}
\label{sub:search space}
The goal of one-shot NAS is to find the best performing \emph{cell}, a fundamental component from which more complex architectures are constructed via stacking.
Following \cite{liu2018darts}, a cell is a directed acyclic graph (DAG) where a node corresponds to the latent representation, and a directed edge transforms predecessor nodes using a given \emph{operation}; common operations used in CNNs include $3 \times 3$ and $5 \times 5$ separable convolutions, $3 \times 3$ max pooling, $3 \times 3$ average pooling, and Identity. Each cell has two input nodes and one output node, and intermediate nodes can only be connected by predecessor nodes including input nodes. Intermediate nodes are wired to two predecessor nodes in CNNs and one predecessor node in RNNs.

\paragraph{Search Strategy and Performance Estimation Strategy.}
\label{sec:search strategy and PES}

Having defined a search space, one-shot NAS approaches employ four steps: (i) train a single ``one-shot'' base model that is capable of predicting the performance of sub-architectures; (ii) \emph{randomly sample} sub-architectures of a trained one-shot model and measure performance over a hold-out validation set of samples; (iii) select the candidate (cell) with best validation performance. (iv) retrain a deeper final architecture using the best cell. Using the one-shot model as the proxy measurements of the candidate architecture corresponds to the performance estimation strategy (first and second step). The search strategy of one-shot architecture search on both \cite{bender2018understanding} and \cite{li2019random} is equivalent to random search.

\section{Background on Boolean Function with Fourier Expansion }
\label{sec:background}

We first mention some basic concepts in Fourier analysis of Boolean functions~\cite{o2014analysis}. The \textit{real-valued Boolean Function} is the function which maps n-bit binary vector (or hypercube) to a real value, $$f : \{-1,1\}^n \rightarrow \mathbb{R}$$ which can be represented in a real multilinear polynomials known as Fourier Basis. The Fourier Basis (or monomial) corresponding to any subsets $S \subseteq [n]$ is defined as following:

\begin{definition}[Fourier Basis~\cite{o2014analysis}]
\label{def:fourier basis}
For $S \subseteq [n]$, the Fourier basis is the parity function $\chi_S: \{-1, 1\}^n \rightarrow \{-1, 1\}$ where $$\chi_S(x) = \prod_{i \in S} x_i$$
\end{definition}

The function $f$ has the unique Fourier expansion with the series of a real multilinear polynomial basis (or Fourier Basis) given by:

\begin{align}
\label{eq:Fourier Series}
f(x) = \sum_{S \subseteq [n]} \hat{f}(S) \chi_S(x)
\end{align}
where
\begin{align}
\label{eq:Fourier Coefficients}
\hat{f}(S) = \mathbb{E}_{x\in\{-1,1\}^n}[f(x)\chi_S(x)]
\end{align}
where $x$ is a random binary string taken with respect to the uniform distribution over the nodes of the n-dimensional hypercube. The essential fact of the Fourier expansion of Boolean function is that the basis of parity function forms an $K$-bounded orthonormal system (BOS) with $K=1$ satisfying two properties:
\begin{align}
\label{eq:BOS}
\langle \chi_S, \chi_T \rangle =
\begin{cases}
    1, & \text{if } S = T\\
    0, & \text{if } S \neq T
\end{cases}
\end{align}
and
\begin{align}
\label{eq:BOS2}
    \sup_{x \in \{-1, 1\}^n} |\chi_S(x)| \leq 1 \:\:\:\: \textrm{for all }S \in [n]
\end{align}

Lastly, Boolean function allows restricting the hypercube domain to subcube with the operation called \emph{restriction}. Restriction operation allows not only to reduce the domain size but also the subfunction of $f$ with Fourier expression if we have prior knowledge of some fixed digits indexed by set $J$:

\begin{definition}[Restriction~\cite{o2014analysis}]
\label{def:restriction}
Let $f: \{-1, 1\}^n \rightarrow \mathbb{R}$, $(J, \overline{J})$ be a partition of $[n]$, and $z \in \{-1, 1\}^{\overline{J}}$. The restriction of $f$ to $J$ using $z$ denoted by $f_{J|z}: \{-1, 1\}^J \rightarrow \mathbb{R}$ is the subfunction of $f$ given by fixing the coordinates in $\overline{J}$ to the bit values $z$.
\end{definition}
\section{Conclusion}
\label{sec:conclusion}

We presented CoNAS, a blend of weight-sharing one-shot models with compressive sensing scheme to approximate the proxy model in a multilinear polynomial expression. By merging two different schemes to NAS, CoNAS is able to outperform the state-of-the-art architecture search methods on image classification and match on language modeling tasks remarkably shortening the search time compare to preexisting NAS methodologies.

Our work also shows how the search space can affect performance while adopting the same search strategy. We have expanded the search space in terms of the number of possible connections to the cell and interesting direction to improve not only CoNAS but also for other preexisting NAS algorithms by expanding the search space by expanding the designing rules on existing search strategies.

\section{Experiments and Results}
\label{sec:result}
We experimented on two different NAS problems: a CNN search for CIFAR-10, and an RNN search for Penn Treebank (PTB). We describe the training details for CIFAR-10 and PTB in Sections~\ref{exp:cifar} and~\ref{exp:PTB} respectively. Our experimental setup for building the final architecture is the same as that reported in DARTS and RSWS.

\subsection{CIFAR-10}
\label{exp:cifar}

\begin{figure}[t]
    \centering
    \subfloat[Normal Cell]{\includegraphics[width=0.93\linewidth]{figures/CoNAS_Normal.pdf}} \\
    \vspace{-5mm}
    \subfloat[Reduce Cell]{\includegraphics[width=0.90\linewidth]{figures/CoNAS_Reduce.pdf}}
    \vspace{5mm}
    \caption{Convolution Cell found from CoNAS}
    \label{fig:cnn cell}
\end{figure}

\paragraph{Architecture Search.} We create a one-shot architecture, similar to RSWS, with a cell containing $N=7$ nodes with two nodes as input and one node as output; our wiring rules between nodes are different and as in Section~\ref{subsec:search space CoNAS}. We used five operations: $3 \times 3$ and $5 \times 5$ separable convolutions, $3 \times 3$ max pooling, $3 \times 3$ average pooling, and Identity.
On CIFAR-10, we divide the 50,000 training set to 25,000 training set and 25,000 validation set as the equivalent experiment setup in \cite{li2019random} and \cite{liu2018darts}. We trained a one-shot model with eight layers and 16 initial channels for $150$ epochs.
All other hyperparameters used in training the one-shot model are the same as in RSWS.

\begin{table}[t]
    \centering
    \begin{threeparttable}
    \caption{\textbf{CIFAR-10 Benchmark: Comparison with hand-designed networks and state-of-the-art NAS methods} (Lower test error is better). The results are grouped as follows: manually designed networks, published NAS algorithms, and our experimental results. The average test error of our experiment used five random seeds. Table entries with "-" indicates that either the field is not applicable or unknown. The methods listed in this table are trained with auxiliary towers and cutout augmentation. The search cost from our experiment is based on NVIDIA Titan X GPU.}
    \label{table: CIFAR10 Comparison Table}
    \begin{tabular}{c c c c c c}
        \toprule
        & \multicolumn{2}{c}{\textbf{Test Error}} & \textbf{Params} & \textbf{Search Cost} & \textbf{Search}\\
        \multicolumn{1}{l}{\textbf{Architecture}} & Best & Average & (M) & GPU days & \textbf{Method}\\ [1.0ex]
        \midrule
        \multicolumn{1}{l}{Shake-Shake~\cite{devries2017improved}} & - & $2.56\pm0.07$ & 26.2 & - & manual\\
        \multicolumn{1}{l}{PyramidNet~\cite{yamada2018shakedrop}} & $2.31$ & - &26 & - & manual\\[0.5ex]
        \midrule
        \multicolumn{1}{l}{ProxylessNAS\textsuperscript{+}~\cite{cai2018proxylessnas}} & $2.08$ & - & 5.7 & 4 & gradient \\
        \multicolumn{1}{l}{NASNet-A~\cite{zoph2018learning}} & - & $2.65$ & 3.3 & 2000 & RL \\
        \multicolumn{1}{l}{AmoebaNet-B~\cite{real2018regularized}} & - & $2.55\pm0.05$ & 2.8 & 3150 & evolution \\
        \multicolumn{1}{l}{GHN\textsuperscript{+}~\cite{zhang2018graph}} & - & $2.84 \pm 0.07$ & 5.7 & 0.84 & hypernetwork \\
        \multicolumn{1}{l}{SNAS~\cite{xie2018snas}} & - & $2.85 \pm 0.02$ & 2.8 & 1.5 & gradient \\
        \multicolumn{1}{l}{ENAS~\cite{pham2018efficient}}& 2.89 & - & 4.6 & 0.45 & RL \\
        \multicolumn{1}{l}{Random search~\cite{liu2018darts}} & - & $3.29 \pm 0.15$ & 3.1 & 4 & random \\
        \multicolumn{1}{l}{DARTs (first order)~\cite{liu2018darts}} & - & $3.00 \pm 0.14$ & 3.3 & 1.5 & gradient \\
        \multicolumn{1}{l}{DARTs (second order)~\cite{liu2018darts}} & - & $2.76\pm0.09$ & 3.3 & 4 & gradient \\
        \multicolumn{1}{l}{ASHA~\cite{li2019random}} & 2.85 & $3.03\pm0.13$ & 2.2 & - & random \\
        \multicolumn{1}{l}{RSWS~\cite{li2019random}} & 2.71 & $2.85\pm0.08$ & 4.3 & 2.7 & random \\
        \multicolumn{1}{l}{DARTs\textsuperscript{\#} (second order)~\cite{li2019random}} &2.62 &$2.78 \pm 0.12$ & 3.3 & 4 & gradient \\
        \midrule
        \multicolumn{1}{l}{DARTs\textsuperscript{\dag} (second order)} & 2.59& $2.78\pm0.13$ & 3.4 & 4 & gradient \\
        \multicolumn{1}{l}{RSWS\textsuperscript{*}} & 2.47 & $2.59\pm0.11$&4.7 & 0.7 & random \\
        \multicolumn{1}{l}{CoNAS} & 2.55 & $2.62\pm0.06$ & 4.8 & 0.7 & CS \\
        \bottomrule
    \end{tabular}
    \begin{tablenotes}
    \small
    \item \textsuperscript{\#} DARTS experimental results from~\cite{li2019random}.
    \item \textsuperscript{\dag} Used DARTS search space with five operations for direct comparisons.
    \item \textsuperscript{*} Used the search space defined in Section~\ref{subsec:search space CoNAS}. Used five operations for direct comparison.
    \item \textsuperscript{+} The search space is not comparable to CoNAS.
    \vspace{-4mm}
    \end{tablenotes}
    \end{threeparttable}

\end{table}

We run CoNAS with multiple sparse recovery stages. Specifically, following \cite{hazan2017hyperparameter}, we used the parameters $s=10$, Fourier basis degree $d=2$, and Lasso coefficient $\lambda=1$ for each stage.
Repeating three to four stages ($t=3,4$) of sparse recovery with restriction returned an architecture encoder $\bm{\alpha^*}$ with sufficient number of edges needed to construct both \emph{normal cell} and \emph{reduce cell} (e.g., see Figure~\ref{fig:cnn cell}). The final architecture found from CoNAS used $t=4$.

\paragraph{Architecture Evaluation.}
\label{sec: architecture evaluation}
We re-trained the final architecture with the learned cell with the equivalent hyper-parameter configurations in DARTS and RSWS to make the direct comparisons. We used NVIDIA Titan X for the training process. Since DARTS and RSWS used somewhat different GPUs (GTX1080-Ti and Tesla P100 respectively), we included not only the experimental results that have been published, but also our results of DARTs and RSWS on a Titan X to compare the search time with equivalent GPU environment in Table~\ref{table: CIFAR10 Comparison Table}. Due to hardware memory constraints, the learned cell from CoNAS could not use the same minibatch size (96) used in DARTs and RSWS; instead, we re-trained the final model with minibatch size (56) \footnote{In our experience, a larger minibatch size enabled training the model quicker than minibatches of smaller size.}.

\begin{table}[t]
  \caption{\textbf{Randomly Wired Model Performance on CIFAR-10}. (Lower test error is better) Trained with auxiliary towers and cutout augmentation.}
  \label{table:randomly wired model}
  \centering
  \begin{tabular}{c c c c}
    \toprule
    Method     & Number of Edges    & Size (M)  & Test Error \\
    \midrule
    Randomly Wired Model (1) & 31 & 5.2 & $3.45\%$  \\
    Randomly Wired Model (2) & 35 & 4.5 & $2.89\%$  \\
    Randomly Wired Model (3) & 32 & 3.1 & $3.52\%$  \\
    \bottomrule
    \vspace{-7mm}
  \end{tabular}
\end{table}

\paragraph{Effect of Increasing Connections in the Cell.} We provide a further supplementary experiment of training results of randomly wired cells which have a similar number of edges to the cell found by CoNAS. We randomly sampled the architecture code $\bm{\alpha}$ with $Binomial(140, 1/4)$ for each digit being $1$ that the expected number of total edges in \emph{normal cell} and \emph{reduce cell} is 35. We randomly selected three architectures as shown in Appendix~\ref{appendix:random wired CNN} and trained with the equivalent environment to the section~\ref{sec: architecture evaluation}. Table~\ref{table:randomly wired model} shows that choosing a larger number of randomly chosen edges is not sufficient to improve the model's performance.

\subsection{Penn Treebank}
\label{exp:PTB}

\paragraph{Architecture Search.} Similar to the setup of DARTS and RSWS, CoNAS explores the cell with the following operations: Tanh, ReLU, Sigmoid, and Identity. We augment the RNN cell with a variation of highway connections suggested by Pham et al~\cite{pham2018efficient}. Layers of depth $l$ in Recurrent Highway Network~\cite{zilly2017recurrent} utilize a nonlinear transformation from its hidden state $h_l$ as follows:
$$
h_l = (1-c_l) \otimes \texttt{activation}(h_{l-1} \cdot W_{l, l-1}^{(h)}) + h_{l-1} \otimes c_l
$$
where "$\otimes$" denotes element-wise multiplication. Since we use an expanded search space, we allow for mulitple operations; in such cases, we replace $\texttt{activation}(h_{l-1} \cdot W_{l, l-1}^{(h)})$ with its sum-pooled version, $\frac{1}{n} \sum_{i}^{n} \texttt{activation\textsuperscript{i}}(h_{l-1} \cdot W_{l, l-1}^{(h)})$

Pre-training the above RNN using weight-sharing showed unstable results since the sub-graphs could have some internal nodes with no connections, leading to exploding gradients; increasing the $p$-parameter in the Bernoulli sampling to enforce connectivity significantly slowed down computations. Hence, we added the additional heuristic of randomly activating an edge to connect the intermediate node if the node does not have any input edge according to its architecture encoder $\alpha$.

After obtaining the one-shot model, we randomly sample the measurements of the sub-graph without the above heuristic used in the training stage. While running CoNAS, two stages of sparse recovery with $s=10$ and $s=5$ respectively found the enough number of edges for the RNN cell. If the final resulting cells has intermediate nodes with a disconnected input, we added ReLU operations from the previous intermediate node. The visualization of the RNN cell found by CoNAS is shown in Appendix~\ref{appendix:RNN CoNAS}.

\begin{table}[t]
    \centering
    \begin{threeparttable}
        \caption{\textbf{PTB Benchmark: Comparision of the state-of-the-art NAS methods and hand-designed networks} (Lower perplexity is better). The results are grouped in following orders: manually designed networks, published NAS algorithms, and our experimental results. The average test error of our experiment used five random seeds. Table entries with "-" indicates either the field is not applicable or unknown. The search cost from our experiment is based on a NVIDIA Titan X GPUs.}
        \label{table: PTB Comparison Table}
        \begin{tabular}{c c c c c c c c}
            \toprule
            & \multicolumn{2}{c}{\textbf{Test Perplexity}} & \textbf{Params} & \textbf{Search Cost} & \textbf{Search}\\
            \multicolumn{1}{l}{\textbf{Architecture}} & Valid & Test & (M) & GPU days & \textbf{Method}\\ [1.0ex]
            \midrule
            \multicolumn{1}{l}{Variational RHN~\cite{zilly2017recurrent}} & 67.9 & 65.4 & 23 & - & manual \\
            \multicolumn{1}{l}{LSTM + DropConnect~\cite{merity2017regularizing}} & 60.0 & 57.3 & 24 & -& manual\\
            \multicolumn{1}{l}{LSTM + Mos~\cite{yang2017breaking}} & 56.5 & 54.4 & 22 & - & manual \\
            \midrule
            \multicolumn{1}{l}{NAS~\cite{zoph2016neural}} & - & 64.0 & 25 & 1e4 & RL \\
            \multicolumn{1}{l}{ENAS\textsuperscript{\dag}~\cite{pham2018efficient}} & - & 56.3 & 24 & 0.5 & RL \\
            \multicolumn{1}{l}{Random search\textsuperscript{\dag}~\cite{liu2018darts}} & 61.8 & 59.4 & 23 & 2 & random \\
            \multicolumn{1}{l}{DARTS (1st order)\textsuperscript{\dag}~\cite{liu2018darts}} & 60.2 & 57.6 & 23 & 0.5 & gradient \\
            \multicolumn{1}{l}{DARTS (2nd order)\textsuperscript{\dag}~\cite{liu2018darts}} & 58.1 & 55.7 & 23 & 1 & gradient \\
            \multicolumn{1}{l}{ASHA\textsuperscript{*}~\cite{li2019random}} & 58.6 & 56.4 & 23 & - & random \\
            \multicolumn{1}{l}{RSWS\textsuperscript{*}~\cite{li2019random}} & 57.8 & 55.5 & 23 & 0.25 & random \\
            \midrule
            \multicolumn{1}{l}{RSWS\textsuperscript{\#, +}} & 60.6 & 57.9 & 23 & 0.25 & random \\
            \multicolumn{1}{l}{CoNAS\textsuperscript{+}} & 59.1 & 56.8 & 23 & 0.25 & CS \\
            \bottomrule
        \end{tabular}
        \begin{tablenotes}
        \small
        \item \textsuperscript{\#} We used the RSWS code with adjusting the search time equivalent to CoNAS.
        \item \textsuperscript{\dag} Used NVIDIA GTX 1080Ti GPU for training/searching.
        \item \textsuperscript{*} Used Tesla P100 GPU for training/searching.
        \item \textsuperscript{+} Used NVIDIA Titan X GPU for training/searching.
        \vspace{-5mm}
        \end{tablenotes}
    \end{threeparttable}
\end{table}

\paragraph{Architecture Evaluation.}
\label{sec: architecture evaluation-PTB}
The PTB results for recurrent architectures are presented in Table~\ref{table: PTB Comparison Table}. We trained the final RNN model with the learned cell with the equivalent hyperparameters in DARTs and RSWS, except with minibatch size = 128 (due to hardware constraints). We also included the experimental results with RSWS methods allocating the equivalent search time with our methods to make the direct comparison with CoNAS. Since the published NAS literature in Table~\ref{table: PTB Comparison Table} uses different GPU hardware (e.g. DARTs and ENAS: NVIDIA GTX 1080Ti, RSWS: Tesla P100), a one-to-one comparison of the search cost value listed in Table~\ref{table: PTB Comparison Table} is not applicable. For example, both CoNAS and RSWS required 0.25 GPU-days, but the latter used a more powerful GPU. We will provide a more precise one-to-one comparison if this manuscript passes peer review.

\subsection{Discussion}

Noticeably, CoNAS achieved improved results on CIFAR-10 in both test error and search cost when compared to the previous state-of-the-art algorithms: DARTs, RSWS, and ENAS. We see that RSWS with an equivalent search space to CoNAS obtained a better architecture (with best test error $2.47\%$); however, the standard deviation of test errors over different pseudorandom was greater than CoNAS. Many previous NAS papers have focused on the search strategy while adopted the same search space tol~\cite{zoph2018learning} and~\cite{liu2018darts}. Our experimental results highlight the importance of considering not only seeking new performance strategies, but also the search space.

On PTB, CoNAS found a better RNN architecture than RSWS using an equivalent search cost. In the overall comparison, DARTs (second order) and RSWS outperformed both valid and test perplexity of CoNAS. We could not include true one-to-one comparisons with other algorithms since they used different search spaces.

\section{Proposed Algorithm: CoNAS}
\label{sec:full algorithm}

\paragraph{Overview.}
Our proposed algorithm, Compressive sensing-based Neural architecture Search (CoNAS), infuses ideas from Boolean Fourier analysis into one-shot NAS.  CoNAS consists of two novel components: an expanded search space, and a more effective search strategy.

\paragraph{Search Space.}
\label{subsec:search space CoNAS}
Following the approach of DARTS~\cite{liu2018darts}, we define a directed acyclic graph (DAG) where all predecessor nodes are connected to every intermediate node with all possible operations. We represent any sub-graph of the DAG using a binary string $\bm{\alpha}$ called the \emph{architecture encoder}. Its length is the total number of edges in the DAG, and a $1$ (resp. $-1$) in $\bm{\alpha}$ indicates an active (resp. inactive) edge.

\begin{figure}[!t]
\centering
    \input{diagram.tex}
    \caption{Diagram adapted from~\cite{bender2018understanding}. The architecture encoder $\bm{\alpha}$ samples the sub-architecture for $N=5$ nodes (two intermediate nodes) with five different operations. Each component in $\bm{\alpha}$ maps to the edges one-to-one in all \emph{Choice} blocks in a cell. If a bit in $\bm{\alpha}$ corresponds to $1$, the edge activates, while $-1$ turns off the edge. Since the CNN search space finds both \emph{normal cell} and \emph{reduce cell}, the length of $\bm{\alpha}$ is equivalent to $(2+3) \cdot 5 \cdot 2 = 50$.}
    \label{fig:diagram}
\end{figure}

Figure~\ref{fig:diagram} gives an example of how the architecture encoder $\bm{\alpha}$ samples the sub-architecture of the fully-connected model in case of a convolutional neural network. The goal of CoNAS is to find the ``best' encoder $\bm{\alpha^*}$, which is "close enough" to the global optimum returning the best validation accuracy by constructing the final model with $\bm{\alpha^*}$ encoded sub-graph.

Since each edge can be switched on and off independently, the proposed search space allows exploring a cell with more diverse connectivity patterns than DARTS~\cite{liu2018darts}. Moreover, the number of possible configurations
 exceeds similar previously proposed search spaces with constrained wiring rules \cite{li2019random,pham2018efficient,real2018regularized,zoph2018learning}. For example, suppose the DAG for a convolutional neural network model has $N=7$ nodes with five operations. The total number of edges is equivalent to $E=(2+3+4+5) \cdot 5 \cdot 2 = 140$; if $1/4$ of the edges are activated, then the number of possible combinations

 is equal to $\binom{140}{35} \approx 1.2 \times 10^{33}$, while DARTs with five operations only supports $((5+1)^{14})^2 \approx 6.1 \times 10^{21}$ configurations.

\paragraph{Search Strategy.}
\label{subsec:search strategy CONAS}
We propose a compressive measuring strategy to approximate the one-shot model with a Fourier-sparse Boolean function. Let $f: \{-1, 1\}^n \rightarrow \mathbb{R}$ map the sub-graph of the one-shot pre-trained model encoded by $\bm{\alpha}$ to its validation performance.  Similar to \cite{hazan2017hyperparameter},
we collect a small number of function evaluations of $f$, and reconstruct the Fourier-sparse function $g \approx f$ via compressive sensing with randomly sampled measurements.
Then, we solve $\argmin_{\bm{\alpha}}g(\bm{\alpha})$ by exhaustive enumeration over all coordinates in its sub-cube $\{-1, 1\}^{\overline{J}}$ where $(J, \overline{J})$ partitions $[n]$.
If the solution of the $\argmin_{\bm{\alpha}}g(\bm{\alpha})$ does not return enough edges to construct the cell (some intermediate nodes are disconnected), we simply perform sparse recovery multiple times, following \cite{hazan2017hyperparameter}. In each iteration, we restrict the approximate function $g$ by by fixing the bit values found in the previous solution, and randomly sample sub-graphs in the remaining edges.

\newcommand{\var}[1]{\text{\texttt{#1}}}
\newcommand{\func}[1]{\text{\textsl{#1}}}

\makeatletter
\newcounter{phase}[algorithm]
\newlength{\phaserulewidth}
\newcommand{\setlength{\phaserulewidth}}{\setlength{\phaserulewidth}}
\newcommand{\phase}[1]{%
  \vspace{-1.25ex}

  \Statex\leavevmode\llap{\rule{\dimexpr\labelwidth+\labelsep}{\phaserulewidth}}\rule{\linewidth}{\phaserulewidth}
  \Statex\strut\refstepcounter{phase}\textbf{Stage~\thephase~--~#1

  \vspace{-1.25ex}\Statex\leavevmode\llap{\rule{\dimexpr\labelwidth+\labelsep}{\phaserulewidth}}\rule{\linewidth}{\phaserulewidth}}
\makeatother

\setlength{\phaserulewidth}{.7pt}

\newlength{\whilewidth}
\settowidth{\whilewidth}{\algorithmicwhile\ }
\algdef{SE}[parWHILE]{parWhile}{EndparWhile}[1]
  {\parbox[t]{\dimexpr\linewidth-\algmargin}{%
     \hangindent\whilewidth\strut\algorithmicwhile\ #1\ \algorithmicdo\strut}}{\algorithmicend\ \algorithmicwhile}%
\algnewcommand{\parState}[1]{\State%
  \parbox[t]{\dimexpr\linewidth-\algmargin}{\strut #1\strut}}

\begin{algorithm}[t]
  \caption{\textsc{CoNAS}}
  \label{alg: CoNAS}
  \begin{algorithmic}[1]
      \State\textbf{Inputs:} Number of one-shot measurements $m$, stage $t$, sparsity $s$, lasso parameter $\lambda$
    \phase{Training the One-Shot Model}
        \Procedure{Model Training}{}
        \While{not converged}
            \State Randomly sample a sub-architecture encoded binary vector $\alpha$
            \State Update weights $w_{\bm{\alpha}}$ by descending $\nabla_{w_{{\bm{\alpha}}}} \mathcal{L}_{train}(w_{{\bm{\alpha}}})$
        \EndWhile
        \EndProcedure
    \phase{Search Strategy}
    \Procedure{One-Shot Model Approximation via Compressive Sensing}{}
        \For{$k \in \{1, \ldots, t\}$}
            \State Collect $\mathbf{y} = (f(\bm{\alpha}_1), f(\bm{\alpha}_2), \ldots, f(\bm{\alpha}_m))^\top$.
            \State Solve $$\mathbf{x}^* = \argmin_{\mathbf{x}} \|\mathbf{y} - \mathbf{A} \mathbf{x}\|_2^2 + \lambda\|\mathbf{x}\|_1$$
            \State Let $x_1^*, x_2^*, \ldots, x_s^*$ be the $s$ absolutely largest coefficients of $\mathbf{x}^*$.  Construct

            $$g(\bm{\alpha}) = \sum_{i=1}^{s} x_i^* \chi_i(\bm{\alpha})$$
            \State Compute minimizer $\mathbf{z} = \argmin_{\bm{\alpha}} g(\bm{\alpha})$ and let $\overline{J}$ the set of indices of $z$.
            \State $f = f_{J|\mathbf{z}}$
        \EndFor
        \State Construct the cell by activating the edge where $z_i = 1$ where $i \in [n]$.
    \EndProcedure
  \end{algorithmic}
\end{algorithm}

\paragraph{Full Algorithm.}
\label{subsec:full algorithm CoNAS}
We now describe CoNAS in detail, with pseudocode shown in Algorithm~\ref{alg: CoNAS}.
We first train a one-shot model with standard backpropagation but only updates the weight corresponding to randomly sampled sub-graph edges for each minibatch.
Then, we randomly sample sub-graphs by generating architecture encoder strings $\bm{\alpha} \in \{-1, 1\}^n$ using a $Bernoulli(p)$ distribution for each bit of $\bm{\alpha}$.

In the second stage, we collect $m$ measurements of randomly sampled sub-architecture performance denoting $\mathbf{y} = (f(\bm{\alpha}_1), f(\bm{\alpha}_2), \ldots, f(\bm{\alpha}_m))^T$.  Next, we construct the \emph{graph-sampling matrix} $\mathbf{A} \in \{-1, 1\}^{m \times |\mathcal{P}_d|}$ with entries
\begin{align}
\label{eq:sampling matrix}
\mathbf{A}_{l, k} = \chi_S(\bm{\alpha}_l), \:\:\:\:\:\:\: l \in [m], k \in [|\mathcal{P}_d|], S \subseteq [n], |S| \leq d
\end{align}
where $d$ is the max degree of monomials in the Fourier expansion. We solve the familiar Lasso problem~\cite{tibshirani1996regression}:
\begin{align}
\label{eq:Lasso}
    \mathbf{x}^* = \argmin_{\mathbf{x} \in \mathbb{R}^{|\mathcal{P}_d|}} \|\mathbf{y} - \mathbf{A} \mathbf{x}\|_2^2 + \lambda\|\mathbf{x}\|_1
\end{align}
to (approximately) recover the global optimizer $\mathbf{x}^*$, the vector containing the Fourier coefficients corresponding to $\mathcal{P}_d$. We define an approximate function $g \approx f$ with Fourier coefficients with the top-$s$ (absolutely) largest coefficients from $\mathbf{x}^*$. We compute the minimizer $\bm{\alpha}^* = \argmin_{\bm{\alpha}} g(\bm{\alpha})$ calculating all possible points in the subcube defined by the support of $g$ (the minimizer computation is feasible if $s$ is small). Multiple stages of sparse recovery (with successive restrictions to previously obtained optimal $\bm{\alpha}^*$) enable to approximate additional monomial terms.  Finally, we obtain the cell to construct the final architecture by activating the edges corresponding to all ${i \in [n]}$ such that $\alpha_i^* = 1$. We include theoretical support for CoNAS in Appendix~\ref{sec:theoretical support}.

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

\textbf{Motivation.} Choosing a suitable neural network architecture for complex prediction tasks such as image classification and language modeling often requires a substantial effort of trial-and-error. Therefore, there has been a growing interest to \emph{automatically learn} (or meta-learn) the architecture of neural networks that can achieve competitive (or better) results over hand-designed architectures. The sub-field of neural architecture search (NAS) addresses the problem of designing competitive architectures with as small a computational budget as possible.

\textbf{Our Contributions.} Numerous approaches for neural architecture search already exist in the literature, each with their own pros and cons: these include {black-box optimization} based on reinforcement learning (RL)~\cite{zoph2016neural}, evolutionary search~\cite{real2018regularized}, and Bayesian optimization~\cite{cao2019learnable,kandasamy2018neural}. Though the algorithmic details vary, most of these NAS methods face the common challenge of evaluating the test/validation performance of a (combinatorially) large number of candidate architecture evaluations.
In a departure from traditional methods, we approach the NAS problem via the lens of \emph{compressive sensing}. The field of compressive sensing (or sparse recovery), introduced by the seminal works of \cite{candes2004robust,donoho2006compressed}, has received significant attention in both ML theory and applications over the last decade, and has influenced the development of numerous advances in nonlinear and combinatorial optimization.

We leverage these advances for the NAS problem. We develop a new NAS method called CoNAS (Compressive sensing-based Neural Architecture Search), which merges ideas from sparse recovery with so-called ``one-shot'' architecture search methods~\cite{bender2018understanding}, described in greater detail below. CoNAS consists of two new innovations: (i) a new \emph{search space} that permits exploration of a large(r) number of diverse candidate architectures, and (ii) a new \emph{search strategy} that borrows ideas from recovery of Boolean functions from their (sparse) Fourier expansions.

Our experiments show that CoNAS is able to design a deep convolutional neural network with test error $2.62 \pm 0.06\%$
on CIFAR-10 classification, outperforming existing state-of-the-art methods such as DARTs~\cite{liu2018darts}, ENAS~\cite{pham2018efficient}, and random search with weight-sharing (RSWS)~\cite{li2019random}. Moreover, CoNAS performance is competitive to NASNet~\cite{zoph2018learning} and AmoebaNet~\cite{real2018regularized},  despite requiring less than one GPU-day of computation. Our experiments on designing recurrent neural networks for language modeling are somewhat short of the state-of-the-art~\cite{zilly2017recurrent}, but we find that CoNAS still finds competitive results with lesser search time than previous NAS approaches. Our results are exactly reproducible (having been trained with fixed pseudorandom seeds), and an implementation of CoNAS will be made publicly available post-peer review.

\begin{figure}
\centering
\input{diagram2.tex}
\medskip
\caption{Overview of CoNAS. A one-shot neural network model $f$ is pre-trained, and an appropriate sub-graph of $f$ is chosen by iteratively applying sparse recovery techniques.
}
\end{figure}

Finally, while our original motivation was to devise a empirically useful NAS method, a nice benefit is that CoNAS can be \emph{theoretically analyzed}, since existing theoretical results for Fourier-sparse Boolean functions can be ported over in order to provide upper bounds for the number of one-shot evaluations required. This, to our knowledge, is one of the first results of their kind in the NAS literature and may be of independent interest. We defer this to the supplementary material.

\textbf{Techniques.} The intuition behind compressive sensing is that if a signal (or function) can be represented via a sparse basis expansion, then it can be recovered (either exactly or approximately) from a small number of randomized measurements. CoNAS leverages this intuition in the context of one-shot architecture search~\cite{bender2018understanding}. In one-shot NAS, instead of evaluating several candidate architectures, a single ``base'' neural network model is pre-trained; a class of sub-networks is identified (called the search space) and the performance of each sub-network is evaluated on a validation set; and the best-performing sub-network is finally selected and fine-tuned.

We model the sub-network selection as a sparse recovery problem. More concretely, consider a function $f$ that maps sub-architectures to a measure of performance (validation loss). We assume that $f$ can be written as a sparse, low-degree polynomial in the (discrete) Fourier basis\footnote{Intuitively, this means that $f$ is well-approximated by a decision tree; for a formal proof, see~\cite{hazan2017hyperparameter}.}. If the sparsity assumption is satisfied, then we claim the function $f$ can be reconstructed using a very small number of measurements (evaluations), thus reducing overall compute time. A key challenge lies in defining a suitable search space; we propose one that is considerably larger than the one used in DARTS or ENAS, allowing us to (putatively) search over a more diverse set of candidate architectures.

\textbf{Prior work}. Due to tight space limits, our review of prior work will be unfortunately brief; we refer to the recent paper by~\cite{li2019random} for a thorough treatment. Early NAS approaches used RL-based controllers~\cite{zoph2018learning}) or evolutionary algorithms \cite{real2018regularized}, and showed competitive performance with manually-designed architectures such as deep ResNets~\cite{he2016deep} and DenseNets~\cite{huang2017densely}. However, these approaches required substantial computational resources, running into thousands of GPU-days.  Subsequent NAS works have focus on boosting search speeds by proposing novel search strategies, such as sequential model-based optimization (SMBO) techniques~\cite{liu2018progressive}, Bayesian optimization \cite{cao2019learnable,kandasamy2018neural}) and Monte Carlo tree search (MCTS)~\cite{negrinho2017deeparchitect}. Other recent NAS approaches include weight-sharing~\cite{pham2018efficient}, hypernetworks~\cite{brock2017smash}, one-shot models~\cite{bender2018understanding}, network transformations \cite{cai2018efficient,elsken2018efficient,jaderberg2017population,jin2018efficient,liu2017hierarchical}, gradient descent \cite{cai2018proxylessnas,liu2018darts,luo2018neural,xie2018snas}, and random exploration \cite{li2018massively,li2019random,sciuto2019evaluating,xie2019exploring}.

\section*{Acknowledgements}

The authors would like to thank Amitangshu Mukherjee, Soumik Sarkar, and Alberto Speranzon for helpful discussions. This work was supported in part by NSF grants CCF-1566281, CAREER CCF-1750920, a GPU gift grant from the NVIDIA Corporation, and a faculty fellowship from the Black and Veatch Foundation.

{{
\bibliographystyle{unsrt}

\section{Theoretical Support for CoNAS}
\label{sec:theoretical support}

Boolean function $f$ with binary $n$-bit input may have up to $2^n$ polynomial basis. The system of linear equations $y = Ax$ with the sampling matrix $A \in \{-1, 1\}^{m \times 2^n}$, measurements $y \in \mathbb{R}^m$, and Fourier coefficient $x \in \mathbb{R}^{2^n}$ is the ill-posed problem not only due to its underdetermined system since $m \ll 2^n$ for large $n$, but also its intractable number of basis. Instead, we assume that the Fourier support of $f$ is $s$-sparse and Boolean Fourier basis constructs the sampling matrix $A$ satisfying BOS (equation~\ref{eq:BOS}) and \emph{Restricted Isometry Property (RIP)}. Stobbe et al.~\cite{stobbe2012learning} applies the Theorem 4.4 (recovery condition of sparse polynomials in general BOS) of Rauhut~\cite{rauhut2010compressive} that the sampling matrix with Boolean Fourier basis satisfies RIP condition with high probability. We restate the Theorem in Rauhut~\cite{rauhut2010compressive} and Stobbe et al.~\cite{stobbe2012learning} with our notational setup.

\begin{theorem}
\label{thm:recovery}
Let $A \in \{-1, 1\}^{m \times p}$ be the sampling matrix (eqation~\ref{eq:sampling matrix}) that satisfies BOS for $K=1$. Assume that the random sampling points $t_1, \ldots, t_m$ are chosen independently and uniformly at random. Then there exist universal constants $C_1,C_2>0$ such that $m \geq \max{(C_1ns\ln^2{(s)}, C_2s\ln{(\delta^{-1})})}$, except for an event of probability no more than $\delta$, the following holds for all $f \in\mathcal{H}_d$. Then with probability at least $1-\delta$ every $s$-sparse vector $x \in \mathbb{R}^{|\mathcal{P}_d|}$ is recovered from the samples
$$
\mathbf{y} = \mathbf{Ax} = \big(\sum_{j=1}^{|\mathcal{P}_d|} x_j\mathbf{A}_{i,j}\big)_{i=1}^{m}
$$
by Lasso (equation~\ref{eq:Lasso}).

For any noise level $\eta \geq 0$ and any noisy vector of measurements $y \in \mathbb{R}^m$ within that noise level: suppose $\mathbf{x}^*$ is the solution of
$$
\min_{\mathbf{z} \in\{-1, 1\}^{|\mathcal{P}_d|}}\|\mathbf{z}\|_1\:\:\:\:subject\:to\:\:\|\mathbf{y} - \mathbf{Az}\| \leq \eta
$$
Then the following bound holds for some universal constants $C_3, C_4$ and Lasso solution with noise, $\mathbf{x}^*$:
$$
\|\mathbf{x}-\mathbf{x}^*\| \leq \frac{C_1 \sigma_s(\mathbf{x})_1}{\sqrt{s}} + C_4\eta
$$
where $\sigma(\cdot)_1$ is defined as the $\ell_1$ error of the best $s$-sparse approximation.
$$
\sigma_s(\mathbf{x})_1 := \inf\{\|\mathbf{x}-\mathbf{z}\|_1, \mathbf{z} \textrm{ is s-sparse}\}
$$
\end{theorem}

Therefore, Theorem~\ref{thm:recovery} guarantees successful approximation of the function using the Boolean Fourier basis given sufficient amount of measurements.