Last Time

• We have spent much time on nonsmooth problems

• Since they are very common in statistical models

• Meanwhile they are considered "hard problems" in convex optimization

• So far, the main tool to overcome nonsmoothness is the proximal operator

• Seen in proximal gradient descent, Douglas-Rachford splitting, Davis-Yin splitting, etc.

• In this lecture we introduce a few other frameworks

Today's Topics

• Alternating direction method of multipliers

• Coordinate descent

• Case studies

Problem

Another popular framework to deal with nonsmooth or constrained problems is the alternating direction method of multipliers (ADMM). It considers problems of the form:

$$\begin{array}{cc}\underset{x,z}{min}& F(x,z):=f\left(x\right)+g\left(z\right)\\ \text{subject to}& Ax+Bz=c\end{array}$$

• $x\in R^n$, $z\in R^m$, $A\in R^{p\times n}$, $B\in R^{p\times m}$, and $c\in R^p$.

• $f\left(x\right)$ and $g\left(z\right)$ are convex functions, possibly nonsmooth

• $Ax+Bz=c$ represents general linear equality constraints

Examples

• Although this problem form is a bit exotic compared with those we have seen in previous lectures

• Many familiar problems can be converted into this form by cleverly defining the $x$ and $z$ variables

• Lasso:

$$\underset{x}{min}\frac{1}{2}\parallel b-Ax{\parallel }^2+\lambda \parallel x{\parallel }_1$$

• $f\left(x\right)=\left(1/2\right)\parallel b-Ax{\parallel }^2$, $g\left(z\right)=\lambda \parallel z{\parallel }_1$

• Linear constraint $x-z=0$

Examples

• Least absolute deviations:

$$\underset{x}{min}\parallel b-Ax{\parallel }_1$$

• $f\left(x\right)=0$, $g\left(z\right)=\parallel z{\parallel }_1$

• Linear constraint $Ax-z=b$

Algorithm

ADMM can be expressed as the following algorithm. Given initial values $z^{\left(0\right)}$ and $u^{\left(0\right)}$, for $k=0,1\dots$, iterate

$$\begin{array}{cc}{x}^{(k+1)}& =\underset{x}{\arg min}f\left(x\right)+\left(\rho /2\right)\parallel Ax+B{z}^{\left(k\right)}-c+u^{\left(k\right)}{\parallel }^2\\ z^{(k+1)}& =\underset{z}{\arg min}g\left(z\right)+\left(\rho /2\right)\parallel A{x}^{(k+1)}+Bz-c+u^{\left(k\right)}{\parallel }^2\\ u^{(k+1)}& =u^{\left(k\right)}+A{x}^{(k+1)}+B{z}^{(k+1)}-c\end{array}$$

$\rho >0$ is an arbitray "inverse step size" parameter.

Some Remarks

• The $x$- and $z$-updates look similar to proximal operators

• However, here we have more general quadratic terms

• This means that the $x$- and $z$-updates are typically more difficult than proximal operators

• On the other hand, ADMM does not need to explicitly project to the linear constraint set

• These findings largely guide how to define the $f$ and $g$ functions

Relation with DRS

Consider using ADMM to solve the problem

$$\underset{x}{min}f\left(x\right)+g\left(x\right)\iff \underset{x,z}{min}f\left(x\right)+g\left(z\right)\text{subject to}x-z=0$$

Then the iteration becomes

$$\begin{array}{cc}{x}^{(k+1)}& ={prox}_{\left(1/\rho \right)f}(z^{\left(k\right)}-u^{\left(k\right)})\\ z^{(k+1)}& ={prox}_{\left(1/\rho \right)g}(x^{(k+1)}+u^{\left(k\right)})\\ u^{(k+1)}& =u^{\left(k\right)}+x^{(k+1)}-z^{(k+1)}\end{array}$$

This is equivalent to DRS, by slightly reparameterizing the variables.

Implementation Details

• Many implementation details can be found in the main reference [1]

• Stopping criterion

• Selection of the $\rho$ parameter

• Closed-form expressions for specific problems

Convergence Property ^[2]

The sequence $(x^{\left(k\right)},z^{\left(k\right)},u^{\left(k\right)})$ satisfy

$$\parallel A(x^{(k+1)}-x^{\left(k\right)}){\parallel }^2+\parallel B(z^{(k+1)}-z^{\left(k\right)}){\parallel }^2+\parallel u^{(k+1)}-u^{\left(k\right)}{\parallel }^2=o\left(1/k\right).$$

