压缩感知重构算法之子空间追踪(SP)

2023-10-07 18:10
文章标签 算法 空间 压缩 重构 sp 追踪 感知 之子

本文主要是介绍压缩感知重构算法之子空间追踪(SP)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

题目：压缩感知重构算法之子空间追踪(SP)

转载自彬彬有礼的专栏

        如果掌握了压缩采样匹配追踪(CoSaMP)后，再去学习子空间追踪(Subspace Pursuit)是一件非常简单的事情，因为它们几乎是完全一样的。

        SP的提出时间比CoSaMP提出时间略晚，首个论文版本是参考文献[1]，后来更新了两次，最后在IEEE Transactions on Information Theory发表[2]。从算法角度来讲，SP与CoSaMP差别非常小，这一点作者也意识到了，在文献[1]首页的左下角就有注释：

在文献[2]第2页提到了SP与CoSaMP的具体不同：

从上面可以知道，SP与CoSaMP主要区别在于“Ineach iteration, in the SP algorithm, only K new candidates are added, while theCoSAMP algorithm adds 2K vectors.”，即SP每次选择K个原子，而CoSaMP则选择2K个原子；这样带来的好处是“This makes the SP algorithm computationally moreefficient,”。

        以下是文献[2]中的给出的SP算法流程：

这个算法流程的初始化(Initialization)其实就是类似于CoSaMP的第1次迭代，注意第(1)步中选择了K个原子：“K indices corresponding to the largest magnitude entries”，在CoSaMP里这里要选择2K个最大的原子，后面的其它流程都一样。这里第(5)步增加了一个停止迭代的条件：当残差经过迭代后却变大了的时候就停止迭代。

        不只是SP作者认识到了自己的算法与CoSaMP的高度相似性，CoSaMP的作者也同样关注到了SP算法，在文献[3]中就提到：

文献[3]是CoSaMP原始提出文献的第2个版本，文献[3]的早期版本[4]是没有提及SP算法的。

        鉴于SP与CoSaMP如此相似，这里不就再单独给出SP的步骤了，参考《压缩感知重构算法之压缩采样匹配追踪(CoSaMP)》，只需将第(2)步中的2K改为K即可。

        引用文献[5]的3.5节中的几句话：“贪婪类算法虽然复杂度低运行速度快，但其重构精度却不如BP类算法，为了寻求复杂度和精度更好地折中，SP算法应运而生”，“SP算法与CoSaMP算法一样其基本思想也是借用回溯的思想，在每步迭代过程中重新估计所有候选者的可信赖性”，“SP算法与CoSaMP算法有着类似的性质与优缺点”。

        子空间追踪代码可实现如下(CS_SP.m)，通过对比可以知道该程序与CoSaMP的代码基本完全一致。本代码未考虑文献[2]中的给出的SP算法流程的第(5)步。代码可参见参考文献[6]中的Demo_CS_SP.m。

[plain] view plain copy
print ?
  1. function [ theta ] = CS_SP( y,A,K )  
  2. %CS_SP Summary of this function goes here  
  3. %Version: 1.0 written by jbb0523 @2015-05-01  
  4. %   Detailed explanation goes here  
  5. %   y = Phi * x  
  6. %   x = Psi * theta  
  7. %   y = Phi*Psi * theta  
  8. %   令 A = Phi*Psi, 则y=A*theta  
  9. %   K is the sparsity level  
  10. %   现在已知y和A，求theta  
  11. %   Reference:Dai W，Milenkovic O．Subspace pursuit for compressive sensing  
  12. %   signal reconstruction[J]．IEEE Transactions on Information Theory，  
  13. %   2009，55(5)：2230-2249.  
  14.     [y_rows,y_columns] = size(y);  
  15.     if y_rows<y_columns  
  16.         y = y';%y should be a column vector  
  17.     end  
  18.     [M,N] = size(A);%传感矩阵A为M*N矩阵  
  19.     theta = zeros(N,1);%用来存储恢复的theta(列向量)  
  20.     Pos_theta = [];%用来迭代过程中存储A被选择的列序号  
  21.     r_n = y;%初始化残差(residual)为y  
  22.     for kk=1:K%最多迭代K次  
  23.         %(1) Identification  
  24.         product = A'*r_n;%传感矩阵A各列与残差的内积  
  25.         [val,pos]=sort(abs(product),'descend');  
  26.         Js = pos(1:K);%选出内积值最大的K列  
  27.         %(2) Support Merger  
  28.         Is = union(Pos_theta,Js);%Pos_theta与Js并集  
  29.         %(3) Estimation  
  30.         %At的行数要大于列数，此为最小二乘的基础(列线性无关)  
  31.         if length(Is)<=M  
  32.             At = A(:,Is);%将A的这几列组成矩阵At  
  33.         else%At的列数大于行数，列必为线性相关的,At'*At将不可逆  
  34.             break;%跳出for循环  
  35.         end  
  36.         %y=At*theta，以下求theta的最小二乘解(Least Square)  
  37.         theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解  
  38.         %(4) Pruning  
  39.         [val,pos]=sort(abs(theta_ls),'descend');  
  40.         %(5) Sample Update  
  41.         Pos_theta = Is(pos(1:K));  
  42.         theta_ls = theta_ls(pos(1:K));  
  43.         %At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空间上的正交投影  
  44.         r_n = y - At(:,pos(1:K))*theta_ls;%更新残差   
  45.         if norm(r_n)<1e-6%Repeat the steps until r=0  
  46.             break;%跳出for循环  
  47.         end  
  48.     end  
  49.     theta(Pos_theta)=theta_ls;%恢复出的theta  
  50. end  
