关于如何解决稀疏OLS - 如何在Matlab中应用`l1`最小化（教育目的）

Posted 2021-05-07

min_{a*x=y} +lambda*norm(hat{a},1)是目标函数，其中a是系数的向量，y表示噪声测量，x是未观察到的输入信号。我知道lasso()函数，但我不喜欢使用内置函数，因为这不会帮助我理解这些步骤。有人可以帮助实施l1规范优化吗？

在数学上，我的模型表示为移动平均（MA）系统：y[k] = a_1*x[k] + a_2*x[k-1] + a_{10}*x[k-9] + n[k]，其中n ~ N(0,sigma^2)是加性高斯白噪声，x是零均值高斯白色过程的单位方差，而a_1,a_2,...,a_10是MA模型的系数，已知是稀疏的。但是，我不知道稀疏系数的位置。

在该模型中，仅3个系数不为零，而其余系数均为零或接近零。进行参数估计的一种方法是构造逆滤波器或者也称为最小化预测误差。

通过逆滤波方法，我可以为MA模型创建一个逆滤波器，表示为：u[k] = x[k]-(hat{a_2}*x[k-1]+ hat{a_3}*x[k-3] + hat{a_{4}}*x[k-4] +ldots+hat{a_{10}}*x[k-9] )。

因此，目标函数变为：J = min_{hat{a}*x=y} +lambda*norm(hat{a},1) wherey是观察到的噪声测量值，而hat{a}*x是干净的。设mathbf{hat{a}} = {[hat{a_1},ldots,hat{a_{10}}]}^T表示估计的系数向量。

我的方法是将目标函数J分成两部分 - 第一部分是OLS估计，它被送入l1最小化程序。 l1最小化的输出给出稀疏系数。这种方法合法吗？如果是这样，我需要帮助什么是Matlab中的l1优化器？

以下是我创建模型的代码片段。但我不知道如何解决目标函数。请帮忙。

%Generate input
N=500;
x=(randn(1,N)*100);
L = 10;
Num_lags = 1:L-1;
a = 1+randn(L,1);
%Data preparation into regressors    
a(rand(L,1)<.9)=0; % 90 of the coefficients are zero
X1 = lagmatrix(x, [0 Num_lags]);

答案

下面的代码可以解决l1-optimization argmin{f(x)} s.t.||x||_1<=t。

编辑：更新了@V中的拼写错误（参见SKM的评论）

clc; clear;
%Generate input data
N=500;
Bnum=10;
X=(randn(N,Bnum)*1000);
true_beta = rand(Bnum,1);
Y=X*true_beta+rand(N,1);

%solve lasso using fminunc
lamda=1;
V = @(x) norm(Y-X*x)^2+lamda*norm(x,1);
options=optimoptions('fminunc','Algorithm','quasi-newton','Display','iter');
xopt = fminunc(V,zeros(Bnum,1),options)

但是，我仍然推荐使用QUADPROG函数的第二篇文章中的代码。它会更快，更准确。

Ref

另一答案

由于lasso paper，我们有：

它是具有线性约束的二次规划，并且存在套索参数t>=O，其控制应用于估计的收缩量。在你的情况下，我们可以假设t=0.1*sum(beta0)。 （beta0是完整的最小二乘估计;根据第3页的结论）

---

编辑：假设我们有一个带有2个参数的套索方程（1）。

1. 方程（1）：argmin{(y(1)-alpha-beta_1*x(1,1)-beta_2*x(2,1))^2+(y(2)-alpha-beta_1*x(1,2)-beta_2*x(2,2))^2} s.t. sum(abs(beta))<=t`
2. 由于hat(alpha)=mean(y)，我们通过定义argmin{(hat(y(1))-beta_1*x(1,1)-beta_2*x(2,1))^2+(hat(y(2))-beta_1*x(1,2)-beta_2*x(2,2))^2}得到了Eq（1）'：s.t. sum(abs(beta))<=t hat(y)=y-mean(y)。
3. （1）Eq（1）'可以重写为argmin{(beta)'H(beta)+f(beta)+C}，它的解决方案等于argmin{(beta)'H(beta)+f(beta)}
4. （2）s.t.部分等于2^(length(beta))约束（Ax=b），A是p元组（例如（-1，-1），（ - 1,1），（1，-1），（1,1） ）在b=(t,t)'的2参数情况下）。

---

然后我们可以通过MATLAB中的quadprog解决它：

• 通过将套索函数转换为二项式来计算矩阵H，f。
• 让delta_i，i = 1,2，...，2P是形式的p元组（+ -1，+ - 1，...，+ - 1），我们可以将我们的约束重写为Ax<b（delta_i*beta<=t） 。
• 通过调用函数quadprog(H,f,A,b)来解决套索估计。

这是我的代码：

clc; clear;
%Generate input
N=500;
Bnum=10;
x=(randn(N,Bnum)*1000);
true_beta = 1+randn(Bnum,1);
y=x*true_beta;

%find the solution of lm & t
%lm solution beta0
ls_beta=(x'*x)(x'*y);
%define t
t=0.1*sum(abs(ls_beta));

%%solving the quadratic programming
%calc H
H=zeros(Bnum);
for ii=1:Bnum
    for ij=1:Bnum
      H(ii,ij)=sum(2*x(:,ii).*x(:,ij));
    end
end
H;
%calc f
f=zeros(Bnum,1);
yhat=y-mean(y);
for ii=1:Bnum
  f(ii)=sum(-2*yhat.*x(:,ii));
end
f;
%calc A
A=zeros(power(2,Bnum),Bnum);
for ii=1:power(2,Bnum)
  v=dec2bin(ii-1,Bnum);
  for ij=1:Bnum
    A(ii,ij)=(str2num(v(ij))-0.5)*2;
  end
end
%calc b
b=ones(power(2,Bnum),1)*t;
%calc quadprog
lasso_beta=quadprog(H,f,A,b);
[true_beta ls_beta lasso_beta]

我们会得到一个答案（[true_beta ls_beta lasso_beta]）：

ans =

    0.5723    0.5723    0.0000
    0.4206    0.4206    0.0000
    1.9260    1.9260    0.1109
    1.0055    1.0055    0.0000
    0.3655    0.3655    0.0000
    1.8583    1.8583    0.1582
    0.5192    0.5192    0.0000
    2.4897    2.4897    0.7242
    0.3706    0.3706   -0.0000
    0.4060    0.4060    0.0000

所有重点都来自套索论文。

希望能帮助到你！