This is one of the most general convergence results on ADMM.

Many other faster rates exist by imposing stronger conditions on $f\left(x\right)$ and $g\left(z\right)$, e.g. smoothness, strong convexity, etc.

Motivation

• One more idea to overcome nonsmoothness is to directly find the optimal point of nonsmooth problems

• Of course, this is unrealistic for general problems

• But in many cases this is possible for univariate functions

Example

Lasso model $$\underset{\beta }{min}\ell \left(\beta \right):=\underset{\beta }{min}\frac{1}{2n}\parallel y-X\beta {\parallel }^2+\lambda \parallel \beta {\parallel }_1$$

• Solving the whole $\beta$ vector is hard

• However, assume we fix all elements of $\beta$ except ${\beta }_i$

• Then we can actually derive the optimal ${\beta }_i$ in closed form

Example

• Recall that in Lecture 7 we have shown

$$f\left(x^{\ast }\right)=\underset{x}{min}f\left(x\right)\iff 0\in \partial f\left(x^{\ast }\right)$$

• View $\tilde{\ell }\left({\beta }_i\right)\equiv \ell \left(\beta \right)$ as a function of ${\beta }_i$, then

$$\partial \tilde{\ell }\left({\beta }_i\right)=n^{-1}\parallel X_i{\parallel }^2{\beta }_i+n^{-1}{X}_i^{\prime }(X_{-i}{\beta }_{-i}-y)+\lambda \partial |{\beta }_i|,$$

• $X_i$ is the $i$-column of $X$

• $X_{-i}$ is the matrix formed by removing $X_i$ from $X$

• ${\beta }_{-i}$ is the vector formed by removing ${\beta }_i$ from $\beta$

Example

• $0\in \partial \tilde{\ell }\left({\beta }_i\right)$ gives the solution

$${\beta }_i=S_{n\lambda /\parallel X_i{\parallel }^2}\left(\frac{X_i^{\prime }(y-X_{-i}{\beta }_{-i})}{\parallel X_i{\parallel }^2}\right)$$

Coordinate Descent

These findings motivate the coordinate descent (CD) algorithm. Consider the problem

$$\underset{x}{min}f\left(x\right):=\underset{x}{min}f(x_1,\dots ,x_n)$$ One common type of CD algorithm proceeds as follows. Given initial value $x^{\left(0\right)}$, for $k=0,1,\dots$, iterate:

$$x_i^{(k+1)}=\underset{x_i}{\arg min}f(x_1^{(k+1)},\dots ,x_{i-1}^{(k+1)},x_i,x_{i+1}^{\left(k\right)},\dots ,x_n^{\left(k\right)})$$ for $i=1,\dots ,n$ within each $k$-outer-iteration.

This is typically called the cyclic CD algorithm.

Variants

Another popular variant is to randomly select a coordinate to update. Given initial value $x^{\left(0\right)}$, for $k=0,1,\dots$, do:

1. Randomly select $i_k$ from $\{1,2,\dots ,n\}$ with equal probabilities

2. Fixing other coordinates and update

$$\begin{array}{cc}{x}_{i_k}^{(k+1)}& =\underset{x_{i_k}}{\arg min}f(x_1^{\left(k\right)},\dots ,x_{i_k-1}^{\left(k\right)},x_{i_k},x_{i_k+1}^{\left(k\right)},\dots ,x_n^{\left(k\right)})\\ x_j^{(k+1)}& =x_j^{\left(k\right)},j\ne i_k\end{array}$$

In the literature, this method is sometimes called randomized CD.

Variants

Many other variants exist:

• Permutation CD: for each $k$, use a random permutation of $\{1,2,\dots ,n\}$ to define the order of coordinate updates

• Block CD: each time update a block of coordinates simultaneously, instead of a single coordinate

• Coordinate proximal gradient descent: use proximal operators to update coordinates

• ...

Convergence Property ^[3]

Unfortunately, CD does not always converge to an optimal point for nonsmooth $f$. Instead, consider the problem

$$f\left(x\right)=g\left(x\right)+\sum _{i=1}^n{h}_i\left(x_i\right),$$

where $g\left(x\right)$ is convex and differentiable, and each $h_i\left(x_i\right)$ is convex, possibly nonsmooth.

[3] shows that with very mild regularity conditions, the CD iterates $x^{\left(k\right)}$ converge to an optimal point of $f\left(x\right)$.

