Some first-orderalgorithms for structured optimizations[1]

最新推荐文章于 2024-11-12 12:29:51 发布

luds97

最新推荐文章于 2024-11-12 12:29:51 发布

阅读量118

点赞数

文章标签：算法 matlab

本文链接：https://blog.csdn.net/luds97/article/details/127014373

版权

本文探讨了压缩感知问题中的Proximal Point Algorithm (PPA)、Relaxed PPA、PC Method I/II及自适应投影收缩(SA-PC)方法，通过实例代码展示了如何在Matlab中实现。作者还对比了这些方法与GPSR-BB、FPC、AS和SpaRSA等先进算法，并提供了详细参数设置。

摘要由CSDN通过智能技术生成

1.Chapter3中的几个方法

1. Proximal Point Algrorithm(PPA) 邻近点算法
经典PPA算法

2. Relaxed PPA-based method 松弛的邻近点算法
松弛的PPA算法
3. PC Method I (Projection and Contraction) 投影收缩法I
投影收缩法1
4. PC method II 投影收缩法II
投影收缩法2
5. Self-Adaptive Projection and Contraction(SA-PC)method 自适应投影收缩法
SAPC Method

2. 要解决的问题：压缩感知问题

在这里插入图片描述

3.代码

1.CS problem

%问题中的矩阵A
A=rand(m,n)*2-1;
for i=1:m
	t=norm(A(i,:));
	A(i,:)=A(i,:)/t;
end

%问题中的向量b
spar=0.5；稀疏度设置为0.5，可以自行调整
xop=sprand(n,1,spar);%函数sprand创建一个稀疏的正态分布矩阵，spar代表稀疏度
[ii,jj,kk]=find(xop);
for i=1:size(ii)
	xop(ii(i))=sign(xop(ii(i))*2-1);
end
b0=A*xop;
b=b0.*(ones(m,1)+0.01*randn(m,1));

2.注意每个算法中的参数设置在这里插入图片描述

3.算法对应的代码

%问题的PC1算法
%Matlab Code for PC Method I
%根据要求r设置为 r=(m/n)*max(eig(A*A'))
function [x,k]=PC1(A,n,m,b,maxit,tau,eps)
% r=max(eig(A*A'))/4;
r=max(eig(A*A'))*m/n;
stopc=1;
x=zeros(n,1);
k=0;
Ax=A*x;
gamma=1.9;
while(stopc>eps && k<maxit)
    k=k+1;
    x0=x;
    Ax0=Ax;
    d=((b-Ax0)'*A)'/r +x0;
    xt=d-max(min(d,tau/r), -tau/r);
    ex=x0-xt;
    Axt=A*xt;
    Aex=Ax0-Axt;
    T1=ex'*ex;
    T2=(Aex'*Aex)/r;
    alpha = T1*gamma/(T1+T2);
    x=x0-ex*alpha;
    Ax=Ax0 - alpha*Aex;
    stopc=norm(x-x0,inf);
end
x
k
end

%问题的松弛的PPA算法
%Matlab Code for Relaxed PPA based Method
%Relax PPA-based method with r=1.02*lamda_max(A*A')
function [x,k]=PPA(A,n,m,b,maxit,tau,eps)
% r=max(eig(A*A'))/4;
r=max(eig(A*A'))*1.02;
stopc=1;
x=zeros(n,1);
k=0;
Ax=A*x;
gamma=1.5;
while(stopc>eps && k<maxit)
    k=k+1;
    x0=x;
    Ax0=Ax;
    d=((b-Ax0)'*A)'/r +x0;
    xt=d-max(min(d,tau/r), -tau/r);
    ex=x0-xt;
    Axt=A*xt;
    Aex=Ax0-Axt;
    x=x0-ex*gamma;
    Ax=Ax0 - gamma*Aex;
    stopc=norm(x-x0,inf);
end
x
k
end

%问题的SA_PC算法
% The matlab code of the SA-PC Method of (3.54)
function [x,k]=SA_PC(A,n,m,b,maxit,tau,eps)
r=1; stopc=1;x=zeros(n,1);
k=0; Ax=A*x; l=1;
while(stopc>eps && k<maxit)
    k=k+1;
    x0=x;
    Ax0=Ax;
    d0=((b-Ax0)'*A)';
    l=l+1;
    d=d0+x0*r;
    x=(d-max(min(d,tau),-tau))/r;
    Ax=A*x;
    l=l+1;
    ex=x0-x;
    Aex=Ax0-Ax;
    %3.60式子
    T1=ex'*ex;
    T2=Aex'*Aex;
    TT=T2/(T1*r);
    while(TT>1.9)
        r=r*TT*1;
        d=d0+x0*r;
        x=(d-max(min(d,tau),-tau))/r;
        Ax=A*x;
        l=l+1;
        ex=x0-x;
        Aex=Ax0-Ax;
        T1=ex'*ex;
        T2=Aex'*Aex;
        TT=T2/(T1*r);
    end
    stopc=norm(x-x0,inf);
    r=T2*0.85/T1;
end
x
k
end

%SA_PC1算法
% The matlab code of the SA-PC Method for problem(3.54) with continuation
function [x,k]=SA_PC1(A,n,m,b,maxit,tau,eps)
r=1; stopc=1;  x=zeros(n,1);
Ax=A*x; k=0; l=1; 
tauf=tau;
rr=exp(log(tau/tauf)/40);
while((stopc>eps && k<maxit )|| tau>tauf)
    k=k+1;
    x0=x;
    Ax0=Ax;
    d0=((b-Ax0)'*A)';
    l=l+1;
    d=d0+x0*r;
    x=(d-max(min(d,tau),-tau))/r;
    Ax=A*x;
    l=l+1;
    ex=x0-x;
    Aex=Ax0-Ax;
    T1=ex'*ex;
    T2=Aex'*Aex;
    TT=T2/(T1*r);
    while(TT>1.9)
            r=r*TT*1;
            d=d0+x0*r;
            x=(d-max(min(d,tau),-tau))/r;
            Ax=A*x;
            l=l+1;
            ex=x0-x;
            Aex=Ax0-Ax;
            T1=ex'*ex;
            T2=Aex'*Aex;
            TT=T2/(T1*r);
     end
        stopc=norm(x-x0,inf);
        r=T2*0.85/T1;
        tau=max(tau/rr,tauf);
end
x
k
end

%main.m
%主程序文件
clc
clear all
%% 设置m,n的大小，设置一下三种
m=1024;
n=4096;
% m=1600;
% n=8192;
% m=2000;
% n=12000; 
%% The matrix A 
A=rand(m,n)*2-1;
for i=1:m
    t=norm(A(i,:));
    A(i,:)=A(i,:)/t;
end
%% The vector b
spar=0.5;
xop=sprand(n,1,spar);
[ii,jj,kk]=find(xop);
for i=1:size(ii)
    xop(ii(i))=sign(xop(ii(i))*2-1);
end
b0=A*xop;
b=b0.*(ones(m,1)+0.01*randn(m,1));
%% paramters
eps=1e-4;
%eps=1e-3;
maxit=1000;
tau=norm(A'*b, inf)*0.1;
%tau=norm(A'*b, inf)*0.1;
%% 调用函数:分别测试
tic
% [x,k]=SA_PC(A,n,m,b,maxit,tau,eps)
% [x,k]=PPA(A,n,m,b,maxit,tau,eps)
[x,k]=SA_PC1(A,n,m,b,maxit,tau,eps)
% [x,k]=PC1(A,n,m,b,maxit,tau,eps);
toc
num=sum(sum(x~=0))%统计其中的非零量