分类目录归档：算法

统计算法的理论与实例之二 EM算法实例：形状约束的非参数混合高斯模型

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

最近比较无聊，干脆把这个系列补完。本文的内容是我表妹之前在芝大读书时的一个小题目，作为EM算法的一个比较复杂的实例来演示是不错的。由于我没找到引文，就不引用了。我不是做统计的，有错误欢迎留言。

\subsubsection{解析推导}
这是一个更为复杂的估计问题实例。要估计的分布是样本空间里一组函数图形对应的分布。具体说来，我们有随机变量 $\{\boldsymbol{X},Y,Z\}$ ,其中 $\boldsymbol{X}\in \mathbf{R}^p$ , $Y\in \mathbf{R}$ , $Z \in \{1,...K\}$ ,其中 $Z$ 是混合高斯模型中的类别标签变量,对于不同的 $Z$ 值, $Y$ 和 $\boldsymbol{X}$ 存在一个不同的潜在函数关系 $m_Z(x)$ ,其中 $m_Z(x)$ 为凸函数。概率模型表述为

(1) $\begin{eqnarray*} p(\boldsymbol{X},Y,Z)&=&\sum_j\delta_{Zj}\pi_j\mathcal{N}(m_j(\boldsymbol{X}),\sigma^2_j) \end{eqnarray*}$

现在假定 $Z$ 不可观测，我们从一系列观测值 $\{\boldsymbol{x}_i,y_i\}$ 中估计出模型的参数 $\boldsymbol{\theta}\equiv \{\boldsymbol{\pi},\boldsymbol{m(x)},\boldsymbol{\Sigma}=\{\sigma_1,...\sigma_{K}\}\}$ 的真实值 $\boldsymbol{\theta}_*$ 。

接下来我们用EM算法来处理这一问题。设 $\mathcal{X}\equiv\{\boldsymbol{x}_{1},...\boldsymbol{x}_N\},\mathcal{Y}\equiv \{y_1,...y_N\}$ 为观测样本, $\mathcal{Z}\equiv\{z_1,...z_N\}$ 为观测样本对应的隐变量。我们首先计算 $p(\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)})$

(2) $\begin{eqnarray*} \gamma_{ij}^{(t)}&\equiv&p(z_i=j\vert \boldsymbol{x}_i, y_i, \boldsymbol{\theta}^{(t)})\\ &=&\frac{p(\boldsymbol{x}_i,y_i,z_i=j \vert \boldsymbol{\theta}^{(t)})}{\sum_k p(\boldsymbol{x}_i,y_i,z_i=k \vert \boldsymbol{\theta}^{(t)})}\\ &=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}} \end{eqnarray*}$

类似于边缘分布和条件分布的因子分解性质，可以证明对任意满足 $q(\mathcal{Z})=\prod_i q_i(z_i)$ 的概率分布 $q(\mathcal{Z})$ ,有

(3) $\begin{eqnarray*} E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]=\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}$

这意味着 $Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})$ 可以分解成在每个样本上各自均值之总和的形式

(4) $\begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &\equiv& E_{\mathcal{Z}\vert \mathcal{X},\mathcal{Y},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y},\mathcal{Z})]\\ &=& \sum_i E_{z_i\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i)] \end{eqnarray*}$

我们把该性质的证明留到后面。

于是

(5) $\begin{eqnarray*} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) &=&\sum_{ij} p(z_i=j\vert \boldsymbol{x}_i,y_i,\boldsymbol{\theta}^{(t)}) \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} \log L(\boldsymbol{\theta};\boldsymbol{x}_i,y_i,z_i=j)\\ &=&\sum_{ij} \gamma_{ij}^{(t)} [\log \pi_j - \log\sigma_j - \frac{1}{2}\log{2\pi} - \frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2}] \end{eqnarray*}$

需要注意的是 $\boldsymbol{m}(\boldsymbol{x})=\{\boldsymbol{m}_1(\boldsymbol{x}),...\boldsymbol{m}_{K}(\boldsymbol{x})\}$ 包含了无穷个参数需要优化，为了约化这一问题，我们将其离散化，取 $m_{ji}=m(\boldsymbol{x}_i)$ 为优化参数，仍然简记为 $\boldsymbol{m}$ .这样， $\boldsymbol{m}(\boldsymbol{x})$ 被约化到用 $K\times N$ 个参数来表示。

接下来我们需要求解的极大化问题为

(6) $\begin{eqnarray*} \max_{\boldsymbol{\theta},\boldsymbol{\beta}} &&Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})\\ \rm{subject\ to} &&\sum_j \pi_j=1\\ &&m_{ji}-m_{jk} \geq \boldsymbol{\beta}_{jk}^{T} (\boldsymbol{x}_i-\boldsymbol{x}_k)\qquad \forall i,j,k \end{eqnarray*}$

这里我们为 $m(x)$ 加入了凸函数的形状约束。约束的几何意义很明确，即\textbf{函数总是位于每一点的切线之上方}。

由于 $Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})$ 中 $\boldsymbol{\pi}$ 的极大化和其他参数相互独立，因而可将问题解耦。问题简化为两个子问题，其中 $\boldsymbol{\pi}$ 参数对应的子问题为

