发表于 | 分类于 | | 阅读次数

简介

主成分分析在信号处理中有着极其重要的作用，这里分析的主成分分析和机器学习中PCA降维略有差别，但是其实两者的实质是一样的。这篇博客会从矩阵理论的角度对主成分分析做一个细致的介绍。

PCA 模型

PCA在矩阵问题中的最初模型是找到秩为k的矩阵$Z$来近似矩阵$A$。优化模型如下：
$$ \text{minimize} \quad ||A-Z||_F^2 \\
\text{subject to} \quad rank(Z) \leq k$$
这里假设矩阵$Z$的大小是$Z\in R^{m \times n}$。由于$Z$的秩是小于等于k的，所以可以令$Z=XY$，其中$X\in R^{m \times k}$，$Y\in R^{k \times n}$(为了简便直接令$rank(Z)=k$)。那么模型转化为$$\text{minimize} \quad ||A-XY||_F^2 = ||\Sigma-U^{T}XYV||_F^2$$，所以可以直接求解得到$XY=U\Sigma_k V^{T}$，然后让$X=U_k\Sigma_k^{1/2}$，$Y=\Sigma_k^{1/2} V_k^{T}$。很显然$X$,$Y$的解不是唯一的。例如它们的取值还可以是$X=U_k\Sigma_k^{1/2}G_k$，$Y=G_k^{-1}\Sigma_k^{1/2} V_k^{T}$。

二次正则化PCA

为了对矩阵$X$和$Y$进行一定的约束，在目标函数中加上了二次正则项。这时求解过程会变得相对复杂些。我们首先看优化模型：
$$\text{minimize} ||A-XY||_F^2+r||X||_F^2+r||Y||_F^2$$
当然，$X\in R^{m \times k}$，$Y\in R^{k \times n}$。目标函数是一个二范数的凸函数，通过直接求导得出最小值。对$X$和$Y$的求导结果分别如下：
$$(XY-A)Y^{T} + rX=0, (XY-A)^{T}X + rY^{T} = 0$$其中的求导过程要用到矩阵迹的求导公式，$$||A||_F^2 = \sqrt{tr(A^{T}A)}$$
$$\frac{\partial}{X}Tr(B^{T}X^{T}CXB)=C^{T}XBB^{T}+CXBB^{T}$$
$$\frac{\partial}{X}Tr(AXB)=A^{T}B^{T} \quad \frac{\partial}{X}Tr(AX^{T}B)=BA$$
化简求导结果，我们可以得到如下等式，其中左边的等式乘以了$X^{T}$，右边等式乘以了$Y$。
$$X^{T}(A-XY)Y^{T}=rX^{T}X \quad Y(A-XY)^{T}X=rYY^{T}$$，所以$XX^{T}=Y^{T}Y$。利用上述结论重写我们的优化条件就是：
$$\left[
\begin{array}{cc}
-rI \quad A\\
A^{T} \quad -rI
\end{array}
\right] \left[
\begin{array}{c}
X\\
Y^{T}
\end{array}
\right] =
\left[
\begin{array}{c}
X\\
Y^{T}
\end{array}
\right] (X^{T}X)
$$
矩阵$X^{T}X$和矩阵$\left[
\begin{array}{cc}
-rI \quad A\\
A^{T} \quad -rI
\end{array}
\right]$有相同的特征值（注意是包含于的关系），并且特征向量张成的空间等于$\left[
\begin{array}{c}
X\\
Y^{T}
\end{array}
\right]$的子空间。假设$A=U\Sigma V$，那么矩阵$\left[
\begin{array}{cc}
-rI \quad A\\
A^{T} \quad -rI
\end{array}
\right]$ 的特征值为$-r\pm \sigma_i$，因为矩阵$X^{T}X$的特征值非负，所以其特征值为$-r+ \sigma_i(\sigma_i \geq r)$，对应的特征向量为
$\left[
\begin{array}{c}
u_i\\
v_i
\end{array}
\right]$。由于特征向量张成的空间等于$\left[
\begin{array}{c}
X\\
Y^{T}
\end{array}
\right]$的子空间。所以可以令$X=U_k(\Sigma_k - rI)^{\frac{1}{2}}$， $Y=(\Sigma_k - rI)^{\frac{1}{2}}V_k$，注意此时的$k$并不是前$k$个奇异值的含义，而是表示在正奇异值的集合中取出$k$个值。当然最后的结论肯定是前$k$个奇异值，只需要把我们的$X$和$Y$重新代入到目标函数中，就可以发现该结论。

总结

