Kaczmarz迭代算法在磁粒子成像（MPI）中的应用，迭代公式推导

Kaczmarz迭代算法背景

Kaczmarz迭代算法作用：最初被用来求解线性方程组，在MPI中被用来重建图像。

假设存在一个线性方程组（此方程可以用来建模MPI中系统矩阵和接收信号之间的关系）：
$$Sc=u$$
其中$S\in C^{M\times N},c\in C^N,u\in C^M$。M为M个频率分量，N对应了MPI中体素的个数。

---

MPI背景： 在MPI中S即为系统矩阵，$c$为空间中粒子分布的浓度，$u$为接收到的电压信号。而MPI重建的就是粒子浓度的空间分布$c$，因此可以通过已知的$S，u$求解$c$。

---

Kaczmarz迭代算法核心公式

​ 和其他文章一样，这里首先给出Kaczmarz算法的核心公式：
$$c^{k+1}=c^k+\lambda \frac{u_i-s_i\cdot c^k}{||s_i|{|}_2^2}{s}_i$$
其中，k表示第k次迭代的结果，$i=kmodM+1$, $\lambda$ 为松弛因子，目的是使算法收敛到加权最小二乘解上，迭代M次后记为一次完整迭代。

​ 简单理解这个公式，就是算法通过正向计算得到一个投影$s_i\cdot c^k$，然后将计算的投影与实际测量值$u_i$相比较，然后使用二者的差值优化参数$c$。迭代进行这一优化过程，直到达到最优解（近似）。其中$c$初始可为任意值。

---

Kaczmarz迭代算法公式推导

​ 首先回到方程（1），可以将方程（1）看作M个方程，每个方程有N个变量，这N个变量对应了空间中的N个维度，那么每个方程就可以被视为N维空间中的一个超平面。求解这M个方程就可以被视为求解M个超平面的交平面（（N-1）维）。那么Kaczmarz算法是如何实现这一求解过程的呢？

​ 算法通过迭代求解的方法，使得每一次迭代的结果都能满足M个方程中的一个，而每次迭代满足的方程与前一次不同（通过随机选择，或顺序选择的方式实现不同方程的选择），这样结果会随着迭代的进行逐渐向最优值收敛。

​ 那么又出现一个问题，为什么满足了其中一个方程就能使得结果逐渐收敛呢？换个问法，它是怎么迭代计算满足该次得方程得呢？Kaczmarz通过把当前解向其他任意一个超平面投影得方式来解决这一问题。

​ 试想一下，若干个超平面若相交，那么一个超平面中的点，向另一个超平面进行投影，新的投影点是不是距离最优解更进一步了呢？想不出来？没关系，下面这张二维空间求解图能帮助你理解。

​ 接下来，举个简单的例子实现核心公式的推导。

​ 若只对空间中的两个相邻点（N=2）求解粒子分布，最高采集到了系统矩阵和粒子信号的二次谐波频率，即M=2。数学上将N维空间超平面的求解，简化为求解二维空间中两（M）条直线的交点。求解过程如下：

​ 假设两个方程在空间中的位置关系如下图所示：从初始的随机点开始，迭代k次之后到达S点，第k+1次迭代由S点向另一条直线映射，投影点为P点。P点所在方程为$s_1\ast c_1+s_2\ast c_2=u_{k+1}$

​ 可知向量$\overrightarrow{OP}=\overrightarrow{OS}+\overrightarrow{SP}$, 其中S为第k次迭代求得的结果（为已知点），记为$C_k$,为了求$\overrightarrow{SP}$,做辅助线 AS 与被投影直线平行,过原点做投影直线的垂线，与两条线分别交于A，B两点，所以有：
$$\begin{gathered} \overrightarrow{SP}=\overrightarrow{AB}, \\ \overrightarrow{AB}=\overrightarrow{OB}-\overrightarrow{OA}, \end{gathered}$$
​ 求解$\overrightarrow{OA},\overrightarrow{OB}$即可得到$\overrightarrow{SP}$；由直线方程可知，$\overrightarrow{AB}$的方向向量为$\overrightarrow{a_{k+1}}=(s_1,s_2)$,单位方向向量为$\vec{e}=\frac{\overrightarrow{a_{k+1}}}{|\overrightarrow{a_{k+1}}|}$，因此，由点到直线的距离方程可得：
$$\begin{gathered} |\overrightarrow{OB}|=\frac{|s_1\ast 0+s_2\ast 0+u_{k+1}|}{\sqrt{s_1^2+s_2^2}}=\frac{u_{k+1}}{|\overrightarrow{a_{k+1}}|} \\ |\overrightarrow{OA}|=\overrightarrow{OS}\cdot \vec{e}=\overrightarrow{OS}\cdot \frac{\overrightarrow{a_{k+1}}}{|\overrightarrow{a_{k+1}}|}=\frac{\overrightarrow{C_k\cdot {\stackrel{⃗}{a}}_{k+1}}}{|\overrightarrow{{\stackrel{⃗}{a}}_{k+1}}|} \end{gathered}$$
所以：
$$\overrightarrow{AB}=\overrightarrow{OB}-\overrightarrow{OA}=\frac{u_{k+1}-C_k\cdot {\vec{a}}_{k+1}}{|\overrightarrow{a_{k+1}}|}\overrightarrow{a_{k+1}}$$
由此可得迭代公式：
$$\begin{gathered} \overrightarrow{OP}=\overrightarrow{OS}+\frac{u_{k+1}-C_k\cdot {\vec{a}}_{k+1}}{|\overrightarrow{a_{k+1}}|}\overrightarrow{a_{k+1}} \\ \because \overrightarrow{OP}=C_{k+1},\overrightarrow{OS}=C_k \\ {C}_{k+1}=C_k+\frac{u_{k+1}-C_k\cdot {\vec{a}}_{k+1}}{|\overrightarrow{a_{k+1}}|}\overrightarrow{a_{k+1}} \end{gathered}$$
​ 观察上述迭代公式，可以发现，迭代公式可以简单的理解为，第k+1次迭代为，第k次的结果加上 使用第k次的结果预测第k+1次时的偏差。

参考：

代数重建算法 - Yoshieraの部屋 (gitee.io)

Kaczmarz算法原理 - 知乎 (zhihu.com)