(7) $\begin{equation*} \left\{ \begin{aligned} \boldsymbol{\pi}^{(t+1)}&=argmax_{\boldsymbol{\pi}} \sum_{ij}\gamma_{ij}^{(t)}\log\pi_j\\ \sum_j \pi_j&=1 \end{aligned} \right. \end{equation*}$

可以求出解析解，即

(8) $\begin{eqnarray*} \pi_j^{(t+1)}&=&\frac{\sum_{i}\gamma_{ij}^{(t)}}{\sum_{ij} \gamma_{ij}^{(t)}}=\frac{\sum_i \gamma_{ij}^{(t)}}{N} \end{eqnarray*}$

$\boldsymbol{m}$ 参数对应的子问题

(9) $\begin{equation*} \left\{ \begin{aligned} \boldsymbol{m}^{(t+1)}&=argmax_{\boldsymbol{m},\boldsymbol{\beta}} \sum_{ij}\gamma_{ij}^{(t)}[- \log\sigma_j - \frac{(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &=\sum_j argmax_{\boldsymbol{m},\boldsymbol{\beta}}[-N\log\sigma_j-\frac{\sum_i \gamma_{ij}^{(t)}(y_i-m_{ji})^2}{2\sigma_j^2}]\\ &m_{jk}-m_{ji} \geq \boldsymbol{\beta}_{ji}^{T} (\boldsymbol{x}_k-\boldsymbol{x}_i)\qquad \forall i,j,k \end{aligned} \right. \end{equation*}$

是个二次规划问题。由于 $\sigma_j,m_{ji}$ 只出现在具有相同 $j$ 值的项中，这意味着按不同的 $j$ 值把项分组，每一组都是独立变化的，因而可以作为一个独立的子问题来最大化。可以发现，\textbf{这里的Gauss混合模型中每个component都对应于一个带约束的最小二乘问题}。这一性质可以用来降低计算复杂度，提高算法效率。

在一些情况下，二次规划中的约束可以大幅度的约减。第一种情况是当 $X$ 为单变量时。不失一般性的，我们假设观测样本已经按照 $X$ 的大小做升序排列，此时可以用二阶导数非负来作为凸性约束。于是约束可约化为

(10) $\begin{eqnarray*} m_{jk}-m_{ji} &\geq& \beta_{ji} (x_{k}-x_i)\qquad \forall i,j,k\in\{i-1,i+1\}\\ \beta_{j,i+1}-\beta_{j,i}&\geq&0 \qquad\forall i\in\{1...N-1\},j \end{eqnarray*}$

这里 $\beta_{ji}$ 代表了第 $j$ 个component在第 $x_i$ 处导数的估计值。需要注意的是，为了得到一点上导数（切线）的估计，我们需要两个约束，约束在左右两侧的函数值都要位于切线上方。

另一种情况则是当 $X$ 为多维的向量，且有 $m(\boldsymbol{X})=\sum_\alpha m^\alpha(X_\alpha)$ 时。此时每个分量 $m^\alpha(X_\alpha)$ 的凸性意味着 $m(\boldsymbol{X})$ 也是凸的，因此约束可写为

(11) $\begin{eqnarray*} m^{\alpha}_{jk}-m^{\alpha}_{ji} &\geq& \beta^{\alpha}_{ji} (x_{k}-x_i)\qquad \forall \alpha,i,j,k\in\{i-1,i+1\}\\ \beta^{\alpha}_{j,i+1}-\beta^{\alpha}_{j,i}&\geq&0 \qquad\forall \alpha,i\in\{1...N-1\},j \end{eqnarray*}$

即对向量函数的每个分量，都满足\eqref{mConstraint1D}所示的约束。

最后我们来证明 $E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]$ 的分解性质，其中 $q(\mathcal{Z})$ 为任意满足 $q(\mathcal{Z})=\prod_i q_i(z_i)$ 的概率分布

(12) $\begin{eqnarray*} &&E_{q(\mathcal{Z})}[\log L(\mathcal{X},\mathcal{Z})]\\ &=&\sum_{\mathcal{Z}} q(\mathcal{Z})\log p(\mathcal{X},\mathcal{Z})\\ &=&\sum_i \sum_{\mathcal{Z}} q(\mathcal{Z}) \log p(\boldsymbol{x}_i,z_i) \end{eqnarray*}$

可见对每个sample，都需要在分布 $q(\mathcal{Z})$ 上求 $\log p(\boldsymbol{x}_i,z_i)$ 的均值。但由于 $\mathcal{Z}$ 中包含大量无关的 $z_k\neq z_i$ ,它们都应该被加和为1,从而只剩下 $q_i(z_i)$ :

(13) $\begin{eqnarray*} &=&\sum_i \sum_{\{z_1,...,z_N\}} \prod_j q_j(z_j) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i}\{\underbrace{\sum_{\{z_1,.\slashed z_i.,z_N\}} \prod_{j\neq i} q_j(z_j)}_{=1}\} q_i(z_i) \log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i \sum_{z_i} q_i(z_i)\log p(\boldsymbol{x}_i,z_i)\\ &=&\sum_i E_{q_i(z_i)}[\log L(\boldsymbol{x}_i,z_i)] \end{eqnarray*}$

\subsubsection{算法实现}
算法实现主要需要注意的是防止数值下溢带来的除零问题。技术上，避免概率连乘积，避免指数，让计算的中间结果尽量在对数域中表示是一些基本技巧。

首先，在计算 $\gamma_{ij}$ 时，可能会出现某个样本点远离所有component的情况，此时分子分母都容易下溢为零，而出现 $\frac{0}{0}$ 的情况。我们做如下处理

(14) $\begin{eqnarray*} \gamma_{ij}^{(t)}&=&\frac{\pi_j \mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}{\sum_k \pi_k \mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\mathcal{N}(m_k(\boldsymbol{x}_i),\sigma^2_k)\vert_{y_i}}{\mathcal{N}(m_j(\boldsymbol{x}_i),\sigma^2_j)\vert_{y_i}}}\\ &=&\frac{1}{\sum_k \frac{\pi_k}{\pi_j} \frac{\sigma_j}{\sigma_k} \exp(-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}+\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})}\\ \end{eqnarray*}$

算法中需要计算对数似然函数的上升率，来判断算法的终止条件是否满足。对数似然函数中有大量的概率连乘积，我们作如下处理来消除之

(15) $\begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y}) &=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(\boldsymbol{x}_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log\sum_{z_i} p(\boldsymbol{x}_i,y_i,z_i\vert\boldsymbol{\theta})\\ &=&\sum_i\log \sum_j \frac{\pi_j}{\sigma_j \sqrt{2\pi}} \exp (-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} \exp (-\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2})\\ && \sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\\ &=&\sum_i\{\log \frac{\pi_k}{\sigma_k \sqrt{2\pi}} -\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}\\ && +\log\sum_j \frac{\pi_j}{\pi_k}\frac{\sigma_k}{\sigma_j} \exp (\frac{(y_i-m_k(\boldsymbol{x}_i))^2}{2\sigma_k^2}-\frac{(y_i-m_j(\boldsymbol{x}_i))^2}{2\sigma_j^2})\} \end{eqnarray*}$

对于二次规划子问题，我们需要将其化为标准形式，从而可以利用现成的二次规划库来求解。

统计算法的理论与实例之一最大似然估计与EM算法

发表评论

去年和朋友做了一些统计方面的东西，一些笔记发在这儿。统计我是外行，贻笑大方之处请指正。

\section{记号约定}
本文中用大写字母 $X,Y,Z...$ 表示随机变量，小写字母 $x,y,z...$ 表示随机变量的观测值，粗体表示向量/矩阵。 $\mathcal{X},\mathcal{Y},\mathcal{Z}...$ 表示向量值随机变量或随机变量 $X,Y,Z...$ 的观测样本。

\section{估计理论}
给随机变量一个待定参数的概率模型，通过一系列样本观测，萃取出概率模型最可能参数的理论。

\subsection{一些数学}
考虑随机变量 $X,Y$ 的概率模型 $p(X,Y)$ ,我们经常需要考虑一系列观测样本 $\mathcal{X}=\{x_1,...x_N\}$ , $\mathcal{Y}=\{y_1,...y_N\}$ 的分布

(1) $\begin{eqnarray*} p(\mathcal{X},\mathcal{Y})=\prod_i p(x_i,y_i) \end{eqnarray*}$

观测样本之间的独立性意味着边缘分布 $p(\mathcal{X})$ （或 $p(\mathcal{Y})$ ），条件分布 $p(\mathcal{X}\vert\mathcal{Y})$ （或 $p(\mathcal{Y}\vert\mathcal{X})$ ）也满足类似的因子分解性质。严格说来，假设 $Y$ 取值于 $\Delta$ ,则 $\mathcal{Y}$ 取值于 $\otimes\Delta^N$ ,我们有

(2) $\begin{eqnarray*} p(\mathcal{X})&=&\sum_{\mathcal{Y}}p(\mathcal{X},\mathcal{Y})\\ &=&\sum_{\mathcal{Y}\in \otimes\Delta^N}\prod_i p(x_i,y_i)\\ &=&\prod_i\sum_{y_i\in\Delta} p(x_i,y_i)\\ &=&\prod_i p(x_i)\\ p(\mathcal{X}\vert\mathcal{Y})&=&\frac{p(\mathcal{X},\mathcal{Y})}{p(\mathcal{Y})}\\ &=&\prod_i \frac{p(x_i,y_i)}{p(y_i)}\\ &=&\prod_i p(x_i\vert y_i) \end{eqnarray*}$