这篇博客从矩阵分析的角度对PCA和二次正则化PCA作了简单的介绍，里面涉及到的矩阵理论知识较多，但也总算写完了，并给出了不太好的理解。博客参考了论文《Generalized Low Rank Models》。若有错误，欢迎讨论。该博客欢迎分享，转载，但请声明出处，严禁抄袭，仅用于学习。

发表于 | 分类于 | | 阅读次数

简介

在凸优化问题的求解方法中有很多具体可行的优化方法，例如梯度下降法、牛顿法、交替方向乘子法（ADMM）等。这篇文章会对求解一般无约束问题的牛顿法做个简单粗略的介绍，让我们明白它的基本原理，知道如何使用它。复杂的证明和讨论不在文中涉及，毕竟不是数学专业，学的不深。

牛顿法原理分析

首先我们从最简单的一维情况开始。它的基本思想是在现有的参考点附近对$f(x)$做二阶泰勒展开，然后找到一个极小点作为下一个估计值。假设$x_k$是当前的参考点，$x_{k+1}$是待求的下一个估计值,优化目标函数是$f(x)$。那么$f(x)$在$x_k$附近的二阶泰勒展开是$\varphi(x) = f(x_k) + f’(x_k)(x-x_k) + \frac{1}{2}f’’(x_k)(x-x_k)^2$,求$\varphi’(x)=0$的解即为下一个估计值$x_{k+1}$。
$$f’’(x-x_k)+f’(x_k)=0$$推出$x_{k+1} = x_k-\frac{f’(x_k)}{f’’(x_k)}$。

将上述一维情况推广到矩阵形式。我们的表达式会变为$$\varphi(\mathbf{x})=f(\mathbf{x_k})+\nabla f(\mathbf{x_k})^{T}(\mathbf{x}-\mathbf{x_k})+\frac{1}{2}(\mathbf{x}-\mathbf{x_k})^{T}\nabla^2 f(\mathbf{x_k})(\mathbf{x}-\mathbf{x_k}),$$其中$\nabla^2 f(\mathbf{x_k})$为Hessian矩阵，是Hermitte的。在这里应注意一点，$f(\mathbf{x_k})$是关于向量$\mathbf{x_k}$的函数，但是函数值是实数。所以要注意改写为泰勒展开式的形式。令$\varphi(\mathbf{x})=0$，得到$$\nabla f(\mathbf{x_k})+\nabla^2 f(\mathbf{x_k})(\mathbf{x}-\mathbf{x_k})=0，$$这里要注意矩阵求导的计算，其中$\frac{1}{2}(\mathbf{x}-\mathbf{x_k})^{T}\nabla^2 f(\mathbf{x_k})(\mathbf{x}-\mathbf{x_k})$求导得到$$\frac{1}{2}[\nabla^2 f(\mathbf{x_k})(\mathbf{x}-\mathbf{x_k})+[\nabla^2 f(\mathbf{x_k})]^{H}(\mathbf{x}-\mathbf{x_k})]，$$根据$\nabla^2 f(\mathbf{x_k})$为Hessian矩阵化简得到上述求导结果。最后求解得$\mathbf{x_{k+1}}=\mathbf{x_k}- \mathbf{H_k}^{-1}\nabla f(\mathbf{x_k})$。可以发现矩阵形式和一维情况很类似。当然，上述推导过程很重要的一点是优化函数在参考点的二阶导数或者Hessian矩阵存在。

阻尼牛顿法简介

原始牛顿法由于迭代公式中没有步长因子，因此是定步长迭代。对于非二次型目标函数，有时会出现函数值上升，即$f(x_{k+1})>f(x_k)$的情况。这说明原始牛顿法不能保证函数值稳定地下降。严重情况下会出现迭代点列${x_k}$发散而导致计算失败。阻尼牛顿法每次迭代方向仍采用$d_k=-\mathbf{H_k}^{-1}\nabla f(\mathbf{x_k})$，但是每次迭代需要沿着此方向作一维搜索以寻求最优的步长因子$\lambda_k$，用优化公式表示是$$\lambda_k=argmin_{\lambda\in R}f(x+\lambda d_k)。$$阻尼牛顿法实质是对步长进行了优化，当然除此之外，还有其他的优化方案。

牛顿法的进一步修正

