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

Question

min_{a*x=y} +lambda*norm(\hat{a},1) 是 objective 函数，其中 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] )。

因此，objective 函数变为：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 表示估计的系数向量。

我的方法是将 objective 函数 J 分成两部分——第一部分是 OLS 估计，它被输入到 l1 最小化例程中。 l1 最小化的输出给出了稀疏系数。做法合法吗？如果是这样，我需要帮助了解什么是 Matlab 中的 l1 优化器？

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

%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

由于 lasso paper，我们有：

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

编辑： 假设我们有一个带 2 个参数的套索方程 (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 `
由于 hat(alpha)=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 通过定义 hat(y)=y-mean(y).
(1) Eq(1) ' 可以改写为 argmin{(beta)'H(beta)+f(beta)+C} 其解等于 argmin{(beta)'H(beta)+f(beta)}
(2) s.t. 部分等于 2^(length(beta)) 约束 (Ax=b)，而 A 是一个 p-tuples（例如（-1 ,-1),(-1,1),(1,-1),(1,1) 在 2 个参数的情况下）与 b=(t,t)'.

然后我们可以在 MATLAB 中通过 quadprog 求解：

计算矩阵H,f将lasso函数转为二项式
设delta_i ，i= 1, 2, ... , 2P 为p-tuples 形式（+-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

p.s所有重点都来自套索纸。

希望对您有所帮助！

Answer 2

下面的代码可以解决 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)

不过，我还是推荐第二个post中的代码，它使用了QUADPROG函数。它会更快更准确。

Ref

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

On how to solve for the sparse OLS - how to apply `l1` minimization in Matlab (educational purpose)

optimization

matlab

machine-learning

time-series

lasso-regression