First Order Methods

This module aims to implement the first order methods for sparse signal recovery problems proposed in [BCandesG11]. The implementation is adapted from TFOCS [BCG12].

First order methods exploit information n values and gradients/subgradients (but not Hessians) of the functions comprising the models under consideration [Bec17].

API

Solvers

fom(smooth_f, prox_h, A, b, x0[, options])	First order methods driver routine
scd(prox_f, conj_neg_h, A, b, mu, x0, z0[, …])	First order solver for smooth conic dual problems driver routine
l1rls(A, b, lambda_, x0[, options])	Solver for l1 regulated least square problem
l1rls_jit(A, b, lambda_, x0[, options])	Solver for l1 regulated least square problem
lasso(A, b, tau, x0[, options])	Solver for LASSO problem
lasso_jit(A, b, tau, x0[, options])	Solver for LASSO problem
owl1rls(A, b, lambda_, x0[, options])	Solver for ordered weighted l1 norm regulated least square problem
owl1rls_jit(A, b, lambda_, x0[, options])	Solver for ordered weighted l1 norm regulated least square problem
dantzig_scd(A, b, delta, mu, x0, z0[, options])	Solver for the (smoothed) Dantzig selector problem using smoothed conic dual formulation

Data types

FomOptions	Options for FOCS driver routine
FomState	State of the FOCS method

Utilities

matrix_affine_func([A, b])

Returns an affine function for a matrix A and vector b

First Order Methods

[BCandesG11] considers problems in the unconstrained form:

(1)$$\text{minimize }\varphi (x)=g(x)+h(x)$$

where:

• Both $g,h:{\mathbb{R}}^n\to (\mathbb{R}\cup \mathrm{\infty })$ are convex functions.

• $g:{\mathbb{R}}^m\to \mathbb{R}$ is a smooth convex function.

• $h:{\mathbb{R}}^n\to \mathbb{R}$ is a non-smooth convex function.

Of particular interest are problems where $g$ takes a specific form:

(2)$$g(x)=f(\mathcal{A}(x)+b)$$

where

• $\mathcal{A}$ is a linear operator from ${\mathbb{R}}^n\to {\mathbb{R}}^m$

• $b$ is a translation vector

• $f:{\mathbb{R}}^m\to \mathbb{R}$ is a smooth convex function

and the function $h$ is a prox-capable convex function (to be described later).

We can then rewrite $\varphi$ as:

(3)$$\text{minimize }\varphi (x)=f(\mathcal{A}(x)+b)+h(x)$$

For a smooth function, its gradient $\mathrm{\nabla }f$ must exist and be easy to compute.

For a prox-capable function, it should have an efficient proximal operator:

(4)$$p_f(x,t)=\underset{z\in {\mathbb{R}}^n}{\text{arg}min}(f(x)+\frac{1}{2t}\Vert z-x{\Vert }_2^2)$$

for any $x\in {\mathbb{R}}^n$ and the step size $t>0$.

See the following sections for details:

• Linear Operators

• Smooth Function Generators

• Proximal Operator Generators

The routine fom() provides the solver for the minimization problem described above.

• Unconstrained smooth minimization problems can be handled by choosing $h(x)=0$. See cr.sparse.opt.prox_zero().

• Convex constraints can be handled by adding their indicator functions as part of $h$.

For solving the minimization problem, an initial solution $x_0$ should be provided. If one is unsure of the initial solution, they can provide $0\in {\mathbb{R}}^n$ as the initial solution.

Smooth Conic Dual Problems

A number of sparse signal recovery problems can be cast as conic optimization problems. Then efficient first order methods can be used to solve the problem. This procedure often involves the following steps:

• Cast the sparse recovery problem as a conic optimization problem

• Augment the objective function with a strongly convex term to make it smooth.

• Formulate the dual of the smooth problem.

• Solve the smooth conic dual problem using a first order method.

• Apply continuation to mitigate the effect of the strongly convex term added to the original objective function.

Our goal is to solve the conic problems of the form [BCandesG11]

(5)$$\begin{array}{r}\begin{array}{rlrl}& \underset{x}{\text{minimize}}& & f(x)\\ & \text{subject to}& & \mathcal{A}(x)+b\in \mathcal{K}\end{array}\end{array}$$

where:

