面对COVID-19, 我们能做些什么?

全球COVID-19死亡病例总数(红),每日新增死亡(蓝)

我们采用死亡数曲线来粗略观察。死亡总数曲线仍然是指数增长态势,约六天翻倍,作为其变化率,每日新增死亡也大致符合这一规律。当前的死亡数是64000+. 由于各国的控制措施和免疫人群的增长,新增病例曲线在未来从指数上压下来是肯定的,但尚不知何时。而由于病毒的高传染性,后面一个长尾的存在几乎是肯定的,我们在很长时间内将不得不保持高度警惕,全世界范围将对此付出巨大经济和社会代价。

如今的态势,我们唯一的希望,恐怕就在科学了。不论是疫苗还是有效药物,都将提供其他任何手段不能比拟的武器。作为一个科技从业者,我们也一直在关注,数学家/物理学家/工程师如何能在这场战役中充分发挥专业优势,为生物、医药工作者提供支援。本文即这样一个收集贴,希望后续中读者们能帮忙充实。

××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××

Folding@Home 项目,为COVID-19相关的蛋白质折叠研究提供CPU/GPU算力

该项目是斯坦福大学主持的,世界上最大的分布式计算项目,旨在为蛋白质折叠研究提供分布式算力,由全世界的玩家们共同捐献CPU/GPU算力。目前COVID-19相关的计算问题已经上线Folding@Home.
出于安全性考虑,客户端除官方下载外没有其他安装方式,这是很安全的一个软件。如果您决定加入,可参考该教程或自行google.

Folding@Home具有分组功能,玩家可以组成自己的战队并肩而战。我们特意创建了一个team,您可以选择加入。
team号:259252
team名:SingularityVision_CN
(到2020.4.18,小组排名已经进入前七千名,欢迎加入)

一点额外的tips.

软件会悄悄在你cpu空闲期间运行,如果您希望暂停之,可点击Pause键(直接退出是不行的,它会自己找时间悄悄启动)。

对于ubuntu下显卡无法使用的情况,可尝试如下步骤:
sudo apt install ocl-icd-opencl-dev
打开Configure -> Slots, 删掉所有slots.
expert, gpu设置为true.
重启,Configure -> Slots添加gpu slots,参数默认即可。

××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××××

更多开放性科研项目,期待你的推荐!

Notes on Representation Theory of Finite Groups

2020这个一言难尽的开局,一个月基本都在刷手机和失眠里度过。后来weibo炸号,心里反而轻松了不少,说出来没用的话,不说也罢。
我们来仰望星空吧。
群表示论是物理系学生能接触到的数学里最美的分支之一。读研时一个晚上突然发现物理学里那些常用的正交函数族、Fourier变换与群表示论间有非常漂亮的联系,大概是读书时印象最深刻的一次顿悟。这里贴出一篇过去的有限群表示论笔记,谨以最美好的数学来对抗这段糟糕的日子。

\section{有限群表示论}
本章试图对有限群表示论给出一个简约清晰的介绍。

有限群表示论的一个基本问题是有限群表示的分类与求解。这将会导致如下经典结果:

\begin{enumerate}
\item 有限群表示是完全可约的,且可酉化;
\item 正则表示中包含了全部不可约表示,且其中全部不可约表示的数目等于其等价类数目,重数等于其维数;
\item Schur引理。不可约模间没有非平凡非同构的模同态;
\item 正交关系的导出;
\item 由特征标正交关系得出的不可约表示判定与可约表示约化。
\end{enumerate}

\subsection{主要概念}

\begin{definition}[群表示] 群G的\textbf{表示}为一个同态\rho:G\to GL(V),其中V为域F上的有限维向量空间,这里若无特殊说明我们取F=C.\dim_{F} V称作\textbf{表示的维数};
\ker \rho平凡,则称该表示是\textbf{忠实的};若对所有g\in G,\rho(g)=\mathcal{E},则称表示为\textbf{平凡的}。
群表示可看作V上的G-模。
\end{definition}