The result holds for block CD and arbitrary order of cycle through coordinates.

Convergence Property ^[4]

There exist many convergence results for different variants of CD based on different assumptions on $f\left(x\right)$. Below is one representative result with weak assumptions.

Assume that $g\left(x\right)$ is coordinate-wise $L_i$-smooth, and define $$R_0=\underset{y}{max}\underset{x^{\ast }\in X^{\ast }}{max}\{{\left(\sum _{i=1}^n{L}_i\right(y_i-x_i^{\ast }{)}^2)}^{1/2}:f(y)\le f(x^{\left(0\right)}\left)\right\},$$ where $X^{\ast }$ denotes the optimal set.

(See next slide)

Convergence Property ^[4]

Consider the randomized CD algorithm. Fix a confidence parameter $0<\rho <1$ and an accuracy level $\epsilon <min\{R_0^2,f(x^{\left(0\right)})-f(x^{\ast }\left)\right\}$. Then with

$$k\ge \frac{2n{R}_0^2}{\epsilon }\log \frac{f\left(x^{\left(0\right)}\right)-f\left(x^{\ast }\right)}{\epsilon \rho },$$ we have $f\left(x^{\left(k\right)}\right)-f\left(x^{\ast }\right)\le \epsilon$ with probability at least $1-\rho$.

Convergence with Strong Convexity ^[4]

For simplicity suppose $L_i\equiv L$, and assume that $g\left(x\right)$ is also strongly convex with parameter $m>0$. Fix a confidence parameter $0<\rho <1$ and an accuracy level $\epsilon >0$. If

$$k\ge \frac{4Ln}{m}\log \frac{f\left(x^{\left(0\right)}\right)-f\left(x^{\ast }\right)}{\epsilon \rho },$$

then $f\left(x^{\left(k\right)}\right)-f\left(x^{\ast }\right)\le \epsilon$ with probability at least $1-\rho$.

Lasso

We use the lasso model $$\underset{\beta }{min}\ell \left(\beta \right):=\underset{\beta }{min}\frac{1}{2n}\parallel y-X\beta {\parallel }^2+\lambda \parallel \beta {\parallel }_1$$ to study the convergence property and computational efficiency of several algorithms we have covered in previous lectures.

• Derive the update rules for several algorithms

• Study the change of loss function versus iteration number

• Study the change of loss function versus running time

Simulate Data Set

The design matrix is a simulated sparse matrix, but we tentatively use its dense version.

library(Matrix)
set.seed(123)
n = 500
p = 1000
s = 30
xsp = rsparsematrix(n, p, density = 0.1)
x = as.matrix(xsp)
beta = c(rnorm(s), numeric(p - s))
y = as.numeric(xsp %*% beta) + 0.1 * rnorm(n)
lam = 0.001
objfn = function(beta, x, y, lam)
{
    0.5 * mean((y - as.numeric(x %*% beta))^2) + lam * sum(abs(beta))
}

Compute Optimal Value

Use the celebrated glmnet package to compute the optimal value.

library(glmnet)
res = glmnet(x, y, lambda = lam, thresh = 1e-15,
             standardize = FALSE, intercept = FALSE)
beta_opt = as.numeric(res$beta)
head(beta_opt)

## [1] -0.3862279  1.4304118 -0.2206914  1.1256892  0.1704705  0.8657187

tail(beta_opt)

## [1]  0.000000000 -0.003938278  0.000000000  0.000000000  0.000000000
## [6]  0.008219107

obj_opt = objfn(beta_opt, x, y, lam)
obj_opt

## [1] 0.02782515

Optimization

A general interface for running optimization algorithms.

run_opt = function(x, y, lam, init_state, update_rule, maxit = 1001, ...)
{
    ...
}

• init_state() and update_rule() are two function arguments to characterize the optimization algorithm

• init_state(x, y, lam, ...) is used to create initial state of the algorithm

• update_rule(beta, iter, x, y, lam, state) returns the updated $\beta$ and state, list(beta, state)

Optimization

A general interface for running optimization algorithms.

