笔记 RASL: Robust Alignment by Sparse and Low-Rank Decomposition for Linearly Correlated Images

来源：互联网发布：温长卿知乎编辑：程序博客网时间：2024/06/05 15:39

Peng, Y., Ganesh, A., Wright, J., Xu, W., & Ma, Y. (2012). RASL: Robust alignment by sparse and low-rank decomposition for linearly correlated images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(11), 2233-2246.
本文是这篇 PAMI 期刊论文的笔记，主要是对文中的理论方法进行展开详解。本人学术水平有限，文中如有错误之处，敬请指正。

摘要：此文研究的问题是同时对齐（配准）一批线性相关的图像，而排除一些破坏（如遮挡）。此文的方法在图像变换域里寻找一个最优集，使得变换之后的图像矩阵可以被分解为一个错误部分的稀疏矩阵，和一个恢复的对齐的低秩的矩阵。此文为了解决这样一个具有挑战性的优化问题，使用了一系列的凸规划，最小化稀疏矩阵的 ℓ1 范数，以及低秩矩阵的核（nuclear）范数的和，这可以使用有效的可扩展的凸优化技术解决。此文的特点就是在之前的研究基础上加入了图像变换对齐（配准），图像识别和分类的效果更好，所以称为 Robust 。

D \circ τ = A + E .

原始图像

D 对齐之后的图像

D∘τ

低秩部分

A 误差部分

如图所示，本文的方法是将一批相似度很高的图像（不仅是人脸）作为输入，首先进行平面的对齐（配准），然后通过本文的 RASL 算法，将它们分解为低秩部分（相似度非常高）和稀疏误差部分（图像之间差异值最大的部分）。

此论文并不是一项开山之作，而是在之前的研究基础之上的改进版本。首先介绍一下这个研究课题组的 leader: Prof. Yi Ma, Senior Member, IEEE, University of Illinois at Urbana-Champaign. 研究成果都展示在这里了，论文和代码都公开了，可以供感兴趣的读者下载，自行研究。
Robust Parameterized Component Analysis 1 算法虽然也拟合一个低阶模型，有效地减少了损坏和遮挡的影响。但是这是一个非凸优化问题，没有理论保证鲁棒性和收敛率。最近提出的 Rank-Sparsity 2， RPCA 3 秩最小化已经证明了它确实有可能通过凸优化方法，有效地恢复低秩矩阵。

要想看懂此文之前，还需要一点预备知识，可以学习 Alternating Direction Method 4，Singular Value Thresholding 5，Dual Method 6 和 Augmented Lagrange Multiplier (ALM) Method 7。这些方法思想都和很简单。

2 图像对齐通过矩阵的秩最小化

2.1 矩阵的秩作为衡量图像的相似度

假设我们有 n 张已经对齐的灰度图像 I01,⋯,I0n∈Rw×h。在大部分情况下，这些已经对齐的图像都是线性相关的（linearly correlated）。如果把图像列向量化 Rw×h→Rm，然后把 n 张图像并列至于一个矩阵中，

A = ˙ [v e c (I 01) | \dots | v e c (I 0 n)] \in R m \times n,

那么这个矩阵

A 应该是低秩的（low-rank）,这个性质非常的普遍。

2.2 建模非对齐作为形变

图像的偏差使得在许多应用中存在困难。上面的低秩矩阵的模型，只要在有一点图像未对齐的情况下，就会失效。文中提到，由于图像的3D结构未知，假设图像的形变只在图像平面上。于是，可以将图像未对齐建模为形变。具体来说，如果 I1 和 I2 是两个未对齐的图像，我们假设存在一种可逆变换 g:R2→R2，使得

I 2 (x, y) = (I 1 \circ τ) (x, y) = I 1 (τ (x, y)) .

在大多数实际应用中，我们把形变建模为有限维度组

G （一组参数表达式）,比如相似性

SE(2)×R+，2D仿射 Aff(2)，平面单应

GL(3)，具体可参考 An Invitation to 3-D Vision 8。
把上面两个模型结合起来：假设

I1,I2,⋯,In 表示

