工程科学学报，第39卷.第8期：1254-1260.2017年8月 Chinese Journal of Engineering,Vol.39,No.8:1254-1260,August 2017 DOI:10.13374/j.issn2095-9389.2017.08.016;http://journals.ustb.edu.cn MMSE准则下基于玻尔兹曼机的快速重构算法 刘玲君2)四，谢中华)，冯久超”，杨萃12) 1)华南理工大学电子与信息学院，广州5106412)国家移动超声探测工程技术研究中心，广州510641 区通信作者，E-mail:liuO@gmail.com 摘要全连接的玻尔兹曼机模型可全面描述稀疏系数间统计依赖关系，但时间复杂度较高.为了提高基于玻尔兹曼机的 贝叶斯匹配追踪算法(BM-BMP)的重构速度和质量，本文提出一种改进算法.第一，将BM-BMP算法的最大后验概率(MAP) 估计评估值分解为上一次迭代的评估值与增量，使得每次迭代仅需计算增量，极大缩短了计算耗时.第二，利用显著最大后 验概率估计值平均的方式，有效近似最小均方误差(MMSE)估计，获得了更小的重构误差.实验结果表明，本文算法比 BM-BMP算法的运行时间平均缩短了73.66%，峰值信噪比(PSNR)值平均提高了0.57dB. 关键词稀疏信号重构；快速贝叶斯匹配追踪；玻尔兹曼机；最小均方误差 分类号TP391 Fast recovery algorithm based on Boltzmann machine and MMSE criterion LIU Ling-jun,XIE Zhong-hua),FENG Jiu-chao,YANG Cui) 1)School of Electronic and Information Engineering,South China University of Technology,Guangzhou 510641,China 2)National Engineering Technology Research Center for Mobile Ultrasonic Detection,Guangzhou 510641,China Corresponding author,E-mail:ljliu@gmail.com ABSTRACT Fully connected Boltzmann machine models can be used to provide a comprehensive description of statistical dependen- cies between sparse coefficients but with high time complexity.To improve the speed and quality of the Boltzmann machine-Bayesian matching pursuit (BM-BMP)method,an improved algorithm was proposed.First,the maximum a posteriori (MAP)estimation of the BM-BMP algorithm is decomposed into its value at the last iteration and an increment;thus,it only needs to calculate the increment in each iteration,which greatly reduces the computational time.Second,by calculating the mean of the significant MAP estimations,an effective approximation is obtained for the minimum mean square error(MMSE)estimation and a smaller reconstruction error is a- chieved.Compared with the BM-BMP,this method reduces the running time on average by 73.66%while improving the peak signal to noise ratio PSNR)by 0.57 dB. KEY WORDS sparse signal reconstruction;fast Bayesian matching pursuit;Boltzmann machine;minimum mean square error MMSE) 稀疏信号重构被广泛应用于信号和图像处理中，y的长度)，其中x是m×1的稀疏系数(m是矢量x的 例如：去噪、盲源分离、波达方向估计、压缩和采样、图长度)，A是n×m的冗余字典（因为m>n,所以A是 像复原等-).假设n×1维信号y=Ar+e(n是矢量冗余字典)，e是高斯白噪声，其均值为零方差为2. 收稿日期：2016-09-12 基金项目：国家自然科学基金资助项目(61327005,61302120)：广东省科技计划资助项目(2017A020214011):中央高校基本科研业务费资助 项目(2017M039)
工程科学学报,第 39 卷,第 8 期:1254鄄鄄1260,2017 年 8 月 Chinese Journal of Engineering, Vol. 39, No. 8: 1254鄄鄄1260, August 2017 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2017. 08. 016; http: / / journals. ustb. edu. cn MMSE 准则下基于玻尔兹曼机的快速重构算法 刘玲君1,2)苣 , 谢中华1) , 冯久超1) , 杨 萃1,2) 1) 华南理工大学电子与信息学院, 广州 510641 2) 国家移动超声探测工程技术研究中心, 广州 510641 苣 通信作者, E鄄mail: ljliu0@ gmail. com 摘 要 全连接的玻尔兹曼机模型可全面描述稀疏系数间统计依赖关系,但时间复杂度较高. 为了提高基于玻尔兹曼机的 贝叶斯匹配追踪算法(BM鄄BMP)的重构速度和质量,本文提出一种改进算法. 第一,将 BM鄄BMP 算法的最大后验概率(MAP) 估计评估值分解为上一次迭代的评估值与增量,使得每次迭代仅需计算增量,极大缩短了计算耗时. 第二,利用显著最大后 验概率估计值平均的方式,有效近似最小均方误差(MMSE) 估计,获得了更小的重构误差. 实验结果表明,本文算法比 BM鄄BMP算法的运行时间平均缩短了 73郾 66% ,峰值信噪比(PSNR)值平均提高了 0郾 57 dB. 关键词 稀疏信号重构; 快速贝叶斯匹配追踪; 玻尔兹曼机; 最小均方误差 分类号 TP391 Fast recovery algorithm based on Boltzmann machine and MMSE criterion LIU Ling鄄jun 1,2)苣 , XIE Zhong鄄hua 1) , FENG Jiu鄄chao 1) , YANG Cui 1,2) 1) School of Electronic and Information Engineering, South China University of Technology, Guangzhou 510641, China 2) National Engineering Technology Research Center for Mobile Ultrasonic Detection, Guangzhou 510641, China 苣 Corresponding author, E鄄mail: ljliu0@ gmail. com ABSTRACT Fully connected Boltzmann machine models can be used to provide a comprehensive description of statistical dependen鄄 cies between sparse coefficients but with high time complexity. To improve the speed and quality of the Boltzmann machine鄄Bayesian matching pursuit (BM鄄BMP) method, an improved algorithm was proposed. First, the maximum a posteriori (MAP) estimation of the BM鄄BMP algorithm is decomposed into its value at the last iteration and an increment; thus, it only needs to calculate the increment in each iteration, which greatly reduces the computational time. Second, by calculating the mean of the significant MAP estimations, an effective approximation is obtained for the minimum mean square error (MMSE) estimation and a smaller reconstruction error is a鄄 chieved. Compared with the BM鄄BMP, this method reduces the running time on average by 73郾 66% while improving the peak signal to noise ratio (PSNR) by 0郾 57 dB. KEY WORDS sparse signal reconstruction; fast Bayesian matching pursuit; Boltzmann machine; minimum mean square error (MMSE) 收稿日期: 2016鄄鄄09鄄鄄12 基金项目: 国家自然科学基金资助项目(61327005, 61302120); 广东省科技计划资助项目(2017A020214011); 中央高校基本科研业务费资助 项目(2017MS039) 稀疏信号重构被广泛应用于信号和图像处理中, 例如:去噪、盲源分离、波达方向估计、压缩和采样、图 像复原等[1鄄鄄2] . 假设 n 伊 1 维信号 y = Ax + e( n 是矢量 y 的长度),其中 x 是 m 伊 1 的稀疏系数(m 是矢量 x 的 长度),A 是 n 伊 m 的冗余字典(因为 m > n,所以 A 是 冗余字典),e 是高斯白噪声,其均值为零方差为 滓 2
刘玲君等：MMSE准则下基于玻尔兹曼机的快速重构算法 ·1255· 稀疏信号重构的任务是获取含噪信号y在冗余字典A 广.文献「121首先通过图像训练库学习BM模型的参 下的最优稀疏逼近，进而恢复出无噪信号A. 数：系数的方差、系数间的相关矩阵和其偏差以及冗余 由于0范数的非凸性，求重构问题的全局最优解 字典：然后根据这些参数使用吉比斯采样和模拟退火 是一个非确定多项式(non-deterministic polynomial, 算法求解关于稀疏模式的MAP估计，进而实现信号重 NP)难问题，因此只能采用近似逼近的算法求解.常 构，但该方法收敛速度慢.为此，结构化贝叶斯追踪 用的近似逼近算法可分为凸松弛方法、贪婪方法和贝 (structured soft bayesian pursuit,SSoBaP)[]和基于玻 叶斯重构方法.其中，凸松弛方法如基追踪(basis 尔兹曼机的贝叶斯匹配追踪(boltzmann machine bayes- pursuit,BP)用凸的、易求解的稀疏度量函数代替非凸 ian matching pursuit,BM-BMP)[)研究了更快速的重 的0范数来逼近原问题，重构效果理想，但复杂度高. 构算法.SSoBaP采用基于平均场近似的变分贝叶斯方 贪婪方法在冗余字典中迭代寻找与当前残差最匹配的 法求解MAP估计，通过迭代的方式构造易求解极值的 原子，近似估计原信号，复杂度较低，但重构效果一般， 概率函数最大化地近似MAP函数，得到MAP估计的 典型的算法有正交匹配追踪(orthogonal matching 一个局部近似.BM-BMP则采用了更加快速的贪婪算 pursuit,OMP))]和子空间追踪(subspace pursuit, 法进行MAP估计求解，在相同的计算复杂度下获得比 SP)).与这两类方法不同，贝叶斯重构利用了信号的 模拟退火算法更好的恢复性能.BM-BMP算法基于贝 统计先验信息，使用概率推断和软判决方法估计原信 叶斯匹配追踪的思想，每次迭代依次激活一个稀疏模 号，获得更好的重构质量，但复杂度上升，典型的算法 式的元素使得当前MAP评估值最大，通过若干次迭代 有贝叶斯压缩感知(bayesian compressive sensing, 后寻找出一个稀疏模式的最优组合.另外，BM-BMP BCS))和随机正交匹配追踪(Rand-OMP)[o].贝叶斯 算法还采用块坐标优化方法实现参数估计和信号重构 压缩感知使用层次先验模型来描述重构信号系数的稀 的同步求解，比预训练的模型参数更匹配信号内容 疏性，通过期望最大化(expectation maximization,EM) 综上所述，提高算法的运算速度和重构质量是贝 方法来估计模型的参数，并由贝叶斯公式求解问题的 叶斯重构方法的两个关键问题.本文在BM-BMP的 最大后验概率(maximum a posteriori,MAP)估计，从而 基础上，利用玻尔兹曼机模型来描述稀疏系数间结构 实现信号重构.Rand-OMP则根据MAP评估式按一 信息，并通过增量计算和MMSE近似估计的方法分别 定概率随机确定贪婪搜索的方向估计信号，并由多组 改善BM-BMP算法的运算速度和重构质量.第一，将 OMP估计值平均构成重构问题的最小均方误差(mii- MAP估计的评估值表示成上一次迭代的评估值加上 mum mean square error,MMSE)估计，并证明了当评判 增值部分，以增量计算替代评估值的计算，从而减少了 准则为最小误差而非最稀疏时，MMSE估计比单个 运行时间：第二，考虑到自然图像的非稀疏性，在噪声 MAP估计具有更小的重构误差. 干扰下，单个MAP估计解往往不稳定，为此，本文利用 传统的贝叶斯重构算法假设稀疏系数彼此独立， 多个显著MAP估计值近似MMSE估计，克服了MMSE 未利用系数间的结构特性，近来研究发现利用结构化 估计需遍历所有解的问题，也改善了BM-BMP算法的 稀疏信息可以加强重构问题的规范化约束，缩小重构 重构质量 信号需要搜索的解空间范围，因此可以更加准确有效 地重构信号，例如：分块正交匹配追踪(Blok-OMP)[] 1基于玻尔兹曼机的稀疏重构算法 利用信号的块稀疏结构，假设信号非零系数成块出现： 1.1玻尔兹曼机模型及MAP估计 基于模型的压缩感知(model-based CS,MB_CS)[)利用 在稀疏信号重构中，信号本身或者它的变换系数 信号的树稀疏结构，假设信号中的系数呈现树形连接 常具有稀疏的先验分布，此外，系数间还会呈现出聚集 的特征.但这两者将系数间的结构看作是确定关系， 的结构特性，如小波树结构和块结构.但对于更复杂 而实际信号中的结构用统计分布来描述更为可靠，如 的信号如非平稳的自然图像而言，变换系数不满足树 基于分块结构的稀疏贝叶斯学习(block sparse bayes- 稀疏或块稀疏的特性.为了描述系数间的任意结构特 ian learning,.BsBL)[s-o]和基于小波树结构的压缩感知 性，本文利用玻尔兹曼分布对稀疏系数间的统计依赖 tree-structured wavelet compressive sensing,TSW- 性进行先验描述 CS)t (1)玻尔兹曼机模型.BM分布定义如下： 使用玻尔兹曼机(Boltzmann machine,BM)模型表 示稀疏系数两两之间的统计依赖性2-)，可描述系数 p(S)=zp(b's+2s"ws） (1) 间的任意结构特性，比块稀疏或树稀疏结构适用性更 其中，稀疏模式S是m维的二进制序列，第i个元素
刘玲君等: MMSE 准则下基于玻尔兹曼机的快速重构算法 稀疏信号重构的任务是获取含噪信号 y 在冗余字典 A 下的最优稀疏逼近 x^,进而恢复出无噪信号 Ax^. 由于 L0 范数的非凸性,求重构问题的全局最优解 是一 个 非 确 定 多 项 式 ( non鄄deterministic polynomial, NP)难问题,因此只能采用近似逼近的算法求解. 常 用的近似逼近算法可分为凸松弛方法、贪婪方法和贝 叶斯重构方法. 其中,凸松弛方法如基 追 踪 ( basis pursuit, BP)用凸的、易求解的稀疏度量函数代替非凸 的 L0 范数来逼近原问题,重构效果理想,但复杂度高. 贪婪方法在冗余字典中迭代寻找与当前残差最匹配的 原子,近似估计原信号,复杂度较低,但重构效果一般, 典型 的 算 法 有 正 交 匹 配 追 踪 ( orthogonal matching pursuit, OMP ) [3] 和 子 空 间 追 踪 ( subspace pursuit, SP) [4] . 与这两类方法不同,贝叶斯重构利用了信号的 统计先验信息,使用概率推断和软判决方法估计原信 号,获得更好的重构质量,但复杂度上升,典型的算法 有 贝 叶 斯 压 缩 感 知 ( bayesian compressive sensing, BCS) [5]和随机正交匹配追踪(Rand鄄鄄OMP) [6] . 贝叶斯 压缩感知使用层次先验模型来描述重构信号系数的稀 疏性,通过期望最大化( expectation maximization,EM) 方法来估计模型的参数,并由贝叶斯公式求解问题的 最大后验概率(maximum a posteriori,MAP)估计,从而 实现信号重构. Rand鄄鄄 OMP 则根据 MAP 评估式按一 定概率随机确定贪婪搜索的方向估计信号,并由多组 OMP 估计值平均构成重构问题的最小均方误差(mini鄄 mum mean square error,MMSE)估计,并证明了当评判 准则为最小误差而非最稀疏时,MMSE 估计比单个 MAP 估计具有更小的重构误差. 传统的贝叶斯重构算法假设稀疏系数彼此独立, 未利用系数间的结构特性,近来研究发现利用结构化 稀疏信息可以加强重构问题的规范化约束,缩小重构 信号需要搜索的解空间范围,因此可以更加准确有效 地重构信号,例如:分块正交匹配追踪(Block鄄OMP) [7] 利用信号的块稀疏结构,假设信号非零系数成块出现; 基于模型的压缩感知(model鄄based CS,MB_CS) [8]利用 信号的树稀疏结构,假设信号中的系数呈现树形连接 的特征. 但这两者将系数间的结构看作是确定关系, 而实际信号中的结构用统计分布来描述更为可靠,如 基于分块结构的稀疏贝叶斯学习( block sparse bayes鄄 ian learning,BSBL) [9鄄鄄10]和基于小波树结构的压缩感知 ( tree鄄structured wavelet compressive sensing, TSW鄄鄄 CS) [11] . 使用玻尔兹曼机(Boltzmann machine,BM)模型表 示稀疏系数两两之间的统计依赖性[12鄄鄄15] ,可描述系数 间的任意结构特性,比块稀疏或树稀疏结构适用性更 广. 文献[12]首先通过图像训练库学习 BM 模型的参 数:系数的方差、系数间的相关矩阵和其偏差以及冗余 字典;然后根据这些参数使用吉比斯采样和模拟退火 算法求解关于稀疏模式的 MAP 估计,进而实现信号重 构,但该方法收敛速度慢. 为此,结构化贝叶斯追踪 (structured soft bayesian pursuit,SSoBaP) [13] 和基于玻 尔兹曼机的贝叶斯匹配追踪( boltzmann machine bayes鄄 ian matching pursuit,BM鄄鄄 BMP) [14] 研究了更快速的重 构算法. SSoBaP 采用基于平均场近似的变分贝叶斯方 法求解 MAP 估计,通过迭代的方式构造易求解极值的 概率函数最大化地近似 MAP 函数,得到 MAP 估计的 一个局部近似. BM鄄鄄BMP 则采用了更加快速的贪婪算 法进行 MAP 估计求解,在相同的计算复杂度下获得比 模拟退火算法更好的恢复性能. BM鄄鄄BMP 算法基于贝 叶斯匹配追踪的思想,每次迭代依次激活一个稀疏模 式的元素使得当前 MAP 评估值最大,通过若干次迭代 后寻找出一个稀疏模式的最优组合. 另外,BM鄄鄄 BMP 算法还采用块坐标优化方法实现参数估计和信号重构 的同步求解,比预训练的模型参数更匹配信号内容. 综上所述,提高算法的运算速度和重构质量是贝 叶斯重构方法的两个关键问题. 本文在 BM鄄鄄 BMP 的 基础上,利用玻尔兹曼机模型来描述稀疏系数间结构 信息,并通过增量计算和 MMSE 近似估计的方法分别 改善 BM鄄鄄BMP 算法的运算速度和重构质量. 第一,将 MAP 估计的评估值表示成上一次迭代的评估值加上 增值部分,以增量计算替代评估值的计算,从而减少了 运行时间;第二,考虑到自然图像的非稀疏性,在噪声 干扰下,单个 MAP 估计解往往不稳定,为此,本文利用 多个显著 MAP 估计值近似 MMSE 估计,克服了 MMSE 估计需遍历所有解的问题,也改善了 BM鄄鄄BMP 算法的 重构质量. 1 基于玻尔兹曼机的稀疏重构算法 1郾 1 玻尔兹曼机模型及 MAP 估计 在稀疏信号重构中,信号本身或者它的变换系数 常具有稀疏的先验分布,此外,系数间还会呈现出聚集 的结构特性,如小波树结构和块结构. 但对于更复杂 的信号如非平稳的自然图像而言,变换系数不满足树 稀疏或块稀疏的特性. 为了描述系数间的任意结构特 性,本文利用玻尔兹曼分布对稀疏系数间的统计依赖 性进行先验描述. (1)玻尔兹曼机模型. BM 分布定义如下: P(S) = 1 Z exp ( b T S + 1 2 S TWS ). (1) 其中,稀疏模式 S 是 m 维的二进制序列,第 i 个元素 ·1255·
·1256· 工程科学学报，第39卷，第8期 S,∈{-1,1}.当S=-1时，表示稀疏系数x,=0,字 其中，Σ，为一个k×k对角矩阵，k是系数支撑s的势， 典A中该系数对应的原子没有被选中.反之，当S,=1 即k是s元素的个数，角元素为（Σ三，）=0 时，表示系数x为非零值，该系数对应的原子被选中. 将系数支撑s的元素作为字典A的列索引，选出 W是沿对角线对称的相关矩阵，用来描述稀疏系数两 的列矢量构成一个矩阵A,·根据噪声服从高斯分布， 两间的依赖性，W表示S和S,之间的统计依赖程度 可以得到如下的高斯似然函数： 偏差b的元素b对应系数x,的统计特性，Z为起归一 化作用的配分函数，b和S分别表示b和S的转置. PCyls.). 1 玻尔兹曼机的图模型表示如图1所示，可以看到由于 (5) 玻尔兹曼机模型的全连接性，可通过相关度矩阵W刻 将(3)~(5)式代入(2)式，可求得MAP估计为： 画系数间的任意结构.W和b具有以下的规律：(a) w=ag晋-之Dy-nda(）+ 当W,>0,值越大表示S,和S,的取值越趋于相同，即对 应的系数x,和x,同时为零或非零值的可能性越大；(b) 2S"Ws +b's. (6) 当b,<0,绝对值越大表示S,=-1(对应系数x:=0)的 其中， 可能性越大.因此，当BM模型中偏差b的元素大部 Φ=A,∑AI+σI (7) 分是值很大的负数，且相关矩阵W的元素很少是值很 由公式(6)，按照文献[14]的方法去除了常数部 大的正数，系数为非零值的概率将会很小，从而保证其 分后可得到第k次迭代MAP估计的评估值Val(s): 稀疏性 Va()'-da()+2Ws+2h. (8) 其中，W为W第i列的转置 当得到估计值S后，系数x的MAP估计值可通过 以下公式计算)： cup=arg maxP(xly,S)≈E[xly,S]=AT中；'y. (9) 其中，E[·]表示求随机变量的数学期望，：为一个 m×m对角矩阵，当j∈s时，对角元素为()w=o5, 图1玻尔兹曼机的图模型表示 否则对角元素为零 Fig.1 Graphical model of the Boltzmann machine 1.2基于玻尔兹曼机的快速贝叶斯匹配追踪算法 (2)MAP估计.假设s是稀疏系数x的支撑，s的 在求解MAP估计时，可采用吉比斯采样和模拟退 元素为非零系数的索引，那么稀疏系数x的非零元素 火算法[]或变分贝叶斯方法)，但这些方法的计算 可以用x,表示，每个非零系数x,服从均值为零方差为 复杂度过高，BM-BMP算法)则利用贪婪搜索的思路 σ的高斯分布.由贝叶斯公式可知系数支撑s的 迭代寻找MAP估计的局部最优解.在每次迭代中， MAP估计值计算式为： BM-BMP算法依次激活稀疏模式S的一个元素（即令 5wr =arg maP(sly)=arg mP(yls)P(s).(2) S由-1变为1)，选出使当前MAP评估值，即Val值 最大的索引i,并将其放入系数支撑集中，即s=s-U 其中，集合2包含了2"个s可能的解，当j∈s时，S,= {i.}.为了进一步降低计算复杂度，本文将Val(s)值 1,否则S=-1. 给定系数支撑s,信号y的条件分布： 分解成上一次迭代的评估值Vl(s-1)和增值△：，则每 次迭代仅需计算增值△，加快算法运算速度.推导过 P(yls)P(yl)P(z.ls) (3) 程如下： 由于稀疏系数x可以看成是玻尔兹曼分布的支撑 (1)按照公式(8)计算Vl(s)值，需先计算中 对高斯分布的采样，则给定支撑s,稀疏系数的非零元 和中.现将第k次迭代的中，分解成上次迭代的值 素x,的条件分布为： Φ加上增值部分，由公式(7)有： 1 重=中-+0a,a (10) 其中a,是A的第i列，由矩阵求逆引理：
工程科学学报,第 39 卷,第 8 期 Si沂{ - 1,1}. 当 Si = - 1 时,表示稀疏系数 xi = 0,字 典 A 中该系数对应的原子没有被选中. 反之,当 Si = 1 时,表示系数 xi为非零值,该系数对应的原子被选中. W 是沿对角线对称的相关矩阵,用来描述稀疏系数两 两间的依赖性,Wij表示 Si和 Sj之间的统计依赖程度. 偏差 b 的元素 bi对应系数 xi的统计特性,Z 为起归一 化作用的配分函数,b T和 S T分别表示 b 和 S 的转置. 玻尔兹曼机的图模型表示如图 1 所示,可以看到由于 玻尔兹曼机模型的全连接性,可通过相关度矩阵 W 刻 画系数间的任意结构. W 和 b 具有以下的规律:( a) 当 Wij > 0,值越大表示 Si和 Sj的取值越趋于相同,即对 应的系数 xi和 xj同时为零或非零值的可能性越大;(b) 当 bi < 0,绝对值越大表示 Si = - 1(对应系数 xi = 0)的 可能性越大. 因此,当 BM 模型中偏差 b 的元素大部 分是值很大的负数,且相关矩阵 W 的元素很少是值很 大的正数,系数为非零值的概率将会很小,从而保证其 稀疏性. 图 1 玻尔兹曼机的图模型表示 Fig. 1 Graphical model of the Boltzmann machine (2)MAP 估计. 假设 s 是稀疏系数 x 的支撑,s 的 元素为非零系数的索引,那么稀疏系数 x 的非零元素 可以用 xs表示,每个非零系数 xi服从均值为零方差为 滓 2 x,i的高斯分布. 由贝叶斯公式可知系数支撑 s 的 MAP 估计值计算式为: ^sMAP = arg max s沂赘 P(s| y) = arg max s沂赘 P(y |s)P(s). (2) 其中,集合 赘 包含了 2 m个 s 可能的解,当 j沂s 时,Sj = 1,否则 Sj = - 1. 给定系数支撑 s,信号 y 的条件分布: P(y |s) = 乙 xs沂Rk P(y | xs,s) P(xs |s)dxs . (3) 由于稀疏系数 x 可以看成是玻尔兹曼分布的支撑 对高斯分布的采样,则给定支撑 s,稀疏系数的非零元 素 xs的条件分布为: P(xs |s) = 1 det(2仔撞s) 1 / 2 exp { - 1 2 x T s 撞 - 1 s xs } . (4) 其中,撞s为一个 k 伊 k 对角矩阵,k 是系数支撑 s 的势, 即 k 是 s 元素的个数,角元素为(撞s) i,i = 滓 2 x,si . 将系数支撑 s 的元素作为字典 A 的列索引,选出 的列矢量构成一个矩阵 As . 根据噪声服从高斯分布, 可以得到如下的高斯似然函数: P(y | xs,s) = 1 (2仔滓 2 ) n / 2 exp { - 1 2滓 2椰y - Asxs椰2 2 } . (5) 将(3) ~ (5)式代入(2)式,可求得 MAP 估计为: ^sMAP = arg max s沂赘 - 1 2 y T椎 - 1 s y - 1 2 ln det (椎s) + 1 2 S TWS + b T S. (6) 其中, 椎s = As撞sA T s + 滓 2 In . (7) 由公式(6),按照文献[14]的方法去除了常数部 分后可得到第 k 次迭代 MAP 估计的评估值 Val(s k ): Val(s k ) = - 1 2 y T椎 - 1 s k y - 1 2 ln det(椎s k) + 2W T i S k + 2bi . (8) 其中,W T i 为 W 第 i 列的转置. 当得到估计值 ^s 后,系数 x 的 MAP 估计值可通过 以下公式计算[14] : x^ MAP = arg max x P(x | y,^s)抑E[x | y,^s] = 撞忆^sA T椎 - 1 ^s y. (9) 其中,E[·] 表示求随机变量的数学期望,撞忆^s 为一个 m 伊 m对角矩阵,当 j沂^s 时,对角元素为(撞忆^s ) j,j = 滓 2 x,^sj , 否则对角元素为零. 1郾 2 基于玻尔兹曼机的快速贝叶斯匹配追踪算法 在求解 MAP 估计时,可采用吉比斯采样和模拟退 火算法[12]或变分贝叶斯方法[13] ,但这些方法的计算 复杂度过高,BM鄄鄄BMP 算法[14]则利用贪婪搜索的思路 迭代寻找 MAP 估计的局部最优解. 在每次迭代中, BM鄄鄄BMP 算法依次激活稀疏模式 S k的一个元素(即令 S k i 由 - 1 变为 1),选出使当前 MAP 评估值,即 Val 值 最大的索引 i*,并将其放入系数支撑集中,即 s k = s k -1胰 {i* }. 为了进一步降低计算复杂度,本文将 Val(s k )值 分解成上一次迭代的评估值 Val(s k - 1 )和增值 驻i,则每 次迭代仅需计算增值 驻i,加快算法运算速度. 推导过 程如下: (1) 按照公式(8) 计算 Val( s k ) 值,需先计算 椎s k 和 椎 - 1 s k . 现将第 k 次迭代的 椎s k分解成上次迭代的值 椎s - 1加上增值部分,由公式(7)有: 椎s k = 椎s k - 1 + 滓 2 x,iaia T i . (10) 其中 ai是 A 的第 i 列,由矩阵求逆引理: ·1256·
刘玲君等：MMSE准则下基于玻尔兹曼机的快速重构算法 ·1257· 中=中+o2β.c,c (11) 计值：MP,而后验概率值P(sIy)可看作是MAP估计 其中， 值的权值，因此MMSE估计可由集合2中的所有2"个 e:=pia (12) MAP估计值加权构成.由于无法获得所有2个s的 B=(1+oac). (13) 解，本文选用重构过程中后验概率值显著的MAP估计 (2)进一步化简c的矩阵运算，将c表示成前k- 值构成MMSE的近似估计值.设第k次迭代的评估式 1次迭代值的累加.由公式(11)以及中，=σ21。，其中 为Vl(s),对应了一个后验概率值P(sIy),该值在 L.表示n×n的单位矩阵，可得： 一定程度上反映了该次迭代重构质量，计算式如下： Φ=H-o∑B0c0cm. P(sly)= exp Val(s") (14) (18) 丁 ∑exp Val(s)} 1 4-1 (15) 如果某次迭代的Val值大于所有迭代Val的最大 -1 值max_Val减去阈值T,则称该次迭代的MAP估计是 (3)由公式(10)~(15)，MAP评估值Val(s)可 显著的，该次迭代的序号将会选入集合⊙中： 分解成上次迭代的评估值Val(s-1)和增值△，两部 =Val(s)max_Val -T,kE[1,H.(19) 分，每次迭代只需计算△，： 其中，H是最大迭代次数.阈值T的取值越小系数x Val(s)=Val(s)+=Val(s)++ 越稀疏，重构出的图像越平滑：T的取值越大系数x越 不稀疏，重构出的图像包含细节和噪声越多.综合考 9cP+2ws+24 (16) 虑两者之间的均衡关系，在实验中该阈值由经验值设 因此，基于玻尔兹曼机的快速贝叶斯匹配追踪 定.由于本文算法采用一次重构过程中多个显著MAP (BM-FBMP)算法的重构流程只需将BM-BMP算法重 估计加权构成MMSE近似估计，比Rand-OMP多次重 构流程中计算MAP评估值的公式(8)换成公式(16)， 构计算MMSE估计节省了时间，保持了运算速度快的 转换后每次迭代MAP估计的计算复杂度从O(m2n2) 优势 降到了O(m).计算复杂度分析如下：在原算法的公 1.4BM-FBMP算法重构流程 式(8)中，第一项的计算复杂度是O(m2),第二项的 BM-FBMP算法采用自适应的信号重构方案[ 计算复杂度是0(n),第三项的计算复杂度是O(n), 该方案无需提前训练模型的参数，可以通过对信号数 第四项的计算复杂度是O(1),其中第一项的计算复 据学习得到参数，迭代求解信号的估计值和参数估计 杂度最高(n<m),同时由于每个原子（共有m个原 值（如图2所示），分成两个步骤：(1)已知参数，恢复 子)都需要计算Val值，所以原算法每次迭代MAP估 信号.根据上次迭代的参数A={W,b,{o}},使 计的计算复杂度为0(m2n2).在新算法的公式(16) 用BM-FBMP算法求本次迭代系数支撑s的MAP估 中，第一项的计算复杂度是0(0)，第二项的计算复杂 计值，进而估计稀疏系数x:(2)已知信号，估计参数. 度是0(n),第三项的计算复杂度是O(n),第四项的 根据上次迭代的s和x估计值，使用最大似然估计 计算复杂度是0(n),第五项的计算复杂度是O(1), (maximum likelihood,ML)求解本次迭代的方差 式中计算复杂度最高为O(n),由于每个原子（共有m {o,:}i,最大伪似然估计(maximum pseudo- 个原子)都需要计算Vl值，所以新算法每次迭代 likelihood,MPL)求参数W和b.图2中，h表示第h次 MAP估计的计算复杂度为O(mn). 循环，ML和MPL估计的具体细节可参照文献[14]. 1.3MMSE的近似估计 本文主要针对自然图像恢复问题.由于自然图像 含唤 住uws广 重构 信号 BM_FBMP 收敛？ 信号 的稀疏度较低，幅值较小的稀疏系数与噪声混叠在一 起，在噪声干扰下，单个MAP估计往往不够稳定.而 ML和MPL估计 在文献[6]中表明，在含噪恢复问题中，MMSE估计比 单个MAP估计具有更小的重构误差，即更高的峰值信 图2基于玻尔兹曼机的快速贝叶斯匹配追踪算法流程图 噪比(peak signal to noise ratio,PSNR)值，为此本文采 Fig.2 Flowchart of the BM-FBMP algorithm 用MMSE估计.MMSE估计的计算式如下： BM-FBMP算法已知模型参数，恢复信号，是基于 E=>P(sly)E[xly,s]. (17) 匹配追踪的思想，用贪婪算法进行MMSE估计近似求 由公式(9)可知，E[xly,s]约等于系数的MAP估 解.重构以空的系数支撑集s°，=☑开始，即m×1的
刘玲君等: MMSE 准则下基于玻尔兹曼机的快速重构算法 椎 - 1 s k = 椎 - 1 s k - 1 + 滓 2 x,i茁i ci c T i . (11) 其中, ci = 椎 - 1 s k - 1 ai . (12) 茁i = (1 + 滓 2 x,ia T i ci) - 1 . (13) (2) 进一步化简 ci的矩阵运算,将 ci表示成前 k - 1 次迭代值的累加. 由公式(11)以及 椎s 0 = 滓 2 In ,其中 In 表示 n 伊 n 的单位矩阵,可得: 椎 - 1 s k = 1 滓 2 In - 滓 2 x,i移 k-1 j = 1 茁 (j) c (j) c (j)T . (14) ci = 1 滓 2 ai - 滓 2 x,i移 k-1 j = 1 茁 (j) c (j) c (j)T ai . (15) (3) 由公式(10) ~ (15),MAP 评估值 Val( s k )可 分解成上次迭代的评估值 Val( s k - 1 ) 和增值 驻i 两部 分,每次迭代只需计算 驻i: Val(s k ) = Val(s k - 1 ) + 驻i = Val(s k - 1 ) + 1 2 ln茁i + 滓 2 x,i 2 茁i(y T ci) 2 + 2W T i S k + 2bi . (16) 因此,基于玻尔兹曼机的快速贝叶斯匹配追踪 (BM鄄鄄FBMP)算法的重构流程只需将 BM鄄鄄BMP 算法重 构流程中计算 MAP 评估值的公式(8)换成公式(16), 转换后每次迭代 MAP 估计的计算复杂度从 O(m 2 n 2 ) 降到了 O(mn). 计算复杂度分析如下:在原算法的公 式(8)中,第一项的计算复杂度是 O(mn 2 ),第二项的 计算复杂度是 O(n 3 ),第三项的计算复杂度是 O( n), 第四项的计算复杂度是 O(1),其中第一项的计算复 杂度最高( n < m),同时由于每个原子(共有 m 个原 子)都需要计算 Val 值,所以原算法每次迭代 MAP 估 计的计算复杂度为 O(m 2 n 2 ). 在新算法的公式(16) 中,第一项的计算复杂度是 O(0) ,第二项的计算复杂 度是 O(n),第三项的计算复杂度是 O( n),第四项的 计算复杂度是 O( n),第五项的计算复杂度是 O(1), 式中计算复杂度最高为 O( n),由于每个原子(共有 m 个原子) 都需要计算 Val 值,所以新算法每次迭代 MAP 估计的计算复杂度为 O(mn). 1郾 3 MMSE 的近似估计 本文主要针对自然图像恢复问题. 由于自然图像 的稀疏度较低,幅值较小的稀疏系数与噪声混叠在一 起,在噪声干扰下,单个 MAP 估计往往不够稳定. 而 在文献[6]中表明,在含噪恢复问题中,MMSE 估计比 单个 MAP 估计具有更小的重构误差,即更高的峰值信 噪比(peak signal to noise ratio, PSNR)值,为此本文采 用 MMSE 估计. MMSE 估计的计算式如下: x^ MMSE = 移s沂赘 P(s| y)E[x | y,s]. (17) 由公式(9)可知,E[x | y,s]约等于系数的 MAP 估 计值 x^ MAP ,而后验概率值 P( s | y) 可看作是 MAP 估计 值的权值,因此 MMSE 估计可由集合 赘 中的所有 2 m个 MAP 估计值加权构成. 由于无法获得所有 2 m个 s 的 解,本文选用重构过程中后验概率值显著的 MAP 估计 值构成 MMSE 的近似估计值. 设第 k 次迭代的评估式 为 Val(s k ),对应了一个后验概率值 P( s k | y),该值在 一定程度上反映了该次迭代重构质量,计算式如下: P(s k | y) = exp{Val(s k )} 移r沂专 exp{Val(s r )} . (18) 如果某次迭代的 Val 值大于所有迭代 Val 的最大 值 max_Val 减去阈值 T,则称该次迭代的 MAP 估计是 显著的,该次迭代的序号将会选入集合 专 中: 专 = {k | Val(s k ) > max_Val - T,k沂[1,H]}. (19) 其中,H 是最大迭代次数. 阈值 T 的取值越小系数 x 越稀疏,重构出的图像越平滑;T 的取值越大系数 x 越 不稀疏,重构出的图像包含细节和噪声越多. 综合考 虑两者之间的均衡关系,在实验中该阈值由经验值设 定. 由于本文算法采用一次重构过程中多个显著 MAP 估计加权构成 MMSE 近似估计,比 Rand鄄鄄OMP 多次重 构计算 MMSE 估计节省了时间,保持了运算速度快的 优势. 1郾 4 BM鄄鄄FBMP 算法重构流程 BM鄄鄄FBMP 算法采用自适应的信号重构方案[14] . 该方案无需提前训练模型的参数,可以通过对信号数 据学习得到参数,迭代求解信号的估计值和参数估计 值(如图 2 所示),分成两个步骤:(1) 已知参数,恢复 信号. 根据上次迭代的参数 撰 = {W,b,{滓 2 x,i } m i = 1 },使 用 BM鄄鄄FBMP 算法求本次迭代系数支撑 s 的 MAP 估 计值,进而估计稀疏系数 x;(2) 已知信号,估计参数. 根据上次迭代的 s 和 x 估计值,使用最大似然估计 (maximum likelihood, ML ) 求 解 本 次 迭 代 的 方 差 {滓 2 x,i } m i = 1 , 最 大 伪 似 然 估 计 ( maximum pseudo鄄鄄 likelihood,MPL)求参数 W 和 b. 图 2 中,h 表示第 h 次 循环,ML 和 MPL 估计的具体细节可参照文献[14]. 图 2 基于玻尔兹曼机的快速贝叶斯匹配追踪算法流程图 Fig. 2 Flowchart of the BM鄄鄄FBMP algorithm BM鄄鄄FBMP 算法已知模型参数,恢复信号,是基于 匹配追踪的思想,用贪婪算法进行 MMSE 估计近似求 解. 重构以空的系数支撑集 s 0 * = 芰开始,即 m 伊 1 的 ·1257·
·1258· 工程科学学报，第39卷，第8期 稀疏模式矢量的元素全为-1(S°，=-1x1).首先，每 的系数方差o,都设为1002，BMP系数为非零值的伯 次将向量S中的一个元素取值由-1变为1，即从m 努利概率p=10/256,而本文算法和BM-BMP的玻尔 个元素中依次取一个到支撑集，得到m个互不相同的 兹曼机模型参数初始值设为W=0,b:=0.5ln(p/(1- 向量S.s=s”UHi计表示取向量S中第i个元素为支 p)),显著MAP选择阈值T=ln(10),最大迭代次数 撑，其中i=1,2,…,m.如果ies,则S=1,意味着系 H=n/2.对比实验包括：(1)不同噪声强度下的重构质 数x为非零值，该系数对应的原子被选中，否则S,= 量；(2)算法运行时间：(3)以Barbara为例的主观质量 -1.根据(16)式计算m个互不相同的向量S对应的 评价. 评估值Vl,如果i,所对应的评估值最大，则令支撑集 s。=s°U{i,}.在更新后的系数支撑的基础上，按照 上述方法，从向量S剩下的元素中挑选使后验概率最 大的一个元素，并将其加入支撑集.使用(19)式计算 显著的MAP估计并使用(I7)式计算EMs:BM-FBMP 重构算法的伪代码如下. Input:y,c2,A和模型参数W,b,{oi- Output:E- s9=☑，S9=-1mx1 k=1. repeat fori s=sUlil: 图320幅测试图像 1,j=i, S[U]= Fig.3 20 test images S[j门，j≠i: (16)式计算Val(s). ☆一OMP end for 日—BMP i.=arg max Val(s) BM-BMP 本文算法 1, j=i., s.=sU{i.},S[j]= S-[Uj门，j≠i,； k=k+1. until k>H或后验概率减小（即P(s1y)<P(s-11 y)). 101520 2530 (19)式计算显著的MAP估计. 噪声强度 (17)式计算ns 图4不同噪声强度下各种算法的PSNR平均值 Return XMMsE Fig.4 Average PSNR of the various algorithms at different noise levels 2实验结果与分析 2.1客观质量比较 为了验证算法的有效性，将本文算法与正交匹配 实验选用5种噪声标准差σ∈{5,10,15,20,25}， 追踪(OMP))、子空间追踪(SP)、贝叶斯匹配追踪 按照文献[14]的流程，所有测试图像分别添加各种强 (BMP)[、基于玻尔兹曼机的贝叶斯匹配追踪 度的噪声得到含噪图像，然后分别用上述各种算法重 (BM-BMP)[4种算法进行比较.实验测试图像采用 构出稀疏系数：，进而恢复出无噪图像A.实验中使 如图3所示的20幅512×512的自然图像，每幅图像 用PSNR作为评价标准，经过100次重复实验，图4给 被划分为不重叠的8×8分块并按列扫描成列矢量. 出了不同噪声强度下各种算法对所有图像重构结果的 实验使用的字典为n×m=64×256的冗余离散余弦 PSNR平均值曲线，其值越高越好.从图4可以看出： 变换(discrete cosine transform,DCT)字典，贝叶斯方法 (1)本文算法的PSNR曲线最高，从图中曲线的对比
工程科学学报,第 39 卷,第 8 期 稀疏模式矢量的元素全为 - 1(S 0 * = - 1 m 伊 1 ). 首先,每 次将向量 S 中的一个元素取值由 - 1 变为 1,即从 m 个元素中依次取一个到支撑集,得到 m 个互不相同的 向量 S. s 1 = s 0 * 胰{i}表示取向量 S 中第 i 个元素为支 撑,其中 i = 1,2,…,m. 如果 i沂s 1 ,则 Si = 1,意味着系 数 xi为非零值,该系数对应的原子被选中,否则Si = - 1. 根据(16)式计算 m 个互不相同的向量 S 对应的 评估值 Val,如果 i* 所对应的评估值最大,则令支撑集 s 1 * = s 0 * 胰{i* }. 在更新后的系数支撑的基础上,按照 上述方法,从向量 S 剩下的元素中挑选使后验概率最 大的一个元素,并将其加入支撑集. 使用(19)式计算 显著的 MAP 估计并使用(17)式计算 x^ MMSE . BM鄄鄄FBMP 重构算法的伪代码如下. Input: y,滓 2 ,A 和模型参数 W,b,{滓 2 x,i} m i = 1 . Output: x^ MMSE . s 0 * = 芰, S 0 * = - 1 m 伊 1 . k = 1. repeat for i埸s k - 1 * , s k = s k - 1 * 胰{i}; S k [j] = 1, j = i, S k - 1 { * [j], j屹i; (16)式计算 Val(s k ). end for i* = arg max i {Val(s k )}; s k * = s k - 1 * 胰{i* },S k * [j] = 1, j = i* , S k - 1 * [j], j屹i* { ; k = k + 1. until k > H 或后验概率减小(即 P( s k | y) < P( s k - 1 | y)). (19)式计算显著的 MAP 估计. (17)式计算 x^ MMSE . Return:x^ MMSE . 2 实验结果与分析 为了验证算法的有效性,将本文算法与正交匹配 追踪(OMP) [3] 、子空间追踪( SP) [4] 、贝叶斯匹配追踪 ( BMP ) [14] 、 基 于 玻 尔 兹 曼 机 的 贝 叶 斯 匹 配 追 踪 (BM鄄鄄BMP) [14] 4 种算法进行比较. 实验测试图像采用 如图 3 所示的 20 幅 512 伊 512 的自然图像,每幅图像 被划分为不重叠的 8 伊 8 分块并按列扫描成列矢量. 实验使用的字典为 n 伊 m = 64 伊 256 的冗余离散余弦 变换(discrete cosine transform, DCT)字典,贝叶斯方法 的系数方差 滓 2 x,i都设为 100 2 ,BMP 系数为非零值的伯 努利概率 p = 10 / 256,而本文算法和 BM鄄鄄 BMP 的玻尔 兹曼机模型参数初始值设为 W = 0,bi = 0郾 5ln( p / (1 - p)),显著 MAP 选择阈值 T = ln(10 4 ),最大迭代次数 H = n / 2. 对比实验包括:(1)不同噪声强度下的重构质 量;(2)算法运行时间;(3)以 Barbara 为例的主观质量 评价. 图 3 20 幅测试图像 Fig. 3 20 test images 图 4 不同噪声强度下各种算法的 PSNR 平均值 Fig. 4 Average PSNR of the various algorithms at different noise levels 2郾 1 客观质量比较 实验选用 5 种噪声标准差 滓沂{5,10,15,20,25}, 按照文献[14]的流程,所有测试图像分别添加各种强 度的噪声得到含噪图像,然后分别用上述各种算法重 构出稀疏系数 x^,进而恢复出无噪图像 Ax^. 实验中使 用 PSNR 作为评价标准,经过 100 次重复实验,图 4 给 出了不同噪声强度下各种算法对所有图像重构结果的 PSNR 平均值曲线,其值越高越好. 从图 4 可以看出: (1) 本文算法的 PSNR 曲线最高,从图中曲线的对比 ·1258·