function [ theta ] = stomp( y,A,S,ts)
[y_rows,y_columns] = size(y);
if y_rows<y_columns
y = y';%y should be a column vector
end
[M,N] = size(A);%传感矩阵A为M*N矩阵
theta = zeros(N,1);%用来存储恢复的theta(列向量)
Pos_theta = [];%用来迭代过程中存储A被选择的列序号
r_n = y;%初始化残差(residual)为y
for ss=1:S%最多迭代S次
product = A'*r_n;%传感矩阵A各列与残差的内积
sigma = norm(r_n)/sqrt(M);
Js = find(abs(product)>ts*sigma);%选出大于阈值的列
Is = union(Pos_theta,Js);%Pos_theta与Js并集
if length(Pos_theta) == length(Is)
if ss==1
theta_ls = 0;%防止第1次就跳出导致theta_ls无定义
end
break;%如果没有新的列被选中则跳出循环
end
if length(Is)<=M
Pos_theta = Is;%更新列序号集合
At = A(:,Pos_theta);%将A的这几列组成矩阵At
else
if ss==1
theta_ls=0;
end
break;
end
theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解
本内容试读结束,登录后可阅读更多
下载后可阅读完整内容,剩余1页未读,立即下载