n 张高度相似但是相互没有对齐的图像。那么，存在一组变换参数

τ1,τ2,⋯,τn，使得

I1∘τ1,⋯,In∘τn 在每一个像素级别上都是对齐的，等价来说，矩阵

D \circ τ = ˙ [v e c (I 01) | \dots | v e c (I 0 n)] \in R m \times n,

有很小的秩，其中

I0i=Ii∘τi, i=1,2,⋯,n。所以，一批的图像对齐问题可以用如下的优化问题形式表示

min A, τ r a n k (A) s . t . D \circ τ = A .

2.3 建模损坏和遮挡作为大、稀疏的误差

因为图像之中会有一些部分被遮挡，或者损坏，但是由于这些误差出现的区域小，频率低，我们把它们建模为稀疏的误差。除此之外，图像中还有常见的噪声，但是此文没有过多考虑，只假设噪声是可忽略的。

用 ei 表示对应图像 Ii 的误差，那么 {Ii∘τi−ei}ni=1 就是已经对齐、并处理了损坏和遮挡的图像。上述的优化形式进而转化为

min A, E, τ r a n k (A) s . t . D \circ τ = A + E, | | E | | 0 \leq k .

其中

E=[vec(e1)|⋯|vec(en)]，

ℓ0 范数计算矩阵中非零元素的个数。用 Lagrangian 形式表达，将最后一个约束加入优化目标函数中，

min A, E, τ r a n k (A) + γ | | E | | 0 s . t . D \circ τ = A + E .

其中

γ 是一个权衡的标量参数。文中将这个问题定义为 Robust Alignment by Sparse and Low-rank decomposition (RASL)。

在实际图像中，一般都会有小的加性噪声，上述问题略作修改

min A, E, τ r a n k (A) + γ | | E | | 0 s . t . | | D \circ τ - A - E | | F \leq ϵ .

其中

ϵ>0 表示噪声水平。

3 迭代凸优化求解

3.1 凸松弛

上面的优化问题虽然很直观，但是矩阵的秩 rank 和 ℓ0 范数最小化是非凸的，求解非常困难（NP-hard）;并且它们都是离散值函数，如果图像并不是真的稀疏，解会不稳定。最近，发现如果矩阵 A 的秩如果小，而且 E 的非零元素值也不大，就可以使用凸形式替代它们。具体的，凸松弛方法使用了核范数（nuclear norm，矩阵的奇异值之和 ||A||∗=∑min{m.n}i=1σi(A)）代替秩 rank(⋅)，用 ℓ1 范数 ∑ij|Eij| 代替 ℓ0 范数，将优化问题转化为

min A, E, τ | | A | | * + λ | | E | | 1 s . t . | | D \circ τ - A - E | | F \leq ϵ .

理论中考虑，参数

λ 设为

C/m−−√，C 为常数，通常设为 1。新的目标函数虽然不光滑，但是已经是连续的、凸的。

3.2 迭代线性化

min A, E, Δ τ | | A | | * + λ | | E | | 1 s . t . D \circ τ + \sum i = 1 n J i Δ τ ϵ i ϵ T i = A + E .

其中

Δτi=Δτϵi 表示取矩阵

Δτ 中的第

i 列；而

ϵiϵTi 表示取其左边矩阵的第

i 列（其他列为 0），但矩阵维度不变。

因为这个线性近似只在局部有效，不能直接求解 Δτ。文中使用迭代方式重复求解凸问题，并不断估计当前的 τ。

RSAL 算法步骤：
INPUT: 图像 I1,⋯,In∈Rw×h，形变参数初始值 τ1,⋯,τn 在一个具体的参数组 G 中，权重系数 λ>0。
WHILE 没有收敛 DO // 其中收敛的条件是目标函数在连续两次迭代的改变量小于预先设定的阈值
Step 1: 计算关于形变参数 τ 的 Jacobian 矩阵：

J i \leftarrow \partial \partial ς (v e c ( I i \circ ς ) | | v e c ( I i \circ ς ) | | 2) ∣ ∣ ∣ ς = τ i, i = 1, \dots, n;

Step 2: 变形并归一化图像矩阵：