\subsection{最大似然估计}
假设随机变量 $X$ 的概率模型为 $p(X\vert\theta)$ ,其中 $\theta$ 为模型参数。我们有一系列独立观测样本 $\mathcal{X}=\{x_1,...x_N\}$ ,欲定出最可能的模型参数 $\theta_{\star}$ .

最大似然估计即假设观测结果是最大可能的结果，从而定出分布参数的估计方法。定义对数似然函数

(3) $\begin{eqnarray*} \log L(\theta;\mathcal{X})&=&\log p(\mathcal{X}\vert\theta)\\ &=&\log\prod_i p(x_i\vert\theta)\\ &=&\sum_{i}\log p(x_i\vert\theta) \end{eqnarray*}$

则模型的最大似然估计为

(4) $\begin{eqnarray*} \theta_{\star}=\max_{\theta} \log L(\theta) \end{eqnarray*}$

注意 $-\log p(x_i,\theta)$ 是第 $i$ 次观测所获取到的信息量，此目标函数一个含义是\textbf{调整参数使得观测到的信息量取极小}。

\subsection{概率分布的描述能力}
估计理论使得我们能够以估计分布的形式得到样本空间中的几何对象。概率分布的描述能力是广泛的，这里举几个例子。

\textbf{简单椭球区域}。可以简单的由Gauss分布 $p(\boldsymbol{X})=\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ 刻画， $\boldsymbol{\mu}$ 给出了位置而 $\boldsymbol{\Sigma}$ 则给出了椭球在各方向上的轴长。

\textbf{多个椭球区域}。由多个Gauss分布的叠加，即\textbf{混合Gauss分布}描述。概率密度分布 $\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\})$ 为

(5) $\begin{eqnarray*} p(\boldsymbol{X})=\mathcal{M}(\{ \boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a,\boldsymbol{\pi}_a\}) \sim \sum_a \pi_a \mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}$

分布描述了样本空间里多个简单椭球区域，每个区域可以代表一个数据类。通常在应用中，我们把观测不到的隐含变量，类别标签 $Z$ 加入到统计量集合中，于是概率密度函数成为

(6) $\begin{eqnarray*} p(\boldsymbol{X},Z)=\sum_a\delta_{Za}\pi_a\mathcal{N}(\boldsymbol{\mu}_a,\boldsymbol{\Sigma}_a) \end{eqnarray*}$

\textbf{函数图形}。函数 $f(x)$ 在样本空间里的图形可以由概率分布

(7) $\begin{eqnarray*} p(X,Y)=\mathcal{N}(f(X),\sigma^2) \end{eqnarray*}$

刻画。给定任意 $X$ , $Y$ 都满足一个以 $f(X)$ 为中心的Gauss分布，因而这里刻画了一个对 $Y=f(X)$ 这一背后本质的数量关系所做的测量，而该测量是带有Gauss噪声的。

\subsection{最小二乘法}
考虑函数族 $f_{\boldsymbol{\alpha}}(x)$ ,我们通过一系列测量数据 $\mathcal{X}\equiv\{x_{1},...x_N\},\mathcal{Y}\equiv\{y_{1},...y_N\}$ 去估计背后真实的函数关系 $f_*(x)$ .
假设观测噪声为Gauss分布，于是观测的概率分布为

(8) $\begin{eqnarray*} p(X,Y)=\mathcal{N}(f_*(X),\sigma^2) \end{eqnarray*}$

构造对数似然函数

(9) $\begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Y})&=&\log p(\mathcal{X},\mathcal{Y}\vert\boldsymbol{\theta})\\ &=&\log\prod_i p(x_i,y_i\vert\boldsymbol{\theta})\\ &=&\sum_{i}\log \mathcal{N}(f_{\boldsymbol{\alpha}}(x_i),\sigma)\\ &=&\sum_{i} (-\log\sigma-\frac{1}{2}\log\pi -\frac{(y_i-f(x_i,\boldsymbol{\alpha}))^2}{2\sigma^2}) \end{eqnarray*}$

于是， $\max_{\theta} \log L(\theta;\mathcal{X},\mathcal{Y})$ 给出

(10) $\begin{eqnarray*} &\max_{\boldsymbol{\alpha},\sigma}&\{-N\log\sigma-\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha},\sigma}&\{N\log\sigma+\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{2\sigma^2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\frac{N}{2}\log\frac{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2}{N}+\frac{N}{2}\}\\ \Longleftrightarrow&\min_{\boldsymbol{\alpha}}&\{\sum_i(y_i-f_{\boldsymbol{\alpha}}(x_i))^2\} \end{eqnarray*}$

其中 $\sigma=\sqrt{\sum_i(y_i-f_{*}(x_i))^2/N}$ .

于是我们得到了最小二乘法。

\subsection{EM算法}
对于一个联合分布，实际应用中并非每一个随机变量都可以被观测到。比如上述混合Gauss分布中，通常类别标签 $Z$ 是无法观测到的。这种不完备观测情况下直接使用最大似然估计会出现计算量过大，实际不可解的问题。于是对于这类估计问题，我们使用\textbf{EM算法}来求解。

考虑混合Gauss分布

(11) $\begin{eqnarray*} p(X,Z\vert\boldsymbol{\theta})=\sum_a\delta_{Za}\pi_a\mathcal{N}(\mu_a,\sigma^2_a) \end{eqnarray*}$

设 $\mathcal{X}\equiv\{x_{1},...x_N\}$ 为观测样本, $\mathcal{Z}\equiv\{z_1,...z_N\}$ 为观测样本对应的隐变量。对数似然函数为

(12) $\begin{eqnarray*} \log L(\boldsymbol{\theta};\mathcal{X})&=&\log\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta})\\ &=&\log\sum_{\mathcal{Z}}\prod_i p(x_i,\mathcal{Z}\vert\boldsymbol{\theta}) \end{eqnarray*}$

其中这一对数似然函数在计算上有诸多不便，首先加和号阻止了对数运算作用在 $p$ 上，再者在 $z$ 变量取值返回较广时计算代价也会很大。因此传统的最大似然估计在此处不适用。

在EM算法中，我们在每一个迭代步都构造一个新的函数 $Q^{(t)}$ ,用 $\max Q^{(t)}$ 来近似 $\log L(\boldsymbol{\theta};\mathcal{X})$ 在当前解 $\boldsymbol{\theta}^{(t)}$ 附近邻域上的变化规律,迭代收敛到最终解

(13) $\begin{eqnarray*} Q^{(t)}(\boldsymbol{\theta})\equiv Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)})&\equiv& E_{\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}}[\log L(\boldsymbol{\theta};\mathcal{X},\mathcal{Z})]\\ \boldsymbol{\theta}^{(t+1)}&=&\mathrm{argmax}_{\boldsymbol{\theta}} Q(\boldsymbol{\theta}\vert\boldsymbol{\theta}^{(t)}) \end{eqnarray*}$

其中

(14) $\begin{eqnarray*} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})= \frac{p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})}{\sum_{\mathcal{Z}} p(\mathcal{X},\mathcal{Z}\vert \boldsymbol{\theta}^{(t)})} \end{eqnarray*}$

该算法的正确性证明留到下一小节处理。

\subsection{EM算法的理论基础}
在目标函数较为复杂的情况下，极大/极小值的寻找通常采用迭代算法来处理。计算目标函数精确的形状是困难的，因而我们需要\textbf{在每一迭代步通过有限的计算，探知目标函数在当前迭代解之邻域上的变化规律，用一个简单的函数来刻画这一局域变化趋势，然后求解这一局部的极小/极大化问题，最后进入下一轮迭代}。这样，一个复杂目标函数的极小/极大问题就被分解成了一系列容易求解的简单问题，通过迭代不断靠近极优解。

