l1约束的最小二乘学习

$\ell_{1}$Constrained Least Squares
In sparse learning, $\ell_{1}$ constrained LS, also known as Lasso Regression, is a common learning method:

$$\min_{\theta} J_{LS}(\theta)\quad s.t. \quad \|\theta \|_{1}\leq R$$ where $$\|\theta\|_{1}=\sum_{j=1}^{b}|\theta_{j}|$$
Generally speaking, the solution of an $\ell_{1}$ constrained LS is located on the axis, that is to say, there are several parameters $\theta_{j}$ equal to zero (sparse).
Then how to solve it? Given the indifferentiable property of the absolute value at the origin, solving an $\ell_{1}$ constrained LS is not so easy as solving the $\ell_{2}$ constrained one. However, we can still apply Lagrange multiplier.
$$\min_{\theta}J(\theta),\quad J(\theta)=J_{LS}(\theta)+\lambda\|\theta\|_{1}$$
Note that
$$|\theta_{j}|\leq\frac{\theta_{j}^{2}}{2c_{j}}+\frac{c_{j}}{2}, \quad c_{j}>0$$ i.e. we can optimize the upper-bound of $J(\theta)$. By iteration, we take the current solution $\tilde{\theta}_{j}\neq 0$ as $c_{j}$ so as to formulate the upper bound constraint:
$$|\theta_{j}|\leq\frac{\theta_{j}^{2}}{2|\tilde{\theta}_{j}|}+\frac{|\tilde{\theta}_{j}|}{2}$$ If $\tilde{\theta}_{j}=0$, we should take $\theta_{j}=0$. When we use general inverse, the inequality above can be referred as:
$$|\theta_{j}|\leq\frac{|\tilde{\theta}_{j}|^{\dagger}}{2}\theta_{j}^{2}+\frac{|\tilde{\theta}_{j}|}{2}$$
Therefore, we can get the following $\ell_{2}$ regularized constrained LS problem formulation:
$$\hat{\theta}=\arg\min_{\theta}\tilde{J}(\theta),\quad \tilde{J}(\theta)=J_{LS}(\theta)+\frac{\lambda}{2}\theta^{T}\tilde{\Theta}^{\dagger}\theta+C$$
where $\tilde{\Theta}= \begin{pmatrix} |\tilde{\theta}_{1}| & & \\ & \ddots & \\ & & |\tilde{\theta}_{b}| \end{pmatrix}$ and $C=\sum_{j=1}^{b} |\tilde{\theta}_{j}|/2$ are independent of $\theta$.

Take the parameterized linear model for example

$$f_{\theta}(x)=\theta^{T}\phi(x)$$
Then, by the use of Lagrange multiplier, we can get
$$\hat{\theta}=\left(\Phi^{T}\Phi+\lambda\tilde{\Theta}^{\dagger}\right)^{-1}\Phi y$$
Renew the estimation $\tilde{\theta}$ as $\tilde{\theta}=\hat{\theta}$, go back to calculate the new $\hat{\theta}$ until $\hat{\theta}$ comes to the required precision.

For simplicity, the whole algorithm goes as follows:

1. Initialize $\theta_{0}$ and $i=1$.
2. Calculate $\Theta_{i}$ using $\theta_{i-1}$.
3. Estimate $\theta_{i}$ using $\Theta_{i}$.
4. $i=i+1$, go back to step 2.

An Example:

n=50; N=1000;
x=linspace(-3,3,n)'; X=linspace(-3,3,N)';
pix=pi*x;
y=sin(pix)./(pix)+0.1*x+0.2*rand(n,1);

hh=2*0.3^2; l=0.1; t0=randn(n,1); x2=x.^2;
k=exp(-(repmat(x2,1,n)+repmat(x2',n,1)-2*(x*x'))/hh);
k2=k^2; ky=k*y;
for o=1:1000
    t=(k2+l*pinv(diag(abs(t0))))\ky;
    if norm(t-t0)<0.001, break, end
    t0=t;
end
K=exp(-(repmat(X.^2,1,n)+repmat(x2',N,1)-2*X*x')/hh);
F=K*t;

figure(1); clf; hold on; axis([-2.8,2.8,-1,1.5]);
plot(X,F,'g-'); plot(x,y,'bp');