\begin{definition}[可约与不可约] 一个表示是\textbf{不可约}的如果它不具有恰当的不变子空间(即不可约子表示);一个表示是\textbf{完全可约的}如果它可以分解为一系列不可约子表示的直和。
\end{definition}

\begin{definition}[正则表示] 考虑G上复函数构成的向量空间\mathbb{C}(G),它的基底为\{\mathrm{e}_{g}\}_{g\in G},其中\mathrm{e}_g表示一个在g上值为1,余皆为0的函数。
G可以通过\mathbb{C}(G)上的线性算子\lambda(g)这样作用在基矢上

(1)   \begin{eqnarray*} \lambda(g)(\mathrm{e}_{h})=\mathrm{e}_{gh} \end{eqnarray*}

则称\lambda(g)G的\textbf{左正则表示}。
\end{definition}

\begin{definition}[G-线性映射] 设GU,V上分别有表示\rho_1:G\to GL(U),\rho_2:G\to GL(V),且存在线性映射\sigma:U\to V,满足

(2)   \begin{eqnarray*} \sigma\circ\rho_1(g)(u)=\rho_2(g)\circ\sigma(u),\forall u\in U,g\in G \end{eqnarray*}

则称\sigma为一个\textbf{G-线性映射},或称\textbf{\sigmaG可交换}.若\sigma可逆,则上式可写为

(3)   \begin{eqnarray*} \rho_1(g)=\sigma^{-1}\circ\rho_2\circ\sigma \end{eqnarray*}

此时称\sigma为\textbf{群表示的同构}。
用模的语言来说,G-线性映射正是一个模同态。
\end{definition}

\subsection{酉表示与可约性}
\begin{definition}[酉表示] 向量空间V上的表示\rho被称作\textbf{酉}的,如果V装备了一个在G作用下不变的厄米内积\Braket{\vert}

(4)   \begin{eqnarray*} \Braket{v_1 \vert v_2}=\Braket{\rho(g)v_1 \vert \rho(g)v_2},\forall v_1,v_2\in V,g \in G \end{eqnarray*}

表示称作\textbf{可酉化的}如果V可以装备这样一个内积。

一个表示是酉的当且仅当同态\rho:G\to GL(V)落在酉群U(V)之中。该条件也可以表述为\rho^{\dagger}(g)=\rho(g^{-1}).
\end{definition}

\begin{lemma}[] 设V为群G的酉表示而WV中的一个不变子空间,则W的正交补W^{\bot}也是一个不变子空间。
\end{lemma}
\begin{proof}
由于G保持酉表示空间中的内积不变,于是对于w^{\prime}\in W^{\bot},它被任意g\in G作用后的结果满足

(5)   \begin{eqnarray*} \Braket{\rho(g)(w) \vert \rho(g)(w^{\prime})}=\Braket{w \vert w^{\prime}}=0,\forall w\in W \end{eqnarray*}

由于\rho(g)(w)\in W,于是\rho(g)(w^{\prime})\in W^{\bot}.
\end{proof}

\begin{theorem}[酉表示的完全可约性] 设\rho:G\to GL(V)为任意群G上的有限维酉表示,则该表示完全可约,并且可约化为一系列不可约酉表示的直和。
\end{theorem}
\begin{proof}
首先,根据上述引理,V可以不断迭代地正交分解为一系列不变子空间的直和,这个分解必然终止于一系列不可约子空间。

然后,V上的厄米内积限制在各不可约子空间上,显然G仍然保持这些子空间中的内积不变,因而V被分解成了一系列不可约酉表示的直和。
\end{proof}

\subsection{有限群表示的完全可约性}
\begin{theorem}[有限群表示的酉性] 有限群的有限维复表示必可酉化。
\end{theorem}
\begin{proof}
只需找到一个G-不变的内积即可。

以下式从\Braket{\vert}定义出一个新的G-不变内积\Braket{\vert}^{\prime}

(6)   \begin{eqnarray*} \Braket{w\vert v}^{\prime}=\frac{1}{|G|}\sum_{g\in G}\Braket{gv \vert gw} \end{eqnarray*}

显然具有G-不变和正定性质,因而是一个G-不变的内积。
\end{proof}