具体说来，不同的算法会为目标函数的局域变化建立不同的模型，比如通过计算导数，用平面来建模该处的目标函数的变化；再如计算到二阶导数，用抛物面来近似目标函数，等等。

EM算法也是这样一种算法。\textbf{它的特殊之处在于，由于目标函数，即概然函数的特殊结构，我们可以找到比一次/二次近似更好的局域变化模型}，从而获得比常规算法更好的性能。这一局域变化模型就是 $Q^{(t)}(\boldsymbol{\theta})$ .

首先我们来寻找 $Q^{(t)}(\boldsymbol{\theta})$ 和 $\log p(\mathcal{X}\vert \boldsymbol{\theta})$ 之间的关系

(15) $\begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&\sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)})\log p(\mathcal{X},\mathcal{Z}\vert\boldsymbol{\theta}) - \sum_{\mathcal{Z}} p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \log p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta})\\ \log p(\mathcal{X}\vert \boldsymbol{\theta})&=&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})+ H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}$

易证 $H(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\geq 0$ ,可见

(16) $\begin{eqnarray*} \log p(\mathcal{X}\vert \boldsymbol{\theta})&\geq&Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) \end{eqnarray*}$

即\textbf{ $Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})$ 构成对数似然函数的下界}。将15减去其在 $\boldsymbol{\theta}^{(t)}$ 处的值，又有

(17) $\begin{eqnarray*} \Delta \log p(\mathcal{X}\vert \boldsymbol{\theta}) - \Delta Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)}) &=& H(\boldsymbol{\theta}\vert \boldsymbol{\theta}^{(t)}) - H(\boldsymbol{\theta}^{(t)}\vert \boldsymbol{\theta}^{(t)})\\ &=&D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))\\ &\geq& 0 \end{eqnarray*}$

即\textbf{每一迭代步中必有 $\Delta \log L \geq \Delta Q$ }.

可见 $Q^{(t)}(\boldsymbol{\theta})$ 可以反映 $\log p(\mathcal{X}\vert \boldsymbol{\theta})$ 在 $\boldsymbol{\theta}^{(t)}$ 处的变化规律。更重要的是， $Q^{(t)}(\boldsymbol{\theta})$ 构成了对数概然函数的一个下界，这样以来我们就可以放心的求解 $\max Q^{(t)}(\boldsymbol{\theta})$ 这一问题来得到一个更好的迭代解，而不用像LM算法中那样担心二阶模型只在一个局域的范围内构成目标函数的有效近似，需要用置信域方法避免大步长的优化。

算法的收敛性由这一点来保证。由于 $D(p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}^{(t)}) \mid p(\mathcal{Z}\vert \mathcal{X},\boldsymbol{\theta}))$ 在 $\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}$ 附近的一阶展开项为零，于是17的一阶展开给出

(18) $\begin{eqnarray*} \frac{\mathrm{d}}{\mathrm{d} \boldsymbol{\theta}} \log p(\mathcal{X}\vert \boldsymbol{\theta})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} &=& \frac{\mathrm{d}}{\mathrm{d}\boldsymbol{\theta}} Q(\boldsymbol{\theta} \vert \boldsymbol{\theta}^{(t)})\vert_{\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}} \end{eqnarray*}$

算法迭代到 $Q$ 不再有明显变化( $\boldsymbol{\theta}^{(t+1)}\approx\boldsymbol{\theta}^{(t)}$ )时，上式右端为零，式18意味着 $\boldsymbol{\theta}=\boldsymbol{\theta}^{(t)}$ 同样也是对数似然函数的极值点。于是 $Q$ 的收敛意味着对数似然函数的收敛，算法必然收敛于一个极优解。

Notes on Conjugate Gradient Method – 2

发表评论

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

打算再将这个主页做一做，先更一波旧笔记。还是那句话，这不是教程，不适合初学。本文中记号大致与前一篇Notes on Conjugate Gradient Method一致。

当 $\mat{A}$ 对称正定时，共轭梯度法是一个很好的选择。

此时的极小化问题其实是一个 $\mat{x}$ 与 $\mat{x}^{*}$ 间的距离极小化问题

(1) $\begin{eqnarray*} &&\min_{\mat{x}}\frac{1}{2}\mat{x}^{T}\mat{A}\mat{x}-\mat{b}^{T}\mat{x}\\ \Leftrightarrow&&\min_{\mat{x}}\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}} \end{eqnarray*}$

以 $\mat{A}$ 为度量，在一个 $\mat{A}$ 正交标架下，上述距离极小化问题成为一个简单的欧式距离极小化问题，只需在依次各 $\mat{A}$ 正交方向上极小化距离即可。
这需要解决以下问题：
\begin{enumerate}
\item 一个搜索方向上如何极小化 $\lVert \mat{x}-\mat{x}^{*} \rVert_{\mat{A}}$ 。
\item 如何高效的构造出与之前搜索过的子空间 $\mat{A}$ 正交的新搜索方向；
\end{enumerate}

第一个问题。由简单的几何直观，距离的极小化等价于

(2) $\begin{eqnarray*} &&\langle\mat{e}^{(n)},D_{n}\rangle_{\mat{A}}=0\\ \Leftrightarrow&&\langle\mat{r}^{(n)},D_{n}\rangle=0 \end{eqnarray*}$

其中 $D_{n}=\mathrm{span}\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}$ 表示 $\mat{x}=\mat{x}^{(n)}$ 时已经搜索过的子空间。

本式 $D_{n}$ 换成 $\mat{d}^{(n)}$ 同样成立。所以，只需搜索到 $\mat{r}^{(n)}$ 正交于搜索方向 $\mat{d}^{(n)}$ 时为止即可达到极小化距离的目的。

第二个问题。我们令 $D_{n}=\mathcal{K}_{m}(\mat{A},\mat{r}^{(0)})$ ,其中 $\mat{r}^{(0)}$ 为初始残差，于是有 $D_{n}=\mathrm{span}\{ \mat{r}^{(0)},\mat{A}D_{n-1} \}$ .这样式2就意味着

(3) $\begin{eqnarray*} \langle\mat{r}^{(n)},D_{n-1}\rangle_{\mat{A}}=0 \end{eqnarray*}$

这带来一个极大的便利，即算法可以通过 $\mat{r}^{(n)}$ 快速构造出与 $D_n$ 子空间 $\mat{A}$ 正交的 $\mat{d}^{(n)}$ ,为此只需要将 $\mat{r}^{(n)}$ 中 $\mat{d}^{(n-1)}$ 方向的分量（ $\mat{A}$ 内积意义上）扣除即可。

这样我们可以得到共轭梯度法的算法流程：每次搜索都搜索到残差方向与搜索方向正交为止，新的搜索方向由当前残差用 $\mat{A}$ 内积与 $\{ \mat{d}^{(0)},\dots,\mat{d}^{(n-1)} \}$ 做正交化得到，实际上只需要扣除 $\mat{d}^{(n-1)}$ 方向的分量即可。
算法描述为

(4) $\begin{eqnarray*} \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \mat{d}^{(n)}&=&\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)} \end{eqnarray*}$

其中 $\alpha$ 由 $\mat{r}^{(n+1)}$ 和 $\mat{d}^{(n)}$ 的正交性给出

