稀疏OLS的求解--如何在Matlab中应用“11”最小化(教育目的)

Question

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)，其中y是观测到的噪声测量值，\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]);

Answer 1

下面的代码可以解决L1优化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函数的第二篇文章中的代码。它会更快更准确。

参考

Answer 2

由于套索纸，我们有：

​

​

这是一个具有线性约束的二次规划，并且有一个lasso参数t>=O，它控制用于估计的收缩量。在您的情况下，我们可以假设t=0.1*sum(beta0)。(beta0是完全最小二乘估计；根据第3页的结论)

编辑：假设我们有一个两个参数的拉索方程(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)，通过定义hat(y)=y-mean(y)，我们得到了Eq(1)‘：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} s.t. sum(abs(beta))<=t。
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是具有b=(t,t)'的p元组(例如(-1，-1)，(-1,1)，(1，-1)，(1，-1)，(1，-1))。

然后，我们就可以在MATLAB中用四象限来解决这个问题了：

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

这是我的代码：

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

所有的重点都来自套索纸。

希望能帮上忙！