\begin{theorem}[有限群表示的完全可约性] 有限群的有限维复表示必完全可约。
\end{theorem}
\begin{proof}
有限群的有限维复表示必可酉化,而酉表示必完全可约。
\end{proof}

下面给出上述定理的第二种证明,以期从不同的角度来理解这一结果。
\begin{lemma}[] 设V为群G在特征不整除|G|的域F上的有限维表示,则每个不变子空间W\subset V都有一个不变的补空间W^{\prime},使得W^{\prime}\oplus W = V.
\end{lemma}
\begin{proof}
\pi:V\to WVW的投影算子,仍然通过在G上平均的手法构造一个新的G-不变投影算子\pi^{\prime}

(7)   \begin{eqnarray*} \pi^{\prime}\equiv \frac{1}{|G|}\sum_{g\in G} g\circ \pi \circ g^{-1} \end{eqnarray*}

容易验证\pi^{\prime}有如下性质
\begin{enumerate}
\item \pi^{\prime} \circ h = h \circ \pi^{\prime},\forall h\in G
\item \pi^{\prime}=\id \quad on \, W
\item \img \pi^{\prime}=W
\end{enumerate}
可见\pi^{\prime}是一个V\to WG-线性映射,或者说一个模同态。

由第一条性质知道W^{\prime\prime}\equiv \ker \pi^{\prime}是一个不变子空间;由第二条性质知道W^{\prime\prime} \bigcap W = 0;由第三条性质得\dim W + \dim W^{\prime\prime} = \dim V,
W^{\prime\prime}\oplus W = V.
\end{proof}

\subsection{Schur引理}
\begin{theorem}[Schur引理] 设\rho_V,\rho_W为群G的有限维不可约表示,设\sigma:V\to W是一个G-线性映射,则有
\begin{enumerate}
\item 若两表示不等价,则\sigma=0;
\item 当V=W,\rho_V=\rho_W时,\sigma=\lambda \id
\end{enumerate}
\end{theorem}\label{SchursLemma}
\begin{proof}
\sigma是一个G-线性映射,因此\ker \sigma\img \sigma均为不变子空间,但考虑到V,W都是不可约表示,不包含任何非平凡子空间,因此只有如下两种可能性
\begin{enumerate}
\item \ker \sigma=0,\img \sigma=W,\sigma为表示同构;
\item \ker \sigma=V,\img \sigma=0,\sigma为零映射。
\end{enumerate}
因此第一条得证。

V=W时,由于\sigmaG对易,因此\sigma的本征子空间E_{\sigma}是不变子空间,又由V的不可约性有E_{\sigma}=V,因此\sigma=\lambda \id.
\end{proof}

\subsection{特征标与第一正交关系}
\begin{theorem}[] 设\rho_V,\rho_W为群G的有限维不可约表示,设\sigma:V\to W是一个线性映射.我们由\sigma构造映射\tilde{\sigma}

(8)   \begin{eqnarray*} \tilde{\sigma}=\frac{1}{|G|}\sum_{G}\rho_W(g) \circ \sigma \circ \rho_V(g)^{-1} \end{eqnarray*}

\tilde{\sigma}满足如下性质
\begin{enumerate}
\item 若两表示不等价,则\tilde{\sigma}=0;
\item 当V=W,\rho_V=\rho_W时,\sigma=\lambda \id,\lambda = \frac{\tr \sigma}{\dim V}
\end{enumerate}
\end{theorem}\label{SchurOrth}
\begin{proof}
易证

(9)   \begin{eqnarray*} \rho_W(g) \circ \tilde{\sigma} \circ \rho_V^{-1} = \tilde{\sigma} \end{eqnarray*}

\rho_W(g) \circ \tilde{\sigma} = \tilde{\sigma}\circ \rho_V,可见\tilde{\sigma}是个G-线性映射,从而由??得出两个断言。
关于\lambda的取值,可以这样计算出来