function [ theta ] = CS_SP( y,A,K )
%CS_SP Summary of this function goes here
%Version: 1.0 written by jbb0523 @2015-05-01
%   Detailed explanation goes here
%   y = Phi * x
%   x = Psi * theta
%	y = Phi*Psi * theta
%   令 A = Phi*Psi, 则y=A*theta
%   K is the sparsity level
%   现在已知y和A，求theta
%   Reference:Dai W，Milenkovic O．Subspace pursuit for compressive sensing
%   signal reconstruction[J]．IEEE Transactions on Information Theory，
%   2009，55(5)：2230-2249.[y_rows,y_columns] = size(y);if y_rows<y_columnsy = y';%y should be a column vectorend[M,N] = size(A);%传感矩阵A为M*N矩阵theta = zeros(N,1);%用来存储恢复的theta(列向量)Pos_theta = [];%用来迭代过程中存储A被选择的列序号r_n = y;%初始化残差(residual)为yfor kk=1:K%最多迭代K次%(1) Identificationproduct = A'*r_n;%传感矩阵A各列与残差的内积[val,pos]=sort(abs(product),'descend');Js = pos(1:K);%选出内积值最大的K列%(2) Support MergerIs = union(Pos_theta,Js);%Pos_theta与Js并集%(3) Estimation%At的行数要大于列数，此为最小二乘的基础(列线性无关)if length(Is)<=MAt = A(:,Is);%将A的这几列组成矩阵Atelse%At的列数大于行数，列必为线性相关的,At'*At将不可逆break;%跳出for循环end%y=At*theta，以下求theta的最小二乘解(Least Square)theta_ls = (At'*At)^(-1)*At'*y;%最小二乘解%(4) Pruning[val,pos]=sort(abs(theta_ls),'descend');%(5) Sample UpdatePos_theta = Is(pos(1:K));theta_ls = theta_ls(pos(1:K));%At(:,pos(1:K))*theta_ls是y在At(:,pos(1:K))列空间上的正交投影r_n = y - At(:,pos(1:K))*theta_ls;%更新残差 if norm(r_n)<1e-6%Repeat the steps until r=0break;%跳出for循环endendtheta(Pos_theta)=theta_ls;%恢复出的theta
end
        鉴于SP与CoSaMP的极其相似性，这里就不再给出单次重构和测量数M与重构成功概率关系曲线绘制例程代码了，只需将CoSaMP中调用CS_CoSaMP函数的部分改为调用CS_SP即可，无须任何其它改动。这里给出对比两种重构算法所绘制的测量数M与重构成功概率关系曲线的例程代码，只有这样才可以看出两种算法的重构性能优劣，以下是在分别运行完SP与CoSaMP的测量数M与重构成功概率关系曲线绘制例程代码的基础上，即已经存储了数据CoSaMPMtoPercentage1000.mat和SPMtoPercentage1000.mat：
