thresholded-Landweber

动机

在研究MRI重建时，需要处理$||Ax-b||_F^2$这样的式子，但是A其实是个算子，包括下采样，NUFFT和coil矩阵，总之非常复杂，希望在处理后面约束项的时候不要A在里面，那么怎么处理呢？
基本思想来自于SPARSE SIGNAL RESTORATION和Learning the Sparsity Basis in Low-rank plus Sparse Model for Dynamic MRI Reconstruction这两篇论文。下面先介绍基础。

MM算法

我们想$\min\ J(x)$，但是直接优化并不容易。这时我们假想有这么一个序列{$x_k$},使得$J(x_{k+1})<J(x_k)$,那么怎么来找这个序列呢？MM方法就是先找个新函数作为$J(x)$的上界，然后最小化这个上界来得到$x_{k+1}$(好像EM算法也是这个思想)。参见下图。

Landwebber 迭代

先考察问题$J(x)=||y-Hx||_2^2$,如果直接求解会得到$x=(H'H)^{-1}H'y$,这里求逆会带来很多的问题，顺便说一句，貌似加入松弛变量可以解决不可逆的问题，不过即使如此，运算量也可能比较大。
那么就可以采用MM算法避免直接求解线性方程组。我们找到上界函数$G(x)=J(x)+non-negative function \ of\ x$；同时在x_{k}这个点上，G(x)=J(x)，因此一个通常的选择是：
$G_k(x)=||y-Hx||_2^2+(x-x_k)'(\alpha-H'H)(x-x_k).$
为了确保非负性，$\alpha>=maxeig(H'H)$.
然后我们最小化G(x)。发现$x=x_k+\frac{1}{\alpha}H'(y-Hx_k)$

ISTA–landweber的应用

这里就不介绍1范数的求解了和ISTA推导了，我们可以直接给出这样的过程：
对于求解$J(x)=||y-Hx||_2^2+\lambda||x||_1$这样的函数，注意到
$G_k(x)=\alpha||x_k+\frac{1}{\alpha}H'(y-Hx_k)-x||_2^2+K$ ,K是与x无关的量。
那么我们可以直接化成求解：
$||x_k+\frac{1}{\alpha}H'(y-Hx_k)-x||_2^2+\frac{\lambda}{\alpha}||x||_1$,采用softsholding,就是
$x_{k+1}=soft(x_k+\frac{1}{\alpha}H'(y-Hx_k),\frac{\lambda}{2\alpha})$