FISTA解决地震反演问题 [WIP]

• 迭代收缩阈值算法（Iterative shrinkage-thresholding algorithm, ISTA）是一种用于信号处理和图像重建的优化算法。

  本文将介绍ISTA算法原理和其处理地震反演问题的实际应用。

优化问题

梯度下降

• 对于一个线性变换问题$\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{b}$，其中$\mathbf{A}$和$\mathbf{y}$已知，$\mathbf{b}$为未知噪音，我们需要求解$\mathbf{x}$。

  本质上这就是一个线性回归问题，从线性回归可知我们可以使用梯度下降解决这个问题。

• 梯度下降会带来新的问题，对于无约束的优化问题

  $$\underset{x}{\min}\{F(x)\equiv f(x)\}.\tag{1}$$

  若实函数$F(x)$在点$a$处可微且有定义，梯度下降总是认为$F(x)$在点$a$沿着梯度相反的方向$-\nabla F(a)$下降最快。

  所以当$f(x)$连续可微时，若存在一个足够小的数值$t>0$使得

  $$x_2 = x_1 -t\nabla F(a).\tag{2}$$

  则有$F(x_1) \ge F(x_2)$。

  梯度下降核心便是通过式2找到序列$\{x_k\}$，使得$F(x_k)\ge F(x_{k-1})$。

  显然，此时初值的选取成了关键，即梯度下降可能陷入局部最优，同时$t_k$的选取（机器学习中的学习率）也是关键，太小会导致迭代太慢，太大会导致无法收敛。

ISTA

• ISTA也是基于梯度下降算法，它的关键在于对$\{x_k\}$的选取。对于式1，当其满足Lipschitz连续条件，即$f(x)$的导数有下界（最小下界称为Lipschitz常数$L(f)$）。此时对于任意$L\ge L(f)$，有

  $$f(\mathbf{x})\le f(\mathbf{y}) + \langle\mathbf{x}-\mathbf{y},\nabla F(\mathbf{y})\rangle + \frac{L}{2}\lVert \mathbf{x}-\mathbf{y}\rVert^2,\tag{3}$$

  对于所有$\mathbf{x,y}\in \mathbb{R}^n$成立。此时在点$\mathbf{x}_k$附近的函数值可近似表示为

  $$\hat{f}(\mathbf{x},\mathbf{x}_k) \doteq f(\mathbf{x}_k) + \langle \nabla f(\mathbf{x}_k), \mathbf{x}-\mathbf{x}_k\rangle + \frac{L}{2}\lVert \mathbf{x}-\mathbf{x}_k\rVert^2.\tag{4}$$

  此时再进行梯度下降时，将点$\mathbf{x}_{k-1}$处的近似函数取最小值的点作为下一次迭代起始点$\mathbf{x}_k$，这便是proximal regularization算法。

• 上述思路只能解决无约束问题，ISTA主要用于解决带惩罚项（噪声）的优化问题，例如对于带范数规范化函数$g(x)$的约束优化问题

  $$\min\{F(\mathbf{x}) \equiv f(\mathbf{x}) + g(\mathbf{x}) : \mathbf{x} \in \mathbb{R}^n\}. \tag{5}$$

  同样我们可以得到在点$y$处$F(\mathbf{x})$的二次近似函数$Q_L$为

  $$Q_L(\mathbf{x}, \mathbf{y}) := f(\mathbf{y}) + \langle\mathbf{x} - \mathbf{y}, \nabla f(\mathbf{y})\rangle + \frac{L}{2}\lVert\mathbf{x}-\mathbf{y}\rVert^2 + g(\mathbf{x}).\tag{6}$$

  忽略常数项$f(\mathbf{y}),\nabla F(\mathbf{y})$，该函数的最小值表示$p_L$可以写为

  $$p_L = \underset{\mathbf{x}}{\argmin}\left\{g(\mathbf{x}) + \frac{L}{2}\lVert \mathbf{x} - (\mathbf{y} - \frac{1}{L}\nabla f(\mathbf{y}))\rVert^2 \right\}. \tag{7}$$

  对于ISTA其迭代步骤即为$\mathbf{x}_k = p_L(\mathbf{x}_{k-1})$。

• ISTA的问题则是Lipschitz常数不一定可知或可计算，所以还衍生出了带回溯的ISTA算法

  

  图1. 带回溯的ISTA算法

  同时我们可知其收敛速度为$O(1/k)$。证明略。

FISTA

• FISTA与ISTA的区别在于迭代步骤中近似函数起始点$\mathbf{y}$的选择，ISTA使用前一次迭代求得的近似函数最小值点$\mathbf{x}_{k-1}$，而FISTA使用另一种方法进行计算

  

  图2. FISTA

• 同时由于Lipschitz常数计算复杂问题，FISTA同样有带回溯版本。

  

  图3. 带回溯的FISATA

  由理论证明FISTA的收敛速度可以达到$O(1/k^2)$。非常快。

地震反演

反射系数反演

• 这里我们以反射系数反演问题为例来演示ISTA的应用。

  地震道可以看为是地震子波和地下反射系数的褶积，即

  $$\mathbf{S} = \mathbf{r}*\mathbf{h} + \mathbf{n},\tag{8}$$

  其中$\mathbf{r}$为反射系数序列（由波阻抗进行计算），$\mathbf{h}$为地震子波构建的脉冲源，$\mathbf{n}$为随机噪音。

• 现在我们给定地震道$\mathbf{S}$和地震子波$\mathbf{h}$，需要求解对应的反射系数$\mathbf{r}$。

  本质上这就是一个反褶积问题，我们求解步骤如下：