[plain] view plain copy
print ?
  1. clear all;close all;clc;  
  2. load CoSaMPMtoPercentage1000;  
  3. PercentageCoSaMP = Percentage;  
  4. load SPMtoPercentage1000;  
  5. PercentageSP = Percentage;  
  6. S1 = ['-ks';'-ko';'-kd';'-kv';'-k*'];  
  7. S2 = ['-rs';'-ro';'-rd';'-rv';'-r*'];  
  8. figure;  
  9. for kk = 1:length(K_set)  
  10.     K = K_set(kk);  
  11.     M_set = 2*K:5:N;  
  12.     L_Mset = length(M_set);  
  13.     plot(M_set,PercentageCoSaMP(kk,1:L_Mset),S1(kk,:));%绘出x的恢复信号  
  14.     hold on;      
  15.     plot(M_set,PercentageSP(kk,1:L_Mset),S2(kk,:));%绘出x的恢复信号  
  16. end  
  17. hold off;  
  18. xlim([0 256]);  
  19. legend('CoSaK=4','SPK=4','CoSaK=12','SPK=12','CoSaK=20',...  
  20.     'SPK=20','CoSaK=28','SPK=28','CoSaK=36','SPK=36');  
  21. xlabel('Number of measurements(M)');  
  22. ylabel('Percentage recovered');  
  23. title('Percentage of input signals recovered correctly(N=256)(Gaussian)');  
clear all;close all;clc;
load CoSaMPMtoPercentage1000;
PercentageCoSaMP = Percentage;
load SPMtoPercentage1000;
PercentageSP = Percentage;
S1 = ['-ks';'-ko';'-kd';'-kv';'-k*'];
S2 = ['-rs';'-ro';'-rd';'-rv';'-r*'];
figure;
for kk = 1:length(K_set)K = K_set(kk);M_set = 2*K:5:N;L_Mset = length(M_set);plot(M_set,PercentageCoSaMP(kk,1:L_Mset),S1(kk,:));%绘出x的恢复信号hold on;    plot(M_set,PercentageSP(kk,1:L_Mset),S2(kk,:));%绘出x的恢复信号
end
hold off;
xlim([0 256]);
legend('CoSaK=4','SPK=4','CoSaK=12','SPK=12','CoSaK=20',...'SPK=20','CoSaK=28','SPK=28','CoSaK=36','SPK=36');
xlabel('Number of measurements(M)');
ylabel('Percentage recovered');
title('Percentage of input signals recovered correctly(N=256)(Gaussian)');

        运行结果如下：

        可以发现在M较小时SP略好于CoSaMP，当M变大时二者重构性能几乎一样。

 参考文献：

[1]Dai W，Milenkovic O. Subspace pursuitfor compressive sensing: Closing the gap between performance and complexity.http://arxiv.org/pdf/0803.0811v1.pdf

[2]Dai W，Milenkovic O．Subspacepursuit for compressive sensing signal reconstruction[J]．IEEETransactions on Information Theory，2009，55(5)：2230-2249.

[3]D. Needell, J.A. Tropp．CoSaMP: Iterative signal recoveryfrom incomplete and inaccurate samples. http://arxiv.org/pdf/0803.2392v2.pdf

[4]D. Needell, J.A. Tropp．CoSaMP: Iterative signal recoveryfrom incomplete and inaccurate samples. http://arxiv.org/pdf/0803.2392v1.pdf

[5]杨真真，杨震，孙林慧.信号压缩重构的正交匹配追踪类算法综述[J]. 信号处理，2013，29(4)：486-496.

[6]Li Zeng. CS_Reconstruction. http://www.pudn.com/downloads518/sourcecode/math/detail2151378.html

这篇关于压缩感知重构算法之子空间追踪(SP)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！

原文地址:https://blog.csdn.net/yuehuihui00/article/details/78327044
本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.chinasem.cn/article/159540。

相关文章

MySQL启动报错：InnoDB表空间丢失问题及解决方法

MySQL启动报错：InnoDB表空间丢失问题及解决方法

《MySQL启动报错：InnoDB表空间丢失问题及解决方法》在启动MySQL时,遇到了InnoDB:Tablespace5975wasnotfound,该错误表明MySQL在启动过程中无法找到指定的s... 目录mysql 启动报错：InnoDB 表空间丢失问题及解决方法错误分析解决方案1. 启用 inno
阅读更多...
SpringBoot实现文件记录日志及日志文件自动归档和压缩