run_opt = function(x, y, lam, init_state, update_rule, maxit = 1001, ...)
{
    p = ncol(x)
    beta_hat = numeric(p)
    losses = c()
    # Initialization
    state = init_state(x, y, lam, ...)
    # Main loop
    for(i in 1:maxit)
    {
        res = update_rule(beta_hat, i - 1, x, y, lam, state)
        beta_hat = res$beta
        state = res$state
        loss = objfn(beta_hat, x, y, lam)
        losses = c(losses, loss)
        if((i - 1) %% 100 == 0)
            cat(sprintf("Iter %d, loss = %.6f\n", i - 1, loss))
    }
    list(beta_hat = beta_hat, losses = losses)
}

Subgradient Method

• Choose $\tilde{\nabla }\ell \left(\beta \right)=\left(1/n\right)X^{\prime }(X\beta -y)+\lambda \cdot sign\left(\beta \right)$

• Update rule: $${\beta }^{(k+1)}={\beta }^{\left(k\right)}-{\alpha }_k\tilde{\nabla }\ell \left({\beta }^{\left(k\right)}\right)$$

• Step size ${\alpha }_k={\alpha }_0/\sqrt{k+1}$

Implementation

init_subgrad = function(x, y, lam, alpha0 = 1)
{
    list(alpha0 = alpha0)
}
update_subgrad = function(beta, k, x, y, lam, state)
{
    n = nrow(x)
    alpha0 = state$alpha0
    alphak = alpha0 / sqrt(k + 1)
    Xtr = crossprod(x, as.numeric(x %*% beta) - y)
    subgrad = as.numeric(Xtr) / n + lam * sign(beta)
    new_beta = beta - alphak * subgrad
    list(beta = new_beta, state = state)
}

Result

res_subgrad1 = run_opt(
    x, y, lam, init_subgrad, update_subgrad, alpha0 = 1)

## Iter 0, loss = 0.770908
## Iter 100, loss = 0.084973
## Iter 200, loss = 0.074552
## Iter 300, loss = 0.070417
## Iter 400, loss = 0.067906
## Iter 500, loss = 0.066074
## Iter 600, loss = 0.064626
## Iter 700, loss = 0.063428
## Iter 800, loss = 0.062394
## Iter 900, loss = 0.061477
## Iter 1000, loss = 0.060661

head(res_subgrad1$beta)

## [1] -2.174641e-01  9.316010e-01 -1.026982e-01  6.878443e-01  1.129936e-05
## [6]  2.784390e-01

tail(res_subgrad1$beta)

## [1] -3.352589e-02 -5.879592e-02 -1.432275e-06 -3.051985e-02  5.442229e-02
## [6]  1.141592e-01

Result

res_subgrad10 = run_opt(
    x, y, lam, init_subgrad, update_subgrad, alpha0 = 10)

## Iter 0, loss = 8.881744
## Iter 100, loss = 0.049232
## Iter 200, loss = 0.040807
## Iter 300, loss = 0.036166
## Iter 400, loss = 0.032993
## Iter 500, loss = 0.030828
## Iter 600, loss = 0.029429
## Iter 700, loss = 0.028636
## Iter 800, loss = 0.028227
## Iter 900, loss = 0.028043
## Iter 1000, loss = 0.027961

head(res_subgrad10$beta)

## [1] -0.3913635  1.4284850 -0.2179803  1.1101840  0.1683580  0.8471952

tail(res_subgrad10$beta)

## [1] -2.979032e-04 -3.842055e-03  2.795830e-04 -9.758347e-05 -1.092878e-04
## [6]  1.166201e-02

Proximal Gradient Descent

• Update rule: $${\beta }^{(k+1)}=S_{\alpha \lambda }({\beta }^{\left(k\right)}-(\alpha /n\left)X^{\prime }\right(X{\beta }^{\left(k\right)}-y\left)\right)$$

• Step size $\alpha =1/{\lambda }_{max}\left(n^{-1}{X}^{\prime }X\right)=n/{\lambda }_{max}\left(X^{\prime }X\right)$

Implementation

library(RSpectra)
# Soft-thresholding function
soft_thresh = function(x, a)
{
    sign(x) * pmax(0, abs(x) - a)
}
init_ista = function(x, y, lam)
{
    # Compute the largest eigenvalue of X'X
    alpha = nrow(x) / svds(x, k = 1, nu = 0, nv = 0)$d^2
    list(alpha = alpha)
}
update_ista = function(beta, k, x, y, lam, state)
{
    n = nrow(x)
    alpha = state$alpha
    Xtr = crossprod(x, as.numeric(x %*% beta) - y)
    grad = as.numeric(Xtr) / n
    new_beta = soft_thresh(beta - alpha * grad, alpha * lam)
    list(beta = new_beta, state = state)
}

