text:
\begin{document}
\title[Monge-Amp\`ere equation]{ Global Monge-Amp\'ere equation with asymptotically periodic data} \subjclass{35J96; 35J67}
\keywords{Monge-Ampere equation, Dirichlet problem, a priori estimate, maximum principle, viscosity solution}
\author{Eduardo V. Teixeira}
\address{Universidade Federal do Cear\'a\\ Departamento de
Matem\'atica\\ Av. Humberto Monte, s/n, Campus do Pici - Bloco 914
\\
Fortaleza-CE, Brazil.  CEP 60.455-760} \email{teixeira@mat.ufc.br}

\author{Lei Zhang}
\address{Department of Mathematics\\
        University of Florida\\
        358 Little Hall P.O.Box 118105\\
        Gainesville FL 32611-8105}
\email{leizhang@ufl.edu}

\date{\today}

\begin{abstract} Let $u$ be a convex solution to $\det(D^2u)=f$ in $\mathbb R^n$ where $f\in C^{1,\alpha}(\mathbb R^n)$ is asymptotically close to a periodic function $f_p$. We prove that the difference between $u$ and a parabola is asymptotically close to a periodic function at infinity, for dimension $n\ge 3$.
\end{abstract}

\maketitle

\section{Introduction}

In this article we study convex, entire viscosity solutions $u \colon \mathbb R^n \to \mathbb{R}$, to the Monge-Amp\'ere equation
\begin{equation}\label{mp-g}
	\det(D^2u)=f(x),\quad \mbox{in }\quad \mathbb R^n.
\end{equation}
The forcing term $f$ is assumed to be positive and asymptotically close to a periodic function at infinity. Our main goal is to establish a classification theorem for such solutions.

\par

Monge-Amp\'ere equation with periodic data can be found in various topics in applied mathematics such as homogenization theory, optimal transportation problems, vorticity arrays, etc. Equation (\ref{mp-g}) also appears in differential geometry, when it is lifted from a Hessian manifold \cite{caf-via}. In spite of the profusion in application, Monge-Amp\'ere equation is well known for its analytical difficulty and it is no exception for equation (\ref{mp-g}) when the right hand side is close to a periodic function. In \cite{caf-li-2} Caffarelli and Li proved that if $f$ is a positive periodic function, $u$ has to be a parabola plus a periodic function with the same periodicity of $f$. This theorem can be viewed as an extension to the classification theorems of J\"orgens \cite{jorgens}, Calabi \cite{calabi}, Pogorelov \cite{p5}, Caffarelli-Li \cite{caf-li-1} for the Monge-Amp\`ere equation $det(D^2u)=1$.

The aim of this article is to establish an optimal perturbation result from the Caffarelli-Li's classification theorem, as to cover forcing terms $f$ that are asymptotically a periodic function at infinity.

In more precise terms, we let
\begin{equation}\label{aa}
 \mathcal{A} \quad \mbox{ be the set of all positive definite and symmetric matrices,}
 \end{equation}
  the assumption on $f$ is as follows:
Let $f_p$ be a positive, $C^{1,\alpha}$ periodic function in $\mathbb R^n$, i.e.:
\begin{align}\label{fp}
 &\exists d_0>0 , \alpha\in (0,1), a_1,...,a_n>0 \mbox{ such that }\\
 & d_0^{-1}\le f_p\le d_0, \quad
\|f_p\|_{C^{1,\alpha}(\mathbb R^n)}\le d_0, \nonumber\\
 & f_p(x+a_i e_i)=f_p(x), \quad \forall x\in \mathbb R^n. \nonumber
\end{align}
where $e_1=(1,0,..,0)$,..., $e_n=(0,...,0,1)$. We assume that
$f\in C^{1,\alpha}(\mathbb R^n)$ is asymptotically close to $f_p$ in the following sense:
\begin{align}\label{af}
&\exists d_1>0 \mbox{ and } \beta>2,\mbox{such that } \\
&d_1^{-1}\le f(x)\le d_1, \,\forall x\in \mathbb R^n, \quad \|f\|_{C^{1,\alpha}(\mathbb R^n)}\le d_1,
\nonumber\\
&|\nabla^j(f^{\frac 1n}-f_p^{\frac 1n})(x)|\le d_1(1+|x|)^{-\beta-j},\quad \forall x\in \mathbb R^n,\quad j=0,1,2,3.\nonumber
\end{align}
\begin{rem} $f,f_p\in C^{1,\alpha}(\mathbb R^n)$ but the difference between $f^{\frac 1n}$ and $f_p^{\frac 1n}$ is more smooth.
\end{rem}
Under the above framework, our main theorem is:
\begin{thm}\label{thm1}
Let $n\ge 3$ and $u\in C^{3,\alpha}(\mathbb R^n)$ be a convex solution to
\begin{equation}\label{mp}
\det(D^2u)=f,\quad \mbox{ in }\quad \mathbb R^n
\end{equation}
where $f$ satisfies (\ref{af}). Then there exist $b\in \mathbb R^n$, $A\in \mathcal{A}$ (defined in (\ref{aa})) with
 $det(A)=\fint_{\Pi_{1\le i\le n}[0,a_i]}f$, and $v\in C^{3,\alpha}(\mathbb R^n)$, which is $a_i-$periodic in the $i$-th variable, such that
\begin{equation}\label{mes}
	\left | u(x)- \left (\frac 12x'Ax+b\cdot x+v(x) \right ) \right |\le C \left ( 1+|x| \right )^{-\sigma},\quad \forall x\in \mathbb R^n
\end{equation}
for some $C(d_0,d_1,n,\beta, a_1,...,a_n)>0$ and $\sigma:=\min\{\beta,n-2\}$.
\end{thm}
Here we note that $\fint_{\Pi_{1\le i\le n}[0,a_i]}f$ is the average of $f$ over the box.

\begin{rem} Theorem \ref{thm1} does not imply corresponding (better) estimates on higher order derivatives because of the oscillation of $D^2u$.
\end{rem}

Caffarelli and Li \cite{caf-li-2} proved that
$$u(x)=\frac 12 x'Ax+b\cdot x+v(x), \quad \mathbb R^n,\quad n\ge 2 $$
if $f=f_p$ in (\ref{mp-g}).
Thus Theorem \ref{thm1} is an extension of the theorem of Caffarelli-Li.

The assumption $\beta>2$ is essentially optimal, as one can  observe from the following example: let $f$ be a radial, smooth, positive function such that $f(r)\equiv 1$ for $r\in [0,1]$ and $f(r)=1+r^{-2}$ for $r>2$. Let
$$u(r)=n^{\frac 1n}\int_0^r\bigg (\int_0^s t^{n-1}f(t)dt\bigg )^{\frac 1n}ds,\quad r=|x|. $$
It is easy to check that $det(D^2u)=f$ in $\mathbb R^n$. Moreover for $n\ge 3$,
$$u(x)=\frac 12|x|^2+O(\log |x|)$$
at infinity, which means by taking $f_p\equiv 1$ the estimate in Theorem \ref{thm1} is violated for $n\ge 4$.

One major difficulty in the study of (\ref{mp-g}) is that the right hand side oscillates wildly when it is scaled. The regularity theory for Monge-Amp\'ere equations with oscillating right hand side is very challenging (see \cite{figalli,savin2} ). In this respect, we found that the strategy implemented by  Caffarelli and Li in  \cite{caf-li-2} is nearly perfect, as we can only simplify a small part of their augument in \cite{caf-li-2}. In turn, our arguments are essentially based on the corresponding steps from \cite{caf-li-2} as well as previous works of Caffarelli and Li \cite{caf-ann1,caf-ann2,caf-li-1}.
 The main difference, though, is that in the proof of Theorem \ref{thm1} one needs to take care of perturbational terms in a sharp manner. In order to handle all the perturbations, we need to make use of intrinsic structures implied by Monge-Amp\'ere equations (such as
(\ref{div-free}) below), estimates on Green's functions by Littman-Stampacchia-Weinberger \cite{littman} and Krylov-Safonov Harnack inequalities, etc.

Theorem \ref{thm1} is closely related to the exterior Dirichlet problem: Given a strictly convex set $D$ and the value of $u$ on $\partial D$, can one solve the Monge-Amp\'ere equation in $\mathbb R^n\setminus D$ if the asymptotic behavior of $u$ at infinity is prescribed? Clearly Theorem \ref{thm1} must be established before such a question can be attacked. We plan to address the exterior Dirichlet problem in a future work. The traditional (interior) Dirichlet problem has been fairly well understood through the contribution of many people (see \cite{alek1, b1, n2,calabi,p5,cheng-yau1,cns,caf-ann2,ivo-1,ivo-2,ivo-3,k2,jian-wang,tw7,urbas,savin2,figalli} and the references therein). If $f$ is equal to a positive constant outside a compact set, Delano\"e \cite{delanoe}, Ferrer-Mart\'inez-Mil\'an \cite{FMM1,FMM2} and Bao-Li \cite{bao} studied the exterior Dirichlet problem for $n=2$, Caffarelli-Li \cite{caf-li-1}) studied the case of $n\ge 3$. If $f$ is a perturbation of a positive constant at infinity, Bao-Li-Zhang \cite{bao} studied the exterior Dirichlet problem in \cite{zhang-mp}.

The organization of this article is as follows: The proof of Theorem \ref{thm1} consists of five steps. First in step one we employ the argument in \cite{caf-ann2,caf-li-1} to show that the growth of the solution of (\ref{mp-g}) is roughly similar to that of a parabola. Then in step two we prove that $D^2u$ is positive definite, which makes (\ref{mp-g}) uniformly elliptic. The key point in this step is to consider a second order incremental of $u$ as a subsolution to an elliptic equation. In step three we prove a pointwise estimate of the second incremental of $u$. The proof of Theorem \ref{thm1} for $n\ge 4$ is placed in step four since all the perturbations in this case are bounded. Finally in step five we prove the case $n=3$, which is a little different because of a logarithmic term. We use the Krylov-Safonov Harnack inequality for linear equations to overcome the difficulties caused by the logarithmic term.

\medskip

\noindent {\bf Acknowledgement:} The authors would like to thank the hospitality of Universidade Federal do Cear\'a and University of Florida, where this work was conducted. ET acknowledges support from CNPq and Funcap. Both authors are grateful to the anonymous referee, who read the whole paper very carefully and gave many excellent suggestions.

\section{Proof of Theorem \ref{thm1}}

Since the Monge-Amp\'ere equation is invariant under affine transformation, we assume $a_1=...=a_n=1$ and $\fint_{[0,1]^n} f=1$.

In step one we prove that $u$ grows like a quadratic polynomial at infinity. First we normalize $u$ to make $u(0)=0$ and $u\ge 0$ in $\mathbb R^n$. Since $f$ is bounded above and below by two positive constants, we use the argument in Caffarelli-Li \cite{caf-li-1}, see also \cite{zhang-mp}. Let
$$\Omega_M=\{x\in \mathbb R^n, \quad u(x)<M\}. $$

Then the following properties hold:
\begin{enumerate}
\item $C^{-1}M^{\frac n2}\le |\Omega_M |\le CM^{\frac n2}$,
\item There exists $A_M(x)=a_M\cdot x+b$ such that $det(a_M)=1$ and $B_R\subset A_M(\Omega_M)\subset B_{nR}$ and
\item $\frac 1C\sqrt{M}\le R\le C\sqrt{M}$,
\item $2nR\ge dist(a_M(\Omega_{\frac M2}),\partial a_M(\Omega_M))\ge \frac RC$
\end{enumerate}
where all the constants $C$ only depend on $d_1$ and $n$.

All the properties listed above are proved in \cite{caf-li-1} only based on the assumption that $f$ is bounded above and below by two positive constants. Here for the convenience of the reader we mention the idea of the proof:
First Caffarelli-Li used the following lemma (Lemma 2.9 in \cite{caf-li-1})

\medskip

\emph{Lemma A (Caffarelli-Li):  Let $e_1=(1,0,..,0)$ and
$$B_{\delta}'=\{(0,x_2,...,x_n);\quad |(0,x_2,...,x_n)|<\delta \}. $$
 Let $K$ be the convex hull of $\bar B_{\delta}'\cup \{re_1\}$, and let $u$ be a nonnegative convex viscosity solutions of
$det(D^2u)\ge \lambda>0$ in the interior of $K$.
Assume
$u\le \beta$ on $\bar B_{\delta}'. $
Then there exists $C(n)>1$ such that
$$\max \{\beta, u(re_1)\} \ge \frac{\lambda^{1/n}\delta^{\frac{2(n-1)}n}r^{2/n}}{C}. $$
}

In other words, $u$ cannot be small in one direction for too long. Lemma A is important since it says Pogorelov's famous example of non-strictly convex solution does not exist if the domain is large. A simple application of Lemma A leads to $\Omega_M\subset B_{CM^{\frac n2}}$. A volume preserving affine transformation can be used to make the image of $\partial \Omega_M$ between too balls with comparable radii. A comparison with a parabola gives $|\Omega_M|\sim M^{\frac n2}$ using only the upper bound and lower bound of $f$.  In \cite{caf-ann1} Caffarelli proved that $u$ must depart from its level set in a non-tangential manner, using this we have
$dist(A_M(\partial \Omega_M),A_M(\Omega_{M/2}))\sim M^{\frac 12}$ where $A_M$ is a volume-preserving affine transformation:
$$A_M(x)=a_Mx+b_M,\quad det(a_M)=1. $$
Using $u(0)=0$ and $u\ge 0$ one can further conclude
\begin{equation}\label{like-ball}
B_{R/C}\subset a_M(\Omega_M)\subset B_{2nR}\quad \mbox{ where } \quad R=M^{\frac 12}.
\end{equation}
Equation (\ref{like-ball}) makes it convenient to define
$$u_M(x)=\frac 1{R^2}u(a_M^{-1}(Rx)),\quad x\in O_M:=\frac 1R a_M(\Omega_M), $$
Using $u(0)=0$ and $u\ge 0$ we have
$$B_{1/C}\subset O_M\subset B_{2n}. $$
In order to analyze the level surfaces of $u_M$, which satisfies
\begin{equation}\label{mp-m}
det(D^2u_M(\cdot))=f(a_M^{-1}(R\cdot)),\quad \mbox{ in }\quad O_M.
\end{equation}
We use the following equation to give a good approximation of $u_M$:
\begin{equation}\label{mp-m2}
\left\{\begin{array}{ll}
det(D^2w_p)=f_p(a_M^{-1}(Rx)),\quad \mbox{ in } O_M, \\
w_p=M/R^2,\quad \mbox{ on }\quad \partial O_M.
\end{array}
\right.
\end{equation}
Let $h=u_M-w_p$ be the difference of $u_M$ and $w_p$, then by the Alexandrov estimate ( see \cite{caf-li-1,zhang-mp}) we have
$$\max_{O_M}(h^-)\le C\bigg (\int_{S^+}det(D^2(u_M-w_p))\bigg )^{\frac 1n}$$
where $h^-$ is the negative part of $h$: ($h=h^+-h^-$) and
$$S^+=\{ x\in O_M;\quad D^2(u_M-w_p)>0\quad \}. $$
On $S^+$ by the concavity of $det^{\frac 1n}$ on positive definite matrices we have
$$det^{\frac 1n}(\frac{D^2u_M}2)\ge \frac 12det^{\frac 1n}(D^2(u_M-w_p))+\frac 12det^{\frac 1n}(D^2w_p). $$
Thus
$$\max_{O_M}h^-\le C\bigg (\int_{O_M}|f^{\frac 1n}(a_M^{-1}(Rx))-f_p^{\frac 1n}(a_M^{-1}(Rx))|^ndx\bigg )^{\frac 1n}. $$
It is easy to see from (\ref{af}) that the right hand side is $O(1/R)$. Next we cite the homogenization theorem of Caffarelli-Li (Theorem 3 of \cite{caf-li-2}): Let $w$ satisfy
$$det(D^2w)=1,\quad \mbox{ in } \quad O_M$$ with $w=M/R^2$ on $\partial O_M$, then
$$|w-w_p|\le CR^{-\delta} $$
for some $\delta (d_0,d_1)>0$. Thus
\begin{equation}\label{close-1}
|u_M(x)-w(x)|\le CR^{-\delta}
\end{equation}
for some $\delta>0$. Set
$$E_M:=\{x;\quad (x-\bar x)'D^2w(\bar x)(x-\bar x)\le 1\} $$
where $\bar x$ is the unique minimum of $w$ (note that Caffarelli \cite{caf-ann1} proved that the minimum point of $w$ is unique) that satisfies $dist(\bar x,\partial O_M)>C(d_1,n)$. By the same argument in \cite{caf-ann2,caf-li-1} there exist $\bar k$ and $C$ depending only on $n$ and $d_1$ such that for $\epsilon<<\delta$, $M=2^{(1+\epsilon)k}$, $2^{k-1}\le M'\le 2^k$, $R\sim M^{\frac 12}$,
$$(\frac{2M'}{R^2}-C2^{-\frac{3\epsilon k}2})^{\frac 12}E_M\subset
\frac 1Ra_M(\Omega_{M'})\subset (\frac{2M'}{R^2}+C2^{-\frac{3\epsilon k}2})^{\frac 12}E_M,\quad \forall k\ge \bar k, $$
which can be translated as
$$\sqrt{2M'}(1-\frac{C}{2^{\epsilon k/2}})E_M\subset a_M(\Omega_{M'})\subset \sqrt{2M'}(1+\frac{C}{2^{\epsilon k/2}})E_M. $$
Let $Q$ be a positive definite matrix satisfying $Q^2=D^2w(\bar x)$, $O$ be an orthogonal matrix that makes $T_k=OQa_M$ upper triangular, then $det(T_k)=1$ and by Proposition 3.4 of \cite{caf-li-1}
$$\|T_k-T\|\le C2^{-\frac{\epsilon k}2} $$
for some matrix $T$. By setting
$$v=u\cdot T$$
we have
$$det(D^2 v(x))=f(Tx) $$
and
$$\sqrt{2M'}(1-\frac{C}{2^{\epsilon k/2}})B_1\subset \{x;\quad v(x)<M'\}\subset \sqrt{2M'}(1+\frac{C}{2^{\epsilon k/2}})B_1 $$
for all $M'\ge 2^{\bar k}$. Here $B_1$ stands for the unit ball (throughout the article we use $B(p,r)$ to denote the ball centered at $p$ with radius $r$. If $p$ is the origin we may use $B_r$). Consequently
$$|v(x)-\frac 12 |x|^2|\le C|x|^{2-\epsilon}. $$
The equation for $v$ is
$$det(D^2 v)=f_v(x) $$
where $f_v(x)=f(T(x))$, correspondingly we let $f_{v,p}(x)=f_p(T(x))$.

\medskip
\noindent{\bf Step two: Uniform Ellipticity}
\medskip

The purpose of this step is to show: There exist $c_1$ and $c_2$ depending only on $d_0,d_1,\beta$, $a_1,...,a_n$ and $n$ such that
\begin{equation}\label{uni-ell}
c_1 I\le D^2 v\le c_2 I.
\end{equation}

First we choose $M>100$ so that for $R=M^{\frac 12}$ and
$$v_R(y)=\frac 1{R^2}v(Ry), $$
$$\Omega_{1,v_R}:=\{y;\quad v_R(y)\le 1\} $$
is very close to $B_{\sqrt{2}}$ in the sense that $\partial \Omega_{1,v_R}\subset B_{\sqrt{2}+\epsilon}-B_{\sqrt{2}-\epsilon}$ for some $\epsilon>0$ small. Applying the standard interior estimate for Monge-Amp\'ere equations (\cite{caf-ann2,jian-wang} ) we have
$$\|D^2v_R\|_{L^{\infty}(\Omega_{1,v_R})}=\|D^2v\|_{L^{\infty}(B_{2})}\le C. $$
In general for $|x|>100$ we shall prove (\ref{uni-ell}) for $D^2v(x)$.
To this end we consider the following second order incremental for $v$:
$$\Delta_e^2v(x):=\frac{v(x+e)+v(x-e)-2v(x)}{\|e\|^2} $$
where $e\in \mathbb R^n$ and $\|e\|$ is its Euclidean norm.
Later we shall always choose $e\in E$ which is defined as
\begin{align}\label{E-def}
E:=&\{a_1 v_1+...+a_nv_n;\quad a_1,...a_n\in \mathbb Z,\\
&\quad f_{v,p}(x+v_i)=f_{v,p}(x), \, \forall x\in \mathbb R^n,\quad i=1,...,n. \} \nonumber
\end{align}

Let $R_1=|x|$ (recall $|x|>10$), $e_x=\frac{x}{R_1}$ and
$$v_{R_1}(y)=\frac 1{R_1^2}v(R_1y),\quad y\in B(e_x,\frac 23). $$
 Then by the closeness result of step one, the sections: $S_{v_{R_1}}(e_x,\frac 14)$ and
$S_{v_{R_1}}(e_x,\frac 18)$ around $e_x$ are very close to the corresponding sections of the parabola $\frac 12 |y|^2$. Here we recall that for a convex, $C^1$ function $v$,
$$S_v(x, h):=\{y;\quad v(y)\le v(x)+\nabla v(x)\cdot (y-x)+h,\,\, \}. $$
In other words, $S_{v_{R_1}}(e_x,\frac 14)$ looks very similar to $B(e_x,1/\sqrt{2})$ and $S_{v_{R_1}}(e_x,\frac 18)$ is very close to $B(e_x,1/2)$.
Moreover the equation for $v_{R_1}$ is
\begin{equation}\label{vr2}
det(D^2v_{R_1}(y))=f_v(R_1y)
\end{equation}
where $f_v(R_1y)=f(T(R_1y))$.
Let $e_{R_1}=e/R_1$, then direct computation shows
\begin{align*}
&\Delta_e^2v(x)=\frac{v(x+e)+v(x-e)-2v(x)}{\|e\|^2},\quad x=R_1y\\
=&\frac{(v(R_1y+e)+v(R_1y-e)-2v(R_1y))}{\|e\|^2}\\
=&\frac{R_1^2(v_{R_1}(y+e_{R_1})+v_{R_1}(y-e_{R_1})-2v_{R_1}(y))}{\|e\|}=\Delta_{e_{R_1}}^2v_{R_1}(y).
\end{align*}
Let
$$w(y)=\frac{v_{R_1}(y+e_{R_1})+v_{R_1}(y-e_{R_1})}2 $$
and $F=det^{\frac 1n}$. Then the concavity of $F$ on positive definite matrices gives
\begin{align}
det^{\frac 1n}(D^2w)(y)&\ge \frac 12 det^{\frac 1n}(D^2 v_{R_1}(y+e_{R_1}))+\frac 12 det^{\frac 1n}(D^2v_{R_1}(y-e_{R_1}))\nonumber\\
&=\frac 12 f_v^{\frac 1n}(R_1y+e)+\frac 12 f_v^{\frac 1n}(R_1y-e). \label{det-con-1}
\end{align}
On the other hand
\begin{align} \label{det-con-2}
F(D^2w)&\le F(D^2 v_{R_1})+F_{ij}(D^2 v_{R_1})\partial_{ij}(w-v_{R_1})\\
&=f_v^{\frac 1n}(R_1y)+F_{ij}(D^2v_{R_1})\partial_{ij}(w-v_{R_1}) \nonumber
\end{align}
where
$$F_{ij}(D^2v_{R_1})=\frac{\partial F(D^2v_{R_1})}{\partial_{ij}v_{R_1}}=\frac 1n (det(D^2v_{R_1}))^{\frac 1n-1}cof_{ij}(D^2v_{R_1}).$$
Thus the combination of (\ref{det-con-1}) and (\ref{det-con-2}) gives
$$a_{ij}\partial_{ij}(\Delta_{e_{R_1}}^2 v_{R_1})\ge E_1$$
where
\begin{align*}
a_{ij}&=cof_{ij}(D^2v_{R_1}), \\
E_1:&=ndet(D^2v_{R_1})^{(n-1)/n}\\
&\cdot (\frac{R_1^2}{2\|e\|^2}(f_v^{\frac 1n}(R_1y+e)+f_v^{\frac 1n}(R_1y-e)-2f_v^{\frac 1n}(R_1y)))
\end{align*}
For applications later we state the following fact: For any smooth $u$
\begin{equation}\label{div-free}
\partial_i(cof(D^2u)_{ij})=0.
\end{equation}

Let $f_{v,p}(y)=f_p(Ty)$, then by choosing $e\in E$ (see (\ref{E-def}))
\begin{equation}\label{per-e}
f_{v,p}^{\frac 1n}(R_1y+e)+f_{v,p}^{\frac 1n}(R_1y-e)-2f_{v,p}^{\frac 1n}(R_1y)=0,\quad \mbox{ for all }y\in \mathbb R^n.
\end{equation}

By (\ref{af}) we see that, for $y\in B(e_x, \frac 14)$ (which implies $2>|y|>1-1/\sqrt{2}$) and $e\in E$ with $\|e\|\le \frac 1{10}R_1$
\begin{align}\label{error-1}
E_1(y)&=n(det(D^2v_{R_1})^{(n-1)/n})\frac{R_1^2}{2\|e\|^2}\{(f_v^{\frac 1n}(R_1y+e)-f_{v,p}^{\frac 1n}(R_1y+e))\nonumber\\
&+(f_v^{\frac 1n}(R_1y-e)-f_{v,p}^{\frac 1n}(R_1y-e))-2(f_v^{\frac 1n}(R_1y)-f_{v,p}^{\frac 1n}(R_1y))\}\nonumber\\
=&O(R_1^{-\beta}).\nonumber
\end{align}

Now we construct a function $h_0$ that solves
$$\left\{\begin{array}{ll}
a_{ij}\partial_{ij}h_0=-E_1,  \quad \mbox{ in  } B(e_x,\frac 14) \\
h_0=0\quad \mbox{ on }\quad \partial B(e_x, \frac 14)
\end{array}
\right.
$$

Then we use the following classical estimate of Aleksandrov (see Page 220-222 of \cite{gilbarg} for a proof):

\emph{Theorem A: Let $\Omega$ be a domain in $\mathbb R^n$ and let $v$ be a solution in $\Omega$ of the equation
$$a^*_{ij}\partial_{ij}v=g $$
such that $v=0$ on $\partial \Omega$ and the coefficient matrix $(a^*_{ij})_{n\times n}$ satisfies
$$c_1\le det(a^*_{ij})\le c_2, \quad \mbox{and}\quad (a^*_{ij})>0, $$
 then
$$|v(x)|\le C(n,c_1,c_2) diam(\Omega) \|g\|_{L^n(\Omega)}, \quad \forall x\in \Omega$$ }

Applying Theorem A to $h_0$ we have
\begin{equation}\label{error-2}
|h_0(y)|\le C(n,d_0)R_1^{-\beta},\quad \mbox{ for }\quad y\in B(e_x,\frac 14).
\end{equation}
\begin{rem} The estimate of Theorem A does not depend on constants of uniform ellipticity.
\end{rem}

Thus $\Delta_{e_{R_1}}^2 v_{R_1}+h_0$ is a super solution:
$$a_{ij}\partial_{ij}(\Delta_{e_{R_1}}^2 v_{R_1}+h_0)\ge 0,\quad \mbox{in}\quad B(e_x,\frac 14). $$
In order to obtain a pointwise estimate for $\Delta_{e_{R_1}}^2 v_{R_1}$ we recall the following Harnack inequality of Caffarelli-Gutierrez in \cite{gute}:

 \medskip

 \emph{Theorem B (Caffarelli-Gutierrez): Let $O$ be a convex set in $\mathbb R^n$ and $B_1\subset O\subset B_n$ ($n\ge 2$). Suppose $\phi\in C^2(O)$ satisfies, for $0<\lambda<\Lambda<\infty$,
 $$ \lambda<det(D^2\phi)\le \Lambda, \mbox{ in } \,\, O, \mbox{ and } \quad \phi=0 \mbox{ on } \quad \partial O. $$
 Then for any $r>s>0$ and any $w\in C^2(O)$ satisfying
 $$a_{ij}\partial_{ij}w\ge 0,\quad w\ge 0,$$
 where $a_{ij}=det(D^2\phi)\phi^{ij}$ and $(\phi^{ij})_{n\times n}=(D^2\phi)^{-1}$, we have
 $$\max_{x\in O, dist(x,\partial O)>r}w\le C\int_{x\in O,dist(x,\partial O)>s}w, $$
 for some $C>0$ depending only on $\lambda, \Lambda, r,s$. }

 \medskip

Using Theorem B for $\Delta_{e_{R_1}}^2 v_{R_1}+h_0$ we have
$$\max_{y\in S(v_{R_1},e_x, \frac 1{16})}(\Delta_{e_{R_1}}^2 v_{R_1}+h_0)\le C(n,d_0)\int_{S(v_{R_1},e_x,\frac 18)}(\Delta_{e_{R_1}}^2
v_{R_1}+h_0) $$
Note that the distance between the two sections above is comparable to $1$. It is also important to point out that the estimates in Theorem B does not depend on the regularity of $v_{R_1}$.
 Before we proceed we cite the following Calculus lemma that can be found in \cite{caf-li-2}:

\medskip

\emph{Lemma B: Let $g\in C^2(-1,1)$ be a strictly convex function, and let $0<|h|\le \epsilon$. Then
$$\Delta_h^2g(x)>0,\, \forall |x|\le 1-2\epsilon, \quad \mbox{ and } \quad \int_{-1+2\epsilon}^{1-2\epsilon}\Delta_h^2g\le \frac{C}{\epsilon} osc_{(-1,1)}g, $$
where $C$ is a universal constant and $osc$ stands for oscillation. }

\medskip

 Let $L$ be a line parallel to $e$, then by Lemma B,
$$\int_{S(v_{R_1},e_x,\frac 18)\cap L}\Delta_{e_{R_1}}^2 v_{R_1}\le C. $$
Consequently, letting $L$ go through all directions, we have
$$\int_{S(v_{R_1},e_x,\frac 18)}\Delta_{e_{R_1}}^2 v_{R_1}\le C. $$
Since $h_0=O(R_1^{-\beta})$  we have proved
\begin{equation}\label{quo-2}
0\le \Delta_e^2v(x)=\Delta_{e_{R_1}}^2 v_{R_1}(y)\le C,\quad \mbox{ for } y\in B(e_x,\frac 1{16}),\,\, x=R_1y.
\end{equation}
Note that $\Delta_e^2v\ge 0$ because $v$ is convex. Then the same argument as in \cite{caf-li-2} can be employed to prove that the level surfaces of $v$ are like balls.

For the convenience of the readers we describe the outline of this argument. Given $x\in \mathbb R^n$ with $|x|>100$ we set
$$\gamma:=\sup_{e\in E, \|e\|\le \frac 1{10}\|x\|}\sup_{y\in \mathbb R^n}\Delta_e^2v(y)$$
 and
$$\bar v(z)=v(z+x)-v(x)-\nabla v(x)z.$$
Clearly $\bar v(0)=0=\min_{\mathbb R^n} \bar v$.
By (\ref{quo-2}) it is easy to see
$$\sup_{B_r}\bar v\le C(n)\gamma r^2,\quad \forall r \in (1, \frac{1}{10}\|x\|).  $$
On the other hand for any $\bar z\in \partial B_r$ we show that for $r$ large (but still only depending on $n$ and $d_0$), $\bar v(\bar z)\ge 1$. Indeed, by (\ref{quo-2}),
\begin{equation}\label{uni-2}
\bar v(\frac{\bar z}2+e)+\bar v(\frac{\bar z}2-e)-2\bar v(\frac{\bar z}2)\le \gamma \|e\|^2
\end{equation}
for all $e\in E$ and $\|e\|\le \frac 1{10} \|x\|$. Hence for $z\in \frac{\bar z}2+(-2,2)^n$, we find $e\in E$ with $\|e\|\le C(n)$ so that $z$ is on the same line with $\frac{\bar z}2+e$ and $\frac{\bar z}2-e$
and is between them. Thus(\ref{uni-2}) implies
$$\bar v(z)\le \bar v(\frac{\bar z}2+e)+\bar v(\frac{\bar z}2-e)\le 2\bar v(\frac{\bar z}2)+C(n)\gamma. $$
Further more, by $\bar v(0)=0$ and the convexity of $\bar v$ we have
$$2\bar v(\frac{\bar z}2)\le \bar v(\bar z). $$
Therefore the following holds:
$$\bar v(z)\le \bar v(\bar z)+C(n)\gamma, \quad z\in \frac{\bar z}2+(-2,2)^n. $$
Consider $$w(z)=\frac{\bar v(\frac{\bar z}2+z)}{\bar v(\bar z)+C(n)\gamma}. $$
Clearly $w$ satisfies $det(D^2w)\ge d_0^{-1}/(\bar v(\bar z)+C(n)\gamma)^n$ and
$$w(z)\le 1, \quad \mbox{ for  }\quad z\in \bar z/2+(-2,2)^n. $$
Applying Lemma A  we have
$$\max \{1,\bar v(\bar z)\}\ge \bigg (\frac{d_0^{-\frac 1n} C(n)}{\bar v(\bar z)+C(n)\gamma} \bigg ) r^{\frac 2n}.  $$
Here we used the fact that $\bar v(0)=0$ and $\bar v(\bar z)$ is the maximum value of $\bar v$ on the line segment connecting $0$ and $\bar z$.
If $\bar v(\bar z)\ge \gamma$, no need to do anything ( here we assume $\gamma>1$ without loss of generality), otherwise we obviously have
$$\max\{1,\bar v(\bar z)\}\ge \frac{d_0^{-\frac 1n}}{C(n)\gamma }r^{\frac 2n}. $$
Choose $r$ large but only depending on $n$,$\gamma$ and $d_0$ we can still make the right hand side of the above greater than $1$. Therefore we have proved $\bar v(z)\ge 1$.
By standard argument the sections of $\bar v$:
$S(\bar v,0,h)$ are similar to balls if $h\sim 1$ and $h<1$.
Then interior estimate of Caffarelli, Jian-Wang \cite{jian-wang} can be employed on $S(\bar v, 0, h)$ to obtain the $C^{2,\alpha}$ norm of $\bar v$. In particular
$|D^2v(x)|=|D^2\bar v(x)|\le C$. Once the upper bound is obtained we also have the lower bounded from the Monge-Amp\'ere equation. (\ref{uni-ell}) is established.

\begin{rem} In order to prove (\ref{uni-ell}) the assumptions on the derivatives of $f^{\frac 1n}-f_p^{\frac 1n}$ are not essential. In fact, as long as
$$|(f^{\frac 1n}-f_p^{\frac 1n})(x)|\le C(1+|x|)^{-\beta}, \quad \mbox{ for }\quad x\in \mathbb R^n, $$
(\ref{uni-ell}) still holds.
\end{rem}

\noindent{\bf Step three. Pointwise estimate of $\Delta_e^2v$}

Let $e\in E$, recall the equation for $\Delta_e^2v$ is
\begin{equation}\label{p3e1}
\tilde a_{ij}(x)\partial_{ij}\Delta_e^2v\ge E_2.
\end{equation}
where
$$ \tilde a_{ij}(x)=cof_{ij}(D^2v(x))  $$ is uniformly elliptic and divergence free,
$$E_2=n(det(D^2v)^{\frac{n-1}n})\frac{f_v^{\frac 1n}(x+e)+f_v^{\frac 1n}(x-e)-2f_v^{\frac 1n}(x)}{\|e\|^2}. $$
By the assumption on $f$,
\begin{align*}
|E_2(x)|&=n(det(D^2v)^{\frac{n-1}n})\\
&\bigg(\frac{(f_v^{\frac 1n} (x+e)-f_{p,v}^{\frac 1n}(x+e))+(f_v^{\frac 1n}(x-e)-f_{v,p}^{\frac 1n}(x-e))}{\|e\|^2}\nonumber\\
&-2\frac{(f_v^{\frac 1n}(x)-f_{v,p}^{\frac 1n}(x))}{\|e\|^2}\bigg ).
\end{align*}
Here we have the following important observation:
\begin{equation}\label{e2e}
|E_2(x)|\le C(1+|x|)^{-2-\beta},\quad \forall e\in E
\end{equation}
where $C$ is independent of $e\in E$ and $x$.
Indeed, let $g(x)=f_v^{\frac 1n}(x)-f_{v,p}^{\frac 1n}(x)$, then
$$|E_2(x)|=(g(x+p)+g(x-p)-2g(x))/\|e\|^2. $$
Without loss of generality we assume that $e=(\|e\|,0,...,0)$. Then
\begin{align*}
&E_2(x)=\frac{(g(x+e)-g(x))-(g(x)-g(x-e))}{\|e\|^2}\\
=&\frac{\int_0^{\|e\|}\partial_1g(x_1+t,x')dt-\int_0^{\|e\|}\partial_1g(x_1-\|e\|+t,x')dt}{\|e\|^2}\\
=&\frac{\int_0^{\|e\|}\int_0^{\|e\|}\partial_{11}g(x_1-\|e\|+t+s,x')dtds}{\|e\|^2}\\
=&\int_0^1\int_0^1\partial_{11}g(x_1-\|e\|+\|e\|t+\|e\|s,x')dtds.
\end{align*}
Using
$$|D^2g(x)|\le C(1+|x|)^{-\beta-2},$$
which comes from (\ref{af}) we have
$$|E_2(x)|\le C\int_0^1\int_0^1\bigg ( \big (x_1+\|e\|(-1+t+s)\big )^2+|x'|^2\bigg )^{-\frac{\beta+2}2}dtds. $$
If $|x'|> \frac 18 |x|$ , (\ref{e2e}) holds obviously. So we only consider when $|x_1|> \frac 12 |x|$.
\begin{align*}
&\int_0^1\int_0^1\bigg (\big (x_1+\|e\|(-1+t+s) \big )^2+|x'|^2\bigg )^{-\frac{\beta+2}2}dtds\\
\le & \int_0^1\int_0^1 \bigg | x_1+\|e\|(-1+t+s) \bigg |^{-\beta -2}dtds \\
=&\|e\|^{-\beta -2}\int_0^1\int_0^1 |L-1+t+s |^{-\beta-2}dtds,\quad L=x_1/\|e\|.
\end{align*}
If $|L|>8$, $|L-1+t+s|>|L|/2$, which gives
\begin{align*}
&\int_0^1\int_0^1\bigg (\big (x_1+\|e\|(-1+t+s) \big )^2+|x'|^2\bigg )^{-\frac{\beta+2}2}dtds\\
 &\le C\|e\|^{-\beta -2}|L|^{-\beta-2}\le C\|x_1\|^{-\beta}\le C|x|^{-\beta-2}.
\end{align*}
where all the constants are absolute constants. Clearly (\ref{e2e}) holds in this case.

If $|L|<8$, easy to see that
$$\int_0^1\int_0^1|L-1+t+s|^{-\beta-2}dtds \le C(\beta). $$
Then
\begin{align*}
&\int_0^1\int_0^1\bigg (\big (x_1+\|e\|(-1+t+s) \big )^2+|x'|^2\bigg )^{-\frac{\beta+2}2}dtds\\
\le &C(\beta)\|e\|^{-\beta-2}\le C(\beta)\|x\|^{-\beta-2},
\end{align*}
since $|L|<8$ implies $\|e\|>|x_1|/8>|x|/16$.
 (\ref{e2e}) is established.

By standard elliptic estimates $f\in C^{1,\alpha}$ yields $\tilde a_{ij}\in C^{1,\alpha}$. Moreover $\tilde a_{ij}$ is divergence free. The Green's function $G(x,y)$ corresponding to
$-\partial_{x_i}(\tilde a_{ij}\partial_{x_j})$ satisfies
$$-\partial_{x_i}(\tilde a_{ij}\partial_{x_j}G(x,y))=\delta_y $$
and
\begin{equation}\label{stand-g}
0\le G(x,y)\le C|y-x|^{2-n}, \quad x,y\in \mathbb R^n
\end{equation}
by a result of
 Littman-Stampacchia-Weinberger in \cite{littman}.

Let $E_2^-\ge 0$ be the negative part of $E_2$: $E_2=E_2^+-E_2^-$ and
$$h(x)=-\int_{\mathbb R^n}G(x,y)E_2^-(y)dy.  $$
Then $h(x)$ satisfies
\begin{equation}\label{p3e2}
\partial_{x_i}(\tilde a_{ij}\partial_{x_j})h=E_2^-,\quad \mbox{in}\quad \mathbb R^n.
\end{equation}
Here we claim that $h$ satisfies
\begin{equation}\label{p3h1}
0\le -h(x)\le C(1+|x|)^{-\beta_1}, \quad \beta_1=\min\{n-2,\beta\}
\end{equation}

The estimate of (\ref{p3h1}) is rather standard, for $x\in \mathbb R^n$, we divide $\mathbb R^n$ into three regions:
\begin{align*}
\Omega_1=\{y;\,\, |y-x|<|x|/2\,\}\\
\Omega_2=\{y;\quad |y|<|x|/2\,\, \}, \\
\Omega_3=\mathbb R^n\setminus (\Omega_1\cup \Omega_2).
\end{align*}

Then one observes $\int_{\Omega_i}G(x,y)E_2^-(y)dy$ for each $\Omega_i$ can be estimated easily. (\ref{p3h1}) is established.
From the definition of $h$ we see that
\begin{equation}\label{delta-vij}
\tilde a_{ij}\partial_{ij}(\Delta_e^2v+h)\ge 0,\quad \mbox{ in }\quad \mathbb R^n.
\end{equation}
The main result of this step is:
\begin{equation}\label{2-e}
\mbox{ Given  }e\in E;\quad \Delta_e^2v(x)\le 1-h(x),\quad \mbox{ for }\quad x\in \mathbb R^n.
\end{equation}

Let $v_{\lambda}(x)=v(\lambda x)/\lambda^2$ and $P(x)=\frac 12 |x|^2$. First we claim that
\begin{equation}\label{2-e1}
D^j(v_{\lambda}-P(x))\to 0,\quad j=0,1,\quad \forall x\in K\subset\subset \mathbb R^n, \mbox{  as  }\lambda\to \infty,
\end{equation}
where $K$ is any fixed compact subset of $\mathbb R^n$.

To see (\ref{2-e1}), by step one
\begin{equation}\label{2-e2}
|v_{\lambda}(x)-P(x)|\le C\lambda^{-\epsilon},\quad \forall x\in K\subset\subset \mathbb R^n.
\end{equation}
and $|D^2v_{\lambda}(x)|\le C$ for $x\in K$ we obtain by Ascoli's theorem that $\partial_lv_{\lambda}(x)$ ($l=1,..,n$) tends to a continuous function. By (\ref{2-e2}) this function has to be $x_l$. Thus (\ref{2-e1}) holds.
\begin{rem} We don't have the estimate of $\|Dv_{\lambda}-x\|_{L^{\infty}(K)}$ but we don't need it.
\end{rem}

The proof of (\ref{2-e}) is as follows.

 Let $\alpha=\sup_{\mathbb R^n}(\Delta^2_ev+h)$ for $e\in E$. By step two $\alpha<\infty$. Let $\hat e=e/\lambda$, then
by (\ref{2-e1}) and  Lebesgue's dominated convergence theorem
\begin{equation}\label{caf-li-e}
\lim_{\lambda\to \infty}\int_{B_1}\Delta_{\hat e}^2v_{\lambda}=\int_{B_1} 1dx=|B_1|.
\end{equation}
Indeed, the integral over $B_1$ can be considered as the collection of integration on segments all in the direction of $\hat e$. Since $Dv_{\lambda}\to DP$ in $C^0$ and $DP$ is smooth,
the following Lemma C in \cite{caf-li-2} implies (\ref{caf-li-e}):

\medskip

\emph{Lemma C (Caffarelli-Li) Let $g_i$ converge to $g$ in $C^1[-1,1]$, $g\in C^2(-1,1)$ and $|h_i|\to 0$ as $i\to \infty$. Then for all $-1<a<b<1$,
$$\lim_{i\to \infty} \int_a^b \Delta_{h_i}^2g_i=g'(b)-g'(a)=\int_a^b g''. $$
}

\medskip

Let
$$h_{\lambda}(x)=h(\lambda x)/\lambda^2. $$
It follows immediately from (\ref{p3h1}) that
$$\lim_{\lambda\to \infty}\int_{B_1}(\Delta_{\hat e}^2v_{\lambda}+h_{\lambda})=|B_1|, $$
which obviously implies $\alpha\ge 1$ because otherwise it is easy to see that the left hand side of the above is less than $|B_1|$.  Our goal is to show that $\alpha=1$. If this is not the case we have $\alpha>1$. Then
$$\lim_{\lambda}\sup(\frac{\alpha+1}2|\{\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}\ge \frac{\alpha+1}2\}\cap B_1|\le \lim_{\lambda\to \infty}\int_{B_1}\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}=|B_1|. $$

Thus
$$\frac{|\{\Delta_{\hat e}^2 v_{\lambda}+h_{\lambda}\ge \frac{\alpha+1}2\}\cap B_1|}{|B_1|}\le 1-\mu $$
for $\mu=\frac 12(1-\frac{2}{\alpha+1})>0$. Here we emphasize that it is important to have $\mu>0$. Equivalently for large $\lambda$
$$\frac{|\{\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}\le \frac{\alpha+1}2\}\cap B_1|}{|B_1|}\ge \mu. $$

Then we cite the following well known result which can be found in \cite{caf-cab} (Lemma 6.5):

\medskip

\emph{Theorem C:
 For $v$ satisfying
$$a^*_{ij}\partial_{ij}v\ge 0, \quad \lambda I\le (a^*_{ij}(x))\le \Lambda I  $$
and
$$v\le 1\quad \mbox{ in } B_1 \quad \mbox{ and } \quad |\{v\le 1-\epsilon\}\cap B_1|\ge \mu |B_1|, $$
then $v\le 1-c(n,\lambda,\Lambda,\epsilon,\mu)$ over $B_{1/2}$
}
\medskip

By (\ref{delta-vij}) we see that
$$ \tilde a_{ij}(\lambda \cdot)\partial_{ij}(\Delta_{\hat e}^2v_{\lambda}+h_{\lambda})\ge 0, \quad \mbox{in }\quad B_1. $$
Applying Theorem C to
$\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}$ we have
\begin{equation}\label{2e-3}
\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}\le \alpha-C,\quad \mbox{in }\quad B_{1/2}.
\end{equation}
Since $\Delta_{\hat e}^2v_{\lambda}+h_{\lambda}$ is sub-harmonic,
$$\alpha=\sup_{\mathbb R^n}\Delta_e^2v+h=\lim_{\lambda\to \infty}\sup_{B_{1/2}}(\Delta_{\hat e}^2v_{\lambda}+h_{\lambda})<\alpha. $$
Then we get a contradiction, (\ref{2-e}) is established.

\medskip

\noindent{\bf Step four: The proof of the Liouville theorem for $n\ge 4$.}

\medskip

By a result of Li \cite{li-cpam}
there exists $\xi\in C^{2,\alpha}(\mathbb R^n)$ such that
$$\left\{\begin{array}{ll}
det(I+D^2\xi)=f_{v,p},\quad \mbox{ in }\quad \mathbb R^n, \\
I+D^2\xi>0,
\end{array}
\right.
$$
where $\xi$ is a periodic function with the same period as that of $f_{v,p}$.
Let $P(x)=\frac 12|x|^2$ and
$$w(x)=v(x)-P(x)-\xi(x). $$
Using $\xi(x+e)=\xi(x)$ for all $x\in \mathbb R^n$ and $e\in E$ , $\Delta_e^2P=1$ and (\ref{2-e}) we have
\begin{equation}\label{al-ca}
\Delta_e^2w+h\le 0, \quad \forall x\in \mathbb R^n, \quad \forall e\in E.
\end{equation}

From the equations for $v$ and $P+\xi$ we have,
$$det(D^2v)-det(D^2(P+\xi))=f_v-f_{v,p}, \quad \mbox{in}\quad \mathbb R^n, $$
which gives
$$
\hat a_{ij}\partial_{ij}w=f_v-f_{v,p}\quad \mbox{in }\quad \mathbb R^n
$$
where
$$\hat a_{ij}(x)=\int_0^1 cof_{ij}(tD^2v+(1-t)D^2(P+\xi))dt $$
is uniformly elliptic and divergence free. Using the Green's function corresponding to $-\partial_i(\hat a_{ij}(x)\partial_j)$ we
find $h_3$ to solve
$$\hat a_{ij}\partial_{ij}h_3=f_{v,p}-f_v,\quad \mbox{in}\quad \mathbb R^n $$
and
$$h_3(x)=O((1+|x|)^{-\beta_1}), \quad x\in \mathbb R^n. $$
Therefore we have
\begin{equation}\label{w-sol}
\hat a_{ij}\partial_{ij}(w+h_3)=0,\quad \mbox{in }\quad \mathbb R^n, \quad n\ge 3.
\end{equation}
Next we let
$$M_k=\sup_{B(0,k)}(w+h_3). $$
Then we claim that
\begin{equation}\label{double-M}
M_k\le 2M_{k/2}+C.
\end{equation}
for $C>0$ independent of $k$.

In order to prove (\ref{double-M}) we consider the equation for $\Delta_e^2w$ where $e\in E$.
In view of (\ref{p3e2}) we have
$$\tilde a_{ij}\partial_{ij}(\Delta_e^2w)\ge E_2,\quad \mbox{ in }\quad \mathbb R^n. $$
where $\tilde a_{ij}(x)=cof_{ij}(D^2v)$ is divergence free and $|E_2(x)|\le C(1+|x|)^{-2-\beta}$.

Just like the estimate of $h$ we can find $h_2$ that satisfies
$$\tilde a_{ij}\partial_{ij}(\Delta_e^2w+h_2)\ge 0 \quad \mbox{ in }\quad \mathbb R^n $$
and
$$|h_2(x)|\le C(1+|x|)^{-\beta_1},\quad \mbox{ in }\quad \mathbb R^n. $$

Let $x_0$ be where $\Delta_e^2 w+h_2$ attains its maximum on $B_k$, clearly $x_0\in \partial B_k$. Let $e=\frac{x_0}2+a$ be chosen so that $e\in E$ and $|a|$ is small:
$|a|\le C$.
Then for $k$ large
$$\Delta_e^2w(x-e)\le Ck^{-\beta_1}$$
which is
\begin{equation}\label{lconv}
w(x_0)\le 2w(x_0-e)-w(x_0-2e)+ O(k^{2-\beta_1}).
\end{equation}
Since $|x_0-2e|$ is in the neighborhood of $0$, we have $|w(x_0-2e)|\le C$
and
$$w(x_0)\le 2\max_{B(0,k/2)}w+C. $$
Hence (\ref{double-M}) is verified.

Set
$$g_k(y)=(w(ky)+h_3(ky))/M_k,\quad |y|\le 1. $$

It follows from (\ref{double-M}) , $w_1(0)=0$ and $h_3=O(1)$ that
$$\max_{B_1}g_k\to 1,\quad \max_{B_{1/2}}g_k\ge \frac 14, \quad g_k(0)\to 0, \mbox{ as } \quad k\to \infty. $$
Then we observe that the Harnack inequality holds for $1-g_k$ in $B_1$ because $1-g_k\ge 0$ in $B_1$ and
$$\hat a_{ij}(k\cdot)\partial_{ij}(1-g_k)=0 \quad \mbox{ in }\quad B_1. $$
The uniform ellipticity of $(\hat a_{ij})$ implies that
\begin{equation}\label{har-f}
\max_K(1-g_k)\le C(K)\min(1-g_k),\quad \forall K\subset\subset B_1.
\end{equation}
Consequently $g_k$ converges in $C^{\alpha}$ norm to $g$ in $B_1$ for some $\alpha>0$ depending only on the upper bound and the lower bound of the eigenvalues of $(\hat a_{ij})$. Since the constant in (\ref{har-f}) does not depend on $k$, we also have
$$\max_K(1-g)\le C(K)\min_K(1-g),\quad \forall K\subset\subset B_1. $$

Next we observe that (\ref{lconv}) also holds for all $e\in E$ with $\|e\|\le \frac 34k$, thus
$$g \mbox{ is concave in } B_1 $$
because the perturbation term disappears as $k\to \infty$.
Let $l$ be a linear function that touches $g$ from above around $0$ in $B_{1/2}$. Then $l-g$ is a convex function that takes its minimum at the origin. Since $l-g$ satisfies Harnack inequality, $l-g\equiv 0$ in $B_1$ (by the $C^{\alpha}$ convergence from $g_k$ to $g$, one can find $l_k\to l$ such that $l_k-g_k\ge 0$ in $B_{3/4}$ and $(l_k-g_k)(0)\to 0$. Thus applying Harnack inequality to $l_k-g_k$ we see that $l-g$ satisfies Harnack inequality too.

Since $\max_{B_1}g\ge \frac 34$ we see that $l=a\cdot x$ for $a\neq 0$, $|a|\ge 3/4$.  In other words
$$\frac{w(ky)}{M_k}-ay=o(1),\quad |y|\le \frac 34. $$
 Let $e\in E$ be a direction that $w(ke)\ge \frac 12 M_k k\|e\|$ for $k$ large. Subtract a linear function from $w$ ($w_1=w-$ the linear function) to make
$$w_1(0)=w_1(e)=0. $$
Clearly
\begin{equation}\label{w1-a}
w_1(ke)\ge \frac 14 M_k k\|e\|
\end{equation}
for $k$ large. Using
$$\Delta_e^2w_1(x)\le c_0|x|^{-\beta_1}, \quad |x|\ge 1 $$
we have
$$w_1(2e)\le 2w_1(e)-w_1(0)+c_0\|e\|^2\|e\|^{-\beta_1}=c_0 \|e\|^{2-\beta_1}. $$
Similarly
\begin{align*}
w_1(4e)&\le 2w_1(2e)-w_1(0)+c_0\|2e\|^2\|2e\|^{-\beta}\\
&\le 2c_0\|e\|^{2-\beta_1}+c_0(2\|e\|)^{2-\beta_1}.
\end{align*}

In general
\begin{align}\label{w1-b}
&w_1(2^{N+1}e)\\
\le &\|e\|^{2-\beta_1}c_0(2^{N}+2^{N-1+(2-\beta_1)}+2^{N-2+2(2-\beta_1)}+...+2^{N(2-\beta_1)}) \nonumber\\
= &\|e\|^{2-\beta}c_0 2^N\frac{1-2^{-(N+1)(\beta_1-1)}}{1-2^{1-\beta_1}}\nonumber
\end{align}
if $\beta_1>1$. Here we recall that $\beta_1=\min\{n-2,\beta\}$. Thus for $n\ge 4$ we obviously have
$$w_1(2^{N+1}e)\le C\|2^{N+1}e\|, $$
which leads to a contradiction (\ref{w1-a}). Therefore we have proved that $M_k$ is bounded if the Monge-Amp\`ere equation is defined in $\mathbb R^n$ for $n\ge 4$.
Since $w+h_3$ is bounded from above, it is easy to see that $w+h_3\equiv C$ by Harnack inequality.
 Theorem \ref{thm1} is established for $n\ge 4$.

 \medskip

 \noindent{\bf Step five: Proof of Theorem \ref{thm1} for $n=3$}

 \medskip

 In $\mathbb R^3$ we have $\beta_1=1$ which means (\ref{w1-b}) becomes
$$w_1(2^{N+1}e)\le CN(2^{N+1}\|e\|) $$
which leads to a contradiction if we assume $M_k\ge (\log k)^2$.
Therefore the same argument has proved that
$$w(x)+h_3(x)\le C(\log (2+|x|))^2, \quad x\in \mathbb R^3. $$

In order to finish the proof we need the following

\begin{lem}\label{har-4}
Let $u$ solve
$$a^*_{ij}\partial_{ij}u=0,\quad \mbox{in }\quad \mathbb R^n $$
where $\lambda I\le (a^*_{ij}(x))\le \Lambda I$ for all $x\in \mathbb R^n$ and
$$|u(x)|\le C(1+|x|)^{\delta}, \quad x\in \mathbb R^n. $$
There exists $\epsilon_0(\lambda,\Lambda)>0$ such that if $\delta\in (0,\epsilon_0)$
$u\equiv constant$.
\end{lem}

\noindent{\bf Proof of Lemma \ref{har-4}}:
For $R>1$ let
$$u_R(y)=u(Ry)/R^2, \quad |y|\le 1. $$
Then $u_R$ satisfies
$$a^*_{ij}(Ry)\partial_{ij}u_R=0,\quad \mbox{ in }\quad B_1 $$
and
$$|u_R(y)|\le CR^{\delta-2}, \quad \mbox{in }\quad B_1. $$
By Krylov-Safonov's \cite{safonov} estimate
\begin{equation}\label{safa-k}
\frac{|u_R(y)-u_R(0)|}{|y|^{\epsilon_0}}\le CR^{\delta}/R^2,\quad \forall y\in B_{1/2}
\end{equation}
where $\epsilon_0>0$ only depends on $\lambda$ and $\Lambda$.
Clearly (\ref{safa-k}) can be written as
$$|u(x)-u(0)|\le C|x|^{\epsilon_0}R^{\delta-\epsilon_0}, \quad \forall x\in B_{R/2}. $$
Fix any $x\in \mathbb R^n$, we let $R\to \infty$, then $u(x)=u(0)$ for $\delta<\epsilon_0$.
Lemma \ref{har-4} is established. $\Box$.

\medskip

Applying Lemma \ref{har-4} to $w+h_3$ we see that $w+h_3\equiv c$ in $\mathbb R^3$. Theorem \ref{thm1} is also proved for the equation defined in $\mathbb R^3$. $\Box$

config: arxiv


text:
\begin{document}
\title{Are good local minima wide in sparse recovery?}

\author{\IEEEauthorblockN{Michael Moeller and Otmar Loffeld}
\IEEEauthorblockA{University of Siegen\\
H\"olderlinstra\ss e 3, 57076 Siegen}

\and
\IEEEauthorblockN{J\"urgen Gall}
\IEEEauthorblockA{University of Bonn\\
Endenicher Allee 19a, 53115 Bonn}

\and
\IEEEauthorblockN{Felix Krahmer}
\IEEEauthorblockA{Technical University of Munich\\
Boltzmannstr. 3, 85748 Garching b. M\"unchen}}

\maketitle

\newcommand{\MM}[1]{{\color{red}\textbf{#1}}}
\begin{abstract}

The idea of compressed sensing is to exploit representations in suitable (overcomplete) dictionaries that allow to recover signals far beyond the Nyquist rate provided that they admit a sparse representation in the respective dictionary. The latter gives rise to the sparse recovery problem of finding the best sparse linear approximation of given data in a given generating system. In this paper we analyze the iterative hard thresholding (IHT) algorithm as one of the most popular greedy methods for solving the sparse recovery problem, and demonstrate that systematically perturbing the IHT algorithm by adding noise to intermediate iterates yields improved results. Further improvements can be obtained by entirely rephrasing the problem as a parametric deep-learning-type of optimization problem. By introducing perturbations via dropout, we demonstrate to significantly outperform the classical IHT algorithm, obtaining $3$ to $6$ times lower average objective errors.
\end{abstract}
\IEEEpeerreviewmaketitle

\section{Introduction}
\subsection{Sparse Recovery}
In many applications, the number of possible measurements is limited by physical or financial constraints on the measuring system. The field of compressed sensing has demonstrated that taking significantly fewer (linear) measurements than the number of unknowns still allows for exact recovery when assuming sparsity in a suitable (overcomplete) dictionary. By combining the measurement matrix and the dictionary into a single linear operator $A\in \mathbb{R}^{m \times n}$, the sparse recovery problem can be stated as the solution of the minimization problem
\begin{align}
\label{eq:sparseRecovery}
\min_u \|Au-f\|^2 \qquad \text{s.t.} \quad |u|_0 \leq s,
\end{align}
for $f \in \mathbb{R}^m$ being measured data, $|u|_0$ denoting the number of nonzero entries in $u$, and $s$ being a known (or estimated) level of sparsity of the desired solution.

Although \eqref{eq:sparseRecovery} is NP-hard in general (see e.g. \cite{npHard}), several greedy algorithms as well as convex relaxation approaches guarantee to find the global minimizer provided that $A$ meets certain regularity conditions such as the restricted isometry property (RIP), and often yield faithful approximate solution far beyond the theoretical guarantees, see \cite{Foucart13} for details.

In general, a popular approach for solving constrained optimization problems of the form
\begin{equation}
\label{eq:constrainedOpti}
\min_u E(u) \qquad \text{s.t.}  \qquad u\in M,
\end{equation}
for some set $M$ and an energy $E$ whose gradient $\nabla E$ is $L$-Lipschitz continuous, is to iteratively minimize a majorizer of $E$, giving rise to the so-called \textit{gradient projection} algorithm
\begin{equation*}
 u^{k+1} \in \arg\min_{u\in M} E(u^k) + \langle \nabla E(u^k), u-u^k \rangle + \frac{1}{2\tau}\|u -u^k\|_2^2,
\end{equation*}

for $\tau \in ]0,\frac{1}{L}[$. Convergence follows from the majorization-minimization framework under weak assumptions, see e.g. \cite{Hunter04,Mairal13,Sun17}.

Interestingly, although the set $M:= \{u\in \mathbb{R}^m~|~ |u|_0 \leq s\}$ is not convex, a (possibly not unique) projection can still be computed efficiently by simply keeping the $s$ largest entries in magnitude and setting the remaining entries to zero. Denoting this \textit{hard thresholding} operation by $H_s$, and applying the gradient projection algorithm to \eqref{eq:sparseRecovery} yields the \textit{iterative hard thresholding} (IHT) algorithm \cite{Blumensath2008,Blumensath2009},
\begin{align}
\label{eq:ihs}
u^{k+1} = H_s(u^k-\tau A^T(Au^k-f)).
\end{align}
Despite its simplicity, the IHT algorithm is one of the most popular sparse recovery algorithms, and is among the greedy algorithms with the least restrictive assumptions for exact recovery (\cite[p. 29]{Foucart13}).

\begin{figure}[t!]
\includegraphics[width=0.49\textwidth]{energyIllustration.png}
\caption{Illustrating the energy landscape of \eqref{eq:sparseRecovery} for $m=2$ and $s=1$: The semi-transparent convex surface represents $\|Au-f\|^2$, which, however, is restricted to a feasible set of one-sparse solutions represented by the union of the orange and blue curves. While the red cross illustrates the global minimum, one can see that there is a local minimum on the orange curve, too. In this paper we investigate possible strategies to avoid getting stuck on such (inferior) local minima.}
\label{fig:energyLandscape}
\end{figure}

In general, we have to expect the energy landscape of \eqref{eq:sparseRecovery} to have many local minima on the feasible set, as we have illustrated for the simple case of $m=2$ and $s=1$ in Figure \ref{fig:energyLandscape}. Therefore, an interesting question is if there are systematic strategies that tend to avoid getting stuck in bad local minima?

Of course, one cannot expect a positive answer to this question for arbitrary nonconvex problems. However, let us assume/conjecture that good local minima tend to be wide for \eqref{eq:sparseRecovery} in many practically relevant settings, e.g. for random measurement matrices. As a consequence of the majorization-minimization framework, the IHT algorithm monotonically decreases the objective and provably converges to a local minimizer, see \cite{Blumensath2008}. If one perturbs the locally optimal solution several times and the IHT algorithm converges to the same locally optimal solution, one can conjecture to be in a wide local minimizer of the energy landscape, which - according to our assumption - has a comparably low energy. Vice versa, perturbing a locally optimal solution in a narrow minimum would likely lead to a different solution when continuing to optimize with IHT.

With the above conjecture and intuition in mind, this manuscript investigates the behavior of two different strategies for obtaining perturbed IHT algorithms - one straight-forward algorithm that adds Gaussian noise to the IHT solution after a fixed number of iterations, and one algorithm that parameterizes the solution in terms of the IHT updates themself and adapts the idea of dropout \cite{dropout} from the world of deep learning.

After summarizing some related work in Section \ref{sec:related}, we present the perturbed IHT methods in Sections \ref{sec:nIHT} and \ref{sec:lIHT}. We demonstrate, in Section \ref{sec:numerics}, that the variants indeed exhibit systematically better performance as compared to the classical IHT  and draw (preliminary) conclusions in Section \ref{sec:conclusions}.

\section{Related Work}
\label{sec:related}
\subsection{Sparse recovery}
As the central problem in the field of compressed sensing, the solution of \eqref{eq:sparseRecovery} has been studied in detail in the literature. While the exact solution to \eqref{eq:sparseRecovery} is known to be NP-hard to compute in general \cite{npHard}, many greedy strategies, e.g. orthogonal matching pursuit (OMP) \cite{omp}, CoSaMP \cite{cosamp}, the IHT algorithm \cite{Blumensath2008} discussed in the introduction, or hard thresholding pursuit (HTP) \cite{htp}, as well as $\ell^1$ minimizing convex relaxation approaches \cite{basisPursuit} are known to yield exact sparse recovery under certain conditions usually related to the restricted isometry property (RIP), see \cite{Foucart13} for an overview. Recent work on partial hard thresholding \cite{partialHT} provides a novel framework including new algorithms and nicely summarizes some recent exact recovery results.

In this work we are less interested in exact recovery guarantees, but rather ask the question how close one can get to global minimizers of the nonconvex problem \eqref{eq:sparseRecovery} in practice. For this purpose we focus on perturbations of the IHT algorithm. Since we not only consider the straight-forward approach of adding noise to the iterates, but also a parametric approach based on deep learning ideas, we briefly summarize the main necessary concepts.

\subsection{Deep learning}
The core idea of deep learning is to use a parametric function $\mathcal{N}(x;\theta)$ that maps the input $x$ to the desired output. The \textit{network} $\mathcal{N}$ ideally has to be parameterized by the \textit{weights} $\theta$ in such a way that all desired outputs, but as few (undesirable) other elements as possible, lie in the range of the network.

The simplest and most generic architectures are \textit{fully connected networks}, for which $\mathcal{N}$ is given by a composition
\begin{equation}
\label{eq:mlp}
\mathcal{N}(x;\theta) = \phi(\sigma(\hdots \phi(\sigma(\phi(x;\theta_1));\theta_2)\hdots);\theta_L)
\end{equation}
of affine linear transfer functions
\begin{equation}
\label{eq:linearTransfer}
 \phi(x;\theta_i) = \theta_i^W x + \theta_i^b ,
\end{equation}
and nonlinearities $\sigma$, e.g. rectified linear units $\sigma(x) = \max(x,0)$. The overall parameters $\theta$ consist of the parameters $\theta_i$ of each \textit{layer}, which themselves typically divide into a weight matrix $\theta_i^W$ and a \textit{bias} $\theta_i^b$.

 Once an appropriate architecture for the network $\mathcal{N}$ has been chosen, the \textit{training} consists of an optimization problem
\begin{equation}
\label{eq:learning}
 \hat{\theta} \in \arg\min_\theta \sum_{i \in \text{training set}} \mathcal{L}(\mathcal{N}(x_i;\theta),y_i) ,
 \end{equation}
in which one tries to determine parameters $\hat \theta$ that yield the lowest \text{loss} $\mathcal{L}$ on the \textit{training data}, typically consisting of pairs $(x_i,y_i)$ of inputs $x_i$ and desired outputs $y_i$ of the network.
We refer the reader to \cite{deepLearningBook} for a more detailed introduction.

To obtain a sufficiently expressive network many researchers have turned to deeply nested functions $\mathcal{N}$ with such a large number of weights $\theta$, that the training data can often be fitted almost perfectly. In order to regularize the training \eqref{eq:learning} of the network and prevent \textit{overfitting}, a very successful strategy is to introduce a \textit{dropout layer} \cite{dropout}, i.e. a random perturbation, into the training. More precisely, a dropout layer sets a certain fraction of the entries of an input vector to zero at random positions, and leaves the others untouched. This strategy has proven to be very efficient in preventing overfitting, and avoiding narrow local minima in the energy landscape of the training \eqref{eq:learning}, which makes it an interesting approach for our purposes.

Finally, a common strategy for learning approaches inspired by variational or partial differential equation based approaches is to roll-out suitable algorithms and treat certain parts of the latter as learnable parameters of the resulting network -- see e.g. \cite{Zheng15,learningRDE, Richard17}
for details.
\subsection{Dithering}
Another area where artificially introducing perturbations has proven useful is in analog-to-digital conversion, in this context this approach is known under the name of dithering. See for example \cite{XJ18} for a recent analysis of its benefits.
\section{Perturbing IHT}
As an additional motivation for a perturbed IHT, we do the following test in 2D, i.e. for $n=m=2$ and $s=1$. We draw a random measurement matrix $A$ with columns of norm two, and random data $f$ with norm one. We run the IHT algorithm from $81\times 81$ different initial positions in $[-1,1]^2$ with step size $\tau = 0.05/\|A\|^2$. Figure \ref{fig:2dExample} illustrates the information such a run contains: One (typically) obtains two minima with corresponding regions in which one has to start to converge to the respective point. The conjecture the perturbed IHT algorithms are based on is that the suboptimal local minimum is closer to the region that converges to the global minimum than the global is to the region which converges to the suboptimal local one. The latter would make a perturbed algorithm more systematically more successful. While we picked a case where our assumption is met for illustration purposes in Figure \ref{fig:2dExample}, we ran the above random setting 1000 times, found 890 settings with two local minima in $[-1,1]$ and computed the average distances discussed above: Indeed the average distance of the global minimum to the local region was 0.254 while the distance of the suboptimal local minimum to the global region was 0.144. Of course, the above test is very simplistic and we had to pick a relatively small $\tau$ in order to arrive local minima more often. Nevertheless, we hope for the intuition to carry over to the significantly more difficult recovery problems in higher dimensions. The following section discusses two variants of perturbed IHT algorithms.

\begin{figure}[t!]
\begin{center}
\includegraphics[trim={0cm 0.5cm 0cm 0.5cm}, clip,width=0.49\textwidth]{2dExampleSimple2.png}
\caption{Running IHT in a simple 2D setting often yields cases in which the region of starting points that converge to the global minimum is significantly closer to the suboptimal local minimum than the global minimum is to the region of points from which IHT converges to the suboptimal local minimum.}
\label{fig:2dExample}
\end{center}
\end{figure}

\subsection{Adding noise}
\label{sec:nIHT}
The first (straight-forward) way to perturb the IHT algorithm is to add noise after a fixed number of iterations. For this purpose, we run the IHT algorithm for 600 iterations before adding zero-mean Gaussian noise with standard deviation $\sigma = 0.025$ to its outcome, which is subsequently used as the input to another 600 IHT iterations. We stop the process after the fifth run of the IHT algorithm. While this procedure can surely be optimized in terms of the hyperparameters, the number of times the system is perturbed, and the tracking of the best solution among all reinitializations, the goal of our work is to see if this simple perturbation strategy allows to demonstrate a systematic improvement over the plain IHT algorithm (which we run for $5\cdot 600=3000$ iterations for the sake of fairness in our numerical experiments). We refer to this variant as the \textit{noisy IHT algorithm}.

\subsection{Optimizing the algorithm with dropout}
\label{sec:lIHT}
An entirely different approach is roll out the IHT algorithm and treat the updates themselves as learnable parameters: Note that $L$ iterations of the IHT algorithm can be written as
\begin{equation}
\label{eq:rolledOutIht}
u^L = H_s(\phi( \hdots H_s(\phi(H_s(\phi(u^0;\theta_1));\theta_2))\hdots ; \theta_L)),
\end{equation}
for $\phi$ denoting the linear transfer functions defined in equation \eqref{eq:linearTransfer}, and the parameters $\theta_i$ taking the specific form
\begin{equation}
\label{eq:paramInit}
\theta_i^W = Id - \tau A^TA, \quad \theta_i^b = \tau A^Tf, \quad \forall i\in \{1,\hdots,L\}.
\end{equation}
In full analogy to \eqref{eq:mlp} one could now treat $u^L =: \mathcal{N}(u^0;\theta)$ as a neural network
and try to optimize
\begin{equation}
\label{eq:optimization}
\min_\theta \|A\mathcal{N}(u^0;\theta) -f \|_2^2.
\end{equation}
Note, however, that using \eqref{eq:paramInit} as an initialization for the weights along with a deep network leads to the initial parametrization yielding a local minimizer such that gradient based weight optimization can never progress. Moreover, each layer of the parametrization \eqref{eq:rolledOutIht} has $n(n+1)$-many free parameters, making it expensive to optimize deep architectures.

We therefore fix the first 3000 iterations to the noisy IHT algorithm and merely use two further IHT iterations as our parametrization. To avoid the problem of starting at a local minimum, we pick up the perturbation idea and introduce an intermediate dropout layer. More precisely, we consider the architecture
\begin{equation}
\label{eq:proposed}
\mathcal{N}(u^0; \theta) = H_s(\phi(H_s(\text{drop}(\phi(u^0;\theta^1);5\%));\theta^2)),
\end{equation}
use the result of the noisy IHT algorithm as $u^0$, initialize according to \eqref{eq:paramInit}, and optimize \eqref{eq:optimization}. The layer $\text{drop}(x; 5\%)$ denotes randomly setting $5\%$ of the entries of $x$ to zero and is only used during the optimization, not for the final prediction. We refer the reader to \cite{dropout} for details on such a \textit{dropout}.

For the optimization itself we utilize a subgradient-descent-type algorithm with momentum as commonly used in learning applications. Similar to the way rectified linear or max-pooling units are optimized in the deep-learning literature, we ignore the set of non-differentiable points in the hard-thresholding operator $H_s$ and use
$$ (\nabla H_s(x))_i = (|(H_s(x))_i|>0)\cdot x_i $$
as a (sub-)derivative. Practically, we use Matlab's deep learning framework, set the momentum to $0.9$, the step size (learning rate) to $10^{-4}$, and run the algorithm for 2000 iterations. Despite the lack of convergence guarantee for such an algorithm, it successfully reduced the objective value in all test cases. We refer to this approach as the \textit{parametric IHT method}.

Note that the parametric IHT method is not a learning based approach to sparse recovery, but merely uses a particular parametrization of the solution space which exploits some concepts borrowed from the deep learning literature. Similar ideas of regularization by reparametrization have recently been made for image reconstruction in \cite{deepImagePrior}, where -- similar to the idea of \eqref{eq:optimization} -- a deep network based architecture is used to parameterize natural images.

It is worth pointing out that the bias of the last linear transfer function is of the same size as the desired solution already, such that \eqref{eq:proposed} represents a significant overparametrization, which on the other hand also appears to be the reason why the method outperforms classical optimization algorithm such as IHT significantly, as we will see in the next section. Similar to other areas of nonconvex optimization, e.g. \cite{Zach14}, increasing the dimensionality of the underlying problem seems to help in avoiding local minima even before entering a setting of convex lifting approaches.

\begin{figure*}[tbh]
\includegraphics[width=0.32\textwidth]{errIht.png}
\includegraphics[width=0.32\textwidth]{errNoisyIht.png}
\includegraphics[width=0.32\textwidth]{errParamIht.png} \\
\centering Mean of $\|Au-f\|^2$ over all 20 runs for each method. \\
$ $ \\
\includegraphics[width=0.32\textwidth]{errHistIht.png}
\includegraphics[width=0.32\textwidth]{errHistNoisyIht.png}
\includegraphics[width=0.32\textwidth]{errHistParamIht.png} \\
\centering Count how often $\|Au-f\|^2>0.03$ for each method within the 20 runs. \\
$ $ \\
\includegraphics[width=0.32\textwidth]{recErrIht.png}
\includegraphics[width=0.32\textwidth]{recErrNoisyIht.png}
\includegraphics[width=0.32\textwidth]{recErrParamIht.png} \\
\centering  Mean of $\frac{\|u-u_{gen}\|^2}{\|u_{gen}\|^2}$ for each method over the 20 runs.\\
\caption{Comparison of the average objective error $\|Au-f\|^2$ (first row), the number of times the objective error exceeded $0.03$ among 20 runs (middle row), and the relative difference $\|u-u_{gen}\|^2/\|u_{gen}\|^2$ to the element that generated the data $f=Au_{gen}$.}
\label{fig:comparison}
\end{figure*}

\section{Numerical results}
\label{sec:numerics}
We investigate the behavior of the IHT, the noisy IHT, and the parametric IHT algorithms by generating matrices $A\in \mathbb{R}^{m\times n}$ with entries sampled from a Gaussian distribution, using a fixed $n=200$, and varying $m\in \{50,60,70,\hdots,120 \}$. We generate data $f=Au_{gen}$ for $u_{gen}$ having $s=\text{round}(\mu \cdot n)$ many nonzero entries also drawn from a Gaussian distribution. We vary the relative sparsity $\mu \in \{0.025,0.05,0.075,\hdots, 0.2\}$, and finally normalize to $\|f\|_2=1$. We then run each of the algorithms with the number $s=|u_{gen}|_0$ of nonzero entries we used to generate the data with, and measure their performance in terms of the objective error $\|Au-f\|^2$. To avoid the issue of possibly having taken too few iterations, we solve one final least-squares problem on the support of the solution returned by each of the algorithms.

The first row of figure \ref{fig:comparison} shows the average objective error over 20 realizations of the IHT, the noisy IHT, and the parametric IHT approaches. The second row shows the count how often (among the 20 overall runs) the objective error exceeded a value of $0.03$, and, for the sake of interpretability, the third row shows the normalized difference to the generating element $\|u-u_{gen}\|^2/\|u_{gen}\|^2$ for each of the three methods.

We conducted similar experiments with measurement matrices $A$ arising from a Bernoulli-sampling, and subsampled discrete cosine transform matrices. Since the results were qualitatively similar, we are merely showing the Gaussian case for the sake of brevity.

As an additional summary, Table \ref{tab:comparison} shows the average objective error of all three methods over all experiments of Figure \ref{fig:comparison}.
As we can see, both perturbations of the IHT algorithm led to significant improvements over the classical IHT with the overall average objective error being improved by a factor of $3.4$ by the noisy IHT and by a factor of $6$ by the parametric IHT. While the number $2n(n+1)$ of free parameters of the parametric IHT approach is notably larger than the $n$ parameters of the original unknown (leading to a runtime of about 7-8 seconds on a CPU in the above test), the resulting optimization with dropout is significantly less prone to getting stuck in local minima and yields surprisingly large improvements.

\newpage
As we can see in Figure \ref{fig:comparison}, a small support of $u_{gen}$ along with a small number of measurements caused relatively large objective errors of the classical IHT, which could be improved significantly by its perturbed variants. Since the histogram graphs (second row in Figure \ref{fig:comparison}) look quite similar to the mean objective errors of the first row, we can conclude that the errors of IHT are not just based on a single unlucky random realization, but yield systematic problems with up to $14$ failures to reduces the objective below $0.03$.

Note that the goal of this work is to investigate the minimization of $\|Au-f\|$ under sparsity constraints rather than the exact recovery of the generating element $u_{gen}$. Nevertheless, we included the normalized distance to $u_{gen}$ for the sake of completeness in the last row of Figure \ref{fig:comparison}. We would, however, like to point out that settings with few measurements and large supports, e.g. $m=50$ and $s=40$ (upper right part of the plots in the last row of Figure \ref{fig:comparison}), naturally imply that we can neither expect $u_{gen}$ to be the unique element to meet $f=Au$, nor the sparsest. This means that although the noisy and parametric IHT gave very high values of $\|u-u_{gen}\|^2/\|u_{gen}\|^2$ in such cases, they do provide desirable solutions with almost no objective errors. Moreover, differing from $u_{gen}$ in only one or two components can already lead to comparably high objective errors in the case of few measurements as we can see in the IHT graphs in the upper left. Nevertheless, in settings of sparse $u_{gen}$ and moderately many measurements, the last row of Figure \ref{fig:comparison} also indicates that the perturbed IHT algorithm have a higher success rate in recovering $u_{gen}$.

\begin{table}[tb]
\caption{Average objective error $\|Au-f\|^2$ of the classical IHT, the noisy IHT, and the parametric IHT algorithms, over all relative sparsities and all number of measurements from Figure \ref{fig:comparison}.}
\label{tab:comparison}
\def\arraystretch{1.5}
\begin{center}
\normalsize
\begin{tabular}{l|ccc}
Method & IHT & Noisy IHT & Param. IHT \\\hline
Avg. $100  \cdot \|Au-f\|^2$&$2.0$ & $0.59$& $0.33$
\end{tabular}
\end{center}
\end{table}

Future research will try to give some theoretical explanations for the behavior we observed in the the numerical studies above. In particular, we will consider the case of noisy measurements and investigate strategies to quantify by what amount the IHT algorithm should ideally be perturbed. Furthermore, an interesting question is the optimal design of the network-like parametrization.

\section{Conclusion}
\label{sec:conclusions}
In this paper we have considered the nonconvex optimization problem of minimizing a quadratic objective subject to a constraint on the maximum number of nonzero entries of the unknown. Based on the assertion that good local minima are wide, we compared the IHT algorithm with a simple modification that perturbs the iterates with noise, as well as with a reparametrization of the solution in a network-like architecture which was perturbed by incorporating dropout. It was shown that the noisy IHT yields significant improvements over the plain IHT algorithm, while the parametric IHT approach involving dropout can further improve the noisy IHT results by a large margin. In summary, our findings support the hypothesis that good local minima are wide for randomly sampled measurement matrices and exact data, and demonstrate that a dropout-based parametrization of the solution space is an excellent approach for exploiting this property.

\bibliographystyle{IEEEtran}
\bibliography{refs}
\end{document}


TITLE: Find the numerically greatest coefficient in binomial expansion of $(3x - 2)^{21}$.
QUESTION [0 upvotes]: My approach :
$$T(r) = \binom{n}{r-1}(3x)^{n-r+1}
          (-2)^{r-1}$$
here $T(r)$ is the $r$th term of the binomial expansion of
$(3x-2)^{21}$.
$$T(r) = \binom{n}{r-1}(3)^{n-r+1}
           (-2)^{r-1}(x)^{n-r+1}$$
here $n = 21$ , so ,
$$T(r) = \binom{21}{r-1}(3)^{22-r}
           (-2)^{r-1}(x)^{n-r+1}$$
So we need $\binom{21}{r-1}(3)^{22-r}(-2)^{r-1}$ to be maximum. I am not able to go further, please help me in this method or provide another method to solve the question in title.

REPLY [1 votes]: Compare the term $3^a2^{21-a}\begin{pmatrix}21\\a\\\end{pmatrix}$ with $3^{a+1}2^{20-a}\begin{pmatrix}21\\a+1\\\end{pmatrix}$. The ratio is $$\frac{3(21-a)}{2(a+1)}.$$
Since this is greater than $1$ for $a\le 12$ the maximum value occurs when $a=12$ i.e. $$3^{12}2^{9}\begin{pmatrix}21\\12\\\end{pmatrix}.$$


text:
\DOC ASM_REAL_ARITH_TAC

\TYPE {ASM_REAL_ARITH_TAC : tactic}

\SYNOPSIS
Attempt to prove goal using basic algebra and linear arithmetic over the reals.

\DESCRIBE
The tactic {ASM_REAL_ARITH_TAC} is the tactic form of {REAL_ARITH}. Roughly
speaking, it will automatically prove any formulas over the reals that are
effectively universally quantified and can be proved valid by algebraic
normalization and linear equational and inequality reasoning. See {REAL_ARITH}
for more information about the algorithm used and its scope. Unlike plain
{REAL_ARITH_TAC}, {ASM_REAL_ARITH_TAC} uses any assumptions that are not
universally quantified as additional hypotheses.

\FAILURE
Fails if the goal is not in the subset solvable by these means, or is not
valid.

\EXAMPLE
This example illustrates how {ASM_REAL_ARITH_TAC} uses assumptions while
{REAL_ARITH_TAC} does not. Of course, this is for illustration only: plain
{REAL_ARITH_TAC} would solve the entire goal before application of {STRIP_TAC}.
{
  # g `!x y z:real. abs(x) <= y ==> abs(x - z) <= abs(y + abs(z))`;;
  val it : goalstack = 1 subgoal (1 total)

  `!x y z. abs x <= y ==> abs (x - z) <= abs (y + abs z)`

  # e(REPEAT STRIP_TAC);;
  val it : goalstack = 1 subgoal (1 total)

    0 [`abs x <= y`]

  `abs (x - z) <= abs (y + abs z)`

  # e REAL_ARITH_TAC;;
  Exception: Failure "linear_ineqs: no contradiction".
  # e ASM_REAL_ARITH_TAC;;
  val it : goalstack = No subgoals
}

\COMMENTS
For nonlinear equational reasoning, use {CONV_TAC REAL_RING} or {CONV_TAC
REAL_FIELD}. For nonlinear inequality reasoning, there are no powerful rules
built into HOL Light, but the additional derived rules defined in
{Examples/sos.ml} and {Rqe/make.ml} may be useful.

\SEEALSO
ARITH_TAC, INT_ARITH_TAC, REAL_ARITH, REAL_ARITH_TAC, REAL_FIELD, REAL_RING.

\ENDDOC

subset_name: curated


TITLE: Proving $\frac1{\sec\phi-\tan\phi} - \frac1{\sec\phi+\tan\phi} =2\tan\phi$
QUESTION [2 upvotes]: Prove
  $$\frac1{\sec\phi-\tan\phi} - \frac1{\sec\phi+\tan\phi} =2\tan\phi$$

I have managed to simplify the lefthand side of the equation to $2\sin\phi /\cos^2\phi$ (maybe incorrect) but cannot seem to finish it off.

REPLY [1 votes]: As $(\sec t+\tan t)(\sec t-\tan t)=1,$
for $\sec t\pm\tan t\ne0\iff\sin t\ne\pm1\iff\cos t\ne0$
$$\dfrac1{\sec t\pm\tan t}=\sec t\mp\tan t$$


TITLE: Formalizing a self referential sentence
QUESTION [1 upvotes]: In The logic of provability, by G. Boolos, we are asked to ponder about this statement:

If this statement is consistent, then you will have an exam tomorrow, but you cannot deduce from this statement that you will have an exam tomorrow.

I tried formalizing it the following way:
Let $p$ stand for the whole statement, and $q$ is going to be true iff I have an exam tomorrow.
Then the paragraph is stating that:
$$
p\leftrightarrow \diamond (p\rightarrow q)\wedge \neg\square(p\rightarrow q)
$$
Solving the fixed point using any of the usual methods will give us that:
$$
GL\vdash \boxdot [p\leftrightarrow A(p,q)]\implies \boxdot[p\leftrightarrow \bot]
$$
Where $A(p,q)$ stands for the right hand side of the first formula.
Therefore, we can conclude by neccesitation and modus ponens that if $p$ and $q$ are suitable chosen statements of $PA$ satisfying the first formula then:
$$
PA\vdash \neg p
$$
Therefore we conclude that the statement is false under every interpretation of $PA$, and therefore it gives no information about $q$.
Is my reasoning correct? Is my conclusion pertinent?

REPLY [1 votes]: Take any consistent theory $T$ satisfying the provability conditions and the fixed point lemma.
Take any sentence $q$ over $T$.
$\def\imp{\rightarrow}$
$\def\eq{\leftrightarrow}$
$\def\box{\square}$
$\def\diam{\lozenge}$
Your equation
By the fixed point lemma there is a sentence $p$ over $T$ such that $T \vdash p \eq \diam( p \imp q ) \land \neg \box( p \imp q )$.
Equivalently $T \vdash p \eq \neg \box( p \land \neg q ) \land \neg \box( \neg p \lor q )$.
Note that $T \vdash \box \neg p \imp \box( p \imp q )$ by (D1) and (D2).
Thus $T \vdash p \imp \neg \box \neg p$ and hence $T \vdash \box \neg p \imp \neg p$.
Therefore by Lob's theorem $T \vdash \neg p$.
Anyway...
Literal question
You translated the sentence wrongly. Literally, $T \vdash p \eq ( \diam p \imp q ) \land \neg \box( p \imp q )$. [$\diam$ before $\imp$!]
Equivalently $T \vdash p \eq ( \box \neg p \lor q ) \land \neg \box( p \imp q )$.
It so happens that the above reasoning works here too, giving $T \vdash \neg p$.
So your first claim is correct though I have no idea whether your reasoning is valid.
But I don't get your reasoning about $q$; how does $T \vdash \neg p$ imply anything about $q$? If you're just saying that the teacher is a liar, then I agree with you.
Probable intended question
The original sentence is ambiguous, and I first interpreted it with "but" binding first before "then".
In this case you want $T \vdash p \eq ( \diam p \imp q \land \neg \box( p \imp q ) )$. [$\diam$ first!]
Equivalently $T \vdash p \eq \box \neg p \lor q \land \neg \box( p \imp q )$.
If $T \vdash \neg p$, then $T \vdash \box \neg p$ by (D1) and hence $T \vdash p$, which is impossible.
Thus the teacher is making a consistent statement over $T$.
At this point it would be erroneous to reason as follows:
  $p$ is consistent. Thus $q$ is true and hence $p \imp q$. So we have deduced from $p$ that $q$ is true.
So this is probably the intended interpretation so that there is this 'paradox'. Once you write in terms of $\box$ it is clear that the error is right at the start, because $T$ is unable to prove $\neg \box \neg p$.