月度归档：2017年03月

全景视频技术的产品化之路

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

甚嚣尘上的VR炒作终于在今年平静了，这大概意味着VR技术开始进入技术成熟度曲线的第三个时期：行业的公众关注进入低谷，人们开始冷静客观评估技术的适用范围和潜力，并逐步发现有效的经营模式。

技术成熟度曲线

VR时代的到来是不可避免的，或者说它已经到来，只是还没有推到大众面前。另外，真正具有想象力和冲击力的新技术乃是紧随其后的AR，这一点可能并不像公众预期的那样。这个时代需要由一系列扎实漂亮的产品撑起（不是概念，不是Demo,不是DIY，是产品），我们这次来谈谈全景摄像机的产品化之路上有哪些曲折和挑战。当然，全景摄像机本身并非仅限于VR应用，我们也要包括安防应用。

安防监控领域

泛泛来说有两种全景视频实时拼接方案，即前端（机内）拼接和后端（PC/手机）拼接。在安防领域也是如此。前端拼接直接由全景摄像机输出拼接完成的全景帧，具有很好的兼容性，可以直接像一台普通IPC一样接入旧有系统；而后端拼接是将全景摄像机看做独立的多路IPC,同时接入监控PC服务器，由PC完成实时拼接和监看。后端拼接的优势在于可以完成极高分辨率（目前我们的后端方案全景监控分辨率最高已经有9600万像素）的全景监控，但兼容性不好，需要将全景拼接SDK嵌入平台软件，不能做到“即插即用”。

从实现上来说，大概有如下几种：FPGA/DSP/CUDA/OpenGL/CPU. 前两种用于前端拼接，FPGA的开发和维护都有较高代价，CUDA和OpenGL方案具有最高的处理能力，CPU方案除非无法选择否则是应该排除的。在前端拼接方案里，还要考虑编码问题，全景帧动辄数千万的分辨率编码并不是一个简单问题。这里我们主要谈我们自己比较熟悉的CUDA/OpenGL方案。

安防监控领域对于全景摄像机有一些特殊需求。对于后端拼接全景，其拼接参数应当保存在设备之中，由设备传给平台软件完成实时监看的初始化流程，而平台软件上则对实时拼接的效率，全景模块与其他设备如球机的互动都有颇多要求，我们简单罗列如下。

拼接参数应该是一个很小(几k到几十k)的文件，方便写入设备及在网络上传输；
灵活的裁剪/融合算法。安防全景细分需求繁多，催生大量不同类型的设备，不同目数，不同镜头，不同摆位都会导致非常不一样的应用场景，算法需要在任何场景下都能够完美融合，不产生图像瑕疵。
全景拼接模块应当支持实时的子码流/主码流切换。平台软件实时监看几十上百路网络摄像机，区分子码流/主码流非常重要，这样在小视图模式下采用子码流，而大视图下自动切换到主码流，既保证了性能又保证了操作体验。
全景拼接模块在子码流输入下能够同时完成多达几十路的实时拼接播放。
应该有多种全景投影模式。除了常见的球面/柱面展开，碗形交互式展开在安防监控领域颇受欢迎，如图所示。各种投影模式之间应该能够实时切换。
ROI局域放大。
与球机联动，全景纵览全局，球机实现局域放大。
全景像素坐标到输入图像像素坐标的正反向投影。

碗形交互式全景

解决以上需求的任务并非平凡。简单一例，子/主码流实时切换中，除非子码流与主码流具有同样的视野，否则无法在不重新初始化算法的前提下完成切换，这要求算法具备瞬间初始化完毕的性能。同样的性能要求也出现在实时全景投影模式切换中。

从一个算法到一个成熟产品的道路是长远的。行业里很多学术型团队最终败在懂算法不懂软件工程，无法将一个Demo级的算法提升成一款结构良好，功能灵活，充分解决行业内需求的算法产品，令人扼腕。

对于前端拼接来说，要支持交互式全景类型如碗形、柱面等，同样也需要将一个全景播放模块嵌入平台软件，此种情况下，上述需求中除1、2、8外都仍然需要满足。

前端拼接技术的一个需要解决的问题是，以条带展开型全景作为全景帧类型做编码传输，如何节省带宽？将一个全景球展开为一个平面图像就如同将柚子皮拍扁在桌面上，总会像全球地图那样产生一个畸变，这是无法避免的。投影类型选择不好可能会导致相当大的畸变，比如球面两极地区一个像素被拉伸成一行像素。更好的投影类型应该是立方体展开

立方体展开

这一展开方式可以将畸变控制到很小的程度，但它一定程度上损失了条带全景图那种一览无余的直观性，需要特殊的全景播放器将它重新贴图到全景球或全景展开平面上才能够还原全景。当然也有其他采用更高级的数学方案设计的展开方式，这里不再提。

全景摄像机与全景直播

