利用吉洪若夫正则化及其西尔韦斯特方程来修复受损图像

---

前言

首先要明白，我们在得到图片时是不可能的都是一张张完美的图片，像是曝光，模糊等如果需要得到一张清晰且没有噪声污染的图片往往需要经历图像修复这一步骤。

---

提示：以下仅代表个人理解
论文出处请点击这里.

一、吉洪诺夫正则化(Tikhonov regularization )是什么？

1.正则化

定义：正则化(regularization)，是指在线性代数理论中，不适定问题通常是由一组线性代数方程定义的，而且这组方程组通常来源于有着很大的条件数的不适定反问题。大条件数意味着舍入误差或其它误差会严重地影响问题的结果。 另外给出一个解释性定义：对于线性方程Ax=b，当解x不存在或者解不唯一时，就是所谓的病态问题(ill-posed problem). 但是在很多时候，我们需要对病态问题求解，那怎么做？ 对于解不存在的情况，解决办法是增加一些条件找一个近似解；对于解不唯一的情况，解决办法是增加一些限制缩小解的范围。这种通过增加条件或限制要求求解病态问题的方法就是正则化方法。 正则化的英文是regularization，即规则化，调整。通过一些调整或者其他办法，使病态问题也能得到唯一解。在这个调整的过程中，使用的技术就是正则化技术，所用的方法就是正则化方法。

2.求解线性方程的标准方法是最小二乘法，即求解min

，对于病态的线性方程，吉洪诺夫提出使用
的方法，这里$\Gamma$叫做吉洪诺夫矩阵（Tikhonov matrix）.

二、西尔维斯特方程（Sylvester equation）是什么？

1.基本概念

西尔维斯特方程（Sylvester equation）是控制理论中的矩阵方程，形式如下：
AX+XD=C
其中A、B及C是已知的矩阵，问题是要找出符合条件的X。其中所有矩阵的系数都是复数。为了要使方程成立，矩阵的行和列需要满足一定条件，A和B都要是方阵，大小分别是n和m，而X和C要是n行m列的矩阵，n和m也可以相等，四个矩阵都是大小相同的方阵。

西尔维斯特方程有唯一解X的充分必要条件是A和-B没有共同的特征值。

AX+XB=C也可以视为是（可能无穷维中）巴拿赫空间中有界算子的方程。此情形下，唯一解X的充份必要条件几乎相同：唯一解X的充份必要条件是A和-B的谱互为不交集。

2.求解方法

对于Sylvester方程：



对于变形后的矩阵方程可以很容易得到解：

上述 vec 是对矩阵列的拼接，X = [X1，X2， ……，Xn]， 则 vecX = [X1，X2， ……，Xn]T = np.reshape（(d*n,1)，order = “F”）。d为X行数，n为X列数。

⨂为克罗内克积。

克罗内克积是两个任意大小矩阵间的运算，表示为 A x B。如果A是一个 m x n 的矩阵，而B是一个 p x q 的矩阵，克罗内克积则是一个 mp x nq 的矩阵。

参考：西尔维斯特方程（Sylvester equation）一般求解方法.
一下是一些克罗克内积的公理

这些公理将被用来解斯尔维斯特方程，便于更好理解所列举出来的

三、krylov方法是什么？

Krylov方法是一种 “降维打击” 手段，有利有弊。其特点一是牺牲了精度换取了速度，二是在没有办法求解大型稀疏矩阵时，他给出了一种办法，虽然不精确。

假设你有一个线性方程组：
Ax=b
其中A是已知矩阵，b是已知向量，x是需要求解的未知向量。

当你有这么个问题需要解决时，一般的思路是直接求A的逆矩阵，然后x就出来了：
x=A^{-1}b

但是，如果A的维度很高，比方说1000*1000的矩阵，那么A就是一个大型矩阵，大型矩阵是很难求逆的，如果A还是一个稀疏矩阵，那就更难求了。这时聪明的Krylov想到了一种方法来替换A^{-1}

其中\beta都是未知标量，m是你来假设的一个值，最大不能超过矩阵的维度，比如这里例子里是1000.
瞧，这么一处理，我们就不用算A^{-1}了。我们只要求出方程里那些\beta的值，就齐活儿了。

（Krylov通过数学上的推导证明了，当m趋近于矩阵维度时（这里是1000），算出来的值就是精确解了。当然很少有人会真的把m提到那个数量级来算，那样就等于新构建了一个大型线形方程组，计算量还是很大。不过这么转换一下也不是没有好处，毕竟从稀疏矩阵变为了非稀疏矩阵，好求一点，没准就能直接求逆了。）
(这里省略了几步，还要用Arnoldi方法做个循环，先留个空，有同学需要我再补上)

解β值要带回第一个公式，得到以下方程：

有细心的同学一看，说不对劲啊。b的维度是1000，那就是有1000个方程，β的数量小于1000. 那不是方程数大于未知数了吗？这种情况应该没法儿求解啊。

对的，这种情况确实没法儿精确求解，只能求近似解。
方程数大于未知数时常用的方法之一是最小二乘法。那么这里可不可以用最小二乘法呢？

一般来说，最小二乘法应用的最重要的条件之一，就是方程须是线性的，最小二乘法一般只用来解线性方程，解非线性的就非常困难，需要进行一些“魔改”，比如基于最小二乘法的Levenberg-Marquardt and trust-region methods，就是matlab里的fsolve函数调用的算法，这里我就不铺开讲了，免得读者分心。我们观察了一下这个方程，正好就是线性的，那么就可以用。

（岔个话，非线性方程组的求解一直是个“老大难”的问题，一般可用的方法只有Newton（牛顿）法，对就是三百年前英国那个牛顿，这么些年一直没啥进步。我们研究Krylov方法，其最重要最广泛的应用，就是可以跟Newton法结合起来，把牛顿法里一般需要手动求解的一个非常复杂的Jacobian矩阵给省去了。创造这一天才结合的科学家将这种耦合算法称作JFNK，就是Jacobian-free Newton Krylov的缩写，意图一目了然，从此科学家们省去了手推Jacobian矩阵的烦恼，人人用了都说好，所以学Krylov算法的同学不顺便学一学JFNK就是“入宝山而空手归”了。）

令

这里r^{(m)}是指当m为m时的残量，所谓残量，就是error，就是我们不想要他存在的一个量。从上面的第一个公式就可以看出来，如果我们最终得出的x完全精确，那么r应该等于0. 于是现在这个问题转变为求一个含有多个自变量的表达式的最小值问题。

含有多个自变量的表达式的最小值问题，可以用最小二乘法来解决。最小二乘法的核心就是以下这些个公式：

（注：这里的r指的是r^{(m)}的平方和）

意思就是在r为最小值的时候，r关于所有变量的偏导都应当为0，这是毫无疑问的。

于是问题转化为了一个求m个方程m个未知数的方程组的问题，而且m通常不大（当然，m是你自己设定的，设那么大不是自找麻烦么）

这种问题就很好解了，一般用前面的x=A^{-1}b方法就可以搞定了。
回顾一下，大概是这样一个流程：
大型稀疏矩阵求逆–>Krylov方法–>线性方程最小二乘问题–>小矩阵求逆
出自krylov方法.

---

总结

以上这些仅仅只是要修复图像的一些算法概念，很多时候图像损坏并不是简单的噪声加成或是单个的模糊矩阵，它涉及到一个非线性的问题且得出来的解还不是唯一的；我们只能通过以上的算法去无限的逼近它的准确值从而得出未损坏前的图像像素
后面会讲到细节一步步实现图像修复