User:Lovebeloved/sandbox

Proximal gradient (forward backward splitting) methods for learning is an area of research in optimization and statistical learning theory which studies algorithms for a general class of convex regularization problems where the regularization penalty may not be differentiable. One such example is $$\ell_1$$ regularization (also known as Lasso) of the form
 * $$\min_{w\in\mathbb{R}^d} \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2+ \lambda \|w\|_1, \quad \text{ where } x_i\in \mathbb{R}^d\text{ and } y_i\in\mathbb{R}.$$

Proximal gradient methods offer a general framework for solving regularization problems from statistical learning theory with penalties that are tailored to a specific problem application. Such customized penalties can help to induce certain structure in problem solutions, such as sparsity (in the case of lasso) or group structure (in the case of group lasso).

Relevant background
Proximal gradient methods are applicable in a wide variety of scenarios for solving convex optimization problems of the form
 * $$ \min_{x\in \mathcal{H}} F(x)+R(x),$$

where $$F$$ is convex and differentiable with Lipschitz continuous gradient, $$ R$$ is a convex, lower semicontinuous function which is possibly nondifferentiable, and $$\mathcal{H}$$ is some set, typically a Hilbert space. The usual criterion of $$ x$$ minimizes $$ F(x)+R(x)$$ if and only if $$ \nabla (F+R)(x) = 0$$ in the convex, differentiable setting is now replaced by
 * $$ 0\in \partial (F+R)(x), $$

where $$\partial \varphi$$ denotes the subdifferential of a real-valued, convex function $$ \varphi$$.

Given a convex function $$\varphi:\mathcal{H} \to \mathbb{R}$$ an important operator to consider is its proximity operator $$\operatorname{prox}_{\varphi}:\mathcal{H}\to\mathcal{H} $$ defined by
 * $$ \operatorname{prox}_{\varphi}(u) = \operatorname{arg}\min_{x\in\mathcal{H}} \varphi(x)+\frac{1}{2}\|u-x\|_2^2,$$

which is well-defined because of the strict convexity of the $$ \ell_2$$ norm. The proximity operator can be seen as a generalization of a projection. We see that the proximity operator is important because $$ x^* $$ is a minimizer to the problem $$ \min_{x\in\mathcal{H}} F(x)+R(x)$$ if and only if
 * $$x^* = \operatorname{prox}_{\gamma R}\left(x^*-\gamma\nabla F(x^*)\right),$$ where $$\gamma>0$$ is any positive real number.

Moreau decomposition
One important technique related to proximal gradient methods is the Moreau decomposition, which decomposes the identity operator as the sum of two proximity operators. Namely, let $$\varphi:\mathcal{X}\to\mathbb{R}$$ be a lower semicontinuous, convex function on a vector space $$\mathcal{X}$$. We define its Fenchel conjugate $$\varphi^*:\mathcal{X}\to\mathbb{R}$$ to be the function
 * $$\varphi^*(u) := \sup_{x\in\mathcal{X}} \langle x,u\rangle - \varphi(x).$$

The general form of Moreau's decomposition states that for any $$x\in\mathcal{X}$$ and any $$\gamma>0$$ that
 * $$x = \operatorname{prox}_{\gamma \varphi}(x) + \gamma\operatorname{prox}_{\varphi^*/\gamma}(x/\gamma),$$

which for $$\gamma=1$$ implies that $$x = \operatorname{prox}_{\varphi}(x)+\operatorname{prox}_{\varphi^*}(x)$$. The Moreau decomposition can be seen to be a generalization of the usual orthogonal decomposition of a vector space, analogous with the fact that proximity operators are generalizations of projections.

In certain situations it may be easier to compute the proximity operator for the conjugate $$\varphi^*$$ instead of the function $$\varphi$$, and therefore the Moreau decomposition can be applied. This is the case for group lasso.

Lasso regularization
Consider the regularized empirical risk minimization problem with square loss and with the $\ell_1$ norm as the regularization penalty:
 * $$\min_{w\in\mathbb{R}^d} \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2+ \lambda \|w\|_1, $$

where $$x_i\in \mathbb{R}^d\text{ and } y_i\in\mathbb{R}.$$ The $$\ell_1$$ regularization problem is sometimes referred to as lasso (least absolute shrinkage and selection operator). Such $$\ell_1$$ regularization problems are interesting because they induce  sparse solutions, that is, solutions $$w$$ to the minimization problem have relatively few nonzero components. Lasso can be seen to be a convex relaxation of the non-convex problem
 * $$\min_{w\in\mathbb{R}^d} \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2+ \lambda \|w\|_0, $$

