\documentclass[11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{epsfig}
\usepackage{bbm}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {\bf MATH888: High-dimensional probability and statistics } \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{Scribe: #4}{Lecture #1}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}

% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex
%\renewcommand{\baselinestretch}{1.25}

\begin{document}

\lecture{36 --- December 1, 2021}{Fall 2021}{Sebastien Roch, UW-Madison}{Thanasis Pittas}

\section{Overview}


In Lecture 34, we saw an estimator for sparse linear regression with strong error guaranties. However, that estimator was computationally inefficient when the sparsity level was unknown. In Lecture 35 we introduced the  sub-Gaussian sequence model and observed that under the assumption ORT (imposing that $\frac{1}{n}\mathbb{X}\mathbb{X}^T=I_p$), that model is essentially equivalent to that of linear regression, up to a transformation of the data and a scaling of the variance. In this lecture, we give an estimator for the sparse sub-Gaussian sequence model that is computationally efficient even if the sparsity level is unknown. This implies an estimator for sparse linear regression that has the same statistical guarantees as that of Lecture 34 as well as improved computational complexity.

\section{Sparse Sub-Gaussian Sequence Model}

\subsection{Model Definition}
Recall the sparse sub-Gaussian Sequence Model: The observation is a vector $y \in \mathbb{R}^p$ which is generated by the process 
\begin{align*}
    y = \theta^* + \xi \;,
\end{align*}
where $\xi \sim \mathrm{subG}_p(\sigma^2/n)$. Following the notation of the previous lectures, we denote by $\|\theta^*\|_0$ the sparsity of $\theta^*$, that is, $\|\theta^*\|_0:= \sum_{i=1}^p \mathbbm{1}\{\theta^*_i \neq 0  \}$.
The goal is to find an estimator $\hat{\theta} = \hat{\theta}(y)$ such that $\|\hat{\theta} - \theta^* \|_2$ is small with high probability.

\subsection{Hard Thresholding Estimator}
Since every entry of  $y$ is a perturbed version of the corresponding entry of $\theta^*$ and since the noise $\xi$ is concentrated, we are motivated to threshold the entries of $y$, hoping that this will distinguish the zero entries of $\theta^*$ from the non-zero entries. We thus define the \emph{hard threshold estimator}
\begin{align} \label{eq:hard-thresholding}
    \hat{\theta}_j^H =
    \begin{cases}
    0,  \quad &|y_j| \leq 2\tau \;,\\
    y_j, \quad &|y_j|> 2\tau\;,
    \end{cases} 
\end{align}
where the threshold $\tau$ remains to be specified. Formally, we show the following (also see Theorem 2.11 in \cite{hds}).
\begin{theorem} \label{thm:hard-thr}
Let $C$ be a large enough constant. If $\tau = \sigma \sqrt{(2C/n) \log(2p/\delta)}$, the estimator $\hat{\theta}^H$ from Equation \eqref{eq:hard-thresholding} satisfies
\begin{align*}
    \| \hat{\theta}^H - \theta^* \|_2^2 \lesssim \frac{\sigma^2}{n} \| \theta^*\|_0 \log\left( \frac{2p}{\delta}  \right)\;,
\end{align*}
with probability at least $1-\delta$.
\end{theorem}

In comparison to the estimator from two lectures ago, we note that $\hat{\theta}^H$ is computationally efficient while achieving the same (optimal) error rate.


\begin{proof}
 We start with the following claim, saying that the selection of $\tau$ is such that all additive errors are small.
 \begin{claim}\label{claim:firstcl}
 Let $\mathcal{A}=\{ |\xi_j| \leq \tau \;\forall j \}$. Then $\mathbb{P}[A^c] \leq \delta$.
 \end{claim}
 \begin{proof}[Proof of Claim \ref{claim:firstcl}]
  By writing the complement of $\mathcal{A}$ as a union and using union bound, we have that
  \begin{align*}
      \mathbb{P}[A^c] = \mathbb{P}\left[ \bigcup_{j=1}^p \{ |\xi_j| > \tau\}  \right]
      \leq  \sum_{j=1}^p \mathbb{P}\left[  \{ |\xi_j| > \tau\}  \right]
      \leq 2 p \exp\left(-\frac{\tau^2}{C 2(\sigma^2/n)} \right) \leq \delta\;,
  \end{align*}
  where for the second  inequality we use that $\xi \sim \mathrm{subG}(\sigma^2/n)$, which means that if $e_j$ denotes the $j$-th vector of the orthonormal basis of $\mathbb{R}^p$, then $\xi_j=\langle e_j,\xi \rangle \sim \mathrm{subG}(\sigma^2/n)$. The last inequality uses $\tau = \sigma \sqrt{(2C/n) \log(2p/\delta)}$.
 \end{proof}
 Next, we note that under the event $\mathcal{A}$, the contribution to the error from all coordinates is small.
 \begin{claim} \label{claim:secondcl}
 Under the event $\mathcal{A}=\{ |\xi_j| \leq \tau \; \forall j \}$, the following hold for all $j \in \{1,\ldots,p\}$:
 \begin{enumerate}
     \item If $\theta_j^*=0$, then $\hat{\theta}_j^H = 0$.
     \item If $\theta_j^* \neq 0$, then $|\hat{\theta}_j^H - \theta_j^*| \leq 3\tau$.
 \end{enumerate}
 \end{claim}
 \begin{proof}[Proof of Claim \ref{claim:secondcl}]
  For the first case, assume that $\theta_j^*=0$. Then, $y_j=\xi_j$ and thus 
  \begin{align*}
      |y_j| = |\xi_j| \leq \tau \leq 2\tau\;,
  \end{align*}
  which means that $\theta_j^H=0$. For the second case,  $\theta_j^* \neq 0$, we examine two sub-cases:
  \begin{enumerate}
      \item If $|y_j|>2\tau$, then $\hat{\theta}_j^H=y_j$ and thus $|\hat{\theta}_j^H - \theta_j^*| = |\xi_j|\leq \tau \leq 3\tau$.
      \item If $|y_j|\leq 2\tau$, then $\hat{\theta}_j^H=0$ and thus $|\hat{\theta}_j^H - \theta_j^*| = | \theta_j^*| = |y_j - \xi_j| \leq 3\tau$ by the triangle inequality.
  \end{enumerate}
 \end{proof}
 
 Having these two claims at hand, we note that $\|\hat{\theta}^H - \theta^* \|_2^2 = \sum_{j=1}^p (\hat{\theta}_j^H - \theta_j^*)^2 \leq 9\tau^2 \|\theta^*\|_0$. Plugging the value of $\tau$ completes the proof of Theorem \ref{thm:hard-thr}.
\end{proof}


\subsection{Soft Thresholding Estimator}
An alternative way of thresholding the values of $y$ is used in the following estimator, for which we get the same guarantee as in Theorem \ref{thm:hard-thr}.
\begin{align}\label{eq:soft}
    \hat{\theta}_j^S =
    \begin{cases}
    0, \quad &|y_j| \leq 2\tau\;,\\
    y_j-2\tau, \quad &y_j> 2\tau \;, \\
    y_j + 2\tau,\quad  &y_j < -2\tau\;,
    \end{cases} 
\end{align}
or,  in a more compact form,
\begin{align*}
    \hat{\theta}_j^S = \left( 1 - \frac{2\tau}{|y_j|} \right)_+ y_j
\end{align*} 
\begin{theorem} \label{thm:soft-thr}
Let $C$ be a sufficiently large constant. If $\tau = \sigma \sqrt{(2C/n) \log(2p/\delta)}$, the estimator $\hat{\theta}^S$ from Equation \eqref{eq:soft} satisfies
\begin{align*}
    \| \hat{\theta}^S - \theta^* \|_2^2 \lesssim \frac{\sigma^2}{n} \| \theta^*\|_0 \log\left( \frac{2p}{\delta}  \right)\;,
\end{align*}
with probability at least $1-\delta$.
\end{theorem}
The proof of this is the same as that of Theorem \ref{thm:hard-thr}, up to a small modification of the proof of Claim \ref{claim:secondcl}. More specifically, the case $\theta^*_j=0$ remains the same. Regarding the other case, we again have two sub-cases: If $|y_j|>2\tau$, then  $\hat{\theta}_j^S = y_j - 2\tau$ and thus $|\hat{\theta}_j^S - \theta_j^*| = | \xi-2\tau|\leq 3\tau$. The second sub-case does not need any changes.






\subsubsection{LASSO estimator}
The reason that we consider the soft thresholding estimator is that it is the solution of an optimization problem of the following nice form. 
\begin{align*}
    \hat{\theta}^S = \arg\min_{\theta \in \mathbb{R}^p} \left\{ \|y-\theta\|_2^2 + 4\tau \|\theta\|_1 \right\}.
\end{align*}
It can be easily checked that this is an equivallent form of the soft thresholding estimator by writing  $\|y-\theta\|_2^2 + 4\tau \|\theta\|_1 = \sum_{i=1}^p \left((y_i-\theta_i)^2 + 4\tau |\theta_i |  \right)$ and solving the optimization problem for every coordinate separately.

To slightly generalize the above, we define the LASSO estimator
\begin{align*}
    \hat{\theta}^{LASSO} = \arg\min_{\theta \in \mathbb{R}^p} \left\{ \frac{1}{n}\|y-\mathbb{X}\theta\|_2^2 + 4\tau \|\theta\|_1 \right\}.
\end{align*}
Based on the discussion of the last lecture, under the assumption ORT, the above estimator is equivalent to the soft thresholding estimator. However, the LASSO estimator is natural to consider even in the absence of that assumption. 

Because of its particular form,  there are algorithms for computing $\hat{\theta}^{LASSO}$ numerically. One such method is coordinate descent:

\begin{enumerate}
    \item Until Convergence:
    \begin{enumerate}
        \item For $j=1,\ldots, p$:
        \begin{enumerate}
            \item $\theta_j \gets R_j\left(1-2\tau/|R_j| \right)_+$, where $R_j = \mathbb{X}_j^T\left(y-\sum_{k\neq j}\theta_k \mathbb{X}_k\right)$.
        \end{enumerate}
    \end{enumerate}
\end{enumerate}

Using standard optimization results, it can be shown that this method converges to the true solution. It remains to show statistical guaranties for $\hat{\theta}^{LASSO}$, which will leave for the next lecture.


\bibliographystyle{alpha}

\begin{thebibliography}{77}

\bibitem{hds}
Philippe Rigollet and Jan-Christian H\"{u}tter,
\emph{18.657: High Dimensional Statistics Lecture Notes}, MIT, 2017.

\end{thebibliography}

\end{document}