为什么以下代码输出的重构概率恒为0?
怎么修改啊!
clear all; close all; clc;
%% 参数配置初始化
CNT = 100; % 对于每组(K,M,N),重复迭代次数
N = 256; % 信号x的长度
Psi = eye(N); % x本身是稀疏的,定义稀疏矩阵为单位阵x=Psi*theta
K_set = [4,12,20,28,36]; % 信号x的稀疏度集合
Percentage = zeros(length(K_set),N); % 存储恢复成功概率
beta = 0.5;
max_iter = 100;
%% 主循环,遍历每组(K,M,N)
tic
for kk = 1:length(K_set)
K = K_set(kk); % 本次稀疏度
M_set = K:5:N; % M没必要全部遍历,每隔5测试一个就可以了
PercentageK = zeros(1,length(M_set)); % 存储此稀疏度K下不同M的恢复成功概率
for mm = 1:length(M_set)
M = M_set(mm); % 本次观测值个数
P = 0;
succ = zeros(1,length(M_set)); % 存储每次实验的成功率
for cnt = 1:CNT % 每个观测值个数均运行CNT次
Index_K = randperm(N);
x = zeros(N,1);
x(Index_K(1:K)) = 5*randn(K,1); % x为K稀疏的,且位置是随机的
Phi = randn(N,N); % 测量矩阵为高斯矩阵
A = Phi .* Psi; % 传感矩阵
y = Phi * x; % 得到观测向量y
[theta, ~] = frank_wolfe_modified(A, y, beta, max_iter, x); % 恢复重构信号theta
x_r = Psi * theta; % x=Psi * theta
% 判断成功率
if norm(x_r - x)/norm(x) < 1e-3
succ(mm) = succ(mm)+1;
end
end
PercentageK(mm) = sum(succ)* 100/ CNT ; % 计算恢复概率
end
Percentage(kk,1:length(M_set)) = PercentageK;
end
toc
%% 绘图
S = ['-ks';'-ko';'-kd';'-kv';'-k*'];
figure;
for kk = 1:length(K_set)
K = K_set(kk);
M_set = K:5:N;
L_Mset = length(M_set);
plot(M_set,Percentage(kk,1:L_Mset),S(kk,:)); % 绘出x的恢复信号
hold on;
end
hold off;
xlim([0 256]);
legend('K=4','K=12','K=20','K=28','K=36');
xlabel('Number of measurements(M)');
ylabel('Percentage recovered');
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');
function [succ_rate, err] = frank_wolfe_modified(A, y, beta, maxiter, x)
n = size(A, 2); % 获取系数个数
c = zeros(n, 1); % 初始化系数向量
h = zeros(n, 1);
E = eye(n);
normr = zeros(maxiter, 1); % 保存每次迭代的残差范数
err = zeros(maxiter, 1); % 保存每次迭代的误差
k = 1; % 初始化迭代次数
succ = false; % 初始化成功标志
while k < maxiter && ~succ
r = y - A * c; % 计算残差
h = A' * r;
[h_new, h_index] = sort(abs(h), 'descend');
l = h_index(1);
s = sign(dot(A(:, l), y - A * c)) * beta * E(:, l);
% 计算步长(可以替换为其他稳健的步长计算方式)
gamma = min(1, 2 / (k + 1)); % 例如使用线性搜索
c = c + gamma * (s - c);
k = k + 1;
normr(k) = norm(r, 2);
err(k) = 0.5 * (norm(y - A * c, 2)^2 - norm(y - A * x, 2)^2);
% 判断是否满足终止准则
% 修改停止准则为残差相对变化量小于一定阈值
if k > 1 && abs(normr(k) - normr(k-1)) / normr(k-1) < 1e-6
succ = true;
end
end
% 返回结果
succ_rate = double(succ);
err = err(1:k); % 截取有效迭代次数的误差
end
```