用于求解线性规划的原始-对偶对数障碍函数法在信号分解和压缩感知中的应用

内点法

\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) PrimalDualMethod
本文中采用由 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γx2+21p2s.t.Ax+σp=bx0
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)minx2s.t.y=Axx0

( 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)minx1s.t.Ax=yx0
上述两个模型都是为了实现信号的自动分解, 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.

©️2020 CSDN 皮肤主题: 游动-白 设计师:上身试试 返回首页