where $$\|w\|_0$$ denotes the $$\ell_0$$ "norm", which is the number of nonzero entries of the vector $$w$$. Sparse solutions are of particular interest in learning theory for interpretability of results: a sparse solution can identify a small number of important factors.

Solving for $$\ell_1$$ proximity operator
For simplicity we restrict our attention to the problem where $$\lambda=1$$. To solve the problem
 * $$\min_{w\in\mathbb{R}^d} \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2+ \|w\|_1, $$

we consider our objective function in two parts: a convex, differentiable term $$F(w) = \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2$$ and a convex function $$R(w) = \|w\|_1$$. Note that $$R$$ is not strictly convex.

Let us compute the proximity operator for $$R(w)$$. First we find an alternative characterization of the proximity operator $$\operatorname{prox}_{R}(x)$$ as follows:

$$ \begin{align} u = \operatorname{prox}_R(x) \iff & 0\in \partial \left(R(u)+\frac{1}{2}\|u-x\|_2^2\right)\\ \iff & 0\in \partial R(u) + u-x\\ \iff & x-u\in \partial R(u). \end{align} $$

For $$R(w) = \|w\|_1$$ it is easy to compute $$\partial R(w)$$: the $$i$$th entry of $$\partial R(w)$$ is precisely


 * $$ \partial |w_i| = \begin{cases}

1,&w_i>0\\ -1,&w_i<0\\ \left[-1,1\right],&w_i = 0. \end{cases}$$

Using the recharacterization of the proximity operator given above, for the choice of $$R(w) = \|w\|_1$$ and $$\gamma>0$$ we have that $$\operatorname{prox}_{\gamma R}(x)$$ is defined entrywise by


 * $$\left(\operatorname{prox}_{\gamma R}(x)\right)_i = \begin{cases}

x_i-\gamma,&x_i>\gamma\\ 0,&|x_i|\leq\gamma\\ x_i+\gamma,&x_i<-\gamma, \end{cases}$$

which is known as the soft thresholding operator $$S_{\gamma}(x)=\operatorname{prox}_{\gamma \|\cdot\|_1}(x)$$.

Fixed point iterative schemes
To finally solve the lasso problem we consider the fixed point equation shown earlier:
 * $$x^* = \operatorname{prox}_{\gamma R}\left(x^*-\gamma\nabla F(x^*)\right).$$

Given that we have computed the form of the proximity operator explicitly, then we can define a standard fixed point iteration procedure. Namely, fix some initial $$w^0\in\mathbb{R}^d$$, and for $$k=1,2,\ldots$$ define
 * $$w^{k+1} = S_{\gamma}\left(w^k - \gamma \nabla F\left(w^k\right)\right).$$

Note here the effective trade-off between the empirical error term $$F(w) $$ and the regularization penalty $$R(w)$$. This fixed point method has decoupled the effect of the two different convex functions which comprise the objective function into a gradient descent step ($$ w^k - \gamma \nabla F\left(w^k\right)$$) and a soft thresholding step (via $$S_\gamma$$).

Convergence of this fixed point scheme is well-studied in the literature and is guaranteed under appropriate choice of step size $$\gamma$$ and loss function (such as the square loss taken here). Accelerated methods were introduced by Nesterov in 1983 which improve the rate of convergence under certain regularity assumptions on $$F$$. Such methods have been studied extensively in previous years. For more general learning problems where the proximity operator cannot be computed explicitly for some regularization term $$R$$, such fixed point schemes can still be carried out using approximations to both the gradient and the proximity operator.

Practical considerations
There have been numerous developments within the past decade in convex optimization techniques which have influenced the application of proximal gradient methods in statistical learning theory. Here we survey a few important topics which can greatly improve practical algorithmic performance of these methods.

Adaptive step size
In the fixed point iteration scheme
 * $$w^{k+1} = \operatorname{prox}_{\gamma R}\left(w^k-\gamma \nabla F\left(w^k\right)\right),$$

one can allow variable step size $$\gamma_k$$ instead of a constant $$\gamma$$. Numerous adaptive step size schemes have been proposed throughout the literature. Applications of these schemes suggest that these can offer substantial improvement in number of iterations required for fixed point convergence.

