内点法
\quad
单纯形法
(
S
i
m
p
l
e
x
M
e
t
h
o
d
)
(Simplex Method)
(SimplexMethod)可以用来求解带约束的线性规划命题
(
L
P
)
(LP)
(LP),与之类似的有效集法
(
A
c
t
i
v
e
S
e
t
M
e
t
h
o
d
)
(Active Set Method)
(ActiveSetMethod)可以用来求解带约束的二次规划
(
Q
P
)
(QP)
(QP),而内点法
(
I
n
t
e
r
i
o
r
P
o
i
n
t
M
e
t
h
o
d
)
(Interior Point Method)
(InteriorPointMethod)则是另一种用于求解带约束的优化命题的方法。而且无论是面对
L
P
LP
LP还是
Q
P
QP
QP,内点法都显示出了相当的极好的性能,例如多项式的算法复杂度。内点法主要有障碍函数法
(
B
a
r
r
i
e
r
M
e
t
h
o
d
)
(Barrier Method)
(BarrierMethod)和原始对偶法
(
P
r
i
m
a
l
−
D
u
a
l
M
e
t
h
o
d
)
(Primal-Dual Method)
(Primal−DualMethod)。
本文中采用由
D
o
n
o
h
o
Donoho
Donoho等人发表在论文
A
t
o
m
i
c
Atomic
Atomic
D
e
c
o
m
p
o
s
i
t
i
o
n
Decomposition
Decomposition
b
y
by
by
B
a
s
i
s
Basis
Basis
P
u
r
s
u
i
t
Pursuit
Pursuit中提出的原始
−
-
−对偶对数障碍函数法。该论文利用原始
−
-
−对偶对数障碍函数法和共轭梯度法成功求解了大规模的优化问题。
原始 − - −对偶对数障碍函数法
\quad
本文中主要利用原始
−
-
−对偶对数障碍函数法来求解下面的一个带扰动的线性规划问题,也可以理解为二次规划问题。
min
x
c
T
x
+
1
2
∥
γ
x
∥
2
+
1
2
∥
p
∥
2
s
.
t
.
A
x
+
σ
p
=
b
x
≥
0
{\rm{ }}\mathop {\min }\limits_x {c^T}x + \frac{1}{2}{\left\| {\gamma x} \right\|^2}{\rm{ + }}\frac{1}{2}{\left\| p \right\|^2} \quad {\rm{ s}}{\rm{.t}}{\rm{. }}Ax + \sigma p = b{\rm{ }}x \ge 0
xmincTx+21∥γx∥2+21∥p∥2s.t.Ax+σp=bx≥0
p
p
p是一个高斯白噪声,
σ
\sigma
σ是扰动的水平,
γ
\gamma
γ是
x
x
x的扰动。论文中
γ
\gamma
γ代表扰动,也可以理解为范数。下面是原始
−
-
−对偶对数障碍函数法算法:
如果对该算法感兴趣,可以浏览最后两篇参考文献,特别是第三篇参考文献给出了更一般的二次规划算法。
信号分解中的应用
上述参数设置不同的值,会有下面两个不同的模型。
(
M
O
F
)
min
∥
x
∥
2
s
.
t
.
y
=
A
x
x
≥
0
(MOF) \quad{\rm{ }}\min {\rm{ }}{\left\| x \right\|_2} \quad {\rm{ s}}{\rm{.t}}{\rm{. }} \quad y = Ax \quad {\rm{ }}x\ge 0
(MOF)min∥x∥2s.t.y=Axx≥0
(
B
P
D
N
)
m
i
n
∥
x
∥
1
s
.
t
.
A
x
=
y
x
≥
0
(BPDN) \quad{\rm{ min}}{\left\| x \right\|_1} \quad {\rm{ s}}{\rm{.t}}{\rm{. }} \quad Ax = y \quad {\rm{ }}x\ge 0
(BPDN)min∥x∥1s.t.Ax=yx≥0
上述两个模型都是为了实现信号的自动分解,
M
O
F
MOF
MOF存在两个问题:
(
1
)
(1)
(1)分解出来的系数不稀疏。在某些情况下,如果信号本身就具有稀疏性,这样分解出来的结果就会非常糟糕。
(
2
)
(2)
(2)
M
O
F
MOF
MOF本身存在分辨率受限的问题。
B
P
D
N
BPDN
BPDN就是把
B
P
BP
BP基追踪应用于去噪,所以命名。相比于
M
O
F
MOF
MOF,
B
P
D
N
BPDN
BPDN把
2
2
2范数变为
1
1
1范数,这就可以保持稀疏性。
但是前者的运算速度会远远快于后者。尽管利用下面的程序,可能会得出完全相反的结论。这主要是利用原始
−
-
−对偶对数障碍函数法算法去实现
M
O
F
MOF
MOF并不是明智的选择。
压缩感知中的应用
\quad 在这里,我么将 B P D N BPDN BPDN模型,也就是 B P BP BP应用于压缩感知 C S CS CS恢复重建过程。我们使用高斯随机矩阵作为测量矩阵 A A A,稀疏度 K = 10 K=10 K=10,原始信号是一个只有 K K K个为 1 1 1,其余为 0 0 0的信号。那么只要 σ = 0 \sigma=0 σ=0, γ = 0 \gamma=0 γ=0, c c c是一个全 1 1 1的列向量即可。
算法程序
function [x,z] = PDLB_LP(c,A,b,delta,gama)
%UNTITLED5 此处显示有关此函数的摘要
% 此处显示详细说明
%A Primal-Dual Log-Barrier LP Algorithm
%set parameters
FeaTol=1e-4;%the feasibility tolerance
PDGapTol=1e-4;%the duality gap tolerance
% delta=0.01;%delta=1,即为BPDN,delta也可以为1e-4
% gama=0.01;%gama,原始变量的扰动
%Initialize
[m,n]=size(A);
x=ones(n,1);%生成全1向量
y=zeros(m,1);%生成全0向量
z=0.5*ones(n,1);
e=ones(n,1);
miu=1;
%Loop
t=c+(gama^2)*x-z-A'*y;
r=b-A*x-(delta^2)*y;
X=diag(x,0);%生成对角线元素为x的对角矩阵
Z=diag(z,0);
v=miu*e-Z*x;
%step a
D=inv(X\Z+(gama^2)*eye(n));%A\b=inv(A)*b
PInfeas=norm(r,2)/(1+norm(x,2));%Primal infeasibility
DInfeas=norm(t,2)/(1+norm(y,2));%Dual infeasibility
DuGap=(z'*x)/(1+norm(x,2)*norm(z,2));
while(PInfeas>FeaTol || DInfeas>FeaTol || DuGap>PDGapTol)
%step b:solve
if (delta>0)
A1=[(D^(1/2))*A';delta*eye(m)];%(n+m,m)
B=[(D^(1/2))*(t-X\v);r/delta];
y1=pinv(A1)*B;%广义逆,最小二次法求delta_y,即y1
else
y1=(A*D*A')\(r+A*D*(t-X\v));
end
x1=D*(A'*y1+X\v-t);
z1=X\(v-Z*x1);
%step c:Calculate the primal and dual step sizes ρp,ρd and update the variables:
index1=find(x1<0);
index2=find(z1<0);
if isempty(index1)==0 && isempty(index2)==0%不是空
x2=x1(index1);%储存负数
z2=z1(index2);
x3=x(index1);%储存x1中负数对应的x值
z3=z(index2);
p1=abs(x3./x2);
p2=abs(z3./z2);
elseif isempty(index1)==0 && isempty(index2)==1
x2=x1(index1);%储存负数
x3=x(index1);%储存x1中负数对应的x值
p1=abs(x3./x2);
p2=p1;
elseif isempty(index1)==1 && isempty(index2)==0
z2=x1(index2);%储存负数
z3=z(index2);%储存x1中负数对应的x值
p2=abs(z3./z2);
p1=p2;
else
p1=min(x./x1);
p2=min(z./z1);
end
rou_p=0.99*min(p1);
rou_d=0.99*min(p2);
x=x+rou_p*x1;
y=y+rou_d*y1;
z=z+rou_d*z1;
miu=(1-min([rou_p,rou_d,0.99]))*miu;
%step a
t=c+(gama^2)*x-z-A'*y;
r=b-A*x-(delta^2)*y;
X=diag(x,0);%生成对角线元素为x的对角矩阵
Z=diag(z,0);
v=miu*e-Z*x;
D=inv(X\Z+(gama^2)*eye(n));%A\b=inv(A)*b
PInfeas=norm(r,2)/(1+norm(x,2));%Primal infeasibility
DInfeas=norm(t,2)/(1+norm(y,2));%Dual infeasibility
DuGap=(z'*x)/(1+norm(x,2)*norm(z,2));
end
end
实验结果
https://blog.csdn.net/dymodi/article/details/46441783
Chen S S , Saunders D M A . Atomic Decomposition by Basis Pursuit[J]. Siam Review, 2001, 43(1):129-159.
P. E. Gill, W. Murray, D. B. Ponceleón, and M. A. Saunders, Solving Reduced KKT Systems in Barrier Methods for Linear and Quadratic Programming, Report SOL 91-7, Stanford University, Stanford, CA, July 1991.