Result

res_ista = run_opt(x, y, lam, init_ista, update_ista)

## Iter 0, loss = 0.494847
## Iter 100, loss = 0.044493
## Iter 200, loss = 0.032321
## Iter 300, loss = 0.028040
## Iter 400, loss = 0.027830
## Iter 500, loss = 0.027825
## Iter 600, loss = 0.027825
## Iter 700, loss = 0.027825
## Iter 800, loss = 0.027825
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_ista$beta)

## [1] -0.3862280  1.4304114 -0.2206905  1.1256887  0.1704708  0.8657158

tail(res_ista$beta)

## [1]  0.000000000 -0.003938134  0.000000000  0.000000000  0.000000000
## [6]  0.008219716

Accelerated Proximal Gradient Descent

• Update rule:

$$\begin{array}{cc}{u}^{(k+1)}& ={\beta }^{\left(k\right)}+\frac{k-1}{k+2}({\beta }^{\left(k\right)}-{\beta }^{(k-1)})\\ {\beta }^{(k+1)}& =S_{\alpha \lambda }(u^{\left(k\right)}-(\alpha /n\left)X^{\prime }\right(X{u}^{\left(k\right)}-y\left)\right)\end{array}$$

• Step size $\alpha =1/{\lambda }_{max}\left(n^{-1}{X}^{\prime }X\right)=n/{\lambda }_{max}\left(X^{\prime }X\right)$

Implementation

init_fista = function(x, y, lam)
{
    p = ncol(x)
    alpha = nrow(x) / svds(x, k = 1, nu = 0, nv = 0)$d^2
    list(alpha = alpha, old_beta = numeric(p), betau = numeric(p))
}
update_fista = function(beta, k, x, y, lam, state)
{
    n = nrow(x)
    alpha = state$alpha
    state$betau = if(k <= 1) beta else
        beta + (k - 1) / (k + 2) * (beta - state$old_beta)
    Xtr = crossprod(x, as.numeric(x %*% state$betau) - y)
    grad = as.numeric(Xtr) / n
    new_beta = soft_thresh(state$betau - alpha * grad, alpha * lam)
    state$old_beta = beta
    list(beta = new_beta, state = state)
}

Result

res_fista = run_opt(x, y, lam, init_fista, update_fista)

## Iter 0, loss = 0.494847
## Iter 100, loss = 0.027830
## Iter 200, loss = 0.027825
## Iter 300, loss = 0.027825
## Iter 400, loss = 0.027825
## Iter 500, loss = 0.027825
## Iter 600, loss = 0.027825
## Iter 700, loss = 0.027825
## Iter 800, loss = 0.027825
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_fista$beta)

## [1] -0.3862278  1.4304116 -0.2206913  1.1256895  0.1704706  0.8657180

tail(res_fista$beta)

## [1]  0.000000000 -0.003938369  0.000000000  0.000000000  0.000000000
## [6]  0.008219289

Douglas-Rachford Splitting

• Let $g\left(\beta \right)=\lambda \parallel \beta {\parallel }_1$, $f\left(\beta \right)=\frac{1}{2n}\parallel y-X\beta {\parallel }^2$. Then

$$\begin{array}{cc}{prox}_{\alpha g}\left(z\right)& =S_{\alpha \lambda }\left(z\right)\\ {prox}_{\alpha h}\left(\beta \right)& =(I+(\alpha /n\left)X^{\prime }X{)}^{-1}\right(\beta +\left(\alpha /n\right)X^{\prime }y)\end{array}$$

• Update rule:

$$\begin{array}{cc}{\beta }^{(k+1)}& =S_{\alpha \lambda }\left(u^{\left(k\right)}\right)\\ u^{(k+1)}& =u^{\left(k\right)}+A^{-1}(2{\beta }^{(k+1)}-u^{\left(k\right)}+(\alpha /n\left)X^{\prime }y\right)-{\beta }^{(k+1)}\end{array}$$ where $A=I+\left(\alpha /n\right)X^{\prime }X$

Implementation