• $x\in {\mathbb{R}}^n$ is the optimization variable

• $f:{\mathbb{R}}^n\to \mathbb{R}$ is a convex objective function which is possibly extended valued but not necessarily smooth

• $\mathcal{A}:{\mathbb{R}}^n\to {\mathbb{R}}^m$ is a linear operator

• $b\in {\mathbb{R}}^m$ is a translation vector

• $\mathcal{K}\subseteq {\mathbb{R}}^m$ is a closed, convex cone

Dual Problem

The Lagrangian associated with the conic problem is given by:

(6)$$\mathcal{L}(x,\lambda )=f(x)-⟨\lambda ,\mathcal{A}(x)+b⟩$$

where $\lambda \in {\mathcal{R}}^m$ is the Lagrange multiplier which must lie in the dual cone ${\mathcal{K}}^{\ast }$.

The dual function is:

(7)$$g(\lambda )=\underset{x\in {\mathbb{R}}^n}{inf}\mathcal{L}(x,\lambda )=-f^{\ast }({\mathcal{A}}^{\ast }(\lambda ))-⟨b,\lambda ⟩$$

• ${\mathcal{A}}^{\ast }:{\mathbb{R}}^m\to {\mathbb{R}}^n$ is the adjoint of the linear operator $\mathcal{A}$

• $f^{\ast }:{\mathbb{R}}^n\to (\mathbb{R}\cup \mathrm{\infty })$ is the convex conjugate of $f$ defined by

(8)$$f^{\ast }(z)=\underset{z}{sup}⟨z,x⟩-f(x)$$

Thus, the dual problem is given by:

(9)$$\begin{array}{r}\begin{array}{rlrl}& \underset{\lambda }{\text{maximize}}& & -f^{\ast }({\mathcal{A}}^{\ast }(\lambda ))-⟨b,\lambda ⟩\\ & \text{subject to}& & \lambda \in {\mathcal{K}}^{\ast }\end{array}\end{array}$$

Smoothing

The dual function $g$ may not be differentiable (or finite) on all of ${\mathcal{K}}^{\ast }$.

We perturb the original problem as follows:

(10)$$\begin{array}{r}\begin{array}{rlrl}& \underset{x}{\text{minimize}}& & f_{\mu }(x)\triangleq f(x)+\mu d(x)\\ & \text{subject to}& & \mathcal{A}(x)+b\in \mathcal{K}\end{array}\end{array}$$

• $\mu >0$ is a fixed smoothing parameter

• $d:{\mathbb{R}}^n\to \mathbb{R}$ is a strongly convex function obeying

(11)$$d(x)\ge d(x_0)+\frac{1}{2}\Vert x-x_0{\Vert }_2^2$$

A specific choice for $d(x)$ is:

(12)$$d(x)=\frac{1}{2}\Vert x-x_0{\Vert }_2^2$$

The new objective function $f_{\mu }$ is strongly convex.

The Lagrangian becomes:

(13)$${\mathcal{L}}_{\mu }(x,\lambda )=f(x)+\mu d(x)-⟨\lambda ,\mathcal{A}(x)+b⟩$$

The dual function becomes:

(14)$$g_{\mu }(\lambda )=\underset{x\in {\mathbb{R}}^n}{inf}{\mathcal{L}}_{\mu }(x,\lambda )=-(f+\mu d{)}^{\ast }({\mathcal{A}}^{\ast }(\lambda ))-⟨b,\lambda ⟩$$

First order methods may be employed to solve this problem with provable performance.

The smoothed conic dual problem is then given by:

(15)$$\begin{array}{r}\begin{array}{rlrl}& \underset{\lambda }{\text{maximize}}& & -(f+\mu d{)}^{\ast }({\mathcal{A}}^{\ast }(\lambda ))-⟨b,\lambda ⟩\\ & \text{subject to}& & \lambda \in {\mathcal{K}}^{\ast }\end{array}\end{array}$$

Composite Forms