Elastic net (mixed norm regularization)
Elastic net regularization offers an alternative to pure $$\ell_1$$ regularization. The problem of lasso ($$\ell_1$$) regularization involves the penalty term $$R(w) = \|w\|_1$$, which is not strictly convex. Hence, solutions to $$\min_w F(w) + R(w),$$ where $$F$$ is some empirical loss function, need not be unique. This is often avoided by the inclusion of an additional strictly convex term, such as an $$\ell_2$$ norm regularization penalty. For example, one can consider the problem
 * $$\min_{w\in\mathbb{R}^d} \frac{1}{n}\sum_{i=1}^n (y_i- \langle w,x_i\rangle)^2+ \lambda \left((1-\mu)\|w\|_1+\mu \|w\|_2\right), $$

where $$x_i\in \mathbb{R}^d\text{ and } y_i\in\mathbb{R}.$$ For $$0<\mu\leq 1$$ the penalty term $$\lambda \left((1-\mu)\|w\|_1+\mu \|w\|_2\right)$$ is now strictly convex, and hence the minimization problem now admits a unique solution. It has been observed that for sufficiently small $$\mu > 0$$, the additional penalty term $$\mu \|w\|_2$$ acts as a preconditioner and can substantially improve convergence while not adversely affecting the sparsity of solutions.

Exploiting group structure
Proximal gradient methods provide a general framework which is applicable to a wide variety of problems in statistical learning theory. Certain problems in learning can often involve data which has additional structure that is known  a priori. In the past several years there have been new developments which incorporate information about group structure to provide methods which are tailored to different applications. Here we survey a few such methods.

Group lasso
Group lasso is a generalization of the lasso method when features are grouped into disjoint blocks. Suppose the features are grouped into blocks $$\{w_1,\ldots,w_G\}$$. Here we take as a regularization penalty


 * $$R(w) =\sum_{g=1}^G \|w_g\|_2,$$

which is the sum of the $$\ell_2$$ norm on corresponding feature vectors for the different groups. A similar proximity operator analysis as above can be used to compute the proximity operator for this penalty. Where the lasso penalty has a proximity operator which is soft thresholding on each individual component, the proximity operator for the group lasso is soft thresholding on each group. For the group $$w_g$$ we have that proximity operator of $$\lambda\gamma\left(\sum_{g=1}^G \|w_g\|_2\right) $$ is given by


 * $$\widetilde{S}_{\lambda\gamma }(w_g) = \begin{cases}

w_g-\lambda\gamma \frac{w_g}{\|w_g\|_2}, & \|w_g\|_2>\lambda\gamma \\ 0, & \|w_g\|_2\leq \lambda\gamma \end{cases}$$

where $$w_g$$ is the $$g$$th group.

In contrast to lasso, the derivation of the proximity operator for group lasso relies on the Moreau decomposition. Here the proximity operator of the conjugate of the group lasso penalty becomes a projection onto the ball of a dual norm.

Other group structures
In contrast to the group lasso problem, where features are grouped into disjoint blocks, it may be the case that grouped features are overlapping or have a nested structure. Such generalizations of group lasso have been considered in a variety of contexts. For overlapping groups one common approach is known as latent group lasso which introduces latent variables to account for overlap. Nested group structures are studied in hierarchical structure prediction and with directed acyclic graphs.

Method for learning matrices
Proximal operator can be extended from vector function to matrix function when the learning object is a matrix instead of a vector.

Learning matrices
A general framework for such setting is the following: consider a training set $$ S=(X_i^t,y_i^t) $$ where $$X_i^t,W \in \mathbb{R}^{D \times T}, y_i^t, \epsilon_i^t \in \mathbb{R}$$ for $$i=1,...,n_t, t=1,...,T$$, and assume the regression model
 * $$ y_i^t= \langle W,X_i^t\rangle_F + \epsilon_i^t $$