chol_solve = function(R, b)
{
    y = backsolve(R, b, transpose = TRUE)  # R'*y = b       
    x = backsolve(R, y)                    # R*x = y
    x
}
init_drs = function(x, y, lam, alpha = 1)
{
    n = nrow(x)
    p = ncol(x)
    scaled_xty = alpha / n * as.numeric(crossprod(x, y))
    R = chol(diag(p) + alpha / n * crossprod(x))
    list(alpha = alpha, sxty = scaled_xty,
         R = as.matrix(R), betau = numeric(p))
}
update_drs = function(beta, k, x, y, lam, state)
{
    n = nrow(x)
    alpha = state$alpha
    new_beta = soft_thresh(state$betau, alpha * lam)
    rhs = 2 * new_beta - state$betau + state$sxty
    state$betau = state$betau + chol_solve(state$R, rhs) - new_beta
    list(beta = new_beta, state = state)
}

Result

res_drs1 = run_opt(x, y, lam, init_drs, update_drs, alpha = 1)

## Iter 0, loss = 1.454710
## Iter 100, loss = 0.053128
## Iter 200, loss = 0.041675
## Iter 300, loss = 0.034461
## Iter 400, loss = 0.029864
## Iter 500, loss = 0.028097
## Iter 600, loss = 0.027848
## Iter 700, loss = 0.027828
## Iter 800, loss = 0.027826
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_drs1$beta)

## [1] -0.3863738  1.4304559 -0.2205224  1.1254740  0.1704870  0.8652928

tail(res_drs1$beta)

## [1]  0.000000000 -0.003902763  0.000000000  0.000000000  0.000000000
## [6]  0.008261868

Result

res_drs10 = run_opt(x, y, lam, init_drs, update_drs, alpha = 10)

## Iter 0, loss = 1.454710
## Iter 100, loss = 0.027825
## Iter 200, loss = 0.027825
## Iter 300, loss = 0.027825
## Iter 400, loss = 0.027825
## Iter 500, loss = 0.027825
## Iter 600, loss = 0.027825
## Iter 700, loss = 0.027825
## Iter 800, loss = 0.027825
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_drs10$beta)

## [1] -0.3862278  1.4304117 -0.2206913  1.1256895  0.1704706  0.8657182

tail(res_drs10$beta)

## [1]  0.000000000 -0.003938377  0.000000000  0.000000000  0.000000000
## [6]  0.008219268

ADMM

• Let $f\left(\beta \right)=\lambda \parallel \beta {\parallel }_1$, $g\left(\beta \right)=\frac{1}{2n}\parallel y-X\beta {\parallel }^2$. Then

• Update rule:

$$\begin{array}{cc}{\beta }^{(k+1)}& =S_{\lambda /\rho }(z^{\left(k\right)}-u^{\left(k\right)})\\ z^{(k+1)}& =A^{-1}({\beta }^{(k+1)}+u^{\left(k\right)}+(n\rho {)}^{-1}{X}^{\prime }y)\\ u^{(k+1)}& =u^{\left(k\right)}+{\beta }^{(k+1)}-z^{(k+1)}\end{array}$$ where $A=I+(n\rho {)}^{-1}{X}^{\prime }X$

• A reparameterized DRS

Implementation

init_admm = function(x, y, lam, rho = 1)
{
    n = nrow(x)
    p = ncol(x)
    scaled_xty = as.numeric(crossprod(x, y)) / n / rho
    R = chol(diag(p) + crossprod(x) / n / rho)
    list(rho = rho, sxty = scaled_xty, R = as.matrix(R),
         betaz = numeric(p), betau = numeric(p))
}
update_admm = function(beta, k, x, y, lam, state)
{
    n = nrow(x)
    rho = state$rho
    new_beta = soft_thresh(state$betaz - state$betau, lam / rho)
    rhs = new_beta + state$betau + state$sxty
    state$betaz = chol_solve(state$R, rhs)
    state$betau = state$betau + new_beta - state$betaz
    list(beta = new_beta, state = state)
}

Result

res_admm1 = run_opt(x, y, lam, init_admm, update_admm, rho = 1)

## Iter 0, loss = 1.454710
## Iter 100, loss = 0.053098
## Iter 200, loss = 0.041664
## Iter 300, loss = 0.034453
## Iter 400, loss = 0.029859
## Iter 500, loss = 0.028096
## Iter 600, loss = 0.027848
## Iter 700, loss = 0.027828
## Iter 800, loss = 0.027826
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_admm1$beta)

## [1] -0.3863737  1.4304559 -0.2205226  1.1254743  0.1704871  0.8652934

tail(res_admm1$beta)

## [1]  0.000000000 -0.003902788  0.000000000  0.000000000  0.000000000
## [6]  0.008261821

Result