Often, it is convenient to split the dual variable $\lambda \in {\mathbb{R}}^m$ as $\lambda =(z,\tau )$ where $z\in {\mathbb{R}}^{m-\overline{m}}$ and $\tau \in {\mathbb{R}}^{\overline{m}}$. For example, if ${\mathcal{K}}^{\ast }={\mathcal{L}}_1^n$, the ${\ell }_1$ norm cone (i.e. the epigraph of the ${\ell }_1$ norm function, then, if we split $\lambda \in {\mathbb{R}}^{n+1}$ as $\lambda =(z,\tau )$ where $z\in {\mathbb{R}}^n$ and $\tau \in \mathbb{R}$, then:

(16)$$\lambda \in {\mathcal{L}}_1^n⟺\Vert z{\Vert }_1\le \tau$$

This is particularly useful, if the smoothed dual function $g_{\mu }(\lambda )=g_{\mu }(z,\tau )$ is linear in $\tau$. Then, we can rewrite $g_{\mu }$ as:

(17)$$g_{\mu }(\lambda )=-g_{\text{sm}}(z)-⟨{\nu }_{\mu },\tau ⟩$$

Since $g_{\mu }$ is a smooth concave function, $g_{\text{sm}}$ will be a smooth convex function.

Problem Instantiations

In the rest of the document, we will discuss how specific sparse signal recovery problems can be modeled and solved using this module as either primal problems or a smooth conic dual problem.

L1 regularized least square problem

We consider the problem:

(18)$$\text{minimize}\frac{1}{2}\Vert \mathcal{A}x-b{\Vert }_2^2+\lambda \Vert x{\Vert }_1$$

Choose:

• $f(x)=\frac{1}{2}\Vert x{\Vert }_2^2$ see cr.sparse.opt.smooth_quad_matrix()

• $h(x)=\Vert \lambda x{\Vert }_1$ see cr.sparse.opt.prox_l1()

• $\mathcal{A}$ as the linear operator

• $-b$ as the translate input

With these choices, it is straight-forward to use fom() to solve the L1RLS problem. This is implemented in the function l1rls().

LASSO

LASSO (least absolute shrinkage and selection operator) s a regression analysis method that performs both variable selection and regularization in order to enhance the prediction accuracy and interpretability of the resulting statistical model.

We consider the problem:

(19)$$\begin{array}{r}\begin{array}{c}\underset{x}{\text{minimize}}\frac{1}{2}\Vert \mathcal{A}x-b{\Vert }_2^2\\ \text{subject to }\Vert x{\Vert }_1\le \tau \end{array}\end{array}$$

Choose:

• $f(x)=\frac{1}{2}\Vert x{\Vert }_2^2$ see cr.sparse.opt.smooth_quad_matrix()

• $h(x)=I_C(x)$ as the indicator function for l1-ball $C=\{x:\Vert x{\Vert }_1\le \tau \}$, see cr.sparse.opt.prox_l1_ball()

• $\mathcal{A}$ as the linear operator

• $-b$ as the translate input

With these choices, it is straight-forward to use fom() to solve the LASSO problem. This is implemented in the function lasso().

Ordered weighted L1 regularized least square problem

We consider the problem:

(20)$$\underset{x\in {\mathbb{R}}^n}{\text{minimize}}\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\sum _{i=1}^n{\lambda }_i|x{|}_{(i)}$$

described in [lBvdBSC13].

Choose:

• $f(x)=\frac{1}{2}\Vert x{\Vert }_2^2$ see cr.sparse.opt.smooth_quad_matrix()

• $h(x)=\sum _{i=1}^n{\lambda }_i|x{|}_{(i)}$ see cr.sparse.opt.prox_owl1()

• $\mathcal{A}$ as the linear operator

• $-b$ as the translate input

With these choices, it is straight-forward to use fom() to solve the ordered weighted L1 regularized least square problem. This is implemented in the function owl1rls().

The Dantzig selector

We wish to recover an unknown vector $x_0\in {\mathbb{R}}^n$ from the data $y\in {\mathbb{R}}^m$ and the model:

(21)$$y=A{x}_0+z$$

where

• $A$ is a known $(m,n)$ sized design matrix

• $z$ is the noise term

• There are fewer observations/measurements than unknowns ($m<n$)

We consider the optimization problem

(22)$$\begin{array}{r}\begin{array}{rlrl}& \underset{x}{\text{minimize}}& & \Vert x{\Vert }_1\\ & \text{subject to}& & \Vert A^{\ast }(y-Ax){\Vert }_{\mathrm{\infty }}\le \delta \end{array}\end{array}$$

• We are minimizing the ${\ell }_1$ norm of the solution (thus promoting sparsity)

• We have an upper bound $\delta$ on the correlation between the residual vector $r=y-Ax$ and the columns/atoms of $A$.

This formulation is known as the Dantzig selector.

The conic form

We can rewrite the Dantzig selector problem as:

(23)$$\begin{array}{r}\begin{array}{rlrl}& \underset{x}{\text{minimize}}& & \Vert x{\Vert }_1\\ & \text{subject to}& & (A^{\ast }(y-Ax),\delta )\in {\mathcal{L}}_{\mathrm{\infty }}^n\end{array}\end{array}$$

where ${\mathcal{L}}_{\mathrm{\infty }}^n$ is the epigraph of the ${\ell }_{\mathrm{\infty }}$-norm. It is a convex cone.

Then the Dantzig selector can be modeled as a standard conic form as follows:

• $f(x)=\Vert x{\Vert }_1$

• $\mathcal{A}(x)=(-A^{\ast }Ax,0)$ is a mapping from ${\mathbb{R}}^n\to {\mathbb{R}}^{n+1}$

• $b=(A^{\ast }y,\delta )$; note that $b\in {\mathbb{R}}^{n+1}$

• $\mathcal{K}={\mathcal{L}}_{\mathrm{\infty }}^n$

Note carefully that:

(24)$$\mathcal{A}(x)+b=(-A^{\ast }Ax,0)+(A^{\ast }y,\delta )=(A^{\ast }(y-Ax),\delta )$$

Thus:

(25)$$\mathcal{A}(x)+b\in \mathcal{K}={\mathcal{L}}_{\mathrm{\infty }}^n⟺\Vert A^{\ast }(y-Ax){\Vert }_{\mathrm{\infty }}\le \delta$$

Dual problem

The dual of the ${\mathcal{L}}_{\mathrm{\infty }}^n$ cone is the ${\mathcal{L}}_1^n$ cone. The dual variable $\lambda$ will lie in the dual cone ${\mathcal{L}}_1^n$. It will be easier to work with defining $\lambda =(z,\tau )$ such that

(26)$$\lambda \in {\mathcal{L}}_1^n⟺\Vert z{\Vert }_1\le \tau$$

The convex conjugate of the l1 norm function $f(x)=\Vert x{\Vert }_1$ is the indicator function for the ${\ell }_{\mathrm{\infty }}$ norm ball.

(27)$$\begin{array}{r}{f}^{\ast }(x)=I_{{\ell }_{\mathrm{\infty }}}(x)=\{\begin{array}{ll}0& \text{if }\Vert x{\Vert }_{\mathrm{\infty }}\le 1\\ \mathrm{\infty }& \text{otherwise}\end{array}\end{array}$$

The adjoint of the linear operator is given by:

(28)$${\mathcal{A}}^{\ast }(\lambda )={\mathcal{A}}^{\ast }((z,\tau ))=-A^{\ast }Az$$

Plugging into the dual problem definition:

(29)$$\begin{array}{r}\begin{array}{c}\begin{array}{rlrl}& \underset{\lambda }{\text{maximize}}& & -f^{\ast }({\mathcal{A}}^{\ast }(\lambda ))-⟨b,\lambda ⟩\\ & \text{subject to}& & \lambda \in {\mathcal{K}}^{\ast }\end{array}\end{array}\end{array}$$

We get:

(30)$$\begin{array}{r}\begin{array}{c}\begin{array}{rlrl}& \underset{z}{\text{maximize}}& & -I_{{\ell }_{\mathrm{\infty }}}(-A^{\ast }Az)-⟨A^{\ast }y,z⟩-\delta \tau \\ & \text{subject to}& & (z,\tau )\in {\mathcal{L}}_1^n\end{array}\end{array}\end{array}$$

For $\delta >0$, the solution $\lambda$ will be on the boundary of the convex cone ${\mathcal{L}}_1^n$, giving us:

(31)$$\Vert z{\Vert }_1=\tau$$

Plugging this back, we get the unconstrained maximization problem:

(32)$$\underset{z}{\text{maximize}}-I_{{\ell }_{\mathrm{\infty }}}(-A^{\ast }Az)-⟨A^{\ast }y,z⟩-\delta \Vert z{\Vert }_1$$