SpringBoot实现文件记录日志及日志文件自动归档和压缩

《SpringBoot实现文件记录日志及日志文件自动归档和压缩》Logback是Java日志框架,通过Logger收集日志并经Appender输出至控制台、文件等,SpringBoot配置logbac... 目录1、什么是Logback2、SpringBoot实现文件记录日志，日志文件自动归档和压缩2.1、
阅读更多...
在Java中基于Geotools对PostGIS数据库的空间查询实践教程

在Java中基于Geotools对PostGIS数据库的空间查询实践教程

《在Java中基于Geotools对PostGIS数据库的空间查询实践教程》本文将深入探讨这一实践,从连接配置到复杂空间查询操作,包括点查询、区域范围查询以及空间关系判断等,全方位展示如何在Java环... 目录前言一、相关技术背景介绍1、评价对象AOI2、数据处理流程二、对AOI空间范围查询实践1、空间查
阅读更多...
MySQL表空间结构详解表空间到段页操作

MySQL表空间结构详解表空间到段页操作

《MySQL表空间结构详解表空间到段页操作》在MySQL架构和存储引擎专题中介绍了使用不同存储引擎创建表时生成的表空间数据文件,在本章节主要介绍使用InnoDB存储引擎创建表时生成的表空间数据文件,对... 目录️‍一、什么是表空间结构1.1 表空间与表空间文件的关系是什么?️‍二、用户数据在表空间中是怎么
阅读更多...
CentOS7增加Swap空间的两种方法

CentOS7增加Swap空间的两种方法

《CentOS7增加Swap空间的两种方法》当服务器物理内存不足时,增加Swap空间可以作为虚拟内存使用,帮助系统处理内存压力,本文给大家介绍了CentOS7增加Swap空间的两种方法：创建新的Swa... 目录在Centos 7上增加Swap空间的方法方法一：创建新的Swap文件（推荐）方法二：调整Sww
阅读更多...
使用雪花算法产生id导致前端精度缺失问题解决方案

使用雪花算法产生id导致前端精度缺失问题解决方案

《使用雪花算法产生id导致前端精度缺失问题解决方案》雪花算法由Twitter提出,设计目的是生成唯一的、递增的ID,下面：本文主要介绍使用雪花算法产生id导致前端精度缺失问题的解决方案,文中通过代... 目录一、问题根源二、解决方案1. 全局配置Jackson序列化规则2. 实体类必须使用Long封装类3.
阅读更多...
Springboot实现推荐系统的协同过滤算法

Springboot实现推荐系统的协同过滤算法

《Springboot实现推荐系统的协同过滤算法》协同过滤算法是一种在推荐系统中广泛使用的算法,用于预测用户对物品（如商品、电影、音乐等）的偏好,从而实现个性化推荐,下面给大家介绍Springboot... 目录前言基本原理 算法分类 计算方法应用场景 代码实现 前言协同过滤算法（Collaborativ
阅读更多...
使用Python实现矢量路径的压缩、解压与可视化

使用Python实现矢量路径的压缩、解压与可视化

《使用Python实现矢量路径的压缩、解压与可视化》在图形设计和Web开发中,矢量路径数据的高效存储与传输至关重要,本文将通过一个Python示例,展示如何将复杂的矢量路径命令序列压缩为JSON格式,... 目录引言核心功能概述1. 路径命令解析2. 路径数据压缩3. 路径数据解压4. 可视化代码实现详解1
阅读更多...
SpringBoot使用GZIP压缩反回数据问题

SpringBoot使用GZIP压缩反回数据问题

《SpringBoot使用GZIP压缩反回数据问题》：本文主要介绍SpringBoot使用GZIP压缩反回数据问题,具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不吝赐教... 目录SpringBoot使用GZIP压缩反回数据1、初识gzip2、gzip是什么，可以干什么？3、Spr
阅读更多...
openCV中KNN算法的实现

openCV中KNN算法的实现

《openCV中KNN算法的实现》KNN算法是一种简单且常用的分类算法,本文主要介绍了openCV中KNN算法的实现,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的... 目录KNN算法流程使用OpenCV实现KNNOpenCV 是一个开源的跨平台计算机视觉库，它提供了各
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2