(5) $\begin{eqnarray*} &&\mat{r}^{(n+1)^T}\mat{d}^{(n)}=0\\ \Leftrightarrow&&(\mat{e}^{(n)}+\alpha_{n}\mat{d}^{(n)})^T\mat{A}\mat{d}^{(n)}=0\\ \Leftrightarrow&&\alpha_{n}=-\frac{\langle \mat{d}^{(n)},\mat{e}^{(n)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}=\frac{\langle \mat{d}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}} \end{eqnarray*}$

而 $\beta$ 则这样计算出

(6) $\begin{eqnarray*} &&\mat{r}^{(n)}=-\mat{A}\mat{e}^{(n)}=\mat{r}^{(n-1)}-\alpha_{n-1}\mat{A}\mat{d}^{(n-1)}\\ \Leftrightarrow&&\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle = -\alpha_{n-1} \langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}\\ \Rightarrow&&\beta_{n,n-1}=-\frac{\langle \mat{r}^{(n)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{1}{\alpha^{(n-1)}}\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n-1)},\mat{d}^{(n-1)} \rangle_{\mat{A}}}\\ &&=\frac{\langle \mat{r}^{(n)},\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{r}^{(n-1)} \rangle} \end{eqnarray*}$

共轭梯度法在加入precondition时会遇到一个问题，即如何保证 $\mat{M}^{-1}\mat{A}$ 仍是对称矩阵。
对于对称正定矩阵 $\mat{M}$ ,存在一个分解 $\mat{M}=\mat{E}\mat{E}^T$ .容易证明 $\mat{M}^{-1}\mat{A}$ 和 $\mat{E}^{-1}\mat{A}\mat{E}^{-T}$ 具有相同的本征值，因而具有同样的precondition效果。于是可以将问题转化为

(7) $\begin{eqnarray*} \mat{E}^{-1}\mat{A}\mat{E}^{-T}\hat{\mat{x}}=\mat{E}^{-1}\mat{b},\hat{\mat{x}}=\mat{E}^{T}\mat{x} \end{eqnarray*}$

再用CG求解。 $\mat{E}^{-1}\mat{A}\mat{E}^{-T}$ 的几何意义是将问题转入一个新的坐标系下求解，这个坐标系下二次型的超椭球被变换的更接近超球了一些。这样要求计算出 $\mat{E}$ ,但实际上我们可以不用计算 $\mat{E}$ 而直接用 $\mat{M}$ 完成计算。为此只需要将原来的迭代公式做坐标变换，写出新坐标系下的迭代公式即可。要注意迭代公式在该变换下并不完全保持形式一致，具体的对应关系如下所示

Rendered by QuickLaTeX.com
有稍许微分几何或相对论背景的同学不难看出这其实是一组在坐标系变换下具有不同变换性质的量之变换式，包含了协/逆变矢量，标量的变换性质，以及如何将当前坐标系下定义，非形式不变的量 $\langle \mat{r},\mat{r}^{\prime} \rangle$ 提升为任意坐标系下均保持形式不变的量 $\langle \mat{r},\mat{M}^{-1}\mat{r}^{\prime} \rangle$ 的操作.

替换后可将 $\mat{E}$ 全部消除掉，剩下仅用 $\mat{M}$ 表达的迭代公式

(8) $\begin{eqnarray*} \mat{d}^{(0)}&=&\mat{M}^{-1}\mat{r}^{(0)}\\ \mat{d}^{(n)}&=&\mat{M}^{-1}\mat{r}^{(n)}-\beta_{n,n-1}\mat{d}^{(n-1)}\\ \mat{x}^{(n+1)}&=&\mat{x}^{(n)}+\alpha^{(n)}\mat{d}^{(n)}\\ \alpha_{n}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{d}^{(n)},\mat{d}^{(n)} \rangle_{\mat{A}}}\\ \beta_{n,n-1}&=&\frac{\langle \mat{r}^{(n)},\mat{M}^{-1}\mat{r}^{(n)} \rangle}{\langle \mat{r}^{(n-1)},\mat{M}^{-1}\mat{r}^{(n-1)} \rangle} \end{eqnarray*}$

Levenberg-Marquardt算法

发表评论

匆匆一更。还没写完，对不起我太懒了T_T

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

Levenberg-Marquardt算法（下文简称LM算法）通常用于非线性最小二乘法的目标函数极小化。这是一个置信域方法（Trust-Region Method），为了防止步长太大而跳到非预期的局部极小值，这类算法自适应的调整步长。

假设 $f(x)$ 为一个向量函数，每个分量都刻画了一个样本和模型预测之间的某种偏差，则非线性最小二乘目标函数为

(1) $\begin{eqnarray*} E&=&\transp{f}f \end{eqnarray*}$

通常基于梯度的数值极小化算法，都是用一点的局域特征（如导数）建立该函数变化趋势的模型，以此来推测一个最可取的下降方向和步长。取一阶导数近似是线性模型，用切平面来近似函数在该点附近的变化趋势；取到二阶导数近似则是二阶模型，用一个抛物面来近似。当然，也可以取到更多阶，模型的几何直观也会越来越复杂，越来越精确的逼近真实的函数变化趋势，但是会付出更大的计算代价。当取到无穷阶时，模型在级数收敛区间内可以做到完全精确的刻画函数，此时这就是模型就是函数的泰勒展开了。

我们把 $f$ 在任一点 $x$ 附近展开 $f(x+\delta x)=f(x)+J(x)\delta x$ ,这样 $E$ 在 $x$ 附近就是

(2) $\begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2\\ &=&|f|^2+2\transp{f}J(x)\delta x + |J(x)\delta x|^2 \end{eqnarray*}$

注意 $J$ 是 $f$ 而非 $E$ 的导数，这里我们得到了 $E$ 的二阶展开，即用一个抛物面来近似目标函数 $E$ 在一点附近的形态。这个二阶模型显然是一个线性最小二乘问题——这正是我们所要的，把一个非线性最小二乘问题转换为一系列线性最小二乘问题来求解。

因此一个简单直接的想法就是，每一迭代步构造这样一个线性最小二乘问题，通过极小化它下降到下一个迭代点，循环直到收敛。但这个方案显然有问题。我们的二阶展开模型只能在一个小邻域上有效，当 $\delta x$ 增大时，被我们忽略掉的高阶项就会越来越明显，最后让二阶模型彻底失效。这样，直接根据二阶模型计算出来的迭代点对目标函数而言可能并非我们想要的，它的函数值可能比上一点还要高（因而目标函数不降反增）；也可能尽管比上一点下降了但是却跳到了另一个我们不想要的凸区域里（因而算法会收敛到离初始猜测解更远的局部极小值点上）；更严重的情况，模型甚至可能是退化的，即一种极端形式的抛物面——就像一本被半卷起的书，此时我们在弯折的方向仍能找到极小值，但在其垂直方向已经无法找到极小值（或者说处处极小）。

为了解决以上问题，我们需要给模型加入一些约束，把迭代步长约束在一个合理的范围内，保证在这个范围内模型足够有效。同时，我们也希望在二阶模型退化的时候，约束能够拯救它，保证算法总能求出一个极小值。这个约束很简单，如下

(3) $\begin{eqnarray*} E(x+\delta x)&=&|f(x)+J(x)\delta x)|^2 + \lambda |C(x)\delta x|^2 \end{eqnarray*}$