这里特意避免了提“VR全景”这一概念，因为严格来说VR全景和普通的全景摄像机并非同一概念，前者要求具有视觉深度感，后者只是个普通的2D曲面，沉浸感不强。但由于普通全景摄像机技术较前者简单，所以目前市面上大都为此类产品。

我的个人观点是，目前全景摄像机难以普及的一个关键是没有标准格式。并不像传统数码相机，全景输出格式杂乱无标准，全景视频播放器无法自动化决定采用何种投影类型播放，使得全景视频成了少数geek一族的玩具。但在真正的行业标准出现之前，让自己的产品对各种不同的输出格式都做好准备不失一个办法，而且不难。

这类产品中低端以前端拼的双鱼眼为主，高端以后端拼的多目摄像机为主，但迄今几乎没有很让人满意的产品出现。

双鱼眼方案的优势在于廉价且可以极小化拼缝。在所有可能的基于拼接算法的方案里，双鱼眼的拼缝是最小的。拼缝大小取决于多个摄像机投影中心的距离，摄像机的投影中心位置大致在sensor中心向后一个焦距远的地方，通常这是个很短的距离。理论上只有各个摄像机的投影中心重合于一点才能够产生出无缝的全景图，但这种情况下相机的体积需要压缩到极限，几乎不可达到，通常只能将尺寸压缩到极小以期更好的拼缝效果。除了双鱼眼方案，它是可以真的做到投影中心重合的。

所以，对于做基于拼接算法的全景摄像机的厂商，一个忠告是，将相机尺寸做小。

全景直播机似乎有很长一段时间卡在很高的软件授权费和拼接服务器价格上，但这是比较奇怪的，因为这一技术并不困难——至少在安防领域，四年前就已经有公司做到了上千万像素的全景监控。像安防领域一样，最高性能且具有很好平台兼容性的方案就是OpenGL方案，现在的显卡处理几千万像素的全景拼接融合如同砍瓜切菜，顺便搞个硬编码做推流是不难的——我们自己的技术在这方面早已验证过。

全景直播机通常并不是多个摄像机拼一块儿这么简单粗暴，它需要解决两个基本问题，一是摄像机之间的帧同步，一是摄像机之间的成像参数同步。前者保证人通过拼缝时不会出现消失又出现这种诡异效果，后者使得全景画面亮度、色彩具有一致性，不出现尖锐的过渡。

但实际上，我们并不真的需要成像参数同步。理论上，多个摄像机各自自动曝光，可以实现HDR（高动态范围）全景，因而目前在硬件上做成像参数同步只是一个过渡方案，将来为了生成HDR效果全景，这一机制是必然要废弃的。

要有更好的拼接质量，可以选择CUDA或OpenCL,它比OpenGL提供更多控制力，使得开发者可以采用更复杂的图像处理算法。我们目前就在基于CUDA尝试HDR全景算法的开发。

VR全景

VR全景是万众期待众望所归，出于不可描述之原因，这一技术被视作新时代的宅男福利。但一定要冷静！因为我们真还有很多技术问题要解决。

实现3D效果，目前主要有基于传统拼接算法拼左右眼全景图（参见我们的文章《DIY 3D全景摄像机》）和光流算法（Google Jump/Facebook Surround360等）两种。

基于拼接算法基本是没有前途的（所以我们直接做成了DIY教程-_-!）。这一方案的死结在于，3D全景中深度感最强的近景，正好是拼缝最大的，而且你不能够通过缩小设备尺寸来解决，因为它至少应该有人的瞳距（~6.2cm）那么大，否则你戴上眼镜后，会发现自己缩小了——周边的一切都大了一遭。

光流算法是目前给出效果最好的，光流刻画了两个图像的像素是如何对应的，算法利用光流来插值计算没有被相机所采集的光线之颜色，从而产生出完全无缝的全景效果。但目前效率不高，关键是光流本身的计算是相当繁重的，而且算法对于每队图像还需要计算正反向两个光流，再考虑上光流在时间轴上的一致性，带来了非常大的计算开销。

实现VR视频采集，本质上是通过有限个相机采集几个点上的物理光线，然后用这些光线来猜测、插值出其他空间位置上任意光线。这在计算机视觉领域早已研究多年（想想黑客帝国里的子弹时间镜头是怎么来的），这个方向叫做”Image-Based Rendering”.

理想自然是通过采集有限个点上的光线就能够计算出一个邻域上的光场。这一定程度上做得到，而且有很好的工作，但付诸应用仍然有距离。

所以，仅就目前的情况来说，基于光流算法来做后期，做高质量近景VR视频是没问题的，但想要直播，还得等等。

统计算法的理论与实例之一最大似然估计与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*}$

Mind Spectrum

Formerly Singularity Vision

月度归档：2017年03月

全景视频技术的产品化之路

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

安防监控领域

全景摄像机与全景直播

VR全景

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

Notes on Conjugate Gradient Method – 2

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