D \circ τ \leftarrow [v e c ( I 1 \circ τ 1 ) | | v e c ( I 1 \circ τ 1 ) | | 2 ∣ ∣ ∣ \dots ∣ ∣ ∣ v e c ( I n \circ τ n ) | | v e c ( I n \circ τ n ) | | 2];

Step 3(内循环): 求解线性凸优化问题；

(A *, E *, Δ τ *) \leftarrow arg min A, E, Δ τ | | A | | * + λ | | E | | 1 s . t . D \circ τ + \sum i = 1 n J i Δ τ ϵ i ϵ T i = A + E;

Step 4: 更新形变参数

τ←τ+Δτ∗；
END WHILE
OUTPUT: 输出

A∗,E∗,τ∗。

3.3 收敛性和最优性证明

涉及的内容太多，这里不做详细解释，可参考原论文。

3.4 増广 Lagrangian 乘子法

算法中主要的计算代价就是 Step 3，这是一个半正定的线性凸优化问题。此文使用快速的一阶方法，増广 Lagrangian 乘子法。首先定义 h(A,E,Δτ)=D∘τ+∑ni=1JiΔτϵiϵTi−A−E，我们构造一个 Lagrangian 函数

L μ (A, E, Δ τ, Y) = | | A | | * + λ | | E | | 1 + < Y, h (A, E, Δ τ) > + μ 2 | | h (A, E, Δ τ) | | 2 F,

其中

Y∈Rm×n 是 Lagrange 乘子矩阵，

<⋅,⋅> 是矩阵内积（

<X,Y>=trace(XTY)），

||⋅||F 是 Frobenius 范数，

μ 是一个标量（非变量），设为一个单调递增的序列，随着迭代，值不断变大。只要合适的选择 Lagrange 乘子矩阵

Y 和足够大的常数

μ，那么増广 Lagrangian 函数可以获得与原问题一致的最小解。増广 Lagrangian 算法通过迭代地最小化求解

(A k + 1, E k + 1, Δ τ k + 1) Y k + 1 = arg min A, E, Δ τ L μ k (A, E, Δ τ, Y k), = Y k + μ k h (A k + 1, E k + 1, Δ τ k + 1) .

直接最小化多个参数的 Lagrangian 函数很困难，那么采用近似求解的方式，即每一次分别最小化一个变量，而固定其他变量，

A k + 1 E k + 1 Δ τ k + 1 = arg min A L μ k (A, E k, Δ τ k, Y k), = arg min E L μ k (A k + 1, E, Δ τ k, Y k), = arg min Δ τ L μ k (A k + 1, E k + 1, Δ τ, Y k) .

迭代的每一步都可以求出闭式解，计算的效率很高。为了更清楚的说明解，还要定义软阈值（soft-thresholding）或收缩（shrinkage）操作，

S α (x) = s i g n (x) \cdot max {| x | - α, 0},

其中

α>0 是设定的一个阈值。这是一个标量函数，对于矩阵或向量操作都是 elementwise 的。结合软阈值函数，给出求解的迭代步骤

RSAL 内循环：
INPUT: (A0,E0,Δτ0)∈Rm×n,Rm×n,Rp×n，权重系数 λ>0。
WHILE 没有收敛 DO // 其中收敛的条件是 h(Ak+1,Ek+1,Δτk+1) 小于预先设定的阈值，或迭代次数达到最大

(U, Σ A k + 1 E k + 1 Δ τ k + 1 Y k + 1, V) = s v d (D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i + 1 μ k Y k - E k); = U S 1 μ k (Σ) V T; = S λ μ k (D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i + 1 μ k Y k - A k + 1); = \sum i = 1 n J † i (A k + 1 + E k + 1 - D \circ τ - 1 μ k Y k) ϵ i ϵ T i; = Y k + μ k h (A k + 1, E k + 1, Δ τ k + 1) .

// * 中间3步的详细的推导在 Appendix *
END WHILE
OUTPUT: 输出最优解

A∗,E∗,τ∗。