上述牛顿法作出了一个很重要的假设前提就是二阶导数或者Hessian矩阵存在，其中还要强调一点的是要保证Hessian矩阵正定。若Hessian矩阵奇异，那么不能确定后继点，简单理解就是对应一维情况时二阶泰勒展开的二次项系数为0，此时函数的泰勒展开是条直线，无法确定后继点。若Hessian矩阵非奇异，也未必正定，牛顿方向可能不是下降方向，简单说是二阶泰勒展开的二次项系数可能为负数，这样在选取后继点时，算法可能会失效。

为了解决Hessian矩阵非正定的情况，是修正$\nabla^2 f(\mathbf{x_k})$，构造出一个对称正定矩阵$\mathbf{G_k}$，从而得到$d_k=-\mathbf{G_k}^{-1}\nabla f(\mathbf{x_k})$，再沿此方向作一维搜索。构造矩阵$\mathbf{G_k}$的方法是令
$$\mathbf{G_k}=\nabla^2 f(\mathbf{x_k})+\epsilon_k \mathbf{I}$$。此时，若$\alpha_k$是$\nabla^2 f(\mathbf{x_k})$的特征值，那么$\mathbf{G_k}$的特征值是$\alpha_k + \epsilon_k$，原因可以结合对称矩阵是可相似对角化的去理解。因此，只要$\epsilon_k$足够大，就能保证$\mathbf{G_k}$的正定性。

其中需要我们注意的是当$\mathbf{x_k}$为鞍点时，有
$$\nabla f(\mathbf{x_k}) = 0, \nabla^2 f(\mathbf{x_k}) \text{不定}$$因此原来的搜索方向不能使用，此时的$d_k$应选取负曲率方向，即$d_k^{T} \nabla^2 f(\mathbf{x_k}) d_k < 0$。这里涉及到负曲率方向的定义，在张贤达的《矩阵分析与应用》中说到，当矩阵$\mathbf{H}$是非线性函数$f(x)$的Hessian矩阵时，称满足$\mathbf{p^{H}Hp} > 0$ 的向量$\mathbf{p}$为函数$f$的正曲率方向，满足$\mathbf{p^{H}Hp} < 0$ 的向量$\mathbf{p}$为函数$f$的负曲率方向。标量$\mathbf{p}^{H}\mathbf{Hp}$称为函数$f$沿着方向$\mathbf{p}$的曲率。沿着负曲率方向进行一维搜索必能使目标函数值下降。

思考

牛顿法的迭代收敛速度很快，这里就要和传统的梯度下降法进行比较一下，梯度下降法只是利用了优化函数的一阶导数，相当于利用了速度的概念，而牛顿法利用了二阶导数，也就是梯度的梯度，相当于考虑了加速度的信息，所以收敛速度更快。更通俗的理解是如果我们的泰勒级数足够多，那么求其倒数为0，只需一步就可以直接求出最优解，当然这是极端情况，例如二次函数。但是，牛顿法涉及到求解Hessian矩阵，所以计算量、存储量很大，此篇博客参考了《最优化理论与算法》陈宝林老师的著作。若有错误，欢迎讨论。该博客欢迎分享，转载，但请声明出处，严禁抄袭，仅用于学习。

发表于 | 分类于 | | 阅读次数

矩阵空间分析

矩阵和线性空间有着密不可分的关系。在矩阵$A\in R^{m \times n}$中存在着四种简单的空间，分别是列空间用$C(A)$表示，即矩阵的列向量组成的空间，进一步说是矩阵的不相关的列向量形成的空间。与列空间正交的空间是其对应的零空间即$N(A^{T})$，零空间是由$A^{T}X=0$的解向量形成的空间。列空间的维度等于矩阵的秩$r(A)$，而$N(A^{T})$零空间的维度是$m-r(A)$。至于零空间$N(A^{T})$和列向量空间正交的原因很简单$\left[
\begin{array}{c}
a_1^{T}\\
a_2^{T}\\
\vdots\\
a_n^{n}
\end{array}
\right]\mathbf{x}=
\left[
\begin{array}{c}
0\\
0\\
\vdots\\
0
\end{array}
\right]$即空间中向量满足正交性。

另外两个空间就是行向量形成的空间$C(A^{T})$和与其对应的零空间$N(A)$。行空间的维度等于矩阵的秩$r(A)$，对应的零空间的维度是$n-r(A)$。它们之间的维度关系可以从矩阵方程解的角度去考虑，因为在基础解系中所含向量的个数等于$n-r(A)$即其中一种零空间的维度。空间之间的正交关系其实也是一种正交互补的关系。

四种空间的图形示意为：
space

投影矩阵

project