(10)   \begin{eqnarray*} \tr \tilde{\sigma}&\xlongequal{1}& |G|^{-1}\sum_G \tr \rho(g)\ \circ \sigma \circ \rho(g)^{-1}\\ &=&|G|^{-1}\sum_G \tr \sigma=\tr\sigma\\ &\xlongequal{2}& \tr \lambda \mathcal{E} = (\dim V)\lambda\\ &\Rightarrow& \lambda=\frac{\tr \sigma}{\dim V} \end{eqnarray*}

\end{proof}

考虑群G上的复函数空间,我们在其中定义任意两个复函数\phi,\psi的内积

(11)   \begin{eqnarray*} \Braket{\phi\vert\psi} \equiv \frac{1}{|G|}\sum_G \bar{\phi}(g)\psi(g) \end{eqnarray*}

\begin{definition}[特征标] 一个表示\rho的\textbf{特征标}\chi_{\rho}:G\to \mathbb{C}定义为

(12)   \begin{eqnarray*} \chi_{\rho}(g)\equiv \tr \rho(g) \end{eqnarray*}

\end{definition}

可以证明特征标是一个完备不变量,可以在同构的意义上对不可约表示进行分类。
\begin{theorem}[特征标的基本性质] 特征标具有以下性质
\begin{enumerate}
\item 特征标是类函数,\chi(g)=\chi(hgh^{-1}),\forall g,h \in G
\item \chi(1)=\dim V
\item \chi(g^{-1})=\bar{\chi}(g)
\item 对于任意两个表示V,W,\chi_{V\oplus W}=\chi_V+\chi_W,\chi_{V\otimes W}=\chi_V * \chi_W
\item 对于对偶表示V^{*},有\chi_{V^{*}}(g)=\chi_V(g^{-1})
\end{enumerate}
\end{theorem}
\begin{proof}
显然。
\end{proof}

\begin{theorem}[第一正交关系] 特征标满足如下正交关系
\begin{enumerate}
\item 若\rho_V是不可约的,则\Vert\chi_{V}\Vert=1;
\item 若\rho_V,\rho_W不可约且不同构,则\Braket{\chi_W\vert \chi_V}=0.
\end{enumerate}
\end{theorem}
\begin{proof}
作如下计算

(13)   \begin{eqnarray*} \Braket{\chi_V\vert \chi_W}&=&\frac{1}{|G|}\sum_G \bar{\chi}_V(g)*\chi_W(g)\\ &=&\frac{1}{|G|}\sum_G \chi_W(g) * \chi_V(g^{-1})\\ &=&\frac{1}{|G|}\sum_{G,i,j} \Braket{w_i\vert \rho_W(g) \vert w_i} \Braket{v_j \vert \rho_V(g^{-1}) \vert v_j}\\ &=&\sum_{i,j} \bra{w_i} \frac{1}{|G|}\sum_{G} [\rho_W(g) \ket{w_i} \bra{v_j} \rho_V(g^{-1})] \ket{v_j}\\ \end{eqnarray*}

由??,令\sigma = \ket{w_i} \bra{v_j},可知V,W不等价时上述内积为零,若V,W等价则有

(14)   \begin{eqnarray*} &=&\sum_{i,j} \frac{1}{\dim V}\Braket{v_i\vert v_j}\\ &=&1 \end{eqnarray*}

\end{proof}

于是我们有如下结论。
\begin{theorem}[] 对于有限群的复表示,我们有
\begin{enumerate}
\item 不可约表示V在表示W中的重数等于\Braket{\chi_V\vert \chi_W};
\item 两表示同构当且仅当它们具有相同的特征标;
\item \Vert \chi \Vert^2总是一个整数,V不可约当且仅当\Vert \chi_V \Vert^2=1;
\end{enumerate}
\end{theorem}
\begin{proof}
考虑表示W的不可约分解,我们最终要证明如下等价关系

(15)   \begin{eqnarray*} W=\oplus_i V_i^{\oplus m_i} \iff \chi_W=\sum_i m_i\chi_{V_i} \end{eqnarray*}

\label{RepRingCharCorrespondence}
证明上式的正向箭头是显然的,由于不可约表示特征标的正交性,不可约表示V_kW中的重数可这样计算