res_admm0.1 = run_opt(x, y, lam, init_admm, update_admm, rho = 0.1)

## Iter 0, loss = 1.454710
## Iter 100, loss = 0.027825
## Iter 200, loss = 0.027825
## Iter 300, loss = 0.027825
## Iter 400, loss = 0.027825
## Iter 500, loss = 0.027825
## Iter 600, loss = 0.027825
## Iter 700, loss = 0.027825
## Iter 800, loss = 0.027825
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_admm0.1$beta)

## [1] -0.3862278  1.4304117 -0.2206913  1.1256895  0.1704706  0.8657182

tail(res_admm0.1$beta)

## [1]  0.000000000 -0.003938377  0.000000000  0.000000000  0.000000000
## [6]  0.008219268

Coordinate Descent

• Each iteration contains a full cycle of ${\beta }_i$ updates

• Update rule:

$${\beta }_i\leftarrow S_{n\lambda /\parallel X_i{\parallel }^2}\left(\frac{X_i^{\prime }(y-X_{-i}{\beta }_{-i})}{\parallel X_i{\parallel }^2}\right)$$

• A "cleverer" scheme is to note that

$${\beta }_i\leftarrow S_{n\lambda /\parallel X_i{\parallel }^2}({\beta }_i+\frac{X_i^{\prime }r}{\parallel X_i{\parallel }^2}),r=y-X\beta$$

• Then update $r\leftarrow r-({\beta }_i^{new}-{\beta }_i^{old})X_i$

Implementation

# Soft-thresholding function
soft_thresh = function(x, a)
{
    sign(x) * pmax(0, abs(x) - a)
}
init_cd = function(x, y, lam)
{
    xi_square = colSums(x^2)
    thresh = nrow(x) * lam / xi_square
    list(xi_square = xi_square, thresh = thresh)
}
update_cd = function(beta, k, x, y, lam, state)
{
    new_beta = beta
    # betai <- S(betai + xi' * r / ||xi||^2, t), r = y - x * beta
    r = y - c(x %*% beta)
    for(i in seq_along(new_beta))
    {
        old = new_beta[i]
        newb = old + sum(x[, i] * r) / state$xi_square[i]
        thresh = state$thresh[i]
        new_beta[i] = soft_thresh(newb, thresh)
        r = r - sum(new_beta[i] - old) * x[, i]
    }
    list(beta = new_beta, state = state)
}

Result

res_cd = run_opt(x, y, lam, init_cd, update_cd)

## Iter 0, loss = 0.046708
## Iter 100, loss = 0.027825
## Iter 200, loss = 0.027825
## Iter 300, loss = 0.027825
## Iter 400, loss = 0.027825
## Iter 500, loss = 0.027825
## Iter 600, loss = 0.027825
## Iter 700, loss = 0.027825
## Iter 800, loss = 0.027825
## Iter 900, loss = 0.027825
## Iter 1000, loss = 0.027825

head(res_cd$beta)

## [1] -0.3862278  1.4304117 -0.2206913  1.1256895  0.1704706  0.8657182

tail(res_cd$beta)

## [1]  0.000000000 -0.003938377  0.000000000  0.000000000  0.000000000
## [6]  0.008219268

Comparison

K = 1000
gdat = data.frame(
    iter = rep(1:K, times = 8),
    loss = c(res_subgrad1$losses[1:K], res_subgrad10$losses[1:K],
             res_ista$losses[1:K], res_fista$losses[1:K],
             res_drs1$losses[1:K], res_drs10$losses[1:K],
             res_admm1$losses[1:K], res_cd$losses[1:K]),
    method = rep(c("Subgrad α=1", "Subgrad α=10", "ISTA", "FISTA",
                   "DRS α=1", "DRS α=10", "ADMM ρ=1", "CD"), each = K)
)

Comparison

library(ggplot2)
ggplot(gdat, aes(x = iter, y = log10(loss - obj_opt))) +
    geom_line(aes(color = method)) +
    scale_color_hue("Method") + xlab("Iteration") + ylab("Log-error") +
    theme_bw()

More Considerations

The previous comparison is neither rigorous nor complete, and more works need to be done:

• The plot only shows error v.s. iteration, but not error v.s. computing time

• The direct factorization of $A=I+\left(\alpha /n\right)X^{\prime }X$ is not optimal (recall that $p>n$)

• Design of algorithms when $X$ is sparse (i.e., operating on the xsp object)

Left as after-class exercises.

References