基追问题源码matlab-CSDN博客

本文链接：https://blog.csdn.net/m0_58786865/article/details/135466669

基追问题源码matlab

命令行
ADLPMM
GEM
PGA_a1
PGA_b1

命令行

%2024/1/8 by xiaoyuhang
clc;clear;
% 生成稀疏矩阵A，真实向量x_ture，起始点x0
m = 1000;n = 1100;       
A = randn(m,n);
x_ture = zeros(n,1);
x_ture(3:8:80)=1;
x_ture(7:8:80)=-1;
b = A*x_ture;
x0=ones(n,1);
eps=1e-6;
iter = 2000;
lambda = zeros(m,1);
[~,D]=eig(A'*A);
beta_min=1/max(diag(D));
%% 基追问题(adlpmm,gem,pga)
% 定义函数
E = @(x) sum(abs(x-x_ture)); % 计算误差
F =@(x) norm(x,1);
%% ADLPMM
rho=0.05;
beta=0.05;
[x_adlpmm,k_adlpmm,t_adlpmm]= ADLPMM_bp(A,b,x0,lambda,iter,eps,rho,beta);
f_adlpmm=[];
for i=1:k_adlpmm
    ff_adlpmm=F(x_adlpmm(:,i));
    f_adlpmm=[f_adlpmm,ff_adlpmm];
end
%% GEM
nu=0.7;
mu=0.3;
[x_gem,k_gem,t_gem]= GEM_bp(A,b,x0,nu,mu,beta,lambda,iter,eps);
conv_gem = E(x_gem);
f_gem=[];
for i=1:k_gem
    ff_gem=F(x_gem(:,i));
    f_gem=[f_gem,ff_gem];
end
%% PGAa1
gama=1.3;
[x_pgaa1,k_pgaa1,t_pgaa1]= PGAa1_bp(A,b,x0,nu,mu,beta,gama,lambda,iter,eps);
conv_pgaa1 = E(x_pgaa1);
f_pgaa1=[];
for i=1:k_pgaa1
    ff_pgaa1=F(x_pgaa1(:,i));
    f_pgaa1=[f_pgaa1,ff_pgaa1];
end
%% PGAb1
[x_pgab1,k_pgab1,t_pgab1]= PGAb1_bp(A,b,x0,nu,mu,beta,gama,lambda,iter,eps);
conv_pgab1 = E(x_pgab1);
f_pgab1=[];
for i=1:k_pgab1
    ff_pgab1=F(x_pgab1(:,i));
    f_pgab1=[f_pgab1,ff_pgab1];
end
%% 画图
figure;
semilogy(f_adlpmm,'--','LineWidth',1,...
    'DisplayName','ADLPMM','Color','black')
hold on
semilogy(f_gem,'-.','LineWidth',1,...
    'DisplayName','GEM','Color','red')
semilogy(f_pgaa1,'--','LineWidth',1,...
    'DisplayName','PGA_{a1}','Color','magenta')
semilogy(f_pgab1,'_.','LineWidth',1,...
    'DisplayName','PGA_{b1}','Color','green')
hold off
xlabel('迭代次数')
ylabel('$f(\widetilde{x}^k)$','Interpreter','latex')
legend('Location','northeast')
title('ADLPMM,GEM,PGA_{a1},PGA_{b1}在基追问题上f值的变化趋势')

ADLPMM

function [x_s,k,timecost] = ADLPMM_bp(A,b,x0,lambda,iter, eps,rho,beta)
% 2024/1/8 by xiaoyuhang
% 用ADMM解决基追问题:
% min  ||x||_1
% s.t. Ax=b
% 输入
%=======================================
% A b x0
% x0......... 起始点
% rho........ 惩罚参数
% eta ....... 常数,确定利普西茨常数L0
% lambda .... 拉格朗日乘子
% eps ....... 停止准则
% iter ...... 迭代次数
%=======================================
% 输出
% x_s ....... adlpmm生成的序列 {xk} 
% k ......... 迭代终止时迭代次数
% timecost .. 迭代终止时花费时间
t0=clock;
x_s = [];
x_old = x0;
lambda_old=lambda;
% 主循环
T = @(tau, x) max(0, x - tau) + min(0, x + tau);
for i = 1:iter
    x_new =  T(beta, x_old - rho*beta*A'*(A*x_old-b+lambda_old/rho));
    lambda_new =lambda_old+rho*(A*x_new-b);
    % 记录
    x_s = [x_s, x_new];
    % 检查停止
    f_x_lambda=T(1,[x_new;lambda_new]-([zeros(1100,1100) -A';A zeros(1000,1000)]*[x_new;lambda_new]+[zeros(1100,1);-b]));
    e = norm([x_new;lambda_new]-f_x_lambda,inf);
    k=i;
    if e < eps|| i==2000
         timecost=etime(clock,t0);
         break
    end
    % 更新
    x_old = x_new;
    lambda_old=lambda_new;
end
end

GEM

function [x_s,k, timecost] = GEM_bp(A,b,x0,nu,mu,beta,lambda,iter, eps)
% 2024/1/8 by xiaoyuhang
% 用GEM解决基追问题:
% min  ||x||_1
% s.t. Ax=b
% min F（x）=1/2*|| Ax - b ||_2^2 + \lambda || x ||_1
%=======================================
% 输入 A b
% x0......... 起始点
% nu......... 常数,确定beta
% eta ....... 常数,确定利普西茨常数L0
% lambda .... l_1范数的正则化参数
% eps ....... 停止准则
% iter ...... 迭代次数
%=======================================
% 输出
% x_s ....... GEM生成的序列 {xk} 
% k ......... 迭代终止时迭代次数
% timecost .. 迭代终止时花费时间
t0=clock;
x_s = [];
x_old = x0;
lambda_old=lambda;
j=0;
% 主循环
T = @(tau, x) max(0, x - tau) + min(0, x + tau);
for i = 1:iter
    x_bar=T(beta, x_old + beta*A'*lambda_old);
    lambda_bar=lambda_old-beta*(A*x_old-b);
    r=beta*norm([A'*(lambda_old-lambda_bar);A*(x_old-x_bar)])/norm([x_old-x_bar;lambda_old-lambda_bar]);
    if r>nu
        beta=(2/3)*beta*min(1,1/r);
        j=j+1;
        continue
    end
    x_new =  T(beta, x_old + beta*A'*lambda_bar);
    lambda_new=lambda_old-beta*(A*x_bar-b);
    if r<=mu
        beta=1.5*beta;
    end
    % 记录
    x_s = [x_s, x_new];
    % 检查停止
    e = norm(x_new-x_old,inf);
    k=i-j;
    if e < eps
        timecost=etime(clock,t0);
        break
    end
    % 更新
    x_old = x_new;
    lambda_old=lambda_new;
end
end

PGA_a1

function [x_s,k,timecost] = PGAa1_bp(A,b,x0,nu,mu,beta,gama,lambda,iter, eps)
% 2024/1/8 by xiaoyuhang
% 用PGAa1 解决基追问题
% min  ||x||_1
% s.t. Ax=b
%=======================================
% 输入 A b
% x0......... 起始点
% nu......... 常数,确定beta
% eta ....... 常数,确定利普西茨常数L0
% lambda .... l_1范数的正则化参数
% eps ....... 停止准则
% iter ...... 迭代次数
%=======================================
% 输出
% x_s ....... pgaa1生成的序列 {xk} 
% k ......... 迭代终止时迭代次数
% timecost .. 迭代终止时花费时间
t0=clock;
x_s = [];
x_old=x0;
lambda_old=lambda;
j=0;
% 主循环
T = @(tau, x) max(0, x - tau) + min(0, x + tau);
for i = 1:iter
    x_bar=T(beta, x_old + beta*A'*lambda_old);
    lambda_bar=lambda_old-beta*(A*x_old-b);
    r=beta*norm([A'*(lambda_old-lambda_bar);A*(x_old-x_bar)])/norm([x_old-x_bar;lambda_old-lambda_bar]);
    if r>nu
        beta=(2/3)*beta*min(1,1/r);
        j=j+1;
        continue
    end
    alpha=norm([x_old-x_bar;lambda_old-lambda_bar])^2/norm([x_old-x_bar+beta*A'*(lambda_old-lambda_bar);lambda_old-lambda_bar-beta*A*(x_old-x_bar)])^2;
    x_new=x_old -gama*alpha*(x_old-x_bar+beta*A'*(lambda_old-lambda_bar));
    lambda_new=lambda_old-gama*alpha*(lambda_old-lambda_bar-beta*A*(x_old-x_bar));
    if r<=mu
        beta=1.5*beta;
    end
    % 记录
    x_s = [x_s, x_new];
    % 检查停止
    e = norm([x_new;lambda_new]-[x_old;lambda_old],inf);
    k=i-j;
    if e < eps || i==2000
        timecost=etime(clock,t0);
        break
    end

    % 更新
    x_old = x_new;
    lambda_old=lambda_new;
end
end

PGA_b1

function [x_s,k,timecost] = PGAb1_bp(A,b,x0,nu,mu,beta,gama,lambda,iter, eps)
% 2024/1/8 by xiaoyuhang
% 用PGAb1 解决基追问题
% min  ||x||_1
% s.t. Ax=b
%=======================================
% 输入 A b
% x0......... 起始点
% nu......... 常数,确定beta
% eta ....... 常数,确定利普西茨常数L0
% lambda .... l_1范数的正则化参数
% eps ....... 停止准则
% iter ...... 迭代次数
%=======================================
% 输出
% x_s ....... pgab1生成的序列 {xk} 
% k ......... 迭代终止时迭代次数
% timecost .. 迭代终止时花费时间
t0=clock;
x_s = [];
x_old=x0;
lambda_old=lambda;
j=0;
% 主循环
T = @(tau, x) max(0, x - tau) + min(0, x + tau);
for i = 1:iter
    x_bar=T(beta, x_old + beta*A'*lambda_old);
    lambda_bar=lambda_old-beta*(A*x_old-b);
    r=beta*norm([A'*(lambda_old-lambda_bar);A*(x_old-x_bar)])/norm([x_old-x_bar;lambda_old-lambda_bar]);
    if r>nu
        beta=(2/3)*beta*min(1,1/r);
        j=j+1;
        continue
    end
    alpha=([x_old-x_bar;lambda_old-lambda_bar]'*[x_old-x_bar+beta*A'*(lambda_old-lambda_bar);lambda_old-lambda_bar-beta*A*(x_old-x_bar)])/norm([x_old-x_bar+beta*A'*(lambda_old-lambda_bar);lambda_old-lambda_bar-beta*A*(x_old-x_bar)])^2;
    x_new=x_old -gama*alpha*(x_old-x_bar+beta*A'*(lambda_old-lambda_bar));
    lambda_new=lambda_old-gama*alpha*(lambda_old-lambda_bar-beta*A*(x_old-x_bar));
    if r<=mu
        beta=1.5*beta;
    end
    % 记录
    x_s = [x_s, x_new];
    % 检查停止
    e = norm([x_new;lambda_new]-[x_old;lambda_old],inf);
    k=i-j;
    if e < eps || i==2000
        timecost=etime(clock,t0);
        break
    end

    % 更新
    x_old = x_new;
    lambda_old=lambda_new;
end
end