\documentclass{article} \usepackage[left=1.25in,top=1.25in,right=1.25in,bottom=1.25in,head=1.25in]{geometry} \usepackage{amsmath,amssymb,amsthm} \usepackage{verbatim,float,url,enumerate} \usepackage{graphicx,subfigure,psfrag} \usepackage{natbib,xcolor,dsfont} \usepackage{microtype} \newtheorem{algorithm}{Algorithm} \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}}} \def\R{\mathbb{R}} \def\E{\mathbb{E}} \def\P{\mathbb{P}} \def\S{\mathbb{S}} \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\dom{\mathrm{dom}} \def\prox{\mathrm{prox}} \begin{document} \title{Homework 4} \author{\Large Convex Optimization 10-725} \date{{\bf Due Friday October 25 at 11:59pm} \\ \bigskip Submit your work as a single PDF on Gradescope. Make sure to prepare your solution to each problem on a separate page. (Gradescope will ask you select the pages which contain the solution to each problem.) \\ \bigskip Total: 84 points (+ 5 bonus points)} \maketitle \section{Newton's method convergence analysis (18 points)} In this problem, we will prove that Newton's method obtains a (local) quadratic rate of convergence under suitable assumptions. \subsection{Univariate setting (5 points)} As a warm up, let's start by looking at the univariate setting, where we are minimizing function $f : \R \to \R$ that is convex and three times continuously differentiable. Assume that $f''(x) \geq C_1 > 0$ and $|f'''(x)| \leq C_2$ for all $x$. Let $x^\star$ be the global minimizer of $f$, and show the Newton's method iterates, under pure step sizes, satisfy $$ |x^{(k)} - x^\star| \leq \frac{C_2}{2C_1} |x^{(k-1)} - x^\star|^2, $$ for all $k=1,2,3,\ldots$. What does this imply about global convergence? Hint: use a Taylor expansion of $f'$. \subsection{General setting (13 points)} Now let's move to the general case, where $f : \R^n \to \R$ is convex and twice continuously differentiable. Assume that $\nabla f$ is Lipschitz with parameter $L$, $f$ is strongly convex with parameter $m>0$, and $\nabla^2 f$ is Lipschitz with parameter $M$: $$ \| \nabla^2 f(x) - \nabla^2 f(y) \|_{\mathrm{op}} \leq M \|x-y\|_2, $$ for all $x,y$, where $\|\cdot\|_{\mathrm{op}}$ is the operator norm. Again let $x^\star$ denote the global minimizer of $f$, and consider Newton's method with backtracking. We will make our lives easier by assuming that we are already in the in the ``pure phase'', where backtracking results in steps of size $t=1$. (Note: you will not need to use the Lipschitz condition on $\nabla f$ in what follows, but this is required for showing that we reach the ``pure phase'' to begin with.) \bigskip \noindent (a, 3 pts) Using that fact that $f$ is strongly convex, prove that for all $x$, $$ f(x) - f(x^\star) \leq \frac{1}{2m} \|\nabla f(x)\|_2^2. $$ Hint: use the quadratic lower bound on $f(y)$ around $f(x)$ that is implied by strong convexity, then minimize both sides over $y$. \bigskip \noindent (b, 6 pts) Show that at an arbitrary iteration of Newton's method, the current and next iterates $x,x^+$ satisfy $$ \frac{M}{2m^2} \|\nabla f(x^+)\|_2 \leq \bigg(\frac{M}{2m^2} \|\nabla f(x)\|_2\bigg)^2. $$ Hint: use a Taylor expansion of $\nabla f$. \bigskip \noindent (c, 4 pts) Suppose that it took us $k_0$ iterations to reach the pure phase. Note that this means by definition that $$ \|\nabla f(x^{(k_0)})\|_2 < \eta \leq \frac{m^2}{M}. $$ By iterating the inequality in part (b) and invoking the result in part (a), show that at any iteration $k>k_0$, we have $$ f(x^{(k)}) - f(x^\star) \leq \frac{2m^3}{M^2} \Big( \frac{1}{2} \Big)^{2^{k-k_0+1}}. $$ \section{GLMs + Newton's method = IRLS (15 points)} In this problem we will study the convexity properties underlying exponential families and generalized linear models (GLMs), and Newton's method applied to perform maximum likelihood. Consider an exponential family density (or mass) function over $y \in D \subseteq \R^n$, of the form \begin{equation} \label{eq:expfam} f(y; \theta) = \exp \big(y^T\theta - b(\theta)\big) f_0(y). \end{equation} Note $\theta \in \R^n$ is called the ``natural parameter''. \bigskip \noindent (a, 2 pts) Prove that the domain of $\theta$, $C = \{ \theta : b(\theta) < \infty \}$, is a convex set. \bigskip \noindent (b, 5 pt) Prove that $b : C \rightarrow \R$ is a convex function. Hint: prove that it is twice differentiable and its Hessian is positive semidefinite everywhere. \bigskip \noindent (c, 1 pt) Suppose that we model $\theta = X \beta$, where $X \in \R^{n\times p}$ is a feature matrix. Note that under this parametrization, the exponential family model in \eqref{eq:expfam} is called a generalized linear model (GLM). Prove that the domain of $\beta$, $B = \{\beta : X\beta \in C\}$, is a convex set. \bigskip \noindent (d, 3 pts) Show that maximizing the log likelihood in \eqref{eq:expfam}, when $\theta=X\beta$, is equivalent to \begin{equation} \label{eq:glm_loglik} \min_\beta \; -y^T X \beta + b(X \beta). \end{equation} Argue that this is a convex optimization problem. What choice of $b$ recovers linear regression? What choice recovers logistic regression? \bigskip \noindent (e, 4 pts) Let $$ \mu_\beta = \nabla b(X\beta) \;\;\;\ \text{and} \;\;\; V_\beta = \nabla^2 b (X\beta). $$ Show that pure Newton's method (with step size $t=1$) applied to \eqref{eq:glm_loglik} is equivalent to the updates $$ \beta^+ = (X^T V_\beta X)^{-1} X^T V_\beta z_\beta, $$ where $z_\beta = X \beta + V_\beta^{-1} (y - \mu_\beta)$. Explain why this is called {\it iteratively reweighted least squares} (IRLS): what optimization problem (for fixed $\beta$) would have $\beta^+$ as its solution? \section{Binary sequence denoising: dual Newton method (23 points)} Recall from Q4 of Homework 3 the binary sequence denoising problem, where we are given $z_i \in \{0,1\}$, $i=1,2,\ldots,n$, and we can think of an underlying logistic model giving us probabilities $$ p_i = \frac{e^{\theta_i}}{1+e^{\theta_i}}, \;\;\; i=1,\ldots,n. $$ We compute estimates \smash{$\hat\theta_i$}, $i=1,\ldots,n$ by solving the logistic fused lasso problem: \begin{equation} \label{eq:log_fl2} \min_\theta \; \sum_{i=1}^n \big(-z_i \theta_i + \log(1+e^{\theta_i}) \big) + \lambda \sum_{i=1}^{n-1} |\theta_i-\theta_{i+1}|. \end{equation} By setting $y_i=2z_i-1$, $i=1,\ldots,n$ and $D \in \R^{(n-1) \times n}$ to be a bidiagonal matrix with lower and upper diagonals $-1$ and 1, the problem \eqref{eq:log_fl2} is equivalent to \begin{equation} \label{eq:log_fl} \min_\theta \; \sum_{i=1}^n \log(1+e^{-y_i\theta_i}) + \lambda\|D\theta\|_1. \end{equation} Recall also that the dual of \eqref{eq:log_fl} is \begin{equation} \label{eq:log_fl_dual} \begin{alignedat}{2} &\min_u \quad &&\sum_{i=1}^n \Big( y_i(D^Tu)_i \log(y_i(D^Tu)_i) + (1-y_i(D^Tu)_i) \log (1-y_i(D^Tu)_i) \Big) \\ &\st \quad && 0 \leq y_i(D^Tu)_i \leq 1, \; i=1,\ldots,n, \; \|u\|_\infty \leq \lambda. \end{alignedat} \end{equation} and the relationship between the solutions in \eqref{eq:log_fl}, \eqref{eq:log_fl_dual} is \begin{equation} \label{eq:log_fl_pd} \theta_i = y_i \log\bigg( \frac{1-y_i(D^Tu)_i}{y_i(D^Tu)_i} \bigg), \; i=1,\ldots,n. \end{equation} Now come the questions. \bigskip \noindent (a, 15 pts) Write down the gradient and Hessian of the ``softened'' dual problem \begin{multline} \label{eq:log_fl_dual2} \min_u \; \sum_{i=1}^n \Big( y_i(D^Tu)_i \log(y_i(D^Tu)_i) + (1-y_i(D^Tu)_i) \log (1-y_i(D^Tu)_i) \Big) \\ - \frac{1}{\tau} \cdot \sum_{i=1}^n \bigg( \log (y_i(D^Tu)_i) + \log (1-y_i(D^Tu)_i)\bigg) - \frac{1}{\tau} \cdot \sum_{i=1}^n \bigg( \log(\lambda - u_i) + \log(u_i + \lambda))\bigg). \end{multline} Implement Newton's method to solve \eqref{eq:log_fl_dual2}, using backtracking to determine the step size at each iteration, and stopping when the difference in criterion value across successive iterations is less than a user-specified tolerance level $\epsilon$. Remember to attach all of your code in an appendix to your homework. Hint: you will need to solve an LP in order to initialize Newton's method at a feasible dual point; recall Q4(c) of Homework 3. After implementing it, run your algorithm on the data {\tt binseq.txt} from Homework 3. Use $\tau=10^3$. Also use $\lambda=20$ as the tuning parameter, $\beta=0.8$ as the contraction factor for backtracking, and $\epsilon=10^{-6}$ as the stopping tolerance. Report how many iterations it took to converge, and compare this to the 50,000 iterations it took gradient descent to reach rough approximate solution in Q4(d) of Homework 3. It should be much, much less! After using the primal-dual relationship \eqref{eq:log_fl_pd} to get a primal solution from your dual solution, compute the estimated probabilities \smash{$\hat{p}_i$}, $i=1,\ldots,n$, and plot them. On the same figure, plot the estimated probabilities from the solution obtained by running proximal gradient descent directly on the primal, from Q4(a) of Homework 3. Do these two look similar? \bigskip \noindent (b, 8 pts) Now (modifying your primal proximal gradient and dual Newton implementations if needed, so that all of the intermediate iterates in the algorithm are saved), rerun primal proximal gradient descent and dual Newton's method, each with $\epsilon=0$, and for 100 iterations. Compute the primal criterion $f(\theta^{(k)})$ at each iterate of $\theta^{(k)}$ of these algorithms. Take the minimum observed criterion value across these two algorithms, subtract $10^{-6}$, and call this $f^\star$. Then plot the criterion gap as a function of iteration number, for both primal proximal gradient and dual Newton, on the same figure. Make sure the y-axis is on a log scale. What do you notice about how fast each algorithm converges? What do you notice about the point at which the dual Newton method converges to? It should noticeably be suboptimal; explain why. Repeat the same analysis (rerunning the algorithms, computing the criterions, computing $f^\star$, producing a figure) when $\lambda=0.02$. Address the same questions; comment on the differences you see to the case $\lambda=20$. {\bf Bonus (2 pts):} explain why the differences between large and small $\lambda$ cases occur. \section{Binary sequence denoising: dual barrier method (28 points)} (a, 10 pts) Implement the barrier method to solve the logistic fused lasso dual \eqref{eq:log_fl_dual}. With each outer iteration of the barrier method, you will solve the barrier problem \eqref{eq:log_fl_dual2} via Newton's method, as you implemented in Q3(a). Your implementation should take an initial value $\tau_0$ for the barrier parameter, a multiplicative update factor $\mu>1$ for the barrier parameter, and an ``outer'' tolerance $\epsilon_{\mathrm{outer}}$ for stopping when the duality gap is smaller than this value. Remember to attach all of your code in an appendix to your homework. Run your implementation on the binseq data ({\tt binseq.txt} from Homework 3), with $\lambda=20$, the ``outer'' parameters set at $\tau_0=5$, $\mu=10$, and $\epsilon_{\mathrm{outer}}=10^{-8}$, and the ``inner'' parameters for Newton's method set as in Q3(a) ($\beta=0.8$ for the contraction factor for backtracking, and $\epsilon_{\mathrm{inner}}=10^{-6}$ for the inner stopping tolerance). Report how many total iterations it took the barrier method to converge (sum of the Newton steps over all outer iterations). After using the primal-dual relationship \eqref{eq:log_fl_pd} to get a primal solution from your dual solution, compute the estimated probabilities \smash{$\hat{p}_i$}, $i=1,\ldots,n$, and plot them. On the same figure, plot the estimated probabilities from the solution obtained by primal proximal gradient descent (from Q4(a) of Homework 3), and dual Newton's method (from from Q3(a) of this homework). Does the dual barrier method solution look closer to the primal proximal gradient solution than the dual Newton's method solution does? Hint: storing $D$ as a sparse and/or structured (banded) matrix should make a big difference in computational speed, both because multiplying by $D$ or $D^T$ should be much faster, and (more importantly) because once things are properly set up, solving a linear system in the Hessian should be much, much faster. \bigskip \noindent (b, 8 pts) Rerun your primal proximal gradient descent implementation on the binseq data, over 80 values of $\lambda$ in between 0.001 and 200, equally spaced on the log scale. Starting from the largest $\lambda$ value to the smallest, run proximal gradient using both warm starts: initializing from the previously computed solution, and using cold starts: initializing from $\theta=0$. For the rest of the parameter values, use the same defaults as in Q4(a) of Homework 3 ($t=1$ as the initial step size before each backtracking loop, $\beta=0.8$ as the contraction factor in backtracking, and $\epsilon=10^{-6}$ as the stopping tolerance), and use a maximum number of 1000 iterations. For each strategy: warm/cold starts, record the total number of iterations (sum the number of backtracking iterations, i.e., prox evaluations, performed by the algorithm) at each value of $\lambda$. Remember that the prox evaluation is more or less the fundamental unit of computation for proximal gradient descent. Also record the final criterion value (upon convergence or termination at 1000 iterations), for later analysis (not in this question). Plot the number of total iterations taken by the algorithm as a function of $\lambda$, overlaying the curves for both warm and cold starts, and with the x-axis on a log scale. Do you see a difference? And more broadly, what do the curves portray about the difficulty of solving the problem \eqref{eq:log_fl} as a function of $\lambda$? \bigskip \noindent (c, 8 pts) Repeat the analysis as in part (b), but using your dual barrier method implementation. Now you will go from the smallest $\lambda$ value to the largest. For warm starts, as before, you will just use the previously computed solution as your initial point (note that in the other direction, the warm starts would not be feasible). For the cold starts, you can simply initialize using a single fixed feasible point $u_0$ obtained by solving the feasibility LP at the smallest $\lambda$ value, once at the start, which will then remain feasible for all larger $\lambda$. For the rest of the parameter values, use the same defaults as in Q4(a) of this homework, and use a maximum number of 100 iterations. For each strategy: warm/cold starts, record the total number of iterations (sum the number of Newton steps) at each value of $\lambda$. Remember that the Newton step is more or less the fundamental unit of computation for the barrier method. Also record the final criterion value (upon convergence or termination at 100 iterations), again for later analysis (not in this question). Plot the number of total iterations taken by the algorithm as a function of $\lambda$, overlaying the curves for both warm and cold starts, and with the x-axis on a log scale. Comment on the difference between the strategies, and any differences to the results for the proximal gradient case in part (b). \bigskip \noindent (d, 2 pts) Plot the difference in final criterion values between primal proximal gradient and dual barrier method, as a function of $\lambda$, from your computations parts (b) and (c) (use the criterions from the warm starts). Which algorithm is more accurate, and are there any noticeable patterns? \bigskip \noindent {\bf Bonus (3 pts):} What is the difference between computational complexity of one call to the prox and one Newton step? Are they of the same order (as a function of $n$)? Explain, and then perform numerical timings to support your arguments. \end{document}