其中 Sα(⋅) 就是为了近似代替优化求解过程中的 ||A||∗+λ||E||1 而加入的；也就是说，在梯度求解 Ak+1,Ek+1,Δτk+1 的过程中，并没有考虑这两项，而是用减去 Ak 的部分（较小的）奇异值 S1/μk(Σ) 和减去 Ek 的较小的元素值 Sλ/μk(Ek)，来代替直接求闭式解表达式。

实验结果表明，

该算法收敛性非常好，速度也比此论文作者之前的方法 Accelerated proximal gradient (APG) 9 快 5~10 倍左右。
尽管増广 Lagrangian 乘子法的收敛性已经在优化的文献中都已经给出证明，但是没有证明它的近似方法（本文采用的迭代优化每一个变量）也能收敛。
困难就在于这里有三项变量在交替最小化，而 Alternating Direction Method of Multipliers 方法交替优化两项变量已经在优化文献中给出了证明。
然而 Yuan X. et al. 10 表明三项变量交替优化过程中收敛在实际上和理论上有差距，而且收敛速度比较慢。

总之，目前要严格证明此论文方法的收敛性还很困难。

3.5 实现细节

μk=ρkμ0，其中 ρ 和 μ0 分别设置为 1.25 和 1.25/||D||spec。RASL 算法的内循环停止的阈值为 10−7，外循环停止的阈值为 10−2。内循环算法中有一点小难点就是要计算 J†，文中使用了 QR 分解，Ji=QiRi∈Rm×p（论文中为 RTi，应该是错了），然后使用正交的 Qi∈Rm×p 代替 Ji 进行计算，相对应的，输出的结果也为 Δτ′i=RiΔτi。因为 Ri∈Rp×p 是可逆的，所以 Δτi=R−1iΔτ′i 可以很轻易地计算得到。虽然这在理论上不影响收敛性，但是在实验中可以获得更稳定的结果。

实验验证

知道了本论文的理论之后，去看实验已经没有什么难处了，本论文的实验有好几个，但是给人印象是没有很多与其它论文的定量对比结果。这里不再一一讲述。

Appendix

这是本人根据对论文的理解，自己推导了 Ak+1,Ek+1,Δτk+1 的迭代公式，首先给出

h (A, E, Δ τ) \partial h ( A , E k , Δ τ k ) \partial A \partial h ( A k + 1 , E , Δ τ k ) \partial E \partial h ( A k + 1 , E k + 1 , Δ τ ) \partial Δ τ = D \circ τ + \sum i = 1 n J i Δ τ ϵ i ϵ T i - A - E, = - 1, = - 1, = \sum i = 1 n (ϵ i ϵ T i) \otimes J T i = \sum i = 1 n J T i .

其中

⊗ 表示 Kronecker 积。

1. 求解 Ak+1，

求当梯度为 0 时，A∗ 的值，（由于 ||A||∗+λ||E||1 非光滑，求梯度时不考虑）

\partial \partial A L μ k (A, E k, Δ τ k, Y k) = \partial \partial A {t r (Y T k h (A, E k, Δ τ k)) + μ k 2 | | h (A, E k, Δ τ k) | | 2 F} = \partial h ( A , E k , Δ τ k ) \partial A Y k + \partial h ( A , E k , Δ τ k ) \partial A μ k h (A, E k, Δ τ k), - Y k - μ k (D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i - A * - E k) = 0, A * = D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i + 1 μ k Y k - E k, (U, Σ *, V) = s v d (A *), A k + 1 = U S 1 μ k (Σ *) V T .

证毕。

2. 求解 Ek+1，

求当梯度为 0 时，E∗ 的值，（由于 ||A||∗+λ||E||1 非光滑，求梯度时不考虑），与求解 Ak+1 基本一致，

\partial \partial E L μ k (A k + 1, E, Δ τ k, Y k) = \partial \partial E {t r (Y T k h (A k + 1, E, Δ τ k)) + μ k 2 | | h (A k + 1, E, Δ τ k) | | 2 F} = \partial h ( A k + 1 , E , Δ τ k ) \partial E Y k + \partial h ( A k + 1 , E , Δ τ k ) \partial E μ k h (A k + 1, E, Δ τ k), - Y k - μ k (D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i - A k + 1 - E *) = 0, E * = D \circ τ + \sum i = 1 n J i Δ τ k ϵ i ϵ T i + 1 μ k Y k - A k + 1, E k + 1 = S λ μ k (E *) .