(16)   \begin{eqnarray*} m_k=\Braket{\chi_{V_k}\vert \chi_W} \end{eqnarray*}

这就证明了第二条.

考虑??的逆向箭头,对于任意不可约表示V_i,都可以通过\Braket{\chi_{V_i}\vert \chi_W}计算出它在W中的重数,从而可以根据特征标重建出W的不可约分解,因而逆向箭头成立,第二条命题也得证。

对于不可约表示V_i,必有\Vert \chi_{V_i} \Vert = 1.可以算出

(17)   \begin{eqnarray*} \Vert \chi_W \Vert^2=\sum_i m_i,m_i\in \mathbb{N} \end{eqnarray*}

可见\Vert \chi \Vert^2必为整数,\Vert \chi \Vert = 1时必然意味着该表示不可约。第三条证毕。
\end{proof}

\subsection{正则表示}

\begin{proposition}[] 对于正则表示\lambda,我们有如下结果
\begin{enumerate}
\item

(18)   \begin{eqnarray*} \chi_{reg}(g)= \begin{cases} |G|,&if\quad g=e\cr 0,& otherwise \end{cases} \end{eqnarray*}

\item 任意不可约表示\rho_V\lambda中的重数等于\dim V;
\item 全部不可约表示的维数平方和等于群的阶,即\sum_V \dim^2_V=|G|.
\end{enumerate}
\end{proposition}
\begin{proof}
第一条: 正则表示下除了单位元e,任何群元都没有非零对角元;

第二条: \Braket{\chi_V \vert \chi_{reg}}=\frac{1}{|G|}\chi_V(e)\chi_{reg}(e)=\chi_V(e)=\dim V;

第三条: 显然。
\end{proof}

现在我们已经证明了不可约表示特征标是一组正交归一向量,接下来我们可以进一步证明它们在类函数空间里是完备的。
证明的基本思路就是假设类函数\phi正交于所有的不可约表示特征标,说明它必为零,从而不存在这样非平凡的\phi

对于任意函数\phi:G\to \mathbb{C}G的任意表示\rho_V,定义如下V的自映射\rho(\phi)

(19)   \begin{eqnarray*} \rho(\phi)\equiv \sum_G \bar{\phi}(g)\rho(g) \end{eqnarray*}

上述自映射有如下性质
\begin{lemma}[] \phi是个类函数当且仅当在任何表示下\rho(\phi)G对易。
\end{lemma}
\begin{proof}
首先假设\phi是个类函数,即\phi(h^{-1}gh)=\phi(g),\forall g,h\in G,则有

(20)   \begin{eqnarray*} \rho(h)\rho(\phi)\rho(h^{-1})=\sum_G \bar{\phi}(g)\rho(hgh^{-1})=\sum_G \bar{\phi}(h^{-1}gh)\rho(g)=\sum_G \bar{\phi}(g)\rho(g)=\rho(\phi) \end{eqnarray*}

然后假设任何表示下\rho(\phi)G对易,取\rho_V为正则表示\lambda,有

(21)   \begin{eqnarray*} \lambda(h)\rho(\phi)(\mathrm{e}_1)&=&\sum_G \bar{\phi}(g)\mathrm{e}_{hg}=\sum_G \bar{\phi}(h^{-1}g)\mathrm{e}_g\\ \rho(\phi)\lambda(h)(\mathrm{e}_1)&=&\sum_G \bar{\phi}(g)\mathrm{e}_{gh}=\sum_G \bar{\phi}(gh^{-1})\mathrm{e}_g \end{eqnarray*}

由对易假设知两式相等,从而有\phi(h^{-1}g)=\phi(gh^{-1}),证毕。
\end{proof}

\phi是类函数且\rho_V不可约时,Schur引理告诉我们\rho(\phi)是个常数映射.我们取迹以计算它

(22)   \begin{eqnarray*} \tr \rho(\phi)&=&\sum_G \bar{\phi}(g)\chi_V(g)=|G|\Braket{\phi \vert \chi_V}\\ \Rightarrow \rho(\phi)&=&\frac{|G|}{\dim V}\cdot\Braket{\phi \vert \chi_V}\cdot\id \end{eqnarray*}