where $$ \langle W,X_i^t\rangle_F = Tr(W'X_i^t) $$ is the Frobenuis (Hilbert-Schmidt) inner product. Typical examples of learning matrices are the special cases of the above general model:

Let $$ (e_t)_t $$ is the canonical basis in $$ \mathbb{R}^T $$,
 * $$ X_i^t =e_t \otimes x_i^t $$    $$\rightarrow  $$    Linear Multi-task learning:           $$ y_{i,t} = \langle W^t, x_i^t\rangle _{\mathbb{R}^D} + \epsilon_{i,t} $$      where $$ W^t \in \mathbb{R}^D $$ is the regression vector for $$ x_i^t \in \mathbb{R}^D $$.
 * $$ X_i^t =e_t \otimes x_i $$    $$\rightarrow  $$    Linear Multivariate Regression:   $$ y_i=W'x_i + \epsilon_i $$                       where $$ y_i, \epsilon_i \in \mathbb{R}^T $$ and $$ x_i \in \mathbb{R}^D $$
 * $$ X_i^t =e_t \otimes e_i' $$    $$\rightarrow  $$    Matrix completion:                        $$ y_{i,t}=W_{i,t} + \epsilon_{i,t} $$                     where $$ (e_i')_i $$ is the canonical basis in $$ \mathbb{R}^D $$

The corresponding penalized regularization scheme for the model is:
 * $$ \min_{W \in \mathcal{H}} \hat {\mathcal{E}} (W) + R(W) $$

where $$\mathcal{H} $$ is the space induced by Hilbert-Schmidt operator, $$\hat{\mathcal{E}}$$ is an empirical error defined by a loss function $$ V: \mathbb{R} \times \mathbb{R} \rightarrow [0, \infty)$$ (e.g., the square loss function), and $$R(W)$$ is a regularizer. As in the vector version of the problem, we can consider $$R(W)=\lambda \| W \|_p$$.

Proximal operator for matrix norm
Different matrix norm will lead to different proximal operator. We discuss the proximal operator associated with two types of matrix norm below, and for simplicity, we assume $$\lambda=1$$.

Entrywise norm
An entrywise $$p$$-norm is defined as:
 * $$ \| W \|_p = (\sum_{i=1}^{D} \sum_{j=1}^{T} \| W_{i,j} \|^p)^\frac{1}{p} $$

Since it simply treats matrix $$W \in \mathbb{R}^{D \times T}$$ as a vector $$\in \mathbb{R}^{DT}$$, the proximal operator is the same as that of a vector function. For example, when $$ p $$ =1, the proximal operator is entrywise soft thresholding.

Schatten norm
The Schatten $$p$$-norm is defined as:
 * $$ \| W \|_p =  (\sum_{j=1}^{m} | \sigma_j |^p)^\frac{1}{p} $$

where $$\sigma_1,...,\sigma_m \ge 0, m = \min(T,D)$$ are the singular values of $$W$$.

Let $$\sigma$$ be the singular value map: $$\mathbb{R}^{D \times T} \rightarrow \mathbb{R}^m$$, the function that takes a matrix and returns a vector of its singular values in nonincreasing order, then
 * $$W=U\textbf{diag}(\sigma(W))V$$

where $$U \in \mathbb{R}^{D \times D}$$ and $$V \in \mathbb{R}^{T \times T}$$ are unitary matrices, and $$\textbf{diag} (\sigma(W)) \in \mathbb{R}^{D \times T}$$.

Schatten $$p$$-norm is unitarily invariant such that $$\| W \|_p= \| \textbf{diag}(\sigma(W))\|_p$$.

It can be further shown that the sub-differential of $$R(W)$$ is:


 * $$\partial (\| W \|_p) =U \textbf{diag}[\partial (\| \sigma (W)\|_p) ]V $$

This implies that:


 * $$ \textbf{prox}_R(W) = U \textbf{diag}[\textbf{prox}_{\|.\|_p}(\sigma(W))]V $$

For example, when $$p$$=1 (nuclear norm), the proximal operator is the singular value thresholding:
 * $$ \textbf{prox}_R(W) = \sum_{i=1}^{m} (\sigma_i-1)_+u_iv_i^T $$

where $$(\sigma_i-1)_+$$ is the soft thresholding operator on $$\sigma_i, u_i$$ is the $$i^{th}$$ row of $$U$$, and $$v_i^T$$ is the $$i^{th} $$ column of $$V$$.

Proximal gradient method for learning matrices
The above results suggest an iterative proximal gradient method for learning matrices:

$$ W_{t+1}=prox_{\lambda R}(W_t-\gamma \nabla \hat{\mathcal{E}}(W_t)), t=0,..., $$

for a given some initialization $$ W_0 $$.