计算成像问题数学基础
成像问题建模
成像问题被表述为一个线性逆问题，这个逆问题的前向过程表示为公式（1）
\[ y=Ax+\epsilon \]
逆问题的目标在于求解 \(x\)。在公式（1）中，\(y\) 是我们可以测量得到的测量值，\(x\) 是我们想要求解的内容，\(A\) 是前向算子，表示施加在原始图像中的物理过程，\(\epsilon\) 是前向过程中引入的噪声值。
成像问题的逆问题的一个重要特征是，前向算子 \(A\) 通常是一个 \(m<n\) 的矩阵，这意味 \(A\) 的逆矩阵不存在，也无法求得唯一的解析解。虽然线性代数中提供了伪逆 \(A^{\dagger}\) 来求解病态的 \(y=Ax\) 问题，但是由于在前向过程中噪声的引入，更常见的方法是使用最优化方法进行建模分析求解上面的问题。
根据《最优化：建模、算法与理论》书中内容1可以注意到，公式（1）的形式满足线性回归模型的形式（虽然定义域和值域不同，但是形式一致）。书中指出，当假设误差是高斯白噪声时，最小二乘法的解是线性回归模型的最大似然解。因此常见的作法是使用最小二乘法来给出病态的 \(y=Ax\) 问题的最佳近似解。最小二乘法最常见的做法是极小化误差的 \(\ell_{2}\) 范数的平方，即
\[ \min_{x\in\mathbb{R}^n}\quad\frac12\|Ax-y\|_2^2. \]
由于在上面这个最小二乘问题中，矩阵 \(A\) 的 \(m < n\) 时，最优解很可能不止一个，且并不是所有的解都是我们想要的，所以通常借助引入先验知识来选出我们需要的特定性质的解。解的先验知识需要对具体问题进行具体分析。常见的，如果我们希望解向量的长度尽可能的短（欧几里得长度尽量短），那么可以加入 \(\ell_{2}\) 范数作为正则项；如果我们希望解中非零元素尽可能多（解是稀疏的），那么可以引入 \(\ell_{1}\) 范数作为正则项。
接下来的分析，我们借助具体的例子进行分析。假设我们希望得到的图像（逆问题的解）是稀疏的，那么可以引入 \(\ell_{1}\) 范数作为正则项。\(\ell_{1}\) 范数定义为将矩阵中的全部元素绝对值相加。当最小化 \(\ell_{1}\) 范数时，可以认为是让尽可能多的元素为0，从而筛选出包含最多零的解。当公式（2）选择\(\ell_{1}\) 范数作为正则项时，此时的问题被称作 LASSO 问题，即
\[ LASSO:\;\;\;\; \min_{x\in\mathbb{R}^n}\quad\frac12\|Ax-b\|_2^2+\mu\|x\|_1 \]
LASSO 问题是一个典型的复合优化问题，可以被抽象地定义为
\[ \min_{x\in\mathbb{R}^n}\quad\psi(x)\xrightarrow{\mathrm{def}}f(x)+h(x) \]
其中 \(f(x)\) 为光滑的可微函数，\(h(x)\) 为不可微函数。对于光滑函数 \(f(x)\)，可以使用最常见且有效的梯度下降法来求解，迭代公式如公式（5）所示。
\[ x^{k+1}=x^k-\alpha_k\nabla f(x^k) \]
但是，对于不可微函数 \(h(x)\)，梯度下降法没有办法使用。因此需要引入新的算法来实现求解。第一个方法是近端梯度下降法（Proximal Gradient Descent）。
近端梯度下降法
类似于一般的梯度下降法，近端梯度下降法给出了迭代公式进行最优化求解，即
\[ y^{k}=x^k-t_k\nabla f(x^k)\\ x^{k+1}=\mathrm{prox}_{t_kh}(y^k)\\ \]
观察公式（6）可以注意到，近端梯度下降法的两步迭代中，第一步是一个一般的梯度下降法，第二步是使用了一个被称为近端算子的运算来实现。对于一个凸函数 \(h(x)\) ，可以将其近端算子定义为
\[ \mathrm{prox}_h(x)=\arg\min_{u\in\mathrm{dom}h}\left\{h(u)+\frac12\|u-x\|^2\right\}. \]
可以看到，邻近算子的目的是求解一个距 \(x\) 不算太远的点2，并使函数值 \(h(x)\) 也相对较小。可以证明的是，给出的定义中优化问题解的存在具有唯一性3。在 LASSO 问题中，为了更方便的求出近端算子，可以借助次梯度来求解。邻近算子与次梯度的关系如公式（8）所示。
\[ u=\mathrm{prox}_h(x)\Longleftrightarrow x-u\in\partial h(u). \]
次梯度被定义为在函数 \(f\) 在某一点 \(x\) 的任意邻域内，对于所有的 \(y\) 都满足 \(f(y) \geq f(x) + g^T (y - x)\) 的向量 \(g\)。这里给出 \(\ell_1\) 范数的次梯度求解。
已知 \(\ell_1\) 范数为
\[ \|x\|_1 = \sum_{i} |x_i| \]
可以注意到到其在 \(x_i=0\) 的时候不可导，此处次梯度为 \([-1, 1]\) 区间内的所有值，写作
\[ \partial\|x\|_1 = \begin{cases} \{1\},&x>0,\\ [-1,1],&x=0,\\ \{-1\},&x<0, \end{cases} \]
因此，当公式（4）中的 \(h(x)\) 为 \(\ell_1\) 范数时，可以根据公式（8）和公式（10）得到
\[ x-u\in t\partial\|u\|_1 = \begin{cases} \{t\},&u>0,\\ [-t,t],&u=0,\\ \{-t\},&u<0, \end{cases} \]
因此，当 \(x > t\) 时，\(u = x - t\)；当 \(x < -t\) 时，\(u = x + t\)；当 \(x \in [-t, t]\) 时，\(u = 0\)，即
\[ u = \mathrm{prox}_h(x) = \text{sign}(x) \max\{|x| - t, 0\} \]
将公式（12）带入到公式（6）中，可以给出LASSO问题的近端梯度下降法具体的迭代公式
\[ y^{k}=x^k-t_kA^\mathrm{T}(Ax^k-b)\\ x^{k+1}=\mathrm{prox}_h(y^k)=\operatorname{sign}(y^k)\max\{|y^k|-t_k\mu,0\} \]
FISTA 算法简介
一个非常简单的 FISTA 算法介绍：
从算法时间复杂度角度上考虑，近端梯度下降法属于 \(\mathcal{O}\left(\frac1k\right)\) 问题。FISTA 算法的目的在于降低算法的时间复杂度到 \(\mathcal{O}\left(\frac{1}{k^2}\right)\)。
实现方法是，先沿着前两步的计算方向计算一个新点，然后在该新点处做一步近似点梯度迭代，即
\[ y^k=\quad x^{k-1}+\frac{k-2}{k+1}(x^{k-1}-x^{k-2})\\ x^k=\quad\mathrm{prox}_{t_kh}(y^k-t_k\nabla f(y^k)) \]