压缩感知重构算法十一

### 压缩感知重构算法概述

压缩感知理论允许信号在远低于奈奎斯特采样率的情况下被有效采集和重建。对于稀疏或可压缩信号，通过少量线性测量即可实现高精度恢复。CoSaMP (Compressed Sampling Matching Pursuit) 是一种高效的迭代贪婪算法，在每次迭代过程中选择最可能属于原信号支撑集的位置并更新估计值[^1]。

该算法具有多项优势特性：

- **快速收敛**：通常只需几次迭代就能达到稳定解
- **精确度高**：能够以接近最优的方式恢复原始信号
- **计算复杂度适中**：相比其他方法更易于硬件实现

### CoSaMP 算法流程描述

以下是简化版的 CoSaMP 工作原理说明：

1. 初始化残差向量 r0=y ，设定最大迭代次数 K 和支持集合大小 S。
2. 对于每一轮迭代 k=1,..K :
- 找到 |g| 中最大的 2S 个位置形成候选支持 T'
- 解最小二倍范数问题求得新近似 x'=argmin ||A_T' z-y||_2 ,其中 A_T' 表示仅保留列索引位于 T' 的子矩阵
- 更新最终支持 T←support(x')[:S],即取前 S 大绝对值分量对应的下标作为新的支持
- 使用 LS 方法重新拟合这些位置上的系数 x[T]=inv(A[:,T]) * y 得到改进后的逼近
- 更新剩余误差 r[k]=y-A*x

此过程持续至满足停止条件为止，如达到预设的最大迭代次数或者连续两次迭代间变化幅度小于阈值。

import numpy as np

def cosamp(A, b, s):
    m, n = A.shape
    
    # Initialization
    residual = b.copy()
    support_set = set()

    max_iterations = 10*s
    for _ in range(max_iterations):

        # Compute the product of measurement matrix and current residuals.
        gradient = abs(np.dot(A.T, residual))
        
        # Identify positions with largest magnitude entries.
        candidate_support = np.argpartition(-gradient, 2*s)[:2*s]

        combined_support = list(support_set.union(set(candidate_support)))
        
        if not combined_support:
            break
            
        # Solve least squares problem on selected columns only.
        ls_solution = np.linalg.lstsq(A[:,combined_support],b,rcond=None)[0]
        
        temp_signal=np.zeros(n)
        temp_signal[list(combined_support)] = ls_solution
        
        # Prune to top-s elements by absolute value.
        pruned_indices = np.argsort(abs(temp_signal))[-s:]
        new_estimate = np.zeros_like(b)
        new_estimate[pruned_indices] = temp_signal[pruned_indices]
        
        support_set=set(pruned_indices)

        # Update residual error term using updated estimate.
        residual=b-np.dot(A,new_estimate)
    
    return new_estimate