问题描述
考虑线性测量 b=Ax+e，其中 b 为 50 维的测量值，A 为 50ⅹ100 维的测量矩阵，x 为 100维的未知稀疏向量且稀疏度为 5，e 为 50 维的测量噪声。从 b 与 A 中恢复 x 的一范数规范化最小二乘模型如下：
m
i
n
1
2
∣
∣
A
x
−
b
∣
∣
2
2
+
p
∣
∣
x
∣
∣
1
min \ \frac{1}{2}||Ax-b||_2^2 + p||x||_1
min 21∣∣Ax−b∣∣22+p∣∣x∣∣1
其中 p 为非负的正则化参数。
在实验中的数据要求：
- 设 x 的真值中的非零元素服从均值为 0 方差为 1 的高斯分布
- A 中的元素服从均值为 0 方差为 1 的高斯分布
- e 中的元素服从均值为 0 方差为 0.1 的高斯分布
请设计下述算法求解该问题：
- 邻近点梯度
- 交替方向乘子法
- 次梯度法
对于每种算法，请给出每步计算结果与真值的距离以及每步计算结果与最优解的距离。此外，请讨论正则化参数 p 对计算结果的影响。
解决思路
- 通过随机生成数据，生成符合题目要求的A，x和e，再根据 b = A x + e b=Ax+e b=Ax+e 求出b。并保存A，x和b，以供后续算法使用的是同一个数据。
- 设计三种算法（邻近点梯度法、交替方向乘子法、次梯度法）解决该问题。
- 记录每一步优化得到的 x k x^k xk ，便于绘制每步计算结果与真值的距离、每步计算结果与最优解的距离。
- 讨论正则化参数 p 对计算结果的影响
- 分析实验结果，对比三种算法的各项性能
生成数据
按照题目要求生成数据并保存，以便每个算法用的数据都是相同的。
要求： X 为100*1的矩阵，稀疏度为5，X 的真值中的非零元素服从均值为 0 方差为 1 的高斯分布。A 为50ⅹ100 维的测量矩阵，A 中的元素服从均值为 0 方差为 1 的高斯分布。e 为 50 维的测量噪声，e 中的元素服从均值为 0 方差为 0.1 的高斯分布。b通过 b=Ax+e 算出。
%% 生成数据
p = 0.01;
xishudu = 5;
A = random('Normal', 0, 1, 50, 100); %生成均值为0，方差为1的50*100矩阵
e = random('Normal', 0, 0.1, 50, 1); %生成均值为0，方差为0.1的50*1矩阵
X = zeros(100,1); %X = sparse(100,1)
index = fix( 100*rand(1,xishudu) );
for i = 1:xishudu
X(index(i)) = randn();
end
b = A*X + e;
save('A.mat','A')
save('b.mat','b')
save('X.mat','X')
算法设计
一、邻近点梯度法
根据邻近点梯度法，若函数有结构非光滑，设
m
i
n
f
0
(
x
)
=
s
(
x
)
+
r
(
x
)
min \ f_0(x)=s(x)+r(x)
min f0(x)=s(x)+r(x)
其中s(x)是光滑部分，可微易求导，在本题中s(x)对应的是 $ \frac{1}{2}||Ax-b||_2^2 $
r(x)是非光滑部分，不可微，但可求得邻近点投影
m
i
n
r
(
x
)
+
1
2
∗
α
∗
∣
∣
x
−
x
^
∣
∣
min \ r(x)+\frac{1}{2*α}*||x-\hat{x}||
min r(x)+2∗α1∗∣∣x−x^∣∣
因此有：
x
k
+
1
2
=
x
k
−
α
⋅
▽
s
(
x
k
)
x^{k+ \frac{1}{2}}=x^k-α\cdot▽s(x^k)
xk+21=xk−α⋅▽s(xk)
x k + 1 = a r g m i n x r ( x ) + 1 2 ∗ α ∣ ∣ x − x k + 1 2 ∣ ∣ 2 x^{k+1}=argmin_x \ r(x)+\frac{1}{2*α}||x-x^{k+\frac{1}{2}}||^2 xk+1=argminx r(x)+2∗α1∣∣x−xk+21∣∣2
在光滑的部分直接求导，非光滑部分求次梯度。若能找到次梯度为0的点，则该点为最优解。
将原问题代入方程，可得到：
x
k
+
1
2
=
x
k
−
2
∗
α
∗
A
T
(
A
x
k
−
b
)
x^{k+ \frac{1}{2}}=x^k-2*α*A^T(Ax^k-b)
xk+21=xk−2∗α∗AT(Axk−b)
x k + 1 = a r g m i n x p ∣ ∣ x ∣ ∣ 1 + 1 2 ∗ α ∣ ∣ x − x k + 1 2 ∣ ∣ 2 x^{k+1}=argmin_x \ p||x||_1+\frac{1}{2*α}||x-x^{k+\frac{1}{2}}||^2 xk+1=argminx p∣∣x∣∣1+2∗α1∣∣x−xk+21∣∣2
而求解argmin的方法可以用软门限法。
因为在计算中无法准确找到次梯度为0的点，我们近似地两步之间的二范数< 1 0 − 5 10^{-5} 10−5的点。
代码
%% 邻近点梯度法
load('A.mat');
load('b.mat');
load('X.mat');
p = 0.01;
alpha = 0.001;
Xsize = size(X);
Xk = zeros(Xsize);
X_dis2real = [];
X_opt_step = zeros(100,30000);
count = 1;
while 1
Xk_2 = Xk - 2*alpha*A'*(A*Xk-b);
Xk_ud = zeros(Xsize);
for i = 1: Xsize
if Xk_2(i) - alpha*p > 0
Xk_ud(i) = Xk_2(i) - alpha*p;
elseif Xk_2(i) + alpha*p < 0
Xk_ud(i) = Xk_2(i) + alpha*p;
end
end
X_dis2real = [X_dis2real norm(Xk_ud - X,2)];
count = count + 1;
X_opt_step(:,count) = Xk_ud;
if norm(Xk_ud - Xk, 2) < 1e-5
Xk = Xk_ud;
break
else
Xk = Xk_ud;
end
end
X_opt = X_opt_step( :, count );
for i = 1:count
X_dis2opt(i) = norm(X_opt_step(:,i) - X_opt, 2);
end
hold on
plot(X_dis2real,'LineWidth',1);
plot(X_dis2opt,'LineWidth',1);
title('Distance');
legend('distance to real', 'distance to opt');
print('Adjacent point gradient','-dpng');
二、交替方向乘子法
根据交替方向乘子法有
m
i
n
f
1
(
x
)
+
f
2
(
y
)
s
.
t
.
A
x
+
B
y
=
0
min \ \ f_1(x)+f_2(y) \\ s.t. \ Ax+By=0
min f1(x)+f2(y)s.t. Ax+By=0
其增广拉格朗日函数为：
L
c
(
x
,
y
,
λ
)
=
f
1
(
x
)
+
f
2
(
y
)
+
<
λ
,
A
x
+
B
y
>
+
c
2
∣
∣
A
x
+
B
y
∣
∣
2
2
L_c(x,y,λ)=f_1(x)+f_2(y)+<λ,Ax+By>+\frac{c}{2}||Ax+By||^2_2
Lc(x,y,λ)=f1(x)+f2(y)+<λ,Ax+By>+2c∣∣Ax+By∣∣22
其中c是一个大于0的常数。
更新
x
k
+
1
x^{k+1}
xk+1 和
y
k
+
1
y^{k+1}
yk+1：
(
x
k
+
1
,
y
k
+
1
)
=
a
r
g
m
i
n
x
,
y
L
c
(
x
,
y
,
λ
k
)
λ
k
+
1
=
λ
k
+
c
(
A
x
k
+
1
+
B
y
k
+
1
)
(x^{k+1},y^{k+1})=argmin_{x,y} \ L_c(x,y,λ^k) \\ λ^{k+1}=λ^k+c(Ax^{k+1}+By^{k+1})
(xk+1,yk+1)=argminx,y Lc(x,y,λk)λk+1=λk+c(Axk+1+Byk+1)
固定
y
k
y^k
yk 和
λ
k
λ^k
λk ，则有
x
k
+
1
=
a
r
g
m
i
n
x
L
c
(
x
,
y
k
,
λ
k
)
y
k
+
1
=
a
r
g
m
i
n
y
L
c
(
x
k
+
1
,
y
,
λ
k
)
λ
k
+
1
=
λ
k
+
c
(
A
x
k
+
1
+
B
y
k
+
1
)
x^{k+1}=argmin_{x} \ L_c(x,y^k,λ^k) \\ y^{k+1}=argmin_{y} \ L_c(x^{k+1},y,λ^k) \\ λ^{k+1}=λ^k+c(Ax^{k+1}+By^{k+1})
xk+1=argminx Lc(x,yk,λk)yk+1=argminy Lc(xk+1,y,λk)λk+1=λk+c(Axk+1+Byk+1)
对于原问题，为了使用交替方向乘子法，引入一个新变量 y = x y=x y=x，所以一致性约束为 x − y = 0 x-y=0 x−y=0。
原问题转换为：
m
i
n
1
2
∣
∣
A
x
−
b
∣
∣
2
2
+
p
∣
∣
y
∣
∣
1
s
.
t
.
x
−
y
=
0
min \ \frac{1}{2}||Ax-b||_2^2 + p||y||_1 \\ s.t. x-y=0
min 21∣∣Ax−b∣∣22+p∣∣y∣∣1s.t.x−y=0
对应的，
f
1
(
x
)
=
1
2
∣
∣
A
x
−
b
∣
∣
2
2
f_1(x)=\frac{1}{2}||Ax-b||_2^2
f1(x)=21∣∣Ax−b∣∣22 ，
f
2
(
y
)
=
p
∣
∣
y
∣
∣
1
f_2(y)=p||y||_1
f2(y)=p∣∣y∣∣1
将原问题代入该算法的方程，则有
L
c
(
x
,
y
,
λ
)
=
1
2
∣
∣
A
x
−
b
∣
∣
2
2
+
p
∣
∣
y
∣
∣
1
+
<
λ
,
x
−
y
>
+
c
2
∣
∣
x
−
y
∣
∣
2
2
L_c(x,y,λ)=\frac{1}{2}||Ax-b||_2^2+p||y||_1+<λ,x-y>+\frac{c}{2}||x-y||^2_2
Lc(x,y,λ)=21∣∣Ax−b∣∣22+p∣∣y∣∣1+<λ,x−y>+2c∣∣x−y∣∣22
x k + 1 = a r g m i n x 1 2 ∣ ∣ A x − b ∣ ∣ 2 2 + p ∣ ∣ y k ∣ ∣ 1 + < λ k , x − y k > + c 2 ∣ ∣ x − y k ∣ ∣ 2 2 x^{k+1}=argmin_{x} \ \frac{1}{2}||Ax-b||_2^2+p||y^k||_1+<λ^k,x-y^k>+\frac{c}{2}||x-y^k||^2_2 xk+1=argminx 21∣∣Ax−b∣∣22+p∣∣yk∣∣1+<λk,x−yk>+2c∣∣x−yk∣∣22
因为x部分是光滑的，可求导。所以可以得到x的梯度
x
k
+
1
=
（
A
T
A
+
c
I
）
−
1
(
A
T
b
+
c
y
k
−
λ
k
)
x^{k+1}=（A^TA+cI）^{-1}(A^Tb+cy^k-λ^k)
xk+1=（ATA+cI）−1(ATb+cyk−λk)
其中
I
I
I 是单位矩阵。
y部分可以求其邻近点投影，这里可以用软门限法来求解。
y
k
+
1
=
a
r
g
m
i
n
y
p
∣
∣
y
∣
∣
1
+
c
2
∣
∣
y
−
x
k
+
1
−
λ
k
c
∣
∣
2
2
λ
k
+
1
=
λ
k
+
c
(
x
k
+
1
−
y
k
+
1
)
y^{k+1}=argmin_{y} \ p||y||_1+\frac{c}{2}||y-x^{k+1}-\frac{λ^k}{c}||^2_2 \\ λ^{k+1}=λ^k+c(x^{k+1}-y^{k+1})
yk+1=argminy p∣∣y∣∣1+2c∣∣y−xk+1−cλk∣∣22λk+1=λk+c(xk+1−yk+1)
代码
初始版本的代码由于每一步都要求矩阵的逆，从探查器看出这一步占用了非常多的时间。所以优化了代码，避免重复运算。另外由于有变量未预分配空间，也导致了它占用了较多的时间，修改过后的版本为变量预分配了空间。
最终版代码如下：
%% 交替方向乘子法(ADMM)
load('A.mat');
load('b.mat');
load('X.mat');
p = 0.01;
c = 0.005;
Xsize = size(X);
Xk = zeros(Xsize);
Yk = zeros(Xsize);
Lk = zeros(Xsize);
I = eye(Xsize(1));
X_dis2real = zeros(1,1000);
X_opt_step = [];
count = 0;
tidai = (A'*A+c*I)^(-1);
A_b = A'*b;
while 1
Xk_1 = tidai * ( A_b + c*Yk - Lk);
Yk_1 = zeros(Xsize);
for i = 1:Xsize
if Xk_1(i) + Lk(i)/c - p/c > 0
Yk_1(i) = Xk_1(i) + Lk(i)/c - p/c;
elseif Xk_1(i) + Lk(i)/c + p/c < 0
Yk_1(i) = Xk_1(i) + Lk(i)/c + p/c;
end
end
Lk_1 = Lk + c*(Xk_1 - Yk_1);
count = count + 1 ;
X_dis2real(count) = norm(Xk_1 - X, 2);
X_opt_step = [X_opt_step norm(Xk_1 - Xk, 2)];
if norm(Xk_1 - Xk, 2) < 1e-5
break
else
Xk = Xk_1;
Yk = Yk_1;
Lk = Lk_1;
end
end
X_opt = X_opt_step( :, count );
for i = 1:count
X_dis2opt(i) = norm(X_opt_step(:,i) - X_opt, 2);
end
hold on
plot(X_dis2real,'LineWidth',1);
plot(X_dis2opt,'LineWidth',1);
title('Distance');
legend('distance to real', 'distance to opt');
print('ADMM','-dpng');
三、次梯度法
根据此梯度法，有
x
k
+
1
=
x
k
−
α
k
∗
g
0
(
x
k
)
x^{k+1}=x^k-α^k*g_0(x^k)
xk+1=xk−αk∗g0(xk)
其中
g
0
(
x
)
g_0(x)
g0(x) 是
f
0
(
x
)
f_0(x)
f0(x) 的次梯度，满足
g
0
(
x
)
∈
∂
f
0
(
x
)
g_0(x)\in \partial{f_0(x)}
g0(x)∈∂f0(x)。
次梯度可以由软门限算法求得。
次梯度法的步长选取有几种方法，如固定步长，不可加但平方可加的变长步长，不可加且平方不可加的变长步长。这里我选择的是变长步长 α k = 0.001 k α^k=\frac{0.001}{k} αk=k0.001 。
代码
%% 次梯度法（subgradient）
load('A.mat');
load('b.mat');
load('X.mat');
p = 0.01;
alpha = 0.001;
count = 0;
Xsize = size(X);
Xk = zeros(Xsize);
alphak = alpha;
g = @(x) A'*(A*x-b) + p * ran_x(x);
X_dis2real = [];
X_dis2opt = [];
X_opt_step = [];
while 1
Xk_1 = Xk - alphak * g(Xk);
alphak = alpha / (count+1);
count = count + 1;
X_dis2real = [X_dis2real norm(Xk_1-X)];
X_opt_step = [X_opt_step Xk_1];
if norm(Xk_1 - Xk,2) < 1e-5
break
else
Xk = Xk_1;
end
end
X_opt = X_opt_step( :,count);
for i = 1:count
X_dis2opt(i) = norm(X_opt_step(:,i) - X_opt, 2);
end
hold on
plot(X_dis2real,'LineWidth',1);
plot(X_dis2opt,'LineWidth',1);
title('Distance');
legend('distance to real', 'distance to opt');
print('Subgradient','-dpng');
使用到的函数ran_x代码如下，该函数是求一范数规范化的次梯度的函数 。
function [x_new] = ran_x(x)
x_new = zeros(size(x));
for i = 1:size(x)
if x(i) == 0
x_new(i) = 2*random('Normal', 0, 1) - 1;
else
x_new(i) = sign(x(i));
end
end
end
扫一扫
1543
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
