http://www.cnblogs.com/daniel-D/p/3222576.html
声明
-
之前虽然听过压缩感知和稀疏表示,实际上昨天才正式着手开始了解,纯属新手,如有错误,敬请指出,共同进步。
-
主要学习资料是 Coursera 上 Duke 大学的公开课——Image and video processing, by Pro.Guillermo Sapiro 第 9 课。
-
由于对图像处理的了解也来自与该课程,没正经儿看过几本图像方面的书籍,有些术语只能用视频中的英文来表达,见谅哈!
1. Denoising 与 MAP
故事从 denoising 说起,话说手头上有一张含有噪音的图片 Lena,如何除去噪音得到好的 clean image 呢?
对于上面的问题,用 x 值表示某个像素的灰度值,我们可以建立这样一个最小化的数学模型:
其中, y 表示已知的观测值,也就是含有噪声的原图, x 表示要恢复成 clean image 的未知值。
模型的第一项的直观作用就是,预测值 x 不要离观测值 y 太远。数学上的解释是, x 的取值概率可以看做是以 y 为均值的高斯分布,即图像带有 Gaussian noise, 第二项是规则化项。由来如下:假设 x 本来是就带有某种先验概率的分布,现在又已知观测值 y, 根据贝叶斯原理, 现在 x 的分布(后验)正比于先验概率分布与高斯分布的乘积。如果先验概率分布也正是指数分布,将乘积取负对数,就可以得到上述在机器学习里非常常见的 MAP 模型。
现在的问题是:最好的先验 (prior) 究竟是什么? G(x) 应该取什么形式? 定义图像信号的最好空间是什么?
在学术界,这方面的工作已经做得非常多,对这个问题的探讨过程可以比喻成类人猿向人类进化的过程:
第一张图, prior 假设 clean image 能量尽量小, x 要尽可能地小。第二张图, prior 认为恢复后的图像要光滑,于是产生了 Laplacian 和 low energy 的结合,朝前进化了一步。第三张图,prior 认为要考虑 edges 是不光滑滴,需要不同情况不同处理…… Sparse and Redundant 是正在讨论的问题,目前是最新的进化版本,而后面也有一些算法,虽然也成功进化成人类,可惜太胖了,行动不便—— computationally expensive and difficult。 Sparse modeling 的先验究竟是什么?要回答这个问题,还需要了解一些基础概念。
2. Sparsity and Lp Norm
-
How to Represent Sparsity
表示一个向量的稀疏程度可以用 Lp norm, 对于 alpha 向量的某一个元素为 x, Lp norm 的计算公式和函数图像如下:
我们希望不管 x 多大,它非零的惩罚是相同的,L0 norm 正好满足这个要求,它表示的意思是数出 alpha 向量中非零的个数。
-
Sparse Modeling of Signal
一张 8×8 的图片,可以表示成 64 维的向量 x ,如何进行稀疏表示?下图中假设 N = 64:
左边矩阵 D 是字典矩阵,由 K 个 N 维的列向量组成。 根据 K 与 N 的关系,又可以划分为:
-
K
>
N: over-complete, 这种情况在稀疏表示里面最常见 -
K = N: complete, 例如傅里叶变换和 DCT 变换都是这种情况
- K
<
N: under-complete
-
中间列向量 alpha 是一个稀疏向量,特点是非零项很少,图中只有三个非零项,代表 D 矩阵对应行向量的线性组合。
最后 x 向量表示恢复后的向量。
atoms 表示 D 的列向量
实际上 DCT 变换也可以看做是一种稀疏表示,它的 D 向量是由固定的且刚好完备的正交基向量组成,并且 alpha 向量也具有一定稀疏性。
对于上图,假设 D 矩阵 K >
N,并且是满秩的,那么对于任意个 N 维的向量 b (图中是 x ),肯定有 Ax = b。现在加入 Lp norm 的约束条件,限制只能用少量的 A 的列向量 (atoms 作为基,向量 b 就被固定在某个 span 内,成为了一个 Lp 优化问题:
用紫色表示平面,用青色表示 norm 取同一个值的球形(等高线),问题如下:在平面 Ax = b 平面内选出 norm 最小的最优解
当 p >= 1时,norm ball和平面的交点有多个。这是一个凸优化问题,可以用拉格朗日乘子来解决这个问题。
当 0 < p < 1 时, norm ball 可行解十分稀疏,是一个非凸优化问题,解决这类问题很难,但是却有很好的稀疏性。
当 p = 0 时, norm ball 上的点除了坐标轴,其他部分无限收缩,与平面的交点在某一个坐标轴上,非零系数只有一个。
回到第一节将的 MAP 模型, Sparse Modeling 模型就是非零系数限制在 L 个之内(意味着解在至多 L 个 atoms 组成的 span 里),尽可能接近平面:
这样,我们用少量的 atoms 组合成真实信号,而 noise cannot be fitted very well, 在投影到低维空间的过程中起到了降噪的作用。
3. Some Issues:
模型可以改成 L0 norm 的形式和其他形式来计算或者求近似吗?
解集 alpha 向量是唯一的吗?我们可以求它的近似吗?如果可以,如何估计近似程度?
应该采用什么样的字典矩阵 D 才能较好地消除噪声?字典 D 如何确定?
1. Uniqueness
假设我们已知字典矩阵 D 和稀疏向量 a, 计算出一个信号 x,即 Da = x, x 存在一个关于 D 的稀疏表示。反过来现在已知前面的 D 和 x,根据 L0 的优化问题,可以归纳为:
的解是唯一的吗?
显然不一定。比如, D 中某些 atoms 恰好相等,或者 column1 = column2 + column3, 以前由 column2 和 column3 现在只用 column1 表示即可。当然也有正面的例子,比如 DCT 变换, 基向量完全正交,解是唯一的。这与 D 中 atoms 的不相关性和数目 K 有关。
2. Sparse Coding
和上面一样,现有字典 D 和带有噪声的信号 y,进行稀疏编码的问题可以表示的 L0 优化问题:
这是一个组合优化问题。假设 alpha 的非零项数目为 L (sparse Level), 先令 L = 1, 每一个列向量尝试一遍,看看是否又满足条件的,共有 K 种组合。如果没有,再令 L = 2, 再次尝试,共有 K(K-1)/2 中组合。还没有满足条件的,则令 L = 3......组合的数目呈指数增长,这是一个 NP-hard 的问题。 实际应用中的 K = 1000, L = 10, 要穷尽所有的排列组合大概需要计算几百万年,因此要采用近似算法, 目前主要有 relaxation methods 和 greedy methods。
-
Relaxation Methods - the Basis Pursuit (BP)
我们知道, L0 norm 可以数出向量中非零 entries 的数目,具有很好的现实意义,但是由于它数学特性(求导等)极差,非常不适合作为一个优化模型中目标函数。在线性分类器中,你可以把误分点的数目作为目标函数,但是没法优化,所以,我们看到的线性分类器的的目标函数一般是 L1 norm(感知器算法), L2 norm(LMS 算法和最小二乘法)以及最大熵(Logistic Regresson)等,也能达到比较好的效果。在上一篇博客中,可以看到 L1 是菱形, L2 是球体,L1 具有更好的稀疏性(解更靠近坐标轴),所以我们采用松弛方法将 L0 norm 转换为 L1 norm:
虽然我们把 count number 变成了 count the magnitude, 但是在某些条件下,上式的解与松弛之前的解等价。上述方法也叫 BP,新定义的问题是一个凸优化问题,解决的方法很多,有 Interior Point methods, Sequential shrinkage for union of ortho-bases, Iterative shrinkage 等。 -
Greedy Methods - Matching Pursuit (MP)
第一步,找到最接近(平行) y 的 atom, 等效与在 alpha 向量上仅取一个非零项,求出最接近的 atom, 保留下来
第二步,计算误差是否满足要求,如果满足,算法停止,否则,计算出残差信号,和第一步类似,找到最接近残差向量的 atom, 保留下来
第三步,调整已选向量的系数,使得 Da 最接近 y,重复第二步 (OMP, Orthogonal Matching Pursuit)
总结一下解决这个问题的算法有:
![](http://images.cnitblog.com/blog/533521/201307/29213932-c80baea7cd26476f9c3cd2a4d6541681.png)
3. Dictionary Learning - K-SVD
字典学习的一个假设是——字典对于一张 good-behaved 的图像具有稀疏表示。因此,选择字典的原则就有能够稀疏地表达信号。有两种方法来设计字典,一种是从已知的变换基中选择,或者可以称为 beyond wavelet 变换,比如 DCT 实际上就是一个稀疏表示(高频部分系数趋向于 0),这种方法很通用,但是不能够 adapted to the signal。第二种方法是字典学习,即通过已有的大量图片数据进行训练和学习。
比如,现在有 P 个信号(张图片)要进行稀疏表示,如何学习一个字典?
![](http://images.cnitblog.com/blog/533521/201307/29214026-e3dcad9088c1425eb17b0ee2e4a08991.png)
上式字典矩阵 D 和 alpha 组成的稀疏表示 A 矩阵都是可变量,目前有几种算法解决这个问题,下面介绍 K-SVD 算法(K-Means的一种变种),idea 非常简单。假设现在有原始信号矩阵 X^T, 该矩阵的每一行表示一个信号或者一张图片, D 矩阵是字典矩阵,右下方是 sparse coding 矩阵,红色的点表示非零项:
![](http://images.cnitblog.com/blog/533521/201307/29214041-23efd2eaab8047f69085c0c3c2175fb5.png)
算法步骤如下:
Step 1: Initialize。在 X^T 矩阵中随机挑选一些行向量(一些原图),填满矩阵 D。( K-means 随机选点初始化)
Step 2: Sparse Coding. 用上一小节的方法(松弛或者贪婪法)进行稀疏编码,Row-by-Row 计算出稀疏矩阵。
Step 3: Dictionary Update. 字典以列向量为基,自然要 Column-by-Column 地更新字典。比如现在更新某一列, 下方对应的红点,根据红点找到对应的信号(图像),然后除掉其他不相关的图像,得到示意图如下:
![](http://images.cnitblog.com/blog/533521/201307/29214109-962cffa67360491a9a7d0eafff4b8206.png)
上图中字典的 atom 对四张图片都有贡献,我们调整 atom 的目的是使得这个贡献更大,从而使稀疏性表示效果更好。当然,一个 atom 只能表示一条直线,三张图片的信号极有可能不在这条直线上,我们要做的是将中间的误差降到最小,这其实就是一个最小二乘(MSE)的问题。具体做法是将最右下角的矩阵进行 SVD 分解(SVD 相关知识可参考之前我写的博客),找出主成分,然后回到 Step2, 迭代。
1. From Local to Global Treatment
图片尺寸有大有小,在 DCT 变换中,我们一般取 8×8 的方块作为一组 64 维的变换信号,在稀疏表示中,我们同样也不能把整张图片作为 X^T 矩阵,而是在大图片中取一定尺寸的 patch (假设是 8×8 的方块)作为一个 signal。对于图片中的所有的 patch (假设 ij 是 patch 的左上角坐标)组成的信号,已知字典 D 和噪声图片 y ,估计公式如下:
y: 带有噪音的图片—— the whole image
x: 要恢复的 clear image
Rij x: 以 i,j 为左上角坐标的 patch, Rij 是从 x 中提取 patch 的 0-1 矩阵
D: 字典 for all the overlapping patches
字典 D 从哪里学习?第一种选择是基于图片的数据库,第二种是直接使用要降噪的图片进行训练。还有一种可能性是:首先基于图片的数据库得到字典 D (off-line),接着来了一张要降噪的图片,我们的做法是新建一个以 D 为初始化的字典,在要处理的图片上再进行迭代(on-line),得到新字典,这个新字典更适合降噪,代价是多一些计算。
2. K-SVD Image Denoising
在上一小节中,我们提出的可能性是 D 也需要根据要降噪的图片进行再适应,所以,图片降噪的公式多了一参数:
有三个变量,处理方法是先固定其中两个,优化一个,然后迭代。从整体上来说,先用 K-SVD 算法得到字典矩阵 D 和系数编码 alpha,保持它们不动,再优化 x:
x 的最优解实际上就是所有包含 x 像素点的 patch 的平均值,比如 patch 的大小是 8×8, 那么包含图片中某一个像素点的 patch 就有 64 个,这个像素点最优解就是取这 64 个patch 对应位置的平均值。当然,你也可以用权重来调节不同位置的 patch 对 pixel 的影响,比如 pixel 在中间的 patch,权重大,pixel 在 patch 边边角角的地方,权重小。
3. Compressed Sensing
前面我们探讨了 sparse represent 的等式,这里主要讲 compressed sensing 的概念,即在稀疏表示的等号两边同时乘以矩阵 Q:
就变成了:
用公式可以表达为:
可以看到,变换后的信号被大大压缩了。在一直 x波浪 和 D波浪 的情况下求 alpha 这个问题和前面 sparse coding 非常类似。一个关键问题是:在什么条件下由已知信号 x波浪 的情况下恢复稀疏表示 alpha?显然,这个问题与矩阵 Q,字典 D 和 alpha 的 sparse level 有关,背后涉及很多数学理论。
4. Structured Sparse Models and GMM
待续...
5. Sparse Modeling and Classification-Activity Recognition
待续...