\begin{theorem}[] 不可约表征标是类函数空间的一组完备基。
\end{theorem}
\begin{proof}
假设类函数\phi正交于所有的不可约表示特征标,则在任何不可约表示或完全可约表示上有\rho(\phi)=0,

考虑正则表示,

(23)   \begin{eqnarray*} \rho(\phi)(\mathrm{e}_1)=\sum_G \phi(g)\mathrm{e}_g=0 \Rightarrow \phi(g)=0 \end{eqnarray*}

证毕。
\end{proof}

\subsection{第二正交关系}
我们考虑第一正交关系,有

(24)   \begin{eqnarray*} &&\frac{1}{|G|}\sum_G \bar{\chi}_r(g)\chi_s(g)=\delta_{rs}\\ &\Rightarrow&\frac{1}{|G|}\sum_C |c|\bar{\chi}_r(c)\chi_s(c)=\delta_{rs}\\ &\Rightarrow&\sum_C \bar{T}_{rc}T_{cs}=\delta_{rs},T_{rc}\equiv \sqrt{\frac{|c|}{|G|}}\chi_r(c) \end{eqnarray*}

可见T是一个酉矩阵,第一正交关系正是说该矩阵的行向量相互正交。我们知道酉矩阵的行与列均正交,于是可以得出列向量的正交关系

(25)   \begin{eqnarray*} \sum_R \frac{|c|}{|G|}\bar{\chi}_r(c)\chi_r(c^{\prime})=\delta_{cc^{\prime}} \end{eqnarray*}

其中R为所有的不可约表示集合。

update.

或许熟悉我们的朋友已经知道,我们大概有两年没有在全景技术方面继续展开新的工作了。这里的主要原因是,我们跳进了另一个坑,算法交易,并且作为一个袖珍团队,很快就被吸干了精力。另一方面,像全景这样的技术,通过技术合作或者技术服务来盈利并非易事,你很难确保知识产权不被侵犯(企业跟你边“合作”边“学习”,事了另起炉灶留你干瞪眼儿这种事儿清华计算机系的老师们都没少遇到),另外整个产品化周期过长,自身团队之外的种种不可控因素过多,都成为要考虑的问题。最初选择这种形式来运营只有一个原因,就是自由和独立。保持团队有限精力不被过多侵蚀,保持做“闲事儿”的权利。如今跳了另一个坑,也无非是这样一个动机。

由于两年间精力都在构筑我们的算法交易平台,我并没有太多时间来考虑全景技术的未来。直到后来看到海思推出了芯片解决方案,我心里一沉,感觉14年时就曾经在担心的事情,终于还是来了。去年我开始考虑将一部分未曾做落地应用,单纯玩儿出来的算法做开源,但并没有足够多时间来整理这些工作,迄今仍未实现。

我们经常嘲笑自己,埋头在家做了好多东西,很多都是工业级强度的东西,都扔在硬盘里睡大觉,就又跑去搞别的新奇玩意儿。这是一个纯技术团队的很大缺陷,而我们直到今天也没(lan)有(de)克服。

所以今天对于我们全景技术的未来做一些说明。

首先,之前的合作我们一直在持续,产品也一直都稳定量产着,有兴趣的客户可以联系深圳华途数字了解产品具体情况。这些高清晰度安防全景产品应用了我们的量产标定技术以及实时拼接融合技术,由我们提供二次开发SDK,算法质量顶级(量产三年来大概只fix了1次较严重bug[在某些型号的机器上,某个全景投影类型和别的类型混淆,已经完全通过软件升级解决],加上我不记得的小bug,总量应该小于5个),您可以信赖。

然后,我们会尽可能花空闲时间整理过去的未公开技术并完成开源。至于工业级代码是否开源,取决于今后的具体情况。

最后,目前已经产品化的量产标定技术、实时拼接融合技术,如果您是安防或VR业内top5的靠谱公司,欢迎来聊聊技术收购(客官我们只卖术不卖身吆)。如果您是海思这样的高富帅芯片厂商,请直接来谈朋友(求你了)!如果你们都不来,我就择吉日都开源掉,我们这样的geek还是适合做开源。

