协方差 & 特征值的交集

Posted on Edited on In Views: Symbols count in article: 6.7k Reading time ≈ 6 mins.

Eigens & Covs

​ 线性代数确实很有趣，但是展开成\(\sum\)的形式就没有趣了。本文意在分析矩阵特征值在某些场合下的应用以及一些线性代数结论（特别是矩阵/向量求导）（烦人的分子分母布局）。

Figure 1. 2D-2D单应匹配

矩阵特征量与法向量估计

​ 对于一个给定矩阵A而言，其特征方程是： \[ \begin{equation}\label{eigen} \begin{array}{l} A\pmb{x}=\lambda\pmb{x}\\ (I-\lambda A)\pmb{x}=\pmb{0} \end{array} \end{equation} \] ​ 以上两个式子定义实际是同一个意思，只不过写法不同。如何理解特征（eigen）量呢？对于特征向量\(\pmb{x}\)，在矩阵A的映射下（线性变换），会导致\(\pmb{x}\)的方向完全不变（可能反向）。而另一方面，空间变换与矩阵运算的联系也有着紧密的联系，比如3 * 3矩阵就能完全表征二维平面上的平移以及旋转。实际上 3 * 3矩阵能表示任何的透视（perspective）变换： \[ \begin{equation} (x',y',z')^T=\begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix} =A\pmb{x} \end{equation} \] ​ 那么对A可以求其特征量。这些特征量表征的是什么呢？为了简单以及便于人类进行空间想象，以下讨论主要集中于二维矩阵与三维矩阵的讨论。

二维直线法向量估计

​ 实际上，在二维的情况下，估计法向量与估计直线是同一个问题。假设有一堆点，如下图所示，存在噪声并且由于实际表面是个曲面，需要选择一个合理的范围，利用范围内点的信息进行法向量估计：

Figure 2. 2D法向量估计

​ 我们需要估计一条直线（法向量形式）： \[ \begin{equation}\label{normal} a^T\pmb{x}-b=0 \end{equation} \] ​ 那么对于n个离散的点\((x_i,y_i)\rightarrow\pmb{x}_i\)，如果能找到参数\(a,b\)使得： \[ \begin{equation}\label{res1} |b-a^T\pmb{x}| \end{equation} \] ​ 这是由需要估计的参数直线\(\eqref{normal}\)决定的，那么可以将\(\eqref{res1}\)化为模平方的形式，这也就是我们需要优化的目标，实际上点的信息已经包含在直线中了，虽然没有写成常见的直线拟合最小二乘表示： \[ \begin{equation}\label{cost} C=\sum_i(b-a^T\pmb{x})^T(b-a^T\pmb{x}) \end{equation} \] ​ 加上约束，可以得到一个带约束优化问题。约束是：\(a^Ta=1\)。也就是说希望法向量是单位向量： \[ \begin{equation} \left\{ \begin{array}{ll}\label{obj} \text{min }\sum_i(b-a^T\pmb{x})^T(b-a^T\pmb{x}) \\ \text{s.t. } a^Ta=1 \end{array} \right.\rightarrow \text{min }\sum_i(b-a^T\pmb{x})^T(b-a^T\pmb{x}) + \lambda(a^Ta-1)\\ \end{equation} \] ​ 也就是又用Lagrange乘子法将有约束问题转化为无约束问题了。对KKT条件进行讨论，等式约束的KKT条件就不用说了。对\(a\)求导，有（这里有必要穿插一下矩阵与向量的求导知识）。如果x是一个 n * 1的列向量，那么一个关于x的标量函数关于x求导，求出将会是一个行向量。进一步地，可以证明（见【Appendix】）： \[ \begin{equation}\label{lemma} \begin{array}{ll} \text{if }\alpha=x^TAx,\text{where }A\text{ is }n \times n \\ \text{then } \frac{\partial\alpha}{\partial x}=x^T(A+A^T) \end{array} \end{equation} \] ​ 那么根据公式\(\eqref{lemma}\)，并且对\(\eqref{obj}\)进行关于a（2 * 1列向量）求导，此处求导的结果是： \[ \begin{equation}\label{der} 2\lambda a^T-2bn\overline{x}^T+\sum_{i=1}^n2a^T(x_ix_i^T) \end{equation} \] ​ 而对b求导实际上也能得到重要的结论：即\(b=a^T\overline{x}\)，那么带入到公式\(\eqref{der}\)中就会有： \[ \begin{equation}\label{der2} 2\lambda a^T+2a^T\sum_{i=1}^n(x_ix_i^T)-(\overline{x}\overline{x}^T)=2\lambda a^T+2a^T\sum_{i=1}^n(x_i-\overline{x})(x_i-\overline{x})^T \end{equation} \] ​ 而\(\eqref{der2}\)的右半部分求和式实际上是点集的协方差，最后也就得到了（忽略负号给\(\lambda\)代来的影响）： \[ \begin{equation}\label{eig} \left(\sum_{i=1}^n(x_i-\overline{x})(x_i-\overline{x})^T\right)a=\lambda a \end{equation} \] ​ 也就是说，\(\lambda\)为点集协方差矩阵的一个特征值，而\(a\)则为对应的特征向量。也就是说，\(a\)（法向量）的解实际上是点集协方差的一个特征向量，但是对于二维问题而言，特征向量有两个（一个最大特征向量，一个最小特征向量）。那么显然，根据\(\eqref{eig}\)以及\(b=a^T\overline{x}\)，带入到目标函数\(\eqref{obj}\)中可以得到： \[ \begin{equation} L=0+\sum_{i=1}^n[a^T(x_i-\overline{x})][a^T(x_i-\overline{x})]^T=n\lambda \end{equation} \] ​ 可以知道，\(\lambda\)应该选择最小特征值。那么\(a\)的解就是最小特征值对应的特征向量。感觉这个过程就十分的精妙，主要是要完全利用KKT条件，一个都不能漏，并且有意识地构造协方差 / 均值以及\(x^Tx / xx^T\)两种形式。不难发现，这种表示形式可以自然推广到三维空间中去，此时求取的不再是直线的法向量，而是平面的法向量。

2D图像配准闭式解

​ 考虑一个2D问题，两幅图像的配准。假设我们通过SIFT算子 + RANSAC优化，得到了准确的特征点配对关系，那么如何求这个单应矩阵呢？

​ 实际上，纯2D-2D匹配并不能称为单应矩阵（单应矩阵是同一物体3D-3D投影变换的表征）

​ 考虑简单的仿射变换，假设仿射变换矩阵为\(H\)，给定source点集与target点集： \[ \begin{equation} P=(p_1,p_2,...p_n),Q=(q_1,q_2,...,q_2),P,Q\in \mathbb{R}^{3\times n} \end{equation} \] ​ 那么假设\(HP=Q\),也即点集\(P\)经过变换\(H\)后直接成了\(Q\)，那么为了使变换后的投影误差最小，考虑最小二乘问题： \[ \begin{equation} L=\Vert HP-Q\Vert^2=tr[(HP-Q)^T(HP-Q)]\Longleftrightarrow tr[(HP-Q)(HP-Q)^T] \end{equation} \] ​ 那么根据矩阵迹的求导（以下性质将在【Appendix】中进行简要的分析）

性质 结论（分母布局） 分子布局
轮换性 \(tr(ABC)=tr(BCA)\) \(tr(ABC)=tr(BCA)\)
线性1 \(\frac{\partial tr(A^TB)}{\partial A}=B\) \(\frac{\partial tr(A^TB)}{\partial A}=B^T\)
线性2 \(\frac{\partial tr(AB)}{\partial A}=B^T\) \(\frac{\partial tr(AB)}{\partial A}=B\)
二次型1 \(\frac{\partial tr(ABA^T)}{\partial A}=A(B+B^T)\) \(\frac{\partial tr(ABA^T)}{\partial A}=(B+B^T)A^T\)
二次型2 \(\frac{\partial tr(A^TBA)}{\partial A}=(B+B^T)A\) \(\frac{\partial tr(A^TBA)}{\partial A}=(B+B^T)A^T\)

​ 可以得到： \[ \begin{equation}\label{H} \frac{\partial \text{ tr}[(HP-Q)(HP-Q)^T]}{\partial H}=2HPP^T-2QP^T=0 \end{equation} \] ​ 则可以得到：\(H=QP^T(PP^T)^{-1}\)是解。其中\((PP^T)^{-1}\)是伪逆，由于其行列式可能为0。可以用SVD分解得到，正常情况则使用LDLT进行分解即可。

Appendix

向量矩阵求导法则以及基本结论证明

​ 这里涉及的都是基本的线性代数（以及高数中的多元函数Jacobian求解），但是自己从来认真推过。向量矩阵求导主要分析其中几个公式：

基本法则

  • 标量对行向量求导，结果是行向量（朴素结论）
  • 标量对列向量求导，结果是列向量（朴素结论）
  • 标量对矩阵，结果是同样大小的矩阵，并且每个元素与求导矩阵对应（朴素结论）

行列求导

​ 行向量对列向量求导，结果显然是个雅可比，由于求导元素是列向量，那么雅可比矩阵是按列组织的： \[ \begin{equation} d(y_1,y_2,...,y_n)/d \begin{pmatrix} x_1\\ x_2\\ \vdots \\ x_m \end{pmatrix}= \begin{pmatrix} \frac{\partial y_1}{\partial x_1} & \frac{\partial y_2}{\partial x_1} & ... & \frac{\partial y_1}{\partial x_1}\\ \frac{\partial y_1}{\partial x_2} & \frac{\partial y_2}{\partial x_2} & ... & \frac{\partial y_1}{\partial x_2} \\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial y_1}{\partial x_m} & \frac{\partial y_2}{\partial x_m} & ... & \frac{\partial y_1}{\partial x_m} \end{pmatrix} \end{equation} \] ​ 则反之，如果是列向量对行向量求导，那么将会是上矩阵的转置。注意，上式我写的是分母布局，之后我将完全使用分子布局，因为比较自然。分子布局与分母布局的区别是：

  • 分子布局：如果f是个列向量，x是个标量，那么\(\frac{\partial f}{\partial x}\)为列向量（自然表达）
  • 分母布局：如果f是个标量，x是个列向量，那么\(\frac{\partial f}{\partial x}\)是个列向量（感觉不太自然嗷）

​ 也就是说，列向量对标量求导是不变形的，矩阵对标量求导也是不变形的（那么自然，行向量对标量求导也不变形，这就是一个比较自然的表达）。接下来讨论的公式都基于几个共同点：

  • 向量\(\pmb{x}\)是列向量（n * 1）
    • 矩阵A是m * n矩阵或是n * n矩阵79

\[ \begin{equation} \frac{\partial A\pmb{x}}{\partial \pmb{x}}=A\tag{prop 1} \end{equation} \]

​ 这个很显然，讨论\(A\pmb{x}\)时只需要将A按行分块，\(A\pmb{x}\)每个元素都是\(\pmb{a}_j\pmb{x}\)（标量），那么一个列向量\(A\pmb{x}\)中的每个标量元素对列向量求导，得到一个行向量： \[ \begin{equation} \pmb{a}_j\pmb{x}=\sum_{i=1}^na_{ji}x_i \end{equation} \] ​ 显然，对这个标量求导得到\(\pmb{a}_j\)本身。 \[ \begin{equation} \frac{\partial \pmb{x}^TA\pmb{x}}{\partial \pmb{x}}=\pmb{x}^T(A+A^T)\tag{prop 2} \end{equation} \] ​ 展开之后易于证明，由于\(\pmb{x}^TA\pmb{x}\)是标量，标量对列向量求导是行向量，可以得到，行内的第j个元素为： \[ \begin{equation} \sum_{i=1}^na_{ji}x_i+\sum_{j=1}^na_{ji}x_j \end{equation} \] ​ 上式前半部分是行向量\(\pmb{x}^T\)与矩阵\(A^T\)的某一列进行相乘，后半部分则与\(A\)的某一列相乘。\(\text{prop 2}\)的证明又是比较容易的，可以根据基本的求和式与矩阵运算性质得到，关于矩阵行列求导操作，就只提这两个。

矩阵迹的求导法则

​ 迹是针对n * n矩阵而言的（对角线之和），假设\(A_{n\times k},B_{k\times n}\)，则\(\text{tr}(AB)\)有： \[ \begin{equation}\label{tr} \text{tr}(AB)=\sum_{i=1}^n\sum_{j=1}^ka_{ij}b_{ji} \end{equation} \] ​ 那么假设\(\text{tr}(AB)\)对A求导，由分子布局，将A看成一个行向量，每个行元素为一个列向量。可以知道，分子布局下，标量对行求导，得到列向量，列向量的元素将由原来行向量内的元素决定。那么由于原来行向量（1 * k）元素是一个个的列向量（n * 1），那么求导后的每个分量是个行向量（1 * n），也即最后的结果是（k * n）的，下面证明这就是矩阵B。

​ 我们假设当前正在对行向量\(A\)的元素\(\pmb{a}_p\)进行求导，\(\pmb{a}_p=(a_{1p},a_{2p},...,a_{np})^T\)是A矩阵的第p列。那么由公式\(\eqref{tr}\)，迹中与\(a_{ip}\)有关的只有： \[ \begin{equation} \sum_{i=1}^na_{ip}b_{pi} \end{equation} \] ​ 对其求导，由于上式是一个标量而\(\pmb{a}_p\)是一个列向量，应该得到一个行向量结果 \[ \begin{equation} \sum_{i=1}^na_{ip}b_{pi}\mathop{\Longrightarrow}^{d\pmb{a}_p}(b_{p1},b_{p2},...,b_{pn}) \end{equation} \] ​ 也就是说，矩阵迹对\(A\)矩阵的第p个列求导，得到B的第p行。那么由此可以得到\(d\text{ tr}(AB)/dA=B\)。（注意，这和我们上面所列表格不一致，表格是分母布局的，分母布局下这里确实是\(B^T\)）。分母布局和分子布局恰是转置关系，那么公式\(\eqref{H}\)在分子布局下应该是不变的（因为减号前后都进行了转置）。