投影矩阵的作用是把一个向量向另一个空间做投影，这个概念在机器学习中的作用就是降维，另外在线性判别分析（LDA）中也用到了投影的概念。上图中$P$就是投影矩阵，$Pb$就是经过投影矩阵$P$作用后的投影，从图中不难发现该投影是正交投影，另外还有斜投影。它们之间有一定的区别，正交投影矩阵满足Hermitte特性。另外所有的投影矩阵满足$P^2=P$,因为投影的再投影是其本身即$P \cdot Pb=Pb$,并且向量$b$是任意的，所以存在等式$P^2=P$。综合起来，正交投影矩阵应该满足两条性质（1）$P^{H}=P$，（2）$P^2=P$。
其实从投影的结果来看实际上是把原来的向量投影到投影矩阵$P$的列空间里，因为向量$b$就是对矩阵列向量的线性组合。

投影矩阵和最小二乘之间的联系

我们在解决线性方程$Ax = b$的时候常常会遇到无解的情况，其实这时候的一种理解是在矩阵$A$的列向量空间中无法找到$x$将向量$b$表示出来，那么我们就可以想到一种求其近似解的思路，即将向量$b$通过投影矩阵P向矩阵$A$的列空间作投影，那么其投影就是可以用列空间中的向量表示出来的。由于投影的做法会有多种，但是为了保证欧式距离$|| b-Pb ||_2$的距离最小，因此作的投影应该是正交投影。我们假设经过上述近似后，原来的求解问题变成了$A\overline{x}=Pb$的近似解。
既然投影$Pb$是落在矩阵$A$的列空间中，那么通过正交性有$A^{T}(b-Pb)=0$，所以有
$A^{T}(b-A\overline{x})=0\\
A^{T}A\overline{x}=A^{T}b$
若$A$的列向量满足不相关或独立性，那么$\overline{x}=(A^{T}A)^{-1}A^{T}b$，并且有投影矩阵$P=A(A^{T}A)^{-1}A^{T}$。上述就是通过投影矩阵来求解线性方程的近似解。

总结

上述内容对矩阵空间、投影矩阵以及投影矩阵和最小二乘之间的联系作了简单介绍。在MIT的矩阵课中有更详细的介绍。若有错误，欢迎讨论。该博客欢迎分享，转载，但请声明出处，严禁抄袭，仅用于学习。

发表于 | 分类于 | | 阅读次数

简介

图像非局部均值滤波的原理利用了噪声的非相关的特性。如下图所示，在一幅图像中，具有相同像素的图像块是很多的，而其中的噪声是不相关的,当然不相关的假设是很强的，这里应该说是噪声的相关性很弱，该假设也是符合一定的实际应用场景的。
Image

原理分析推导

我们假设无噪声像素块的值为$f(x,y)$，加性噪声为$n(x,y)$,那么添加噪声后的像素块的值为$g(x,y)=f(x,y)+n(x,y)$。我们把多个相似的像素块进行叠加后取均值得到$\overline{g}(x,y)=\frac{1}{K}\sum_{i=1}^{K}g_i(x,y)$。此处应注意相似是个模糊的概念，以量化的观点看是两个相似块间的距离大小，这个距离可以是欧式距离等，为了对原理进行粗略的说明，下面的推导中把相似的概念强化了。应用概率的观点，我们对叠加处理取平均后的像素块求期望得到，$$E[{\overline{g}(x,y)]}=\frac{1}{K}\sum_{i=1}^{K}(E[f_i(x,y)]+E[n(x,y)])$$.
由于像素块的相似性，上式中右边式子的第一项化简为$f(x,y)$，第二项由于假设噪声的均值为0，所以$$E[\overline{g}(x,y)]=f(x,y)$$.另外由于噪声的非相关性，我们可以得到$g(x,y)$的方差计算公式为$${\sigma_{\overline{g}}}^2={\sigma_f}^2+\frac{1}{K^2}({\sigma_{\eta_1}}^2+{\sigma_{\eta_2}}^2+\dots+{\sigma_{\eta_{K}}}^2)$$
所以，最后推导出的方差为：
$${\sigma_{\overline{g}}}^2=\frac{K}{K^2}{\sigma_{\eta}}^2=\frac{1}{K}{\sigma_{\eta}}^2$$.
其中${\sigma_f}^2$为无噪声的像素块的方差值为0.

由此可见，通过非局部均值的方法可以进行图像的去噪。

总结

本文是在阅读论文“Non-Local Means Denoising”后编写的，其中的假设推倒参考了圣经《数字图像处理》,若有错误，欢迎讨论。该博客欢迎分享，转载，但请声明出处，严禁抄袭，仅用于学习。