我们加入了一个二次约束项。单独看该项，它是一个中心在 $x$ 点，开口向上的抛物面，或者说一个二次势阱。这个势阱像一个箍，把 $\delta x$ 箍在 $x$ 点附近的小邻域中。通过调节 $\lambda$ 可以增强或减弱它对 $\delta x$ 的约束：当 $\lambda$ 充分大时，约束项成为主要贡献项，这时模型总有一个极小值，且步长较小；当 $\lambda$ 较小时，约束项几乎可以忽略，算法以大步长快速往极小值接近。

再来看矩阵 $C(x)$ . 这个矩阵通常被设置为对角的。此时，

(4) $\begin{eqnarray*} |C(x)\delta x|^2 &=& \transp{\delta x}\transp{C}C\delta x\\ &=&\sum_{i} \lambda_i \delta x_i^2 \end{eqnarray*}$

$\transp{C}C$ 的本征值控制着每个方向上抛物面上升的快慢。可见，如果一个方向相应的本征值被设置的较小，那么该方向的约束项上升缓慢，约束较弱，步长跨度较大；否则反之。现实中使用的模型，目标函数对于各个参数的敏感度不同，有的参数略微调整就可以导致目标函数的强烈变化，有的则反之。 $C(x)$ 这个参数矩阵使得我们可以针对各个不同的参数设置不同程度的步长约束，从而能够构造出一个稳定的算法。

（未完待补）

图像拼接算法原理 2

1条回复

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

2. 曲面投影

图6. 近 $90^{\circ}$ 时产生越来越大的畸变

通常简单的图像拼接技术，就是如上节所示的基本原理，找出一张大概处于中间位置的图像，然后利用单应性变换把其他图像变换到该中心图像的视角下，再做一些后续的曝光补偿、图像融合等处理即可。但是这一技术有相当大的局限性，最简单的例子，不能直接用它拼出 $360^{\circ}$ 的全景图。

为什么呢？让我们来考虑图6中所示情形。可见，三条光线与 $P^{\prime}$ 相交的三点，原本是近乎等间距均匀分布的，而当它们映射到 $P$ 平面上后，间距却产生了巨大的差异。表现在图像上， $P^{\prime}$ 上的图像变换到 $P$ 上后，会产生相当大的拉伸畸变。当图中两相机的投影平面越来越趋于垂直时，这个畸变越来越大，以至于 $P^{\prime}$ 上普通的一点可能会被映射到 $P$ 平面的无穷远点。这时，这种简单的单应性拼接方案就彻底崩溃了。

图7. 曲面投影解决大视场角下的投影畸变问题

怎样得到更宽视场角下的拼接图？解决方法很简单。以上所出现的这种畸变源自于我们将点投影到一个平面上，设想将 $P$ 掰弯，或者直接弯曲成一个圆柱面，成为图7所示的样子，那原本被投影到 $P$ 平面无穷远处的点就被拉回来了。我们在圆柱面上选取一个足够均匀的坐标系，把坐标对应到像素坐标，就可以得到一个全景图了。

当然，针对不同的应用，我们还可以选取不同的投影曲面，比如选取球面用于 $360^{\circ} * 180^{\circ}$ 的球面全景图，甚至也可以选择一个立方体作为投影曲面。

3. 后续处理

至此全景拼接的几何原理就大致说完了，虽然我们还没有给出数学表达。为了先居高临下的了解整个拼接流程，我们不妨把后续处理的梗概也在此一说。

实际应用中为了创建出完美的全景图，有很多的问题需要考虑。最典型的问题有两个，一个是如何解决不同照片中曝光不一致的问题；一个是如何在拼接缝处完美平滑的融合两张图像的问题。

第一个由曝光补偿算法来解决，大体思路是估计两张图间的曝光差异，然后进行补偿。此处不多说。

第二个问题也有众多解决方案，最为著名的大概就属Multi-Band融合算法。该算法虽然八十年代就已提出，但其效果至今仍让人赞叹。在通常图像间失配程度不大的情况下，Multi-Band可以达到肉眼几乎不可分辨的融合效果。其原理也不复杂，下面略微一提。

融合两张图像，最直接的方案是在两张图像的重合区域用一个平滑渐变的权重对二者加权叠加。该方法的效果并不理想，关键原因是我们无法兼顾拼缝附近的局域细节和大尺度上两张图片的宏观特征（如光照）。当我们希望局域细节能够完好拼接时，需要用较小的平滑渐变区；而当我们希望要宏观上平滑过渡时，又想要较大的渐变区域。这二者似乎不可调和。

但事实上并非如此。Multi-Band的成功之处就是在于它同时兼顾两种需求，当融合宏观特征时，采用一个大的平滑渐变区；融合局域细节时，则采用小的平滑渐变区。那如何才能把这两种情况分开处理呢？很简单，把图像分解为不同频带的分量之加和，图像的宏观特征在它的低频分量图里，而局域特征在高频分量图里。

所以，Multi-Band算法的过程大致就是：把图像按照频率高低展开成一个金字塔，然后高低频分量各自按照不同的方式平滑加权并叠加，最后把各频带分量重新加和，得到最终的融合结果。

该算法融合效果虽好，但对于计算量要求较大，它需要创建多座金字塔并对金字塔进行各种运算，图像像素较高时，在CPU上要达到实时基本无望。当然，GPU上情况就不一样了，我们自己就实现了实时的Multi-Band融合算法，效果很好。

这一系列文章主要以拼接的几何原理为主。下一节开始用数学建模前两节所述的投影模型。

（未完待续）

Notes on Conjugate Gradient Method

发表评论

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

这是一个关于共轭梯度法的笔记。请大家注意的是，这是个笔记，并不是一个教程，因此少不了跳跃和欠解释的地方。对CG方法了解不多的同学请移步这里。

线性方程组和极小化问题

一个关于对称矩阵 $A$ 的线性方程组 $Ax=b$ 等价于求解如下极小值问题：

$f(x) = \frac{1}{2} \transp{x} A x - \transp{b} x$

这很容易说明，我们微分目标函数得

(1) $\begin{eqnarray*} \td f &=& \frac{1}{2}(\td\transp{x}Ax+\transp{x}A\td x)-\transp{b} \td x\\ &=& \transp{(Ax - b)} \td x \end{eqnarray*}$

所以 $\td f=0$ 意味着 $Ax=b$ .

设 $x^*$ 为问题的解， $e$ 为偏离极小值点的位移，即

(2) $\begin{eqnarray*} Ax^*&=&b\\ e&=&\Delta x=x-x^* \end{eqnarray*}$

我们定义残差或负梯度 $r=b-Ax$ ，容易看出， $r$ 正比于偏离梯度零点 $x^*$ 的位移 $e$ :

$r=-Ae$

最速下降法和共轭梯度法简述

最速下降法：搜索方向为本轮迭代初始点的梯度方向，搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法：搜索方向为本轮迭代初始点的梯度方向（残差方向）用A-内积做正交化，扣除掉之前所有搜索方向的分量给出，搜索到梯度与搜索方向正交的位置开始下一轮迭代。

共轭梯度法的理解

首先考虑一个简单情形，即A=I时。此时f退化，且x的不同分量解耦，f的极小值问题分解为各分量上的极小值问题。我们只需要在各分量方向上找到极小值，问题就得解了。具体的说，我们从任意一个初始点出发，在方向 $d_0$ 上寻找极小值，然后从这个极小值出发继续在 $d_1$ 寻找 $d_1$ 方向的极小，以此类推最后到达全局极小值。