至于这个网站,以后会自我放飞成为一个geek的后花园。

统计算法的理论与实例之二 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*}

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

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

版权声明:原创作品,欢迎转载,但转载请以超链接形式注明文章来源(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方案。

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

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

碗形交互式全景

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

从一个算法到一个成熟产品的道路是长远的。行业里很多学术型团队最终败在懂算法不懂软件工程,无法将一个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算法

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

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

\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*}

DIY 3D全景摄像机

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

目前市面上的VR全景多是二维,没有深度感,若真想拥有身临其境般的体验,深度信息必不可少。诺基亚的OZO,Google的Jump,Facebook开源的Surround 360,都是为3D全景而设计。OZO设备8个鱼眼售价高达三十多万,Surround 360搭载的是Point Grey的相机,硬件成本二十多万,Jump也要搭载十几个GoPro,硬件成本少说也要几万,普通玩家真心想玩也要思量的。如果不追求那么高大上,其实自己就可以DIY出一台3D全景相机。

组件:

千兆交换机+网络摄像头模组+180°鱼眼镜头+线材

结构:

要实现深度感,结构是关键。

相机可以理解为光线采集设备,采集到的光线与成像平面的交点即像点。

通常的二维全景要求采集到的所有光线汇于中心点,即视点,以视点为中心的球面或圆柱面为成像面,所有光线交于成像面形成全景图,如下图(a)。二维全景相机要求所有相机共中心摆位,即所有相机的光轴相交于视点。

3D全景即左右眼各对应一个全景图。

两只眼睛分别对应两个不同的视点位置,当转头360度时,两只眼睛转过的轨迹即一个以瞳距为直径的圆,称之为Viewing Circle,3D全景要求采集到的所有光线相切于此圆。左右眼采集到的光线分别与成像面相交形成左眼全景图和右眼全景图, 如下图(b)(c)。

两只眼睛所在的视点位置投影出的图像称为一个立体对,如图(b)左眼光线1和图(c)右眼光线1即可看成一个立体对,同理左眼中的光线2,3,4,5等与右眼中的2,3,4,5等分别构成不同的立体对。

该图引自文献:Stereo Panorama with a Single Camera.

至于相机摆位一般有两种方案,如下图:

该图引自文献:Jump: Virtual Reality Video.

切向摆位如上图(a):每个相机的光轴相切于Viewing Circle,此种方案一半的相机用于左眼全景(图(a)中绿色相机),而另一半的相机用于右眼全景(图(a)中红色相机)。

径向摆位如上图(b): 每个相机的光轴沿Viewing Circle半径方向,此种方案不区分左右眼相机,每个相机都对两眼的全景图有贡献,因此对每个相机水平视场角有更高要求:R越小,要求每个相机的水平视场角越大。

一般R设计比较大时采用径向摆位,R较小时采用切向摆位。

切向摆位最简单的结构设计即正多边形,每条边上放置两个camera,其sensor中心的距离设为瞳距。如果镜头的视场角足够大,可以设计一个正三角形,用六个camera来实现3D全景。本文介绍的是正四边形八个camera的方案,用Solidworks设计一个简单的支架,预留出上camera的安装孔位。

结构设计及3D打印: 

组装:

效果:

原始视频截图

上下3D格式

VR眼镜观看3D效果

目前能实现3D全景的技术无外乎几种:

  • 拼接方案。左右眼的视野分别做拼接融合以达到3D全景的效果,这种方案最简单,可以实时化,其缺点是拼缝难消除。
  • 光流方案。Google Jump以及Facebook开源的算法即此方案,能很好的消除拼缝,但实时化比较困难,适合做后期处理。
  • 光场重建。用有限个相机重建光场,给出真正的3D效果,使用户拥有更多活动自由度,这是终极的VR视频采集方案,但即便在理论上也有很多困难之处。

