Split Bregman迭代介绍及MATLAB实现

Bregman距离

定义一个凸函数$E$中点$u$与点$v$的Bregman距离为
$$D_E^p\left(u,v\right)=E\left(u\right)-E\left(v\right)-<p,u-v>$$
其中，$p$为$E$在点$v$的次梯度。显然，可以得到类似于距离的性质
非负性：$D_E^p\left(u,v\right)\ge 0$
当$u\le w\le v$时，有 $D_E^p\left(u,v\right)\ge D_E^p\left(w,v\right)$

Bregman迭代

考虑如下约束的优化问题
$$\underset{u}{\min }E\left(u\right)\text{+}\lambda H\left(u\right)$$
其中，$E$，$H$均为凸函数，且 $\underset{u}{\min }H\left(u\right)=0$
可以迭代求解上述问题
$$\begin{array}{rl}& u^{k+1}=\underset{u}{\min }{D}_E^p\left(u,u^k\right)\text{+}\lambda H\left(u\right)\\ & =\underset{u}{\min }E\left(u\right)-<p^k,u-u^k>\text{+}\lambda H\left(u\right)\end{array}$$
假设$H$可微，且$0\in \partial \left(D_E^p\left(u,u^k\right)+\lambda H\left(u\right)\right)$表示在$u^{k+1}$处的次梯度，由于$p^{k+1}\in \partial E\left(u^{k+1}\right)$，那么可以得到
$$p^{k+1}=p^k-\nabla H\left(u^{k+1}\right)$$
事实上，在一般条件的约束下，随着$k$的增加，$H\left(u^k\right)$是趋近于0的。当然在$E$，$H$均为凸函数且$H$可微的条件下还能得到$H\left(u\right)$是单调递减且收敛的。
下面用该方法来解决实际的问题。对于以下约束问题
$$\{\begin{array}{rl}& \underset{u}{\min }E\left(u\right)\\ & s.t.Au=b\end{array}$$
可以将其转化为
$$\underset{u}{\min }E\left(u\right)\text{+}\frac{\lambda }{2}{\Vert Au-b\Vert }_2^2$$
利用Bregman迭代可以求解，迭代如下：
$$\begin{array}{rl}& u^{k+1}=\underset{u}{\min }{D}_E^p\left(u,u^k\right)\text{+}\frac{\lambda }{2}{\Vert Au-b\Vert }_2^2\\ & =\underset{u}{\min }E\left(u\right)-<p^k,u-u^k>\text{+}\frac{\lambda }{2}{\Vert Au-b\Vert }_2^2\end{array}$$
$$P^{k+1}=p^k-\lambda A^T\left(A{u}^{k+1}-b\right)$$
其中关于$P^{k+1}$的推导如下
$$\begin{array}{rl}& \nabla H\left(u^{k+1}\right)\text{=}\nabla \frac{\lambda }{2}{\Vert A{u}^{k+1}-b\Vert }_2^2\\ & =\nabla \frac{\lambda }{2}\left[{\left(A{u}^{k+1}-b\right)}^H\left(A{u}^{k+1}-b\right)\right]\\ & =\frac{\lambda }{2}\nabla tr\left[{\left(A{u}^{k+1}-b\right)}^H\left(A{u}^{k+1}-b\right)\right]\end{array}$$
而
$$\begin{array}{rl}& dtr\left[{\left(A{u}^{k+1}-b\right)}^H\left(A{u}^{k+1}-b\right)\right]=tr\left[{\left(Ad{u}^{k+1}\right)}^H\left(A{u}^{k+1}-b\right)+{\left(A{u}^{k+1}-b\right)}^{H}Ad{u}^{k+1}\right]\\ & =tr\left[2{\left(A{u}^{k+1}-b\right)}^{H}Ad{u}^{k+1}\right]\end{array}$$
因此
$$\nabla \left[{\left(A{u}^{k+1}-b\right)}^H\left(A{u}^{k+1}-b\right)\right]\text{=}2A^T\left(A{u}^{k+1}-b\right)$$
即$P^{k+1}=p^k-\lambda A^T\left(A{u}^{k+1}-b\right)$
当$A$为线性时，上述的迭代过程可以简化为
$$\{\begin{array}{rl}& u^{k+1}=\underset{u}{\min }E\left(u\right)+\frac{\lambda }{2}{\Vert Au-b^k\Vert }_2^2\\ & b^{k+1}=b^k+b-A{u}^{k+1}\end{array}$$
Bregman迭代有以下优点：
收敛速度比较快，尤其对于$E\left(u\right)$为$l_1$范数正则化。同时，相对于延拓法来说，Bregman迭代中的$\lambda$是一常数，也就说可以选择合适的$\lambda$的值使得迭代过程快速收敛。

Split Bregman迭代

考虑如下问题
$$\{\begin{array}{rl}& \underset{u,d}{\min }\left|d\right|+H\left(u\right)\\ & s.t.d=\Phi \left(u\right)\end{array}$$
虽然看似$d=\Phi \left(u\right)$这一步是多余的，但是在实际中却能起到很重要的作用。
为了解决该问题，可以将其转化为
$$\underset{u,d}{\min }\left|d\right|+H\left(u\right)\text{+}\frac{\lambda }{2}{\Vert d-\Phi \left(u\right)\Vert }_2^2$$
令$E\left(u,d\right)=\left|d\right|+H\left(u\right)$，$H\left(u,d\right)=d-\Phi \left(u\right)$，那么该表达式就变成了之前介绍的Bregman迭代的更普遍的形式，可以得到
$$\begin{array}{rl}& \left(u^{k+1},d^{k+1}\right)=\underset{u,d}{\min }{D}_E^p\left(u,u^k,d,d^k\right)+\frac{\lambda }{2}{\Vert d-\Phi \left(u\right)\Vert }_2^2\\ & =\underset{u,d}{\min }E\left(u,d\right)-<p_u^k,u-u^k>-<p_d^k,d-d^k>+\frac{\lambda }{2}{\Vert d-\Phi \left(u\right)\Vert }_2^2\end{array}$$
$$p_u^{k+1}=p_u^k-\lambda {\left(\nabla \Phi \right)}^T\left(\Phi u^{k+1}-d^{k+1}\right)$$
$$p_d^{k+1}=p_d^k-\lambda \left(d^{k+1}-\Phi u^{k+1}\right)$$
那么可以简化为
$$\left(u^{k+1},d^{k+1}\right)\text{=}\underset{u,d}{\min }\left|d\right|+H\left(u\right)\text{+}\frac{\lambda }{2}{\Vert d-\Phi \left(u\right)-b^k\Vert }_2^2$$
$$b^{k+1}=b^k+\left(\Phi \left(u^{k+1}\right)-d^{k+1}\right)$$
显然，为了实现上述迭代过程的计算，必须要解决
$$\left(u^{k+1},d^{k+1}\right)\text{=}\underset{u,d}{\min }\left|d\right|+H\left(u\right)\text{+}\frac{\lambda }{2}{\Vert d-\Phi \left(u\right)-b^k\Vert }_2^2$$
那么可以将其分为两步迭代（这也是为什么称其为Split Bregman迭代的原因）
第一步：
$$u^{k+1}=\underset{u}{\min }H\left(u\right)+\frac{\lambda }{2}{\Vert d^k-\Phi \left(u\right)-b^k\Vert }_2^2$$
第二步：
$$d^{k+1}=\underset{d}{\min }\left|d\right|+\frac{\lambda }{2}{\Vert d^k-\Phi \left(u^{k+1}\right)-b^k\Vert }_2^2$$
在第二步中，可以利用收缩算子来进行简化并且提高运算的速度
$$d^{k+1}=shrink\left(\Phi \left(u\right)+b^k,1/\lambda \right)$$
其中， $shrink\left(x,y\right)=\frac{x}{|x|}\ast \max \left(\left|x\right|-y,0\right)$
Split Bregman迭代的具体过程为
While ${\Vert u^k-u^{k-1}\Vert }_2^2>threshold$
for $n=1$ to $N$
$$u^{k+1}=\underset{u}{\min }H\left(u\right)+\frac{\lambda }{2}{\Vert d^k-\Phi \left(u\right)-b^k\Vert }_2^2$$
$$d^{k+1}=\underset{d}{\min }\left|d\right|+\frac{\lambda }{2}{\Vert d^k-\Phi \left(u^{k+1}\right)-b^k\Vert }_2^2$$
end
$b^{k+1}=b^k+\left(\Phi \left(u^{k+1}\right)-d^{k+1}\right)$
end
经验上$N=1$可以使内循环收敛

Split Bregman稀疏信号重构

考虑信号$f$是稀疏信号$u$经过线性算子$A$后得到的，那么可以通过利用$l_1$范数正则化从 $f$中恢复出信号$u$，即
$$\hat{u}=\arg \min {\Vert u\Vert }_1+\frac{\mu }{2}{\Vert Au-f\Vert }_2^2$$
令$\Phi =I$，$H\left(u\right)=\frac{\mu }{2}{\Vert Au-f\Vert }_2^2$，那么可以使用Split Bregman迭代进行求解，迭代如下：
$$u^{k+1}=\underset{u}{\min }\frac{\mu }{2}{\Vert Au-f\Vert }_2^2+\frac{\lambda }{2}{\Vert d^k-u-b^k\Vert }_2^2$$
$$d^{k+1}=\underset{d}{\min }{\Vert d\Vert }_1+\frac{\lambda }{2}{\Vert d^k-u^{k+1}-b^k\Vert }_2^2$$
$$b^{k+1}=b^k+u^{k+1}-d^{k+1}$$
为了解决$u^{k+1}$，通过对其求导，令倒数为0，可以得到
$$\partial \left(\frac{\mu }{2}{\Vert Au-f\Vert }_2^2+\frac{\lambda }{2}{\Vert d^k-u-b^k\Vert }_2^2\right)=0$$
即
$$\mu A^{\ast }\left(A{u}^{k+1}-f\right)-\lambda \left(d^k-u^{k+1}-b^k\right)=0$$
因此
$$u^{k+1}={\left(\mu A^{\ast }A+\lambda I\right)}^{-1}\left(\mu A^{\ast }f+\lambda \left(d^k-b^k\right)\right)$$
而$d^{k+1}$可以通过收缩算子进行求解
$$d^{k+1}=shrink\left(u^{k+1}+b^k,1/\lambda \right)$$
因此，利用Split Bregman进行稀疏信号重构的过程为
初始化：$d^0=b^0=0$
While ${\Vert u^k-u^{k-1}\Vert }_2^2>threshold$
$$u^{k+1}={\left(\mu A^{\ast }A+\lambda I\right)}^{-1}\left(\mu A^{\ast }f+\lambda \left(d^k-b^k\right)\right)$$
$$d^{k+1}=shrink\left(u^{k+1}+b^k,1/\lambda \right)$$
$$b^{k+1}=b^k+u^{k+1}-d^{k+1}$$
$$k=k+1$$
end

仿真参数

参数名称	参数值
稀疏度	12
信号长度	256
$\lambda$	1
$\mu$	10

仿真结果

从图中可以看出，经过Split Bregman重构后的信号基本与原始稀疏信号一致，达到了较好的效果。同时也可以通过重构信号与原始信号之间的误差来评价该方法的性能。
代码如下：
主函数

clear;
close all;
clc;
M = 256;            %观测值个数
N = 256;            %信号x的长度
K = 12;             %信号x的稀疏度
lambda=1;           %正则化参数
mu=10;
x = zeros(N,1);        
Index_K = randperm(N);    
x(Index_K(1:K)) = 5*randn(K,1);     %构造稀疏信号
A=randn(M,M);                            %观测矩阵
y=A*x;
x_re=L1_SplitBregmanIteration(y,A,mu,lambda);   %利用Split Bregman进行稀疏信号重构
figure(1);
plot(x_re,'k.-');
hold on;
plot(x,'r');
legend('重构信号','原始信号')

L1_SplitBregmanIteration.m

function u=L1_SplitBregmanIteration(f,A,mu,lambda)
%======================================
%
% L1 Split Bregman Iteration
% This function compute the solution
% u = arg min |u|+0.5*mu||Au-f||_2^2
%
% by using the Split Bregman Iteration
%
% In this version we consider only
% the case where A is a square matrix
%
% f: measured data
% A: some linear operator in its matrix
% form
% mu: regularization coefficient
% lambda: "spliting" regularization
% coefficient
% Author: Jerome Gilles
% Institution: UCLA - Math Department
% email: jegilles@math.ucla.edu
%
% Note: typically mu=10, lambda=1,
% eps=0.01 work well
%======================================
N=size(f,1);
d=zeros(N,1);
b=zeros(N,1);
u=zeros(N,1);
Z=zeros(N,1);
Ft=mu*A'*f;
IV=inv(mu*A'*A+lambda*eye(N));
up=ones(N,1);
while ((u-up)'*(u-up))>eps
    up=u; %store the previous iteration
    u=IV*(Ft+lambda*(d-b)); %update u
    tmp=u+b;
    d=sign(tmp).*max(Z,abs(tmp)-1/lambda); %update d
    b=tmp-d; %update b
end

参考文献
[1] Goldstein T, Osher S. The split Bregman method for L1-regularized problems[J]. SIAM journal on imaging sciences, 2009, 2(2): 323-343.
[2]Gilles J, Osher S. Bregman implementation of Meyer’s G-norm for cartoon+ textures decomposition[J]. UCLA Cam Report, 2011: 11-73.
[3] Bush J. Bregman algorithms[J]. Senior Thesis. University of California, Santa Barbara, 2011.
[4] Yin W, Osher S, Goldfarb D, et al. Bregman iterative algorithms for ${\ell }_1$-minimization with applications to compressed sensing[J]. SIAM Journal on Imaging sciences, 2008, 1(1): 143-168.