A为正定对称矩阵的一般情形，其实只是上述情形中的x表述在一个新坐标系下的结果（下文详细解释）。为了求解这种情形，我们考察上述方法的求解轨迹在新坐标系里的对应。

首先，A=I时的各个正交搜索方向在新坐标系里对应着一组A-正交的方向，因此我们需要设法找出这样一组正交方向，然后各自求解这些方向上的极小值。

求解特定方向上极小值的方法不变，极小点就是这个方向上梯度正交于搜索方向时的点（原因很简单，这说明搜索方向上梯度分量为零）。

最后一点，如何找到这样一组A-正交方向？

一般方法是找一组现成的完备基，进行A-正交化。问题是计算量比较大。

记迭代到第i步时，所搜索过的这样一组A-正交的搜索方向单位矢量为 $\{d_0,d_1,...,d_{i-1}\}$ ,它们撑起的子空间

(3) $\begin{eqnarray*} D_i=\mathrm{span}\{d_0,d_1,...,d_{i-1}\} \end{eqnarray*}$

即第 $i$ 步前已经搜索过的子空间。

共轭梯度法的一个关键之处在于它可以从每次迭代的 $r_i$ 中快速构造出新的搜索方向 $d_i$ ，而不需要对 $r_i$ 进行完整的正交化手续（即不需依次扣除 $r_i$ 中 $d_0$ , $d_1$ … $d_{i-1}$ 方向的分量）。原因在于 $r_i$ 自然的和 $D_{i-1}$ A-正交，因此构造 $d_i$ 只需要从 $r_{i}$ 中扣除掉 $d_{i-1}$ 分量即可。

因此，共轭梯度算法有两处关键，一处为通过利用A-正交关系来把搜索问题转换到一个特殊坐标系，使得各搜索方向上的极小值问题解耦；一是利用 $r_i$ 和 $D_{i-1}$ 是A-正交的这一规律来充分简化A-正交搜索方向的构造。

接下来我们详细论证这个算法。

第一个关键点，坐标变换和子问题解耦

设

(4) $\begin{eqnarray*} A&=&\transp{T}T\\ y&=&T x\\ b^{\prime}&=&T^{-T} b\\ \end{eqnarray*}$

其中 $T^{-T}$ 表示矩阵 $T$ 求逆后做转置，很拗口，不过这个不重要。

于是我们有

(5) $\begin{eqnarray*} f(x) &=& \frac{1}{2} \transp{y} y - b^{\prime T} y\\ &=& \frac{1}{2}\sum_{n}y_n^2 - b^{\prime}_n y_n \equiv \frac{1}{2}\sum_{n} P_n(y_n) \end{eqnarray*}$

可见 $f$ 已经被分解为一系列独立的子问题 $P_n(y_n)$ ， $T$ 正是将问题解耦所用的坐标变换。

在新的 $\{y_n\}$ 坐标系中，整个极小值问题可以通过在一组正交方向上依次寻找极小点来解决。但求解 $T$ 需要对 $A$ 做分解，并非一个计算上经济的解决方法。这一个变换的真正意义在于，我们可以通过它来找到求解轨迹在 $\{x_n\}$ 坐标系中的对应。

$T$ 是线性变换，因此 $\{y_n\}$ 中的一组直线仍被变换为 $\{x_n\}$ 中的一组直线。 $\{y_n\}$ 坐标系中，我们依次沿着一组正交方向求解极小值，因此只需要找到这组方向在 $\{x_n\}$ 坐标系中的对应（一组未必正交的方向），就可以等价的直接在这组方向之上寻找极小值，而不需要经过坐标变换。

考虑这组正交方向单位矢 $\{y_n\}$

(6) $\begin{eqnarray*} &&\transp{y_i} y_j = \delta_{ij}\\ &\Leftrightarrow& \transp{(T x_i)} T x_j=\delta_{ij}\\ &\Leftrightarrow& \transp{x_i} \transp{T}T x_j = \delta_{ij}\\ &\Leftrightarrow& \transp{x_i} A x_j = \delta_{ij} \end{eqnarray*}$

可见， $\{y_n\}$ 这一组正交的方向在变换到 $\{x_n\}$ 坐标系中后，虽然不是 $\{x_n\}$ 中通常意义的正交矢量组，却可以在A-内积 $\langle \bullet,\bullet \rangle_A$ 意义下成为一组A-正交矢量。这里的A-内积，指的是以对称矩阵 $A$ 定义的内积运算

(7) $\begin{eqnarray*} \langle a,b \rangle_A \equiv \transp{a} A b \end{eqnarray*}$

需要注意的是，上述结论是可逆的，如果 $\langle x_i,x_j \rangle_A = 0$ ,同样也可以得出， $x_i,x_j$ 对应着 $\{y_n\}$ 坐标系中两个正交的方向。因此，可以直接在 $\{x_n\}$ 坐标系下，寻找一组A-正交的方向，然后在这一组方向上依次极小化目标函数，就可以找到整个问题的极小值。

第二个关键点， $D_i$ 子空间结构和正交关系

算法的更新规则为

(8) $\begin{eqnarray*} e_{i+1}&=&e_i + \alpha_i d_i\\ r_{i+1}&=&-A(e_i+\alpha_i d_i)=r_i-\alpha_i A d_i\\ d_{i}&=&r_i + \sum_{j<i} \beta_{ij} d_{j}\\ \beta_{ij}&=& - \frac{\transp{r_i}d_j}{\transp{d_j}Ad_j}\\ &=& - \frac{\langle r_i,d_j \rangle }{\langle d_j,d_j \rangle _A} \end{eqnarray*}$

其中 $\alpha_i$ 为第 $i$ 步搜索方向上的位移长度，它的大小可以通过 $r_{i+1}$ 应该和 $d_i$ 正交这一要求定出

(9) $\begin{eqnarray*} &&\transp{r_{i+1}}d_i=0\\ &\Rightarrow&\transp{(e_i+\alpha_i d_i)} A d_i =0\\ &\Rightarrow&\alpha_i = - \frac{\langle d_i,e_i \rangle_A }{\langle d_i,d_i \rangle_A}= \frac{\langle d_i,r_i \rangle}{\langle d_i,d_i \rangle_A} \end{eqnarray*}$

算法将参数空间划分成了一个层级结构，下面我们来看看这个结构的性质。

根据上述更新规则，可以发现 $D_i$ 可以由另外几个矢量组撑起。首先，由 $d_i$ 的更新规则可以发现， $d_i$ 由 $r_i$ 和 $d_{j<i}$ 线性组合而成，另一方面，我们有初始条件 $d_0=r_0$ ，归纳而知， $d_i$ 可以由 $\{r_{j\leq i}\}$ 线性组合出来。于是我们有

(10) $\begin{eqnarray*} D_i=\mathrm{span}\{r_0,r_1,...,r_{i-1}\} \end{eqnarray*}$

接着我们考虑 $r_i$ 的更新规则。可以看到， $D_{i+1}=\mathrm{span}\{D_{i},Ad_i\}$ .于是我们可以从 $D_1=\{d_0\}$ 或 $D_1=\{r_0\}$ 构造出 $D_i$

(11) $\begin{eqnarray*} D_i&=&\mathrm{span}\{d_0,Ad_0,A^2d_0,...,A^{i-1}d_0\} = \mathrm{span}\{d_0,A D_{i-1}\}\\ &=&\mathrm{span}\{r_0,Ar_0,A^2r_0,...,A^{i-1}r_0\} = \mathrm{span}\{r_0,A D_{i-1}\} \end{eqnarray*}$

