AMP算法实战用Python手把手复现Donoho经典论文中的稀疏信号恢复稀疏信号恢复是信号处理领域的核心问题之一而近似消息传递(AMP)算法因其高效性和理论保证成为该领域的重要工具。本文将带您用Python一步步实现Donoho经典论文中的AMP算法完全避开复杂的数学推导专注于代码实现、调参和可视化。1. 环境准备与数据生成在开始AMP算法实现前我们需要准备好Python环境和生成符合要求的模拟数据。AMP算法处理的核心问题是给定欠采样观测向量y和感知矩阵A如何恢复原始稀疏信号x。首先安装必要的库pip install numpy matplotlib scikit-learn接着生成模拟数据。我们创建一个长度为N1000的稀疏信号其中只有50个非零元素import numpy as np import matplotlib.pyplot as plt # 参数设置 N 1000 # 信号维度 M 300 # 观测维度 K 50 # 稀疏度 sigma 0.1 # 噪声标准差 # 生成稀疏信号 x_true np.zeros(N) nonzero_indices np.random.choice(N, K, replaceFalse) x_true[nonzero_indices] np.random.randn(K) # 生成感知矩阵高斯随机矩阵 A np.random.randn(M, N) / np.sqrt(M) # 生成含噪声的观测 y A x_true sigma * np.random.randn(M)关键点说明感知矩阵A需要满足RIP(Restricted Isometry Property)条件高斯随机矩阵是常用选择观测噪声通常假设为高斯白噪声稀疏信号的非零位置随机分布幅值可以是高斯分布2. AMP算法核心实现AMP算法的核心在于迭代更新信号估计和残差。算法通过软阈值函数逐步修正信号估计同时考虑感知矩阵的反馈效应。2.1 软阈值函数实现软阈值函数是AMP中的关键非线性操作用于促进解的稀疏性def soft_threshold(x, threshold): 软阈值函数 :param x: 输入信号 :param threshold: 阈值 :return: 阈值处理后的信号 return np.sign(x) * np.maximum(np.abs(x) - threshold, 0)为了高效计算我们还需要软阈值函数的导数用于状态演化def soft_threshold_derivative(x, threshold): 软阈值函数的导数 return (np.abs(x) threshold).astype(float)2.2 AMP迭代过程完整的AMP算法实现如下def AMP(y, A, max_iter100, tol1e-5, verboseFalse): AMP算法实现 :param y: 观测向量 :param A: 感知矩阵 :param max_iter: 最大迭代次数 :param tol: 收敛阈值 :param verbose: 是否打印迭代信息 :return: 恢复的信号 M, N A.shape x np.zeros(N) # 初始信号估计 z y.copy() # 初始残差 tau np.var(y) # 初始噪声方差估计 history {x: [], tau: [], mse: []} for t in range(max_iter): # 计算有效观测 r x A.T z # 软阈值处理 x_new soft_threshold(r, tau) # 计算残差更新 div np.mean(soft_threshold_derivative(r, tau)) z y - A x_new z * div # 更新噪声方差估计 tau np.sqrt(np.mean(z**2)) # 记录历史 history[x].append(x_new.copy()) history[tau].append(tau) mse np.mean((x_new - x_true)**2) history[mse].append(mse) # 检查收敛 if np.linalg.norm(x_new - x) / np.linalg.norm(x_new) tol: if verbose: print(f收敛于迭代 {t1}) break x x_new return x, history算法参数说明参数描述典型值max_iter最大迭代次数50-100tol收敛阈值1e-5tau_init初始阈值np.var(y)3. 算法可视化与性能分析理解AMP算法行为的最佳方式是通过可视化观察其迭代过程。我们将从多个角度分析算法性能。3.1 信号恢复过程可视化def plot_recovery_process(history, x_true): plt.figure(figsize(15, 10)) # 绘制MSE变化 plt.subplot(2, 2, 1) plt.plot(history[mse], o-) plt.xlabel(迭代次数) plt.ylabel(MSE) plt.title(均方误差变化) plt.grid(True) # 绘制tau变化 plt.subplot(2, 2, 2) plt.plot(history[tau], o-) plt.xlabel(迭代次数) plt.ylabel(Tau) plt.title(阈值参数变化) plt.grid(True) # 绘制最终恢复信号 plt.subplot(2, 2, 3) plt.stem(x_true, markerfmtbo, label真实信号) plt.stem(history[x][-1], markerfmtro, label恢复信号) plt.xlabel(索引) plt.ylabel(幅值) plt.title(信号恢复对比) plt.legend() plt.grid(True) # 绘制残差变化 residuals [np.linalg.norm(y - A x) for x in history[x]] plt.subplot(2, 2, 4) plt.plot(residuals, o-) plt.xlabel(迭代次数) plt.ylabel(残差范数) plt.title(观测残差变化) plt.grid(True) plt.tight_layout() plt.show()3.2 与Lasso的对比为了评估AMP性能我们将其与Scikit-learn中的Lasso实现进行对比from sklearn.linear_model import Lasso def compare_with_lasso(y, A, x_true): # AMP算法 x_amp, history AMP(y, A, verboseTrue) mse_amp np.mean((x_amp - x_true)**2) # Lasso算法 # 通过交叉验证选择最佳alpha lasso Lasso(alpha0.1, max_iter10000) lasso.fit(A, y) x_lasso lasso.coef_ mse_lasso np.mean((x_lasso - x_true)**2) # 绘制对比图 plt.figure(figsize(12, 6)) plt.stem(x_true, markerfmtbo, label真实信号) plt.stem(x_amp, markerfmtro, labelAMP恢复) plt.stem(x_lasso, markerfmtgx, labelLasso恢复) plt.legend() plt.title(fAMP MSE: {mse_amp:.2e}, Lasso MSE: {mse_lasso:.2e}) plt.grid(True) plt.show() return mse_amp, mse_lasso性能对比要点AMP通常比Lasso收敛更快AMP对噪声更鲁棒AMP不需要精细调节正则化参数AMP提供状态演化理论预测性能4. 参数调优与实用技巧在实际应用中AMP算法的性能很大程度上取决于参数设置和实现细节。以下是几个关键调优点4.1 初始阈值选择初始阈值τ^0的选择至关重要。实践中发现以下策略有效# 改进的初始阈值选择 tau_init np.median(np.abs(A.T y)) / np.sqrt(np.log(N))提示初始阈值应与噪声水平匹配太大会导致过度稀疏太小则收敛慢4.2 迭代停止条件除了相对变化阈值还可以考虑# 改进的停止条件 if (np.mean((y - A x_new)**2) 1.1 * M * sigma**2 and np.linalg.norm(x_new - x) / np.linalg.norm(x_new) tol): break4.3 阻尼技术为防止振荡可以引入阻尼因子# 添加阻尼的AMP更新 damping 0.1 x_new (1 - damping) * soft_threshold(r, tau) damping * x5. 扩展到实际应用场景AMP算法不仅限于理论模型还可以应用于多种实际问题5.1 图像修复def image_inpainting(image, mask): 使用AMP进行图像修复 :param image: 输入图像(灰度) :param mask: 采样掩模(1表示保留0表示缺失) :return: 修复后的图像 # 将图像转换为向量 y image[mask 1] M len(y) N image.size # 构建感知矩阵(采样算子) A np.zeros((M, N)) pos 0 for i in range(image.shape[0]): for j in range(image.shape[1]): if mask[i, j] 1: A[pos, i * image.shape[1] j] 1 pos 1 # 使用AMP恢复 x_rec, _ AMP(y, A) # 恢复图像 recovered x_rec.reshape(image.shape) return recovered5.2 传感器网络数据重构在传感器网络中AMP可用于从部分观测中恢复完整数据def sensor_data_recovery(partial_data, sensor_locations, grid_size): 传感器网络数据恢复 :param partial_data: 部分观测数据 :param sensor_locations: 传感器位置列表 :param grid_size: 完整网格大小 (nx, ny) :return: 完整网格数据 nx, ny grid_size N nx * ny M len(partial_data) # 构建感知矩阵 A np.zeros((M, N)) for i, (x, y) in enumerate(sensor_locations): idx x * ny y A[i, idx] 1 # AMP恢复 x_rec, _ AMP(partial_data, A) return x_rec.reshape((nx, ny))6. 高级主题与性能优化对于需要更高性能的场景我们可以从以下几个方面优化AMP实现6.1 矩阵自由实现对于大型问题显式存储A可能不可行。我们可以用线性算子代替class MatrixFreeOperator: def __init__(self, A_shape, forward_fn, adjoint_fn): self.shape A_shape self.forward forward_fn self.adjoint adjoint_fn def __matmul__(self, x): return self.forward(x) property def T(self): return MatrixFreeOperator( (self.shape[1], self.shape[0]), self.adjoint, self.forward ) # 使用示例 def forward(x): return np.dot(A_dense, x) def adjoint(x): return np.dot(A_dense.T, x) A_op MatrixFreeOperator((M, N), forward, adjoint)6.2 并行化实现利用多核CPU或GPU加速AMPfrom numba import jit, prange jit(nopythonTrue, parallelTrue) def parallel_AMP(y, A, max_iter100): M, N A.shape x np.zeros(N) z y.copy() tau np.var(y) for t in range(max_iter): r x np.dot(A.T, z) x_new np.zeros_like(x) for i in prange(N): if r[i] tau: x_new[i] r[i] - tau elif r[i] -tau: x_new[i] r[i] tau else: x_new[i] 0.0 div np.mean(np.abs(r) tau) z y - np.dot(A, x_new) z * div tau np.sqrt(np.mean(z**2)) if np.sum((x_new - x)**2) / np.sum(x_new**2) 1e-5: break x x_new return x6.3 自适应阈值调整基于噪声水平的自适应阈值策略def adaptive_AMP(y, A, max_iter100, sigma_estNone): M, N A.shape x np.zeros(N) z y.copy() if sigma_est is None: tau np.median(np.abs(A.T y)) / np.sqrt(np.log(N)) else: tau sigma_est * np.sqrt(2 * np.log(N)) for t in range(max_iter): # 更新信号估计 r x A.T z x_new soft_threshold(r, tau) # 更新残差 div np.mean(soft_threshold_derivative(r, tau)) z y - A x_new z * div # 估计噪声水平并调整阈值 sigma_est np.sqrt(np.mean(z**2)) tau sigma_est * np.sqrt(2 * np.log(N)) # 检查收敛 if np.linalg.norm(x_new - x) / np.linalg.norm(x_new) 1e-5: break x x_new return x在实际项目中AMP算法的表现往往取决于具体问题的特性。通过合理调整参数和优化实现可以使其在各种稀疏恢复任务中达到最佳性能。
AMP算法实战：用Python手把手复现Donoho经典论文中的稀疏信号恢复
发布时间：2026/5/28 18:25:47
AMP算法实战用Python手把手复现Donoho经典论文中的稀疏信号恢复稀疏信号恢复是信号处理领域的核心问题之一而近似消息传递(AMP)算法因其高效性和理论保证成为该领域的重要工具。本文将带您用Python一步步实现Donoho经典论文中的AMP算法完全避开复杂的数学推导专注于代码实现、调参和可视化。1. 环境准备与数据生成在开始AMP算法实现前我们需要准备好Python环境和生成符合要求的模拟数据。AMP算法处理的核心问题是给定欠采样观测向量y和感知矩阵A如何恢复原始稀疏信号x。首先安装必要的库pip install numpy matplotlib scikit-learn接着生成模拟数据。我们创建一个长度为N1000的稀疏信号其中只有50个非零元素import numpy as np import matplotlib.pyplot as plt # 参数设置 N 1000 # 信号维度 M 300 # 观测维度 K 50 # 稀疏度 sigma 0.1 # 噪声标准差 # 生成稀疏信号 x_true np.zeros(N) nonzero_indices np.random.choice(N, K, replaceFalse) x_true[nonzero_indices] np.random.randn(K) # 生成感知矩阵高斯随机矩阵 A np.random.randn(M, N) / np.sqrt(M) # 生成含噪声的观测 y A x_true sigma * np.random.randn(M)关键点说明感知矩阵A需要满足RIP(Restricted Isometry Property)条件高斯随机矩阵是常用选择观测噪声通常假设为高斯白噪声稀疏信号的非零位置随机分布幅值可以是高斯分布2. AMP算法核心实现AMP算法的核心在于迭代更新信号估计和残差。算法通过软阈值函数逐步修正信号估计同时考虑感知矩阵的反馈效应。2.1 软阈值函数实现软阈值函数是AMP中的关键非线性操作用于促进解的稀疏性def soft_threshold(x, threshold): 软阈值函数 :param x: 输入信号 :param threshold: 阈值 :return: 阈值处理后的信号 return np.sign(x) * np.maximum(np.abs(x) - threshold, 0)为了高效计算我们还需要软阈值函数的导数用于状态演化def soft_threshold_derivative(x, threshold): 软阈值函数的导数 return (np.abs(x) threshold).astype(float)2.2 AMP迭代过程完整的AMP算法实现如下def AMP(y, A, max_iter100, tol1e-5, verboseFalse): AMP算法实现 :param y: 观测向量 :param A: 感知矩阵 :param max_iter: 最大迭代次数 :param tol: 收敛阈值 :param verbose: 是否打印迭代信息 :return: 恢复的信号 M, N A.shape x np.zeros(N) # 初始信号估计 z y.copy() # 初始残差 tau np.var(y) # 初始噪声方差估计 history {x: [], tau: [], mse: []} for t in range(max_iter): # 计算有效观测 r x A.T z # 软阈值处理 x_new soft_threshold(r, tau) # 计算残差更新 div np.mean(soft_threshold_derivative(r, tau)) z y - A x_new z * div # 更新噪声方差估计 tau np.sqrt(np.mean(z**2)) # 记录历史 history[x].append(x_new.copy()) history[tau].append(tau) mse np.mean((x_new - x_true)**2) history[mse].append(mse) # 检查收敛 if np.linalg.norm(x_new - x) / np.linalg.norm(x_new) tol: if verbose: print(f收敛于迭代 {t1}) break x x_new return x, history算法参数说明参数描述典型值max_iter最大迭代次数50-100tol收敛阈值1e-5tau_init初始阈值np.var(y)3. 算法可视化与性能分析理解AMP算法行为的最佳方式是通过可视化观察其迭代过程。我们将从多个角度分析算法性能。3.1 信号恢复过程可视化def plot_recovery_process(history, x_true): plt.figure(figsize(15, 10)) # 绘制MSE变化 plt.subplot(2, 2, 1) plt.plot(history[mse], o-) plt.xlabel(迭代次数) plt.ylabel(MSE) plt.title(均方误差变化) plt.grid(True) # 绘制tau变化 plt.subplot(2, 2, 2) plt.plot(history[tau], o-) plt.xlabel(迭代次数) plt.ylabel(Tau) plt.title(阈值参数变化) plt.grid(True) # 绘制最终恢复信号 plt.subplot(2, 2, 3) plt.stem(x_true, markerfmtbo, label真实信号) plt.stem(history[x][-1], markerfmtro, label恢复信号) plt.xlabel(索引) plt.ylabel(幅值) plt.title(信号恢复对比) plt.legend() plt.grid(True) # 绘制残差变化 residuals [np.linalg.norm(y - A x) for x in history[x]] plt.subplot(2, 2, 4) plt.plot(residuals, o-) plt.xlabel(迭代次数) plt.ylabel(残差范数) plt.title(观测残差变化) plt.grid(True) plt.tight_layout() plt.show()3.2 与Lasso的对比为了评估AMP性能我们将其与Scikit-learn中的Lasso实现进行对比from sklearn.linear_model import Lasso def compare_with_lasso(y, A, x_true): # AMP算法 x_amp, history AMP(y, A, verboseTrue) mse_amp np.mean((x_amp - x_true)**2) # Lasso算法 # 通过交叉验证选择最佳alpha lasso Lasso(alpha0.1, max_iter10000) lasso.fit(A, y) x_lasso lasso.coef_ mse_lasso np.mean((x_lasso - x_true)**2) # 绘制对比图 plt.figure(figsize(12, 6)) plt.stem(x_true, markerfmtbo, label真实信号) plt.stem(x_amp, markerfmtro, labelAMP恢复) plt.stem(x_lasso, markerfmtgx, labelLasso恢复) plt.legend() plt.title(fAMP MSE: {mse_amp:.2e}, Lasso MSE: {mse_lasso:.2e}) plt.grid(True) plt.show() return mse_amp, mse_lasso性能对比要点AMP通常比Lasso收敛更快AMP对噪声更鲁棒AMP不需要精细调节正则化参数AMP提供状态演化理论预测性能4. 参数调优与实用技巧在实际应用中AMP算法的性能很大程度上取决于参数设置和实现细节。以下是几个关键调优点4.1 初始阈值选择初始阈值τ^0的选择至关重要。实践中发现以下策略有效# 改进的初始阈值选择 tau_init np.median(np.abs(A.T y)) / np.sqrt(np.log(N))提示初始阈值应与噪声水平匹配太大会导致过度稀疏太小则收敛慢4.2 迭代停止条件除了相对变化阈值还可以考虑# 改进的停止条件 if (np.mean((y - A x_new)**2) 1.1 * M * sigma**2 and np.linalg.norm(x_new - x) / np.linalg.norm(x_new) tol): break4.3 阻尼技术为防止振荡可以引入阻尼因子# 添加阻尼的AMP更新 damping 0.1 x_new (1 - damping) * soft_threshold(r, tau) damping * x5. 扩展到实际应用场景AMP算法不仅限于理论模型还可以应用于多种实际问题5.1 图像修复def image_inpainting(image, mask): 使用AMP进行图像修复 :param image: 输入图像(灰度) :param mask: 采样掩模(1表示保留0表示缺失) :return: 修复后的图像 # 将图像转换为向量 y image[mask 1] M len(y) N image.size # 构建感知矩阵(采样算子) A np.zeros((M, N)) pos 0 for i in range(image.shape[0]): for j in range(image.shape[1]): if mask[i, j] 1: A[pos, i * image.shape[1] j] 1 pos 1 # 使用AMP恢复 x_rec, _ AMP(y, A) # 恢复图像 recovered x_rec.reshape(image.shape) return recovered5.2 传感器网络数据重构在传感器网络中AMP可用于从部分观测中恢复完整数据def sensor_data_recovery(partial_data, sensor_locations, grid_size): 传感器网络数据恢复 :param partial_data: 部分观测数据 :param sensor_locations: 传感器位置列表 :param grid_size: 完整网格大小 (nx, ny) :return: 完整网格数据 nx, ny grid_size N nx * ny M len(partial_data) # 构建感知矩阵 A np.zeros((M, N)) for i, (x, y) in enumerate(sensor_locations): idx x * ny y A[i, idx] 1 # AMP恢复 x_rec, _ AMP(partial_data, A) return x_rec.reshape((nx, ny))6. 高级主题与性能优化对于需要更高性能的场景我们可以从以下几个方面优化AMP实现6.1 矩阵自由实现对于大型问题显式存储A可能不可行。我们可以用线性算子代替class MatrixFreeOperator: def __init__(self, A_shape, forward_fn, adjoint_fn): self.shape A_shape self.forward forward_fn self.adjoint adjoint_fn def __matmul__(self, x): return self.forward(x) property def T(self): return MatrixFreeOperator( (self.shape[1], self.shape[0]), self.adjoint, self.forward ) # 使用示例 def forward(x): return np.dot(A_dense, x) def adjoint(x): return np.dot(A_dense.T, x) A_op MatrixFreeOperator((M, N), forward, adjoint)6.2 并行化实现利用多核CPU或GPU加速AMPfrom numba import jit, prange jit(nopythonTrue, parallelTrue) def parallel_AMP(y, A, max_iter100): M, N A.shape x np.zeros(N) z y.copy() tau np.var(y) for t in range(max_iter): r x np.dot(A.T, z) x_new np.zeros_like(x) for i in prange(N): if r[i] tau: x_new[i] r[i] - tau elif r[i] -tau: x_new[i] r[i] tau else: x_new[i] 0.0 div np.mean(np.abs(r) tau) z y - np.dot(A, x_new) z * div tau np.sqrt(np.mean(z**2)) if np.sum((x_new - x)**2) / np.sum(x_new**2) 1e-5: break x x_new return x6.3 自适应阈值调整基于噪声水平的自适应阈值策略def adaptive_AMP(y, A, max_iter100, sigma_estNone): M, N A.shape x np.zeros(N) z y.copy() if sigma_est is None: tau np.median(np.abs(A.T y)) / np.sqrt(np.log(N)) else: tau sigma_est * np.sqrt(2 * np.log(N)) for t in range(max_iter): # 更新信号估计 r x A.T z x_new soft_threshold(r, tau) # 更新残差 div np.mean(soft_threshold_derivative(r, tau)) z y - A x_new z * div # 估计噪声水平并调整阈值 sigma_est np.sqrt(np.mean(z**2)) tau sigma_est * np.sqrt(2 * np.log(N)) # 检查收敛 if np.linalg.norm(x_new - x) / np.linalg.norm(x_new) 1e-5: break x x_new return x在实际项目中AMP算法的表现往往取决于具体问题的特性。通过合理调整参数和优化实现可以使其在各种稀疏恢复任务中达到最佳性能。