证毕。

3. 求解 Δτk+1，

求当梯度为 0 时，Δτ∗ 的值，

\partial \partial Δ τ L μ k (A k + 1, E k + 1, Δ τ, Y k) = \partial \partial Δ τ {t r (Y T k h (A k + 1, E k + 1, Δ τ)) + μ k 2 | | h (A k + 1, E k + 1, Δ τ) | | 2 F} = \partial h ( A k + 1 , E k + 1 , Δ τ ) \partial Δ τ Y k + \partial h ( A k + 1 , E k + 1 , Δ τ ) \partial Δ τ μ k h (A k + 1, E k + 1, Δ τ), (\sum i = 1 n J T i) Y k + (\sum i = 1 n J T i) μ k (D \circ τ + \sum i = 1 n J i Δ τ * ϵ i ϵ T i - A k + 1 - E k + 1) = 0, 1 μ k Y k + D \circ τ + \sum i = 1 n J i Δ τ * ϵ i ϵ T i - A k + 1 - E k + 1 = 0, \sum i = 1 n J i Δ τ * ϵ i ϵ T i = A k + 1 + E k + 1 - 1 μ k Y k - D \circ τ, J i Δ τ * i ϵ T i = (A k + 1 + E k + 1 - 1 μ k Y k - D \circ τ) ϵ i ϵ T i, Δ τ * i ϵ T i = J † i (A k + 1 + E k + 1 - 1 μ k Y k - D \circ τ) ϵ i ϵ T i, Δ τ * = \sum i = 1 n Δ τ * i ϵ T i = \sum i = 1 n J † i (A k + 1 + E k + 1 - 1 μ k Y k - D \circ τ) ϵ i ϵ T i, Δ τ k + 1 = Δ τ * = \sum i = 1 n J † i (A k + 1 + E k + 1 - D \circ τ - 1 μ k Y k) ϵ i ϵ T i .

证毕。

E. Cande`s, X. Li, Y. Ma, and J. Wright, “Robust Principal Component Analysis?” J. ACM, vol. 58, no. 3, pp. 1-37, 2011. ↩
V. Chandrasekaran, S. Sanghavi, P. Parrilo, and A. Willsky,“Rank-Sparsity Incoherence for Matrix Decomposition,” SIAM J. Optimization, vol. 21, no. 2, pp. 572-596, 2011. ↩
E. Cande`s, X. Li, Y. Ma, and J. Wright, “Robust Principal Component Analysis?” J. ACM, vol. 58, no. 3, pp. 1-37, 2011. ↩
Yuan, Xiaoming, and Junfeng Yang. “Sparse and low-rank matrix decomposition via alternating direction methods.” preprint 12 (2009). ↩
E. Cande`s, J. Cai, and T. Shen, “A Singular Value Thresholding Algorithm for Matrix Completion,” SIAM J. Optimization, vol. 20, no. 4, pp. 1956-1982, 2010. ↩
Lin, Z., Ganesh, A., Wright, J., Wu, L., Chen, M., & Ma, Y. (2009). Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix. Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP), 61. ↩
Z. Lin, M. Chen, L. Wu, and Y. Ma, “The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices,” Technical Report UILU-ENG-09-2215, Univ. of Illinois at Urbana-Champaign, 2009. ↩
Y. Ma, S. Soatto, J. Kosecka´, and S.S. Sastry, An Invitation to 3-D Vision. Springer, 2004. ↩
K. Toh and S. Yun, “An Accelerated Proximal Gradient Algorithms for Nuclear Norm Regularized Least Squares Problems,”Pacific J. Optimization, vol. 6, pp. 615-640, 2010. ↩
X. Yuan and M. Tao, “Recovering Low-Rank and Sparse Components of Matrices from Incomplete and Noisy Observations,”SIAM J. Optimization, vol. 21, no. 1, pp. 57-81, 2011. ↩

1 0