本文中所用的实时3D全景拼接软件是奇点视觉实时全景拼接方案,演示中为2700万像素实时拼接融合效果。我们目前可以轻易的实现超高分辨率实时拼接和直播,但仍有以下几个问题:

  1. 硬件使用了安防用网络摄像头模组,不具备同步曝光功能,因此全景图明暗不均较严重;
  2. 近景拼缝明显。3D全景要求设备直径至少等于人眼瞳距,但过大的直径容易导致更严重的拼缝。这是基于拼接方案做3D全景的一个本质困难。

但至少到目前,拼接方案仍然是唯一能做到低开销高分辨率直播的。我们正在研发新一代的实时光流/光场重建算法,希望能解决高质量近景VR直播的问题。

奇点视觉是一个致力于计算机视觉技术的研发和产品化的团队,专注于算法,为了发挥我们的优势,没有去做产品,而是给有能力做产品的公司提供技术解决方案。过去两年,我们专注于安防全景,毫不夸张的说,安防全景技术我们做到了世界顶级水平,已经由客户厂商实现产品化并销往海内外。目前我们致力于把全景技术迁移到VR应用中,并作技术升级,希望能够为更高质量的VR内容贡献一份力量,敬请关注奇点视觉。

全景视频相关需求征集

首先说一下我们目前的情况。到今天,我们的全景技术已经非常完善的产品化,性能、质量和稳定性都趋于收敛,可以不夸张的说,这一技术目前处于世界上领先的地位。

我们有目前最高的实时拼接融合性能,在普通笔记本上可以做到数千万像素的实时拼接。在这一实时拼接模块的驱动下,实现超高清全景直播只需要一台普通中端笔记本。基于这一模块,我们在安防领域实现了4000w像素实时全景监控,且在笔记本上,子码流模式下可以同时进行多达30个全景实时浏览。该模块其余特征包括:

  1. 子/主码流动态切换机制。小窗口时用子码流,放大观看用主码流,后台实时切换,非常适合安防监控领域。
  2. 多种全景类型。包括但不限于交互式球面/半球/柱面,球面180/360展开,柱面180/360展开,小世界,双目上下/左右格式的3D全景等等。
  3. 可以将全景图实时播放,也可以实时回传给图像分析算法。
  4. 可以将全景像素坐标、输入视频像素坐标、物理空间实际方向互相转换。在安防领域可以以此实现全景/球机联动,全景/全景联动实现局域放大等;也可以籍此实现一些增强现实效果。

我们有一套强大的量产标定算法,支持对任意多路,任意多不同类型的摄像机进行最优拼接参数标定,支持无显著重合区情况下的标定,拼缝达到理论极小值,固定步骤可以得到稳定标定结果,返工率极低

以上全部技术完全自主知识产权。不依赖于PTGUI/OpenCV等第三方实现。也正因为如此,我们对技术中全部细节都具有掌控力,不会因为一个问题出在第三方代码中而束手无策。

******************************************************************

然后说一下此文的真正目的,即收集各行业对于全景视频的需求

VR这个风口吹起来了漫天的塑料袋,但没有看到太多扎实漂亮的产品。我们希望做一款这样的产品。所以我们希望收集您对于全景视频最迫切的渴求,以此决定接下来的方向。您只需简略说明在您这里全景的用途,以及对它的一些特殊要求即可,可以在本站留言,也可发邮件:

planckscale1729@163.com

QQ是397692433, 但不太经常上,可能没法及时相应您:-P

开卖了

Shenzhen Civic Center

Shenzhen Civic Center

我们和客户合作研发的安防类多目全景网络摄像机今年上半年已经完成产品化,开卖了。

目前的产品包括4目180度,5目半球,8目半球,以及8目车载全景等,目前最高分辨率为2400W像素,30fps,很快将会支持到4000W像素(update: 目前已经支持到9600w像素)。

产品支持多种全景类型(沉浸式半球/各种投影类型的条带展开),浏览时多码流实时切换,对PC端硬件要求低,最低支持到OpenGL 2.1.

欲知更多细节,欢迎来电。

1500w像素全景监控演示视频

全景视频技术安防行业应用案例

5 6 13 14