\documentclass{article}

\usepackage[left=1.25in,top=1.25in,right=1.25in,bottom=1.25in,head=1.25in]{geometry}
\usepackage{amsfonts,amsmath,amssymb,amsthm}
\usepackage{verbatim,float,url,enumerate}
\usepackage{graphicx,subfigure,psfrag}
\usepackage{natbib}
\usepackage{environ}
\usepackage{listings}
\usepackage[dvipsnames]{xcolor}
\usepackage{hyperref}
\usepackage{pifont}       
\usepackage{bbm}
\usepackage[ruled,linesnumbered]{algorithm2e}

\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newtheorem{corollary}{Corollary}

\theoremstyle{remark}
\newtheorem{remark}{Remark}
\theoremstyle{definition}
\newtheorem{definition}{Definition}

\newcommand{\argmin}{\mathop{\mathrm{argmin}}}
\newcommand{\argmax}{\mathop{\mathrm{argmax}}}
\newcommand{\minimize}{\mathop{\mathrm{minimize}}}
\newcommand{\maximize}{\mathop{\mathrm{maximize}}}
\newcommand{\st}{\mathop{\mathrm{subject\,\,to}}}
\newcommand{\dist}{\mathop{\mathrm{dist}}}
\newcommand{\minimizewrt}[1]{\underset{#1}{\minimize}}   
\newcommand{\maximizewrt}[1]{\underset{#1}{\maximize}}  
\newcommand{\subjectto}{\mbox{subject to}}                       
\newcommand{\ones}{\mathrm 1}                                          
\newcommand{\diag}{\mathop{\rm diag}}                               

\newcommand{\reals}{\mathbb R}
\newcommand{\RR}{\mathbb{R}}
\newcommand{\EE}{\mathbb{E}}
\newcommand{\SSS}{\mathbb{S}}
\newcommand{\prox}{\operatorname{prox}}
\newcommand{\dom}{\operatorname{dom}}
\setlength\parindent{0pt}

\def\R{\mathbb{R}}
\def\E{\mathbb{E}}
\def\P{\mathbb{P}}
\def\Cov{\mathrm{Cov}}
\def\Var{\mathrm{Var}}
\def\half{\frac{1}{2}}
\def\sign{\mathrm{sign}}
\def\supp{\mathrm{supp}}
\def\th{\mathrm{th}}
\def\tr{\mathrm{tr}}
\def\dim{\mathrm{dim}}
\def\hbeta{\hat{\beta}}

\lstset{frame=single,
  aboveskip=3mm,
  belowskip=3mm,
  showstringspaces=false,
  columns=flexible,
  basicstyle={\small\ttfamily},
  numbers=none,
  numberstyle=\tiny\color{gray},
  keywordstyle=\color{blue},
  commentstyle=\color{dkgreen},
  stringstyle=\color{mauve},
  breaklines=true,
  breakatwhitespace=true,
  tabsize=3,
  captionpos=b
}


\usepackage{color}     
\newcommand{\matr}[1]{\begin{bmatrix} #1 \end{bmatrix}}    

\begin{document}

\title{Homework 5}

\author{\Large Convex Optimization 10-725/36-725}
\date{{\bf Due Tuesday November 22 at 5:30pm \\
submitted to Christoph Dann in Gates 8013} \\
(Remember to a submit separate writeup for each \\ 
problem, with your name at the top)\\
\bigskip
\textcolor{blue}{\bf Out of the 4 questions, you are supposed to choose any 3
to solve. You may hand in a solution to all four questions for practice and
receive feedback on all questions but only the scores of 3 questions count
towards your grade. If you do decide to hand in all four problems, please
indicate at the top of each question whether that question should count towards
your grade. If you give no such indication, the \emph{best score will be
discarded (= worst 3 out of 4)}!} \\
\bigskip Total: 75 points\\
\textcolor{red}{v1.3}} 

\maketitle 



\section{Derivation of DFP update (25 points) [Christoph]}

\begin{enumerate}

\item[(a, 5pts)] Assume $y,s\in \R^n$ are non-zero vectors and $B\in \R^{n\times n}$ is a symmetric matrix.  Show that the solution to the Frobenius norm minimization problem
\[
\begin{array}{rl}
\displaystyle\min_{B^+} & \|B^+ - B\|_F^2\\
\st & (B^+)^T = B^+ \\
& B^+ s = y
\end{array}
\]
is
\[
B^+ = B + \frac{(y-Bs)s^T}{s^Ts} + \frac{s(y-Bs)^T}{s^Ts} - 
 \frac{(y-Bs)^Ts}{(s^Ts)^2} ss^T.
\]
{\em Hint:}  show that the KKT conditions hold for a suitable choice of multipliers for the contraints.  Use a low-rank matrix for the first one.


\item[(b, 5pts)] Show via a counterexample that $B^+$ may not be positive definite even if $B$ is symmetric positive definite and $y^Ts > 0$.

{\em Hint:} observe that $B^+$ can be written as 
    \begin{align*}
        B^+ = \left(I - \frac{s s^\top}{s^\top s}\right) B \left(I - \frac{s s^\top}{s^\top s}\right) 
        + \frac{y s^\top}{s^\top s}
        + \frac{s y^\top}{s^\top s}
        - \frac{y^\top s}{(s^\top s)^2} s s^\top.
    \end{align*}

\item[(c, 5pts)] Assume $y,s\in \R^n$ are such that $y^Ts > 0$. Prove that there exists a symmetric and positive definite matrix $M \in \R^{n\times n}$ such that 
\[
Ms= y.
\]
In particular, there exists a non-singular matrix $W \in \R^{n\times n}$ such that $WW^Ts = y$.

        {\em Hint:} \textcolor{red}{use a matrix $M$ that is the sum of the identity and low-rank matrices.}

\item[(d, 5pts)] Assume $y,s\in \R^n$ are such that $y^Ts > 0$ and $B\in \R^{n\times n}$ is a symmetric matrix.  Let $W \in \R^{n\times n}$ be a non-singular matrix such that $WW^Ts = y$.  Show that the solution to the weighted Frobenius norm minimization problem
\[
\begin{array}{rl}
\displaystyle\min_{B^+} & \|W^{-1}(B^+ - B)W^{-T}\|_F^2\\
\st & (B^+)^T = B^+ \\
& B^+ s = y
\end{array}
\]
is the DFP update
\[
B^+ = B + \frac{(y-Bs)y^T}{y^Ts} + \frac{y(y-Bs)^T}{y^Ts} - 
 \frac{(y-Bs)^Ts}{(y^Ts)^2} yy^T = \left(I - \frac{ys^T}{y^T s} \right) B \left(I - \frac{sy^T}{y^T s} \right) + \frac{yy^T}{y^T s}.
\]

{\em Hint:} use part (a) and a suitable change of variables.

\item[(e, 5pts)] Show that if $B$ in the previous question is positive definite then so is $B^+$.

\end{enumerate}

\section{Quasi-Newton Methods in Practice (25 points) [Mariya \& Christoph]}

\begin{enumerate}

    \item[(a, 7pts)]
    In this part, we'll use a quasi-Newton line-search algorithm to minimize the quadratic function:
    \[
    f(x)=\frac{1}{2}x^{T}Qx-b^{T}x,
    \]
    where $Q \in \mathbb{R}^{n\times n}$ is symmetric positive definite. The quasi-Newton line-search
    algorithm can be stated as:   
    
    \begin{algorithm}
     Initialize $x^{0}\in \mathbb{R}^{n}, k=0, B_{0}=I$\;
     \While{$\|\nabla f(x^{k})\|> \epsilon)$ and $(k \leq k_{max})$}{
      $p^{k}:=-(B^{k})^{-1}\nabla f(x^{k})$\;
      $s^{k}:=t_{k}p^{k},${\textit \ where }$t_{k}=\arg\!\min
    \{f(x^{k}+tp^{k}):t \geq 0\}$\;
      $x^{k+1}:=x^{k}+s^{k}$\;
      Update $B^{k+1}$ using $B^{k},\;s^{k},$ and $y^{k}=\nabla
    f(x^{k+1})-\nabla f(x^{k})$\;
      $k:=k+1$\;
     }
    \end{algorithm}   
    
    \begin{enumerate}
        \item[i.]
        Implement this algorithm for both the DFP and BFGS updates (you
    may want to update $H_{k}=B_{k}^{-1}$ instead of $B_{k}$). Use
    $\epsilon=10^{-8}$ and $k_{max} = 100$ quasi-Newton main
    iterations. Initialize $x_0$ as a vector of $n$ normally
    distributed random numbers. Choose $b=0_n$ and construct a random
    positive definite $Q$ according to the following steps: 
        \begin{enumerate}
            \item[1.] Construct a positive definite diagonal matrix $D$ with appropriate entries so that its condition number $\kappa(D) = 10$.
            \item[2.] Construct a random orthogonal matrix $P$.
            \item[3.] Set $Q=P^{T}DP$.
        \end{enumerate}
        \item[ii.]  
        Evaluate your implementation for both updates using the following:  
        \begin{enumerate}
            \item[1.] Three different random matrices $Q \in \mathbb{R}^{2\times 2}$.
            \item[2.] Three different random matrices $Q \in \mathbb{R}^{3\times 3}$.
            \item[3.] Three different random matrices $Q \in \mathbb{R}^{5\times 5}$.
        \end{enumerate}
        For each random Q, report the normalized (divided by the number of elements) Frobenius norm of the difference between the last $B_k$ (or $H_k^{-1}$) and $Q$. Report the final objective value and the number of iterations needed to achieve that value.    
        \item[iii.]  Do the experiments above suggest a general convergence property for the iterates in the special case of convex quadratic functions with exact line search? Can you propose a concrete conjecture?  (Just state a one sentence conjecture.  No need to prove it.)
    \end{enumerate}

    \item[(b, 8pts)]
    In this part, we'll implement a quasi-Newton backtracking line-search algorithm to minimize the Rosenbrock function:
    \[
    f(x_1,x_2) = 100(x_2 - x_1^2)^2 + (1-x_1)^2.
    \]
    This algorithm is very similar to the one above, but instead of exact line-search to determine the step-size $t_k$, we use backtracking. For backtracking, we use the following condition to determine the step-size at each iteration $k$ of the quasi-Newton method: 
    
    \begin{algorithm}
     $t_k = 1$\;
     \While{$f(x^{(k)}+t_k p^{(k)}) > f(x^{(k)}) + 0.0001 t_k\nabla f(x^{(k)})^T p^{(k)}$}{
    $t_k:=0.9t_k$
     }
    \end{algorithm}
     
    \begin{enumerate}
        \item [i.]
        Implement the quasi-Newton backtracking line-search algorithm with $4$ different updates for $B_k$:
    \begin{itemize}
        \item $B^k = I$ for all $k$ (i.e., steepest descent direction).
        \item DFP update.
        \item  BFGS update. 
        \item $B^k = \nabla^2 f(x^{(k)})$ (i.e., Newton direction).
    \end{itemize}
    
        Use $x^{(0)} = \left[ \begin{array}{c} -1.2 \\ 1 \end{array}\right]$ as the initial point. In the second and third case take $B^0 = H^0 = I$. Limit the number of main quasi-Newton iterations to $100$. 
    
        \item[ii.]
         Plot $f^{(k)}$ versus $k$ on a semi-log scale (i.e. where the y-axis is in log scale), where $f^{(k)}$ denotes the objective value obtained by the quasi-Newton method with each of the $4$ updates at each quasi-Newton iteration $k$. Plot in the same graph and include a legend.
    
         Now, we would like to visualize the solution paths (i.e. $x^k$ for each iteration $k$) of the different updates. First, plot the 2D contours of the Rosenbrock function with sufficient number of level sets. In MATLAB, such a plot can be obtained in the following way:
         \begin{verbatim}
             [X,Y] = meshgrid(-2:0.05:2);
             Z = evalRosenbrock([X;Y]);
             contour(X,Y,Z,50)
         \end{verbatim}
         where evalRosenbrock is a user-created function that evaluates the Rosenbrock function.
         Next, plot the solution paths for the $4$ updates on top of the Rosenbrock contour. Plotting without interpolating between different values works best for visualization. Plot the final value of the iterate $x^k$ in larger font so that it is easily visible. Include a legend. 
         
         \item[iii.] Based on the resulting plots, discuss the performance of each of the updates in minimizing the Rosenbrock function.    
    \end{enumerate} 

\item[(c, 10pts)]
In this question, we will implement the limited-memory version of BFGS and apply it to $\ell_2$-penalized logistic regression.
\begin{enumerate}
    \item[i.] Implement the two-loop recursive version of LBFGS
        \begin{itemize}
            \item As line search use a basic Armijo rule of the form 
                \begin{algorithm}
                     $t_k = 1$\;
                     \While{$f(x^{(k)}+t_k p^{(k)}) > f(x^{(k)}) + c_1 t_k\nabla f(x^{(k)})^T p^{(k)}$}{
                    $t_k:=\eta t_k$
                     }
                \end{algorithm}
                with parameters $c_1$ and $\eta$.
          \item For computing the step direction use the two-loop recursive
matrix-vector multiplication from slide 23 of the \textit{Quasi-Newton Lecture
slides}. Use the suggested choice $H^{0,k}$ on this slide. Make sure your
algorithm has run- and memory-complexity $O(mn)$ where $m$ is the length of the
history considered by LBFGS.
          \item Stop the algorithm once at least one of the following three conditions are satisfied:
              \begin{description}
                  \item[Gradient condition:] The gradient satisfies for each component $i \in \{1, \dots n\}$ the condition
             \begin{align}
                \frac{\left|\nabla f(x^{(k)})_i\right| \max\{|x^{(k)}_i|, 1\}}{\max \{f(x^{(k)}), 1\}}  < \epsilon_g
            \end{align}
                    where $\epsilon_g > 0$ is an algorithm parameter.
                  \item[Step condition:] The difference in iterates $s^{(k)} = x^{(k+1)} - x^{(k)}$ satisfies for each component $i \in \{1, \dots n\}$ the condition
\begin{align}
    \frac{\left|s^{(k)}_i\right|}{\max\{\left| x^{(k)}_{\textcolor{red}{i}} \right|, 1\}} < 4 \epsilon_f
\end{align}
where $\epsilon_f$ is machine precision of double floating point numbers.
\item[Iteration limit:] The number of iterations is at least $k \geq k_{\max}$
where $k_{\max}$ is a parameter of the algorithm.
             \end{description}
        \item Try to use as few function and gradient evaluations of the objective function in your code as possible.    
  \end{itemize}

    \item[ii.] Apply your LBFGS implementation to the Rosenbrock function from part
(b) with $c_1 = 10^{-4}, \eta = 0.9, k_{\max} = 10^4$ and $\epsilon_g =
10^{-8}$. 
 Use the same
initial $x^{(0)}$ as in part (b).

Report (1) the number of iterations does your algorithm needs, (2)
its runtime in seconds, (3) how often the Rosenbrock function evaluated, (4)
how often its gradient is evaluated and (5) the final objective value for each
of the following settings of $m \in \{3, 4, 5, 6, 7, 8, 9, 10, 20\}$.

How do the number of function and gradient evaluations compare against BFGS from part (b)?

\item[iii.] We will now use our LBFGS implementation to find a classifier for the age
prediction task from Problem 4 of Homework 2.  Instead of $\ell_1$ logistic
regression, we now want to optimize the $\ell_2$-penalized logistic regression
loss
\begin{align}
    f(\beta) = -y^\top \tilde X \beta + \sum_{i=1}^{n}\log(1 +
    \exp\{(\tilde X\beta)_i\}) + \lambda \|\beta\|_2^2
\end{align}
where $\tilde X \in \reals^{v \times n}$ are the movie ratings $X \in \reals^{v \times (n-1)}$ with an
additional column appended with entries $1$, $y \in \{0, 1\}^{v}$ are the age classes and
$\beta$ are the logistic regression parameters we want to optimize.
Load $X$ and $y$ from the training set in \texttt{Q4c\_movies.zip} available on
the course webpage and use your LBFGS implementation with parameters $c_1 =
10^{-4}, \eta = 0.9, k_{\max} = 10^4, m=7, \epsilon_g = 10^{-8}$ to find the
minimizer of $f(\beta)$ with $\lambda = 100$. Use $\beta^{(0)} = \mathbf 0$ as initial iterate.

Report the final objective value, the number of iterations, the number of
objective function evaluations, the number of gradient evaluations as well as
the total runtime of the optimization procedure in seconds.

\item[iv.] Use the solution obtained in the previous part to make
predictions on the test set, available in \texttt{moviesTest.mat} in \texttt{Q4c\_movies.zip}. What is the classification error?

\item[v.] Apply your BFGS implementation from part (b) to the same logistic regression problem with identical settings. 
Report the final objective value, the number of iterations, the number of
objective function evaluations, the number of gradient evaluations as well as
the total runtime of the optimization procedure in seconds.
How does LBFGS and BFGS compare on this problem?

\end{enumerate}


\end{enumerate}

\textbf{Attach all relevant code to the end of this problem.}

\section{Row Subset Selection for Linear Regression (25 points) [Han]}

In this question we will consider the optimal subset selection problem in the setting of linear regression. A typical linear regression model is:
$$y = X\beta + \varepsilon$$
where $X\in\RR^{n\times p}$ is the design matrix, $y\in\RR^n$ is the response, $\varepsilon\sim\mathcal{N}(0, \sigma^2 I_n)$ is a vector of i.i.d.\ Gaussian noise with variance $\sigma^2$, and $\beta$ is a fixed but unknown $p$-dimensional coefficient vector. Given the design matrix $X$ and the response $y$, the goal is to estimate the regression model $\beta$. In this problem we focus on the low-dimensional setting where $n > p$ and $\text{rank}(X) = p$. An effective estimator of $\beta$ in the low-dimensional setting is the \emph{ordinary least square} (OLS) estimator, given by:
$$\hat{\beta} = \argmin_{\beta\in\RR^p} ||y - X\beta||_2^2 = (X^TX)^{-1}X^Ty$$
A simple calculation shows that the mean squared error associated with the OLS estimator is:
$$\EE_\varepsilon||\beta - \hat{\beta}||_2^2 = \sigma^2\tr\left((X^TX)^{-1}\right).$$

In many cases where $n\gg p$, it is favorable to select a representative subset of rows in $X$ so that the subsampled regression still accurately estimates $\beta$. However, for a given $k$, with $p\leq k\leq n$, determining the optimal $k$-subset of rows that minimizes the mean squared error leads to a combinatorial optimization problem, which is computationally intractable to solve exactly. To this end, we propose the following convex relaxation of the $k$-subset selection problem to be solved:
\begin{equation}
\begin{aligned}
& \text{minimize}_{\pi\in\RR^n} && \tr\left((X^T\diag(\pi)X)^{-1}\right) \\
& \text{subject to} && \pi \geq 0, 1^T\pi = 1
\end{aligned}
\label{equ:convex}
\end{equation}
\textcolor{red}{where we define the objective function to be $\infty$ whenever $X^T\diag(\pi)X$ is not invertible.}

\begin{enumerate}[(a)]
    \item[(a, 5pts)]   Show that (\ref{equ:convex}) is a convex optimization problem. (\emph{Hint}: There are two ways to go about this: (i) show directly that the above objective is convex, or (ii) reformulate the above problem into one that is of ``canonical'' convex form, i.e., one of the LP, QP, SOCP, SDP, Conic problem classes.)
    \item[(b, 5pts)]   Derive the gradient of the objective function (as a function of $\pi$).
    \item[(c, 5pts)]   List an algorithm that you know of from the course that applies to \eqref{equ:convex}, and describe in as much detail as possible the updates. If an update involves a projection operator, then provide a precise description/pseudocode as to how exactly this projection will be computed. 
    \item[(d, 10pts)]   Implement your proposed algorithm to solve the optimization problem \eqref{equ:convex} on the dataset \texttt{madelon.txt}, available on the course website. Report your optimal objective function value, and the number of nonzero components in the optimal $\pi^*$. \textcolor{red}{To answer this question, you can treat all the values less than $10^{-6}$ to be 0.} (\emph{Hint}: In your implementation, make sure that after each update iteration, the new $\pi$ is in the feasible region and has at least $p$ nonzero components).
    \item [(e, 10 pts, Bonus)]  Describe any other algorithm(s), inspired by what you learned in the course or otherwise, that you believe to be efficient for solving \eqref{equ:convex}.  For full bonus points, implement and run on the data from part (d).

\end{enumerate}
\textbf{Attach all relevant code to the end of this problem.}

\section{Quantile regression and ADMM (25pts) [Justin and Alnur]}

\newcommand{\Prob}{\mathop{\mathrm{Prob}}}
\textit{Quantile regression} is a way to estimate the conditional $\alpha$-quantile, for some $\alpha \in [0,1]$ fixed by the implementer, of some response variable $Y$ given some predictors $X_1,\ldots, X_p$.  (Just to remind you: the $\alpha$-quantile of $Y|X_1,\ldots,X_p$ is given by $\inf_t \left\{ t : \Prob(Y \leq t | X_1,\ldots,X_p) \geq \alpha \right\}$.)  It turns out that we can estimate the $\alpha$-quantile by solving the following (convex) optimization problem:
\begin{equation}
\begin{array}{ll}
\minimizewrt{\beta \in \reals^{p}} & \sum_{i=1}^n  \ell_{i}^{(\alpha)} \left( y - X \beta \right),
\end{array}
\label{eq:aa:qr3}
\end{equation}
where $y \in \reals^n$ is a response vector, $X \in \reals^{n \times p}$ is a data matrix, and the map $\ell^{(\alpha)} : \reals^{n} \to \reals^{n}$ is the \textit{quantile loss}, defined as
\[
\ell^{(\alpha)}_{i}(a) = \max \left\{ \alpha a_{i}, (\alpha-1) a_{i} \right\}, \quad i=1,\ldots,n,
\]
for some $a \in \reals^n$.  (It is not too hard to show why this is the case, but you can just assume it is true for this question.)  In Figure \ref{fig:aa:qr}, we plot some estimated $\alpha$-quantiles, for a few values of $\alpha$, obtained by solving problem \eqref{eq:aa:qr3} on some synthetic data. \\
\begin{figure}[h!]
\centering
\includegraphics[scale=0.5]{quants.pdf}
\caption{Various quantiles estimated from synthetic data.}
\label{fig:aa:qr}
\end{figure}

\newcommand*{\vertbar}{\rule[-1ex]{0.5pt}{2.5ex}}
We often want to estimate multiple quantiles simultaneously, but would like them to not intersect: e.g., in Figure \ref{fig:aa:qr}, we would like to \textit{avoid} the situation where the 0.1-quantile line (maybe eventually) lies above the 0.9-quantile line, since this doesn't make any sense.  To do so, we can instead solve the following (again convex) optimization problem, referred to as the \textit{multiple quantile regression problem with non-crossing constraints}:
\begin{equation}
\begin{array}{ll}
\minimizewrt{B \in \reals^{p \times r}} & \sum_{i=1}^n \sum_{j=1}^r \ell_{ij}^{(\mathcal{A})} \left( y \ones^T - X B \right) \\
\subjectto & X B D^T \geq 0.
\end{array}
\label{eq:aa:qr}
\end{equation}
Here, $r$ is the number of $\alpha$'s (i.e., quantile levels); $\mathcal{A} = \{ \alpha_1, \ldots, \alpha_r \}$ is a set of user-defined quantile levels; $y,X$ are as before; $\ones$ is the $r$-dimensional all-ones vector; and $B \in \reals^{p \times r}$ is now a parameter \textit{matrix}, structured as follows:
\[
B = 
\left[
\begin{array}{ccc}
\vertbar & & \vertbar \\
\beta_1 & \cdots & \beta_r \\
\vertbar & & \vertbar
\end{array}
\right],
\]
for $\beta_i \in \reals^{p}, \; i=1,\ldots,r$.  Here, we have also extended the definition of the quantile loss: it is now a map $\ell^{(\mathcal{A})} : \reals^{n \times r} \to \reals^{n \times r}$, defined as
\[
\ell^{(\mathcal{A})}_{ij}(A) = \max \left\{ \alpha_j A_{ij}, (\alpha_j-1) A_{ij} \right\}, \quad i=1,\ldots,n, \; j=1,\ldots,r.
\]
Lastly, $D \in \reals^{(r-1) \times r}$ is the first-order finite differencing operator, given by
\[
D = 
\left[
\begin{array}{cccccc}
-1 & 1 & 0 & \cdots & 0 & 0 \\
0 & -1 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & -1 & 1
\end{array}
\right].
\]
Hence, you can convince yourself that the constraints in
\eqref{eq:aa:qr} prevent the estimated quantiles from crossing.

Phew!  Now, finally, for the problem.

\begin{enumerate}
\item[(a, 5pts)] Derive the prox operator for the quantile loss $\ell_{ij}^{(\mathcal{A})}$ (defined above). 
In other words, write down an expression for $\prox_{\lambda \ell_{ij}^{(\mathcal{A})}} (A)$, where $\lambda > 0$ is a user-defined constant, and $A \in \reals^{m \times n}$ is the point at which the prox operator is evaluated.

\item[(b, 5pts)] We can put problem \eqref{eq:aa:qr} into ``ADMM form'' as follows:
\begin{equation}
\begin{array}{ll}
\minimizewrt{\substack{B \in \reals^{p \times r}, \\ Z_1, Z_2 \in \reals^{n \times r}, \\ Z_3 \in \reals^{n \times (r-1)}}} & \sum_{i=1}^n \sum_{j=1}^r \ell_{ij}^{(\mathcal{A})} \left( Z_1 \right) + I_+(Z_3) \\
\subjectto & Z_1 = y \ones^T - X B, \quad Z_2 = X B, \quad Z_3 = Z_2 D^T
\end{array}
\label{eq:aa:qr2}
\end{equation}
where we introduced the variables $Z_1, Z_2, Z_3$, and $I_+(\cdot)$ is the projection map onto the nonnegative orthant. \\

The augmented Lagrangian associated with problem \eqref{eq:aa:qr2} is then (in ``scaled'' form)
\begin{multline*}
L_\rho(B, Z_1, Z_2, Z_3, U_1, U_2, U_3) = \sum_{i=1}^n \sum_{j=1}^r
\ell_{ij}^{(\mathcal{A})} \left( Z_1 \right) + I_+(Z_3) +
\frac{\rho}{2} \Big( \| y \ones^T - X B - Z_1 + U_1 \|_F^2 +{} \\ 
  \| X B - Z_2 + U_2 \|_F^2 + \| Z_2 D^T - Z_3 +U_3 \|_F^2 - \| U_1 \|_F^2
  - \| U_2 \|_F^2 - \| U_3 \|_F^2 \Big),  
\end{multline*}
where $U_1, U_2, U_3$ are dual variables. \\

Write down the ADMM updates for problem \eqref{eq:aa:qr2}. You can
assume that the quantity $X^T X$ is invertible. 

\item[(c, 7 pts)] Implement your ADMM algorithm from part (b) for the quantiles
  $\mathcal{A} = \{0.1, 0.5, 0.9\}$, using synthetic data from 
  \href{www.stat.cmu.edu/~ryantibs/convexopt/hw4q4.zip}{hw4q4.zip} which contain
  a predictor matrix \texttt{X}, response vector \texttt{y}, and in csv
  files. Note that the predictor matrix in this example data is $1000
  \times 2$, with the 2nd column being the column of all 1s, hence 
  providing an intercept term in the quantile regression.

  Set initial values for all primal and dual variables as zero matrices
  of appropriate size, use an augmented Lagrangian parameter of
  $\rho=1$, and run ADMM for 50 steps. Report your optimal criterion
  value after 50 steps, and plot the optimal criterion value across
  the iterations. Recreate the plot similar to Figure \ref{fig:aa:qr},
  plotting data as points and overlaying as lines the
  \textit{conditional quantile} 
  predictions for the quantiles in $\mathcal{A}$. (Hint: recall that $B$ is a
  matrix whose $r$ columns are coefficient estimates for each quantile
  level $\alpha_1,\cdots,\alpha_r$, so that $X_{new}^TB$ at a new data
  point $X_\text{new} = (x_\text{new}, 1)$ is the fitted estimates for
  each quantile in $\mathcal{A}$.)

\item[(d, 8 pts)] An alternative (more efficient) way to parameterize
  problem \eqref{eq:aa:qr} in ``ADMM form'' is as follows: 
  \begin{equation}
    \begin{array}{ll}
      \minimizewrt{\substack{B \in \reals^{p \times r}, \\ Z_1, Z_2 \in \reals^{n \times r}}} & \sum_{i=1}^n \sum_{j=1}^r \ell_{ij}^{(\mathcal{A})} \left( Z_1 \right) + I_+(Z_2 D^T) \\
      \subjectto & Z_1 = y \ones^T - X B, \quad Z_2 = X B.
    \end{array}
    \label{eq:aa:qr4}
  \end{equation}

  The augmented Lagrangian associated with problem \eqref{eq:aa:qr4} is then
  \begin{multline*}
    L_\rho(B, Z_1, Z_2, U_1, U_2) = \sum_{i=1}^n \sum_{j=1}^r
    \ell_{ij}^{(\mathcal{A})} \left( Z_1 \right) + I_+(Z_2D^T) \\
    \quad + \frac{\rho}{2} \Big( \| y \ones^T - X B - Z_1 + U_1 
      \|_F^2 + \| X B - Z_2 + U_2 \|_F^2 - \| U_1 \|_F^2 - \| U_2
      \|_F^2 \Big).  
  \end{multline*}

  Write down and implement the ADMM updates for problem
  \eqref{eq:aa:qr4}.  Explain why the $Z_2$ update reduces to solving
  $n$ isotonic regressions.  Give the name of an algorithm (by looking
  around in the literature) for isotonic regression that runs in 
  linear-time (or very nearly linear-time).  Most programming
  languages (Python, Matlab, R) should provide access to an
  implementation of such an algorithm.

  Then, run the ADMM steps for \eqref{eq:aa:qr4} on the data from part
  (c), and reproduce the same things requested in the part.  Compare
  the optimal criterion curves between the implementation in part (c)
  and the one from the current implementation---is there a noticeable
  difference? 
\end{enumerate}

\textbf{Attach all relevant code to the end of this problem.}

\end{document}