接下来我们再考虑各子空间和向量之间的正交关系。首先， $e_i$ 在 $D_i$ 的补空间里。这点是明显的，因为每一个优化步都会优化掉 $e_j$ 在相应 $d_j$ 方向上的分量，因而剩余的 $e_i$ 就只位于 $D_i$ 的补空间里。于是我们有

(12) $\begin{eqnarray*} e_i = \sum_{j \geq i} \delta_j d_j \end{eqnarray*}$

这样 $e$ 和 $d$ 之间的A-正交关系就很明显了。由于 $\langle d_i,d_j \rangle_A = \delta_{ij}$ ，而 $e_i$ 只包含 $d_j>i$ 分量，因而

(13) $\begin{eqnarray*} \langle e_i,d_{j<i} \rangle_A =0 \end{eqnarray*}$

即 $e_i$ 和子空间 $D_i$ 也是A-正交的。这意味着 $r_i$ 正交于 $D_i$

(14) $\begin{eqnarray*} &&\transp{e_i} A d_{j<i}=0\\ &\Rightarrow&\transp{(A e_i)} d_{j<i}=0\\ &\Rightarrow&\langle r_i,d_{j<i}\rangle =0 \end{eqnarray*}$

到这里就可以解释为什么我们有 $\langle r_i,d_{j<i-1} \rangle_A =0$ 了：

(15) $\begin{eqnarray*} &&\langle r_i,D_i \rangle =0\\ &\Rightarrow&\langle r_i,\mathrm{span}\{d_0,A D_{i-1}\} \rangle =0\\ &\Rightarrow&\langle r_i,A D_{i-1} \rangle =0\\ &\Rightarrow&\langle r_i,D_{i-1} \rangle_A =0 \end{eqnarray*}$

至此我们完成了第二个关键点的论证，且更清楚的了解了参数空间，尤其是 $D_i$ 子空间的结构。

图像拼接算法原理 1

2条回复

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

图1,2,3. 两张图片的拼接

图像拼接是计算机视觉中一个有趣的领域，它把来自多个不同视角相机的图像变换到同一视角下，无缝拼接成一张宽视野图像（比如360度全景图，甚至360度*180度的球面全景）。上图所示，即为两张图像的拼接，结果基本是完美的。

需要注意的是，由于相机各自的指向角度不一样，因此两图片中来自同样场景的部分并不能够通过平移图像而完全重合。比如上图中那栋带有岭南风格拱顶的房子，它的屋檐在左图中要比右图中更水平些，如果试图通过平移对齐像素来拼接两幅图，结果必然不自然（比如在屋檐处出现拐角）。两图要做到完美叠合，不是一个平移变换能做到的。当然，肯定存在这样一个变换，那它是什么呢？

事实上这里边的数学原理（射影几何）很古老，什么情况下不同位置、视角的相机可以变换到同一视角下拼接起来也是久为人知的，只不过能够利用计算机来进行大规模自动拼接，还是近些年才成熟起来的事。

那到底什么情况下图像可以拼接，如何拼接呢？不妨先摆出结论吧：在两种情况下图像可拼接，一是各相机几何中心重合；二是各相机位置任意，但场景是一个平面。一种特殊情况，场景为远景时，可以近似的等价于一个平面场景，从而也是可拼的。拼接的方法很简单，在满足上述两种情况之一时，存在一个单应性变换（Homography），能够将一个相机的图像变换到另一个相机的视角下从而可以进行拼接。

下面我们用一套直观的语言来讲解这个变换，以及更复杂些的拼接方法。

1. Homography

图 4. 单应性变换

如图，计算机视觉中，我们用一个投影中心点和一张图像平面来刻画一个理想的相机。点与平面的距离为焦距，任意场景点所成像点用这样一种简单的方式来决定：将它与投影中心连线，连线与图像平面的交点即像点位置。图中给出了两个共中心的相机，它们有不同的指向（因而也有不同的图像平面 $P$ 和 $P^{\prime}$ ），同一场景点 $S$ 在两个平面上的像点由紫线与两平面的交点 $D$ 、 $D^{\prime}$ 给出。到这里我们就直观的得到了这样一个变换： $P$ 上像点到 $P^{\prime}$ 上对应像点的一个映射。这个映射可以把两张不同视角下拍摄的照片变换到同一视角下，从而能够完美拼接两张图像。这个映射正是我们要说的单应性变换。

到这里我们可以回头来解释另一个问题，为什么要求相机共中心？不妨反过来考虑，如果两相机不共中心，会出现什么样的情况？如图所示，

图 5. 单应性的破坏

我们让两相机的投影中心 $O$ 、 $O^{\prime}$ 相离。这时可以看到， $S_1$ 、 $S_2$ 两点原本在相机 $O$ 上投影到同一像点，这时在相机 $O^{\prime}$ 上却投影到了不同像点。这意味着相机 $O^{\prime}$ 看到了在 $O$ 视角下被遮挡而看不到的内容，此时我们无法把 $O^{\prime}$ 看到的图像变换到 $O$ 的视角下。考虑一个极端情形可以进一步解释这个问题：如果两个相机分别拍摄同一物体的正面和背面，那它们所看到的内容无论如何也不能变换到同一视角之下，或者说，我们根本找不到这样一个新的相机，在它的视角之下可以看到 $O$ 、 $O^{\prime}$ 二者所看到内容的总和。这就解答了我们之前的问题：非共中心放置时，各相机看到的内容太丰富以至于不能变换到同一视角下。读者可以检查，在相机共中心时，一个相机中被投影到同一像点的场景点，总会在另一相机中也被投影为同一像点。

事实上，在相机摆位任意时，我们看到的信息是如此丰富，以至于可以尝试重建出场景点的三维坐标。现代的三维重建算法可以利用不同位置、视角下拍摄的图像重建出物体的三维模型，这是后话。

当然，读者仍然可能想起，除掉共中心情形，我们还说过，相机摆位任意，但场景为一个平面时也是存在单应变换的。我们在这里不再给出这种情况的直观解释，读者可以自己思考。

（未完待续）

Mind Spectrum

Formerly Singularity Vision

分类目录归档：算法

统计算法的理论与实例之二 EM算法实例：形状约束的非参数混合高斯模型

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

统计算法的理论与实例之一最大似然估计与EM算法

Notes on Conjugate Gradient Method – 2

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

Levenberg-Marquardt算法

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

图像拼接算法原理 2

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

2. 曲面投影

3. 后续处理

Notes on Conjugate Gradient Method

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

线性方程组和极小化问题

最速下降法和共轭梯度法简述

共轭梯度法的理解

第一个关键点，坐标变换和子问题解耦

第二个关键点， $D_i$ 子空间结构和正交关系

图像拼接算法原理 1

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

1. Homography

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

2. 曲面投影

3. 后续处理

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

线性方程组和极小化问题

最速下降法和共轭梯度法简述

共轭梯度法的理解

第一个关键点，坐标变换和子问题解耦

第二个关键点，子空间结构和正交关系

版权声明：原创作品，欢迎转载，但转载请以超链接形式注明文章来源(planckscale.info)、作者信息和本声明，否则将追究法律责任。

前置广告：多路视频实时全景拼接Demo可见我们的视频主页。

0. 引言

1. Homography

第二个关键点， $D_i$ 子空间结构和正交关系