TeXLive + TeXstudio 配置与使用

LaTex 是学术论文写作的标准工具,其排版好,可控制性强同时格式调整容易书写更快。TeXLive + TeXstudio 是一个比较常见的搭配(当然 TeXmaker 也是可以的,并没有太大区别,TeXstudio 高级功能更多一点),这里面简单总结下其安装配置过程。

0 系统环境

Ubuntu 16.04

1 安装步骤

在 Ubuntu 上安装这一组合非常容易仅需两行命令即可(分别安装 TeXLive 包和 TeXstudio 包):

2 中文设置

打开 TeXstudio,菜单 Options -> Configure TeXstudio 修改 Build -> Default Compiler 为 XeLaTex:

同时对于中文文档需要在开头引入 ctex 包,例如:

3 编译运行

点击 图标输入以下内容:

点击 图标进行编译。

如果编译无误则会在右边正确显示:

PS:如果没有显示预览窗口可以点击 图标打开。

4 保存 PDF

其实并没有一个菜单进行导出 PDF 的操作(或者我没有找到),想要保存刚才 Tex 文件编译生成的 PDF,只需在运行 Build & View 后在 Tex 文件所在目录找到同名的 PDF 文件即可。例如上面的文件保存成 hello.tex ,编译的文件就是 hello.pdf:
hello

5 常用快捷键

TeXstudio 有点像编译工具,有些快捷键也很类似 Visual Studio:

功能快捷键说明
Build & ViewF5
CompileF6
ViewF7
CommentCtrl + T
UnCommentCtrl + U
SaveCtrl + S
Inline MathCtrl + Shift + M$ ... $
Display MathAlt + Shift + M\[ ... \]
Numbered EquationCtrl + Shift + N\begin{equation}
...
\end{equation}

参考地址

[1] 官网:http://www.texstudio.org/

论文笔记:Deep ChArUco: Dark ChArUco Marker Pose Estimation

ChArUco 在相机标定中比较常用,同时在一些 AR 应用中也有涉及。这篇文章是针对 ChArUcho 这个矩阵二维码使用 Deep Learning 方法进行姿态估计,取得了比传统方法更佳鲁棒的结果(似乎可以为标志检测之类)。

文章主要贡献在于:

1、两个网络 ChArUcoNet 和 RefineNet:前者用于定位 ChArUco 的角点坐标,后者用于对角点坐标进行亚像素修正
2、使用仿真数据进行自动标注与训练的方式

1 传统方法

ChArUcho 示意图如下:

ChArUco 就是棋盘格和 ArUco 码的结合,整体是一个 5×5 的期盼形状,总共有 16 个 ArUco 码代替了棋盘格中的白色方格,其中的 ID 分别为 0-15。

传统方法j基于图像处理,识别黑白色块并进行定位,这里不再详述。

2 Deep ChArUco

类似于 SuperPoint (文章作者也是 SuperPoint 作者),Deep ChArUco 希望使用关键点监测的方式完成一整套 ChArUco 定位和姿态估计的流程。与很多 End-to-End 方法不同,Deep ChArUco 方法主要通过网络回归关键点坐标,然后通过 PnP 解算姿态,相对来说更加鲁棒和可信。

算法整体结构如下:

图像首先通过 ChArUcoNet 进行粗定位获得关键点,然后通过 RefineNet 进行亚像素精度修正,最后通过 PnP 解算姿态。

2.1 关键点定义

使用 Deep Learning 学习关键点,首先要定义出关键点。文中罗列出三种关键点的定义方法:

对于 a) 边缘点可能比较容易收到干扰;对于 b) 中间的点其实没有明显的特征可以学习;最终作者选择了 c)作为定义的关键点。

2.2 ChArUcoNet

ChArUcoNet 部分,作者采用了类似 SuperPoint 的网络结构,其中 Encoder 部分是一个类似 VGG 的结构。详细可以参照 SuperPoint 的代码:

降采样的部分主要依靠 MaxPooling 实现,总计经过了三个 MaxPooling 操作因此最终输出的 feature map 是原始大小的 1/8。该网络总计有两个输出:

1)\(\mathcal{X} \in \mathbb{R}^{H_{c} \times W_{c} \times 65}\) 是 Detector Head,其中输出的 65 个维度分别代表该 feature 所对应的 8×8 感受野是 keypoint 的概率,额外一个输出表示该 8×8 感受野没有 keypoint。

2)\(\mathcal{C} \in \mathbb{R}^{H_{c} \times W_{c} \times (N_c + 1)}\) 是 Class Head,其中输出的 (N_c) 表示 ChArUco 中的ID,在本文中 \(N_c = 16\),额外一个输出表示该 8×8 感受野没有 keypoint。

显然 ChArUcoNet 参数量比较少,是一个非常适合实时运行的网络,在作者的测试中 320×240 分辨率的图像作为输入,网络在 1080 显卡上可以跑到 100fps。

2.3 RefineNet

由于我们这里要做的是比较精确的定位操作,仅有证书像素级别的关键点定位是非常不够的,通常在 SLAM 或者 AR 中,我们都会采用 subpixel 这样的亚像素操作来进一步提升关键点定位的精度,在深度学习的方法中也不例外。这里作者设计了一个 RefineNet 用于回归每个关键点的亚像素精度。

RefineNet 中输入为关键点周围一个 24×24 的 patch,输出仍然为一个是否为关键点坐标的分类,由于每个关键点位置本身代表了一个 8×8 的原始图像块,在亚像素精度中,我们将每个像素再对应 2×2 的区块,也就是精确到了 0.5 像素精度,这样总计就变成了一个 16×16=4096 的分类输出(softmax)。

可以说整体设计上是比较简单的,规避了常见的上采样输出 heatmap 的方式,改成了多个分类问题,loss function 应该都是交叉熵。总体参数量得到了很好的压缩。

3 训练数据

对于 ChArUcoNet 部分,作者使用了实际采集的数据集(预计真值应该是传统方法或者标注获得)同时加入了很多的 data augmentation。具体参数如下:

对于 RefineNet 部分,作者还是和 SuperPoint 一样采用了仿真数据,这里面没有谈细节,但是应该是对于任意几何形状取亚像素的真值。作者给出了一个示例可见一斑:

4 实验

很明显由于训练数据的大量增强,深度网络是可以比传统方法取得更好的适应能力的,特别是在比较恶劣情况下,作者做了一些实验,证明了方法的鲁棒性。

抗模糊实验

低曝光实验

其他还做了一些例如速度、RefineNet 和 conerSubPix 对比不一而足。

整体上很明显取得了相比传统方法好很多的效果。

个人小结

这篇文章虽然是对应 ChArUco 但实际上只要任何可以定义关键点的 Marker 方法都是可以同样方法来操作的。由于没有采用很多 heatmap 方式输出关键点坐标的方法,整体上网络参数较小非常适合实时。另一个创新是使用了一个分类问题的 RefineNet 网络进行亚像素精度的校准。

猜测整篇文章都采用分类来解决 keypoint 问题的原因一个是分类问题更容易训练,一个是数据也更好标注。

感觉可以对于自定义的标志牌识别使用这一网络进行学习。

相关材料

论文笔记:Deep Closest Point: Learning Representations for Point Cloud Registration

DCP 是一篇基于 Deep Learning 来解决 ICP 问题的,其中 Deep Learning 部分主要用于做匹配,后端仍然沿用 SVD 的方法。比很多 MLP 直接出 Pose 的合理,也取得了更好的效果。在与传统方法例如 Go-ICP 以及深度学习方法 PointNetLK 的对比中,都取得了一定的优势。

1 经典 ICP 问题

这一部分在之前的论文笔记中已经有了比较详细的阐述,参见:使用 SVD 方法求解 ICP 问题

2 Deep Closest Point

参照经典的使用 SVD 求解 ICP 问题的流程,我们采用神经网络提取特征,并且使用注意力机制,最终使用一个可导 SVD layer 进行求解(在 PyTorch 和 Tensorflow 都提供了这样的 Layer)。整体网络结构框架如下图所示:

2.1 特征提取

我们定义点云 \(\mathcal{X}=\{\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{i},\ldots,\boldsymbol{x}_{N}\}\subset\mathbb{R}^{3}\) 和 \(\mathcal{Y}=\{\boldsymbol{y}_{1},\ldots,\boldsymbol{y}_{j},\ldots,\boldsymbol{y}_{M}\}\subset\mathbb{R}^{3}\) 是需要匹配的两对点云,为简便直接考虑 M=N 的情况。定义在网络最后一层之前(第 L 层)提取的局部特征分别为 \(\mathcal{F}_{\mathcal{X}}=\{\boldsymbol{x}_{1}^{L},\boldsymbol{x}_{2}^{L},…,\boldsymbol{x}_{i}^{L},…,\boldsymbol{x}_{N}^{L}\}\) 和 \(\mathcal{F}_{\mathcal{Y}}=\{\boldsymbol{y}_{1}^{L},\boldsymbol{y}_{2}^{L},…,\boldsymbol{y}_{i}^{L},…,\boldsymbol{y}_{N}^{L}\}\)。

在点云特征提取部分,作者分析了两个网络,一种是 PointNet,一种是 DGCNN(当然 DCP 作者也是 DGCNN 作者自然要分析自己的)。这一部分参见:论文笔记:Dynamic Graph CNN for Learning on Point Clouds,已经对二者的区别进行了分析。简单说就是 PointNet 对于每个点进行独立的特征提取没有考虑点和周围点的关系,而 DGCNN 使用 knn 把邻域的信息包含了进来,更符合特征提取的要求,也会取得更好的匹配结果。

2.2 注意力机制

本文的 Embeding 和 Attentation 都是借鉴了 NLP 领域的一些概念。

Embeding

Embeding 可以简单理解为将一个数据映射到另一个空间形成新的表达。本文将 2.1 中 \(\mathcal{F}_{\mathcal{X}}\) 和 \(\mathcal{F}_{\mathcal{Y}}\) 看做相对于原始点云的 Embeding,其实也就是每个点相对位置的特征。

Attention

本文的注意力模型是学习一个函数:\(\phi:\mathbb{R}^{N\times P}\times\mathbb{R}^{N\times P}\to\mathbb{R}^{N\times P}\),其中 N 是代表点的数量,P 是代表 Embeding 也就是特征的维度。也就是说通过函数 对原始点云的特征进行加权转换后求得一个新的 Embeding,这个 Embding 是由 \(\mathcal{F}_{\mathcal{X}}\) 和 \(\mathcal{F}_{\mathcal{Y}}\) 获得,可以理解为是基于两个点云的匹配关系对 Embeding 进行加入 Attention 的重新生成新的 Embeding,其表达式如下:

\(\begin{aligned} \Phi_{\mathcal{X}} &=\mathcal{F}_{\mathcal{X}}+\phi\left(\mathcal{F}_{\mathcal{X}}, \mathcal{F}_{\mathcal{Y}}\right) \\ \Phi_{\mathcal{Y}} &=\mathcal{F}_{\mathcal{Y}}+\phi\left(\mathcal{F}_{\mathcal{Y}}, \mathcal{F}_{\mathcal{X}}\right) \end{aligned} \tag{1}\)

这里面,我们将函数 \(\phi\) 看做是一个残差,通过由 \(\mathcal{F}_{\mathcal{X}}\) 和 \(\mathcal{F}_{\mathcal{Y}}\) 学习到的关系,对于原有的 Embeding 进行的改变。这里的逻辑是,通过这样的映射 \(\mathcal{F}_{\mathcal{X}}\mapsto\Phi_{\mathcal{X}}\),\(\Phi_{\mathcal{X}}\) 不仅包含了原始点云 \(\mathcal{X}\) 的特征,同时也包含了目标点云 \(\mathcal{Y}\) 中结构的知识,因此对于我们所要进行的点云配准任务是更加有利的。反之同理。

这里面的函数 \(\phi\) 本文借鉴 NLP 中的 Transformer 网络来实现,因为点云匹配很像是一个 sep2seq 的问题,因而可以用这样的方式来处理。Transformer 的结构如下图所示:

2.3 匹配过程

相对于传统的寻找匹配点对的方式(这一过程不可导),本文采用的是 soft map 方式,可以理解 x 一个点不是和 Y 中某一个的点匹配,而是由 Y 生成一个 x 匹配的目标点的位置,与此位置实现关联。这一点可以近似理解为是差值,这一差值是由 Y 中所有点加权平均得到,我们所学习的 “匹配” 就是求这样的权重向量(soft map)。对于每一个点 \(\boldsymbol{x}_{i}\in\mathcal{X}\),其相对于点集 \(\mathcal{Y}\) 的 soft map 公式如下:

\(m\left(\boldsymbol{x}_{i}, \mathcal{Y}\right)=\operatorname{softmax}\left(\Phi_{\mathcal{Y}} \Phi_{\boldsymbol{x}_{i}}^{\top}\right) \tag{2}\)

其中 \(\Phi_{\mathcal{Y}}\in\mathbb{R}^{N\times P}\) 表示由 \(\mathcal{Y}\) 生成的 Embeding,\(\Phi_{\boldsymbol{x}_{i}}\) 代表 \(\Phi_{\boldsymbol{x}}\) 的第 i 行。

2.4 SVD模块

与《Least-Squares Rigid Motion Using SVD》文章相同,本文最终也是使用 SVD 求解 R、t,匹配的点是使用上述 soft map 重新生成的点云,生成函数如下:

\(\hat{\boldsymbol{y}}_{i}=Y^{\top} m\left(\boldsymbol{x}_{i}, \mathcal{Y}\right) \in \mathbb{R}^{3} \tag{3}\)

这里定义 \(Y\in\mathbb{R}^{N\times 3}\) 是 \(\mathcal{Y}\) 所有点坐标排列组成的矩阵。这样转换为新的匹配点对 \(\boldsymbol{x}_{i}\mapsto\hat{\boldsymbol{y}}_{i}\) 最终就可以用 SVD 方法求解 R 、t了。

2.5 Loss

误差函数设计比较简单,因为是回归 Pose,所以采用的是典型的 R 和 t 回归的方式:

\(\mathrm{Loss}=\left\|\boldsymbol{R}_{\mathcal{X} \mathcal{Y}}^{\top} \boldsymbol{R}_{\mathcal{X} \mathcal{Y}}^{g}-I\right\|^{2}+\left\|\boldsymbol{t}_{\mathcal{X} \mathcal{Y}}-\boldsymbol{t}_{\mathcal{X} \mathcal{Y}}^{g}\right\|^{2}+\lambda\|\theta\|^{2} \tag{4}\)

其中 g 代表真值(Ground Truth),\(\theta\) 是正则项。

PS:感觉误差函数没有精心设计,一个是 R 没有采用李代数,一个是 R 和 t 没有分别赋予权重,不知道是因为效果足够还是在这里并没有必要。

3 实验

实验部分作者做得还是比较详细的,比如比较 PointNet 和 DGCNN 效果,MLP 与 SVD 求解姿态效果,都能够证明本文方法的有效性。这里简单做一些摘要,详细的实验方法和数据建议阅读原文。

作者将是否使用 Attention 机制的 DCP 算法分别称作 DCPv1DCPv2,下面的实验数据中均有所体现。

3.1 迁移能力 & 鲁棒性

为了证明算法迁移能力作者在 ModelNet40 数据集上做了以下3组实验:

1)在一部分测试集中进行训练一部分未见过的数据集中进行测试(Table 1)
2)在一部分分类中进行训练在另一部分分类中进行测试(Table 2)
3)在测试数据中加入高斯噪声(Table 3)

实验效果均表明 DCP 可以取得很好的效果。

3.2 DCP 后使用 ICP 精调

如下图所示,在一些 Model 上虽然 DCP 可以给出非常接近于真值的 Pose,但是如果再加入 ICP 在此基础上精调的话取得的效果几乎等同于真值(我觉得这里最好是用基于优化的 ICP)。

3.3 PointNet 与 DGCNN 对比

在特征提取部分作者推荐他自己的 DGCNN,并与广为使用的 PointNet 做了对比试验,结果表明 DGCNN 具有明显的优势。

3.4 SVD 与 MLP 对比

很多基于 Deep Learning 方法的 Pose 回归直接使用 MLP,其实这样很不科学,但是似乎也很 work。作者也对 MLP 和 SVD 方法进行了对比,结果表明更加科学的 SVD 明显具有更好的精度(所以不要搞什么 MLP 回归 Pose 的东西了)。

3.5 Embedding 维度

对于 Embedding 维度作者也进行了对比,结果表明对于 DCPv1 维度从 512 提升到 1024 效果有更明显提升,但是 DCPv2 提升不明显。

个人小结

这篇文章从 SVD 求解 ICP 问题的细节原理出发,针对 ICP 问题定制了网络模型,分析非常细致而且也做了比较充足的实验,比如比较 PointNet 和 DGCNN 效果,MLP 与 SVD 求解姿态效果,都能够证明本文方法的有效性。

在我们的复现过程中,实际上重新生成的点云形状看起来和目标点云很不像,然而最终结果至少在 3D 上还是不错的,这一点 Transformer 起到的作用是不是和文章中所述还是比较迷。

相关材料

原文PDF:

官方源码:
https://github.com/WangYueFt/dcp

使用 SVD 方法求解 ICP 问题

本文是结合《Least-Squares Rigid Motion Using SVD》和《Least-Squares Fitting of Two 3-D Point Sets》两篇文章写的一个总结,里面有一些是自己的理解不一定正确。

1 问题定义

假设我们有两个点云集合 \(\mathcal{P}=\left\{\mathbf{p}_{1}, \mathbf{p}_{2}, \ldots, \mathbf{p}_{n}\right\}\) 和 \(\mathcal{Q}=\left\{\mathbf{q}_{1}, \mathbf{q}_{2}, \ldots, \mathbf{q}_{n}\right\}\),则我们定义的 ICP 问题是通过最小化点对之间距离获得相应的 Pose:

\((R, \mathbf{t})=\underset{R \in S O(d), \mathbf{t} \in \mathbb{R}^{d}}{\operatorname{argmin}} \sum_{i=1}^{n} w_{i}\left\|\left(R \mathbf{p}_{i}+\mathbf{t}\right)-\mathbf{q}_{i}\right\|^{2}\tag{1}\)

其中 \(w_i\) 代表每个点的权重。R 和 t 是我们所要求的旋转矩阵和平移向量。

2 计算平移

我们分两步求解旋转和平移,首先求解 t。固定 R,我们要优化的误差函数如下 \(F(\mathbf{t})=\sum_{i=1}^{n} w_{i}\left\|\left(R \mathbf{p}_{i}+\mathbf{t}\right)-\mathbf{q}_{i}\right\|^{2}\),为求最小值可以对 F 求导得到:

\(\begin{aligned} 0 &=\frac{\partial F}{\partial \mathbf{t}}=\sum_{i=1}^{n} 2 w_{i}\left(R \mathbf{p}_{i}+\mathbf{t}-\mathbf{q}_{i}\right)=\\ &=2 \mathrm{t}\left(\sum_{i=1}^{n} w_{i}\right)+2 R\left(\sum_{i=1}^{n} w_{i} \mathbf{p}_{i}\right)-2 \sum_{i=1}^{n} w_{i} \mathbf{q}_{i} \end{aligned}\tag{2}\)
记:

\(\overline{\mathbf{p}}=\frac{\sum_{i=1}^{n} w_{i} \mathbf{p}_{i}}{\sum_{i=1}^{n} w_{i}}, \quad \overline{\mathbf{q}}=\frac{\sum_{i=1}^{n} w_{i} \mathbf{q}_{i}}{\sum_{i=1}^{n} w_{i}}\tag{3}\)
也就是加权平均的质心,并再次带回 (2) 式可以得到:

\(\mathbf{t}=\overline{\mathbf{q}}-R \overline{\mathbf{p}}\tag{4}\)
那么可以知道,t 可以由加权的质心得到。再将 t 带回目标函数 F 可得:

\(\begin{aligned} \sum_{i=1}^{n} w_{i}\left\|\left(R \mathbf{p}_{i}+\mathbf{t}\right)-\mathbf{q}_{i}\right\|^{2} &=\sum_{i=1}^{n} w_{i}\left\|R \mathbf{p}_{i}+\overline{\mathbf{q}}-R \overline{\mathbf{p}}-\mathbf{q}_{i}\right\|^{2}=\\ &=\sum_{i=1}^{n} w_{i}\left\|R\left(\mathbf{p}_{i}-\overline{\mathbf{p}}\right)-\left(\mathbf{q}_{i}-\overline{\mathbf{q}}\right)\right\|^{2} \end{aligned}\tag{5}\)
再记:

\(\mathbf{x}_{i} :=\mathbf{p}_{i}-\overline{\mathbf{p}}, \quad \mathbf{y}_{i} :=\mathbf{q}_{i}-\overline{\mathbf{q}}\tag{6}\)
我们可以直接求质心后对每个点进行归一化再求解。有了 t 之后,我们要求解的 R 如下:

\(R=\underset{R \in S O(d)}{\operatorname{argmin}} \sum_{i=1}^{n} w_{i}\left\|R \mathbf{x}_{i}-\mathbf{y}_{i}\right\|^{2}\tag{7}\)
其中,d 表示 x 和 y 的维度,通常情况下对于三维点云就是 d = 3,对于二维点云就是 d = 2。

3 计算旋转

首先简化下 (7) 式,先看权重后面的部分,则公式推导如下:

\(\begin{aligned}\left\|R \mathbf{x}_{i}-\mathbf{y}_{i}\right\|^{2} &=\left(R \mathbf{x}_{i}-\mathbf{y}_{i}\right)^{\top}\left(R \mathbf{x}_{i}-\mathbf{y}_{i}\right)=\left(\mathbf{x}_{i}^{\top} R^{\top}-\mathbf{y}_{i}^{\top}\right)\left(R \mathbf{x}_{i}-\mathbf{y}_{i}\right)=\\ &=\mathbf{x}_{i}^{\top} R^{\top} R \mathbf{x}_{i}-\mathbf{y}_{i}^{\top} R \mathbf{x}_{i}-\mathbf{x}_{i}^{\top} R^{\top} \mathbf{y}_{i}+\mathbf{y}_{i}^{\top} \mathbf{y}_{i}=\\ &=\mathbf{x}_{i}^{\top} \mathbf{x}_{i}-\mathbf{y}_{i}^{\top} R \mathbf{x}_{i}-\mathbf{x}_{i}^{\top} R^{\top} \mathbf{y}_{i}+\mathbf{y}_{i}^{\top} \mathbf{y}_{i} \end{aligned}\tag{8}\)
上述推导中最后一步解释下:

首先,\(R^{\top} R=I\),因此有 \(\mathbf{x}_{i}^{\top} R^{\top} R \mathbf{x}_{i} = \mathbf{x}_{i}^{\top} \mathbf{x}_{i}\);

其次,由于 \(\mathbf{x}_{i}^{\top} R^{\top} \mathbf{y}_{i}\) 是一个标量(这是很显然的,\(\mathbf{x}_{i}^{\top}\) 是 \(1 \times d\) 的向量,\(R^{\top}\) 是 \(d \times d\) 矩阵,\(\mathbf{y}_{i}\) 是 \(d \times 1\) 向量)。对于标量,我们有性质 \(a=a^{\top}\),因此有:

\(\mathbf{x}_{i}^{\top} R^{\top} \mathbf{y}_{i}=\left(\mathbf{x}_{i}^{\top} R^{\top} \mathbf{y}_{i}\right)^{\top}=\mathbf{y}_{i}^{\top} R \mathbf{x}_{i}\tag{9}\)
带回公式 (8) 可以得到:

\(\left\|R \mathbf{x}_{i}-\mathbf{y}_{i}\right\|^{2}=\mathbf{x}_{i}^{\top} \mathbf{x}_{i}-2 \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}+\mathbf{y}_{i}^{\top} \mathbf{y}_{i}\tag{10}\)
整理公式 (7) 可以得到:

\(\begin{aligned} & \operatorname{argmin}_{R \in S O(d)} \sum_{i=1}^{n} w_{i}\left\|R \mathbf{x}_{i}-\mathbf{y}_{i}\right\|^{2}=\underset{R \in S O(d)}{\operatorname{argmin}} \sum_{i=1}^{n} w_{i}\left(\mathbf{x}_{i}^{\top} \mathbf{x}_{i}-2 \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}+\mathbf{y}_{i}^{\top} \mathbf{y}_{i}\right)=\\=& \operatorname{argmin}_{R \in S O(d)}\left(\sum_{i=1}^{n} w_{i} \mathbf{x}_{i}^{\top} \mathbf{x}_{i}-2 \sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}+\sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} \mathbf{y}_{i}\right)=\\=& \underset{R \in S O(d)}{\operatorname{argmin}}\left(-2 \sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}\right) \end{aligned}\tag{11}\)

最后一个的等号显而易见,因为 \(\sum_{i=1}^{n} w_{i} \mathbf{x}_{i}^{\top} \mathbf{x}_{i}\) 和 \(\sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} \mathbf{y}_{i}\) 两项都与 R 无关。再进一步去掉负号改成求最大值则有:

\(\underset{R \in S O(d)}{\operatorname{argmin}}\left(-2 \sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}\right)=\underset{R \in S O(d)}{\operatorname{argmax}} \sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}\tag{12}\)
根据矩阵迹的性质,我们有以下结论:

\(\sum_{i=1}^{n} w_{i} \mathbf{y}_{i}^{\top} R \mathbf{x}_{i}=\operatorname{tr}\left(W Y^{\top} R X\right)\tag{13}\)
其中:\(W=\left[\begin{array}{ccccc}{w_{1}} & {} & {} & {} & {} \\ {} & {w_{2}} & {} & {} & {} \\ {} & {} & {} & {\ddots} & {} \\ {} & {} & {} & {} & {w_{n}}\end{array}\right]\),\(Y^{\top} = \left[\begin{array}{c}{-\mathbf{y}_{1}^{\top}-} \\ {-\mathbf{y}_{2}^{\top}-} \\ {\vdots} \\ {-\mathbf{y}_{n}^{\top}-}\end{array}\right]\),\(X=\begin{bmatrix}| & | & & |\\x_1 & x_2 & \cdots & x_n\\ | & | & & |\end{bmatrix}\)。

上述公式推导如下:

\(\begin{align*}W Y^{\top} R X &= \left[\begin{array}{ccccc}{w_{1}} & {} & {} & {} & {} \\ {} & {w_{2}} & {} & {} & {} \\ {} & {} & {} & {\ddots} & {} \\ {} & {} & {} & {} & {w_{n}}\end{array}\right]
\left[\begin{array}{c}{-\mathbf{y}_{1}^{\top}-} \\ {-\mathbf{y}_{2}^{\top}-} \\ {\vdots} \\ {-\mathbf{y}_{n}^{\top}-}\end{array}\right]
\left[\begin{array}{ccccc}{} & {} & {} \\ {} & {R} & {} \\ {} & {} & {} \end{array}\right]
\begin{bmatrix}| & | & & |\\x_1 & x_2 & \cdots & x_n\\ | & | & & |\end{bmatrix}\\
&= \left[\begin{array}{c}{-w_{1}\mathbf{y}_{1}^{\top}-} \\ {-w_{2}\mathbf{y}_{2}^{\top}-} \\ {\vdots} \\ {-w_{n}\mathbf{y}_{n}^{\top}-}\end{array}\right]
\begin{bmatrix}| & | & & |\\Rx_1 & Rx_2 & \cdots & Rx_n\\ | & | & & |\end{bmatrix}\\
&=
\left[\begin{array}{cccc}{w_{1} \mathbf{y}_{1}^{\top} R \mathbf{x}_{1}} & {} & {} & {*} \\ {} & {w_{2} \mathbf{y}_{2}^{\top} R \mathbf{x}_{2}} & {} \\ {} & {} & {\ddots} & {} \\ {*} & {} & {} & {w_{n} \mathbf{y}_{n}^{\top} R \mathbf{x}_{n}}\end{array}\right] \end{align*}\tag{14}\)

因此很明显可以得出公式 (13) 中的结论。

之所以用这样表示是因为用矩阵的迹有很多性质可以利用,根据迹的性质我们有对称性:

\(\operatorname{tr}(A B)=\operatorname{tr}(B A)\tag{15}\)
因此我们有:

\(\operatorname{tr}\left(W Y^{\top} R X\right)=\operatorname{tr}\left(\left(W Y^{\top}\right)(R X)\right)=\operatorname{tr}\left(R X W Y^{\top}\right)\tag{16}\)
定义“协方差”矩阵 \(S=X W Y^{\top}\),对 S 进行 SVD 分解:

\(S=U \Sigma V^{\top}\tag{17}\)
将上述公式带入 (16) 同时再根据性质 (15),则有:

\(\operatorname{tr}\left(R X W Y^{\top}\right)=\operatorname{tr}(R S)=\operatorname{tr}\left(R U \Sigma V^{\top}\right)=\operatorname{tr}\left(\Sigma V^{\top} R U\right)\tag{18}\)
因为其中 V、R 和 U 都是正交矩阵。因此 \(M=V^{\top} R U\) 也是正交矩阵。根据正交矩阵的定义其行、列向量均为正交的单位向量,因此对于矩阵 M 中每个列向量 \(\mathbf{m_j}\) 均有 \(\mathbf{m}_{j}^{\top} \mathbf{m}_{j}=1\)。可以推出 M 中每一个元素的绝对值均有 \(\left | \mathbf{m_{ij}}\right |\leqslant 1\)。推导如下:

\(1=\mathbf{m}_{j}^{\top} \mathbf{m}_{j}=\sum_{i=1}^{d} m_{i j}^{2} \Rightarrow m_{i j}^{2} \leq 1 \Rightarrow\left|m_{i j}\right| \leq 1\tag{19}\)
根据 SVD 分解的性质 \(\Sigma =\left(\begin{array}{cccc}{\sigma_{1}} & {} & {} & {} \\ {} & {\sigma_{2}} & {} & {} \\ {} & {} & {\ddots} & {} \\ {} & {} & {} & {\sigma_{d}}\end{array}\right)\),其中对角线的元素 \(\sigma_{1}, \sigma_{2}, \ldots, \sigma_{d} \geq 0\),因此有:

\(\operatorname{tr}(\Sigma M)=\left(\begin{array}{cccc}{\sigma_{1}} & {} & {} & {} \\ {} & {\sigma_{2}} & {} & {} \\ {} & {} & {\ddots} & {} \\ {} & {} & {} & {\sigma_{d}}\end{array}\right)\left(\begin{array}{c}{m_{11} m_{12} \ldots m_{1 d}} \\ {m_{21} m_{22} \ldots m_{2 d}} \\ {\vdots} \\ {m_{d 1} m_{d 2} \ldots m_{d d}}\end{array}\right)=\sum_{i=1}^{d} \sigma_{i} m_{i i} \leq \sum_{i=1}^{d} \sigma_{i}\tag{20}\)
上式等号成立也就是取得最大值的唯一条件是每一个 \(\mathbf{m}_{i i}=1\),由于 M 是正交矩阵,其行列向量均为正交的单位向量,因此只能是 M 为单位矩阵:\(M = I\),得到:

\(I=M=V^{\top} R U \Rightarrow V=R U \Rightarrow R=V U^{\top}\tag{21}\)
这就是 ICP 中求解时计算 \(R=V U^{\top}\) 的推导过程。

4 反射修正

4.1 反射变换

首先理解下反射变换(也叫镜面反射),根据维基百科定义,在数学上反射是把一个物体变换成它的镜像的映射。要反射一个平面图形,需要“镜子”是一条直线(反射轴),对于三维空间中的反射就要使用平面作为镜子。

反射变换同样是一个正交矩阵,显而易见它满足如下性质:

  • 镜面反射是正交变换。
  • 镜面反射的逆变换为镜面反射。
  • 任意一个正交变换都可以表示成若干个镜面反射的乘积。

除此之外,旋转变换和反射变换具有如下特性(对于2D和3D同样成立):

  • 旋转矩阵和反射矩阵都是正交矩阵
  • 旋转矩阵的行列式值为 +1,反射矩阵的行列值为 -1
  • 旋转矩阵 \(R(\theta)\) 的逆矩阵为 \(R(-\theta)\),反射矩阵的逆矩阵为其本身。或者也可以记为,旋转矩阵:\(R^{\top}R=I\),反射矩阵:\(R’R’=I\)。
  • 旋转矩阵和反射矩阵可以相互转换:
    \(\begin{aligned}\operatorname{Ref}(\theta) \operatorname{Ref}(\phi) &=\operatorname{Rot}(2(\theta-\phi)) \\ \operatorname{Rot}(\theta) \operatorname{Rot}(\phi) &=\operatorname{Rot}(\theta+\phi) \\ \operatorname{Rot}(\theta) \operatorname{Ref}(\phi) &=\operatorname{Ref}(\phi+\theta / 2) \\ \operatorname{Ref}(\phi) \operatorname{Rot}(\theta) &=\operatorname{Ref}(\phi-\theta / 2) \end{aligned}\tag{22}\)

最简单的反射变换是沿某个轴/面的镜像,例如相对于 \(Z=0\) 平面的镜像变换为:\(\begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\end{pmatrix}\)。则根据上面的最后一条性质,我们一定可以把任意一个镜面变换拆成一个关于 \(Z=0\) 的镜面变换与一个旋转变换。也就是变成:

\(R’=\begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\end{pmatrix}R\tag{23}\)

4.2 反射修正

通过之前的推导,我们用 SVD 求解的 R 一定是一个正交矩阵,但是并不是所有的正交矩阵都是旋转矩阵,也可能是一个反射矩阵(或者说包含了反射变换的矩阵)。因此我们还需要对所求得的 R 进行行列式判断,判断方法也很简单:

  • 如果 \(\operatorname{det}\left(V U^{\top}\right)=-1\),则所求的 R 包含了镜像;
  • 如果 \(\operatorname{det}\left(V U^{\top}\right)=1\),则所求的 R 是我们所求的旋转矩阵。

这一部分在《Least-Squares Fitting of Two 3-D Point Sets》里面讲得更详细一点。假设 \(\mathcal{P}=\left\{\mathbf{p}_{1}, \mathbf{p}_{2}, \ldots, \mathbf{p}_{n}\right\}\) 和 \(\mathcal{Q}=\left\{\mathbf{q}_{1}, \mathbf{q}_{2}, \ldots, \mathbf{q}_{n}\right\}\) 是可以完美重合的两团点云,两个点云有如下三种情况需要讨论:

1)点云 \(\mathcal{Q}\) 不共面

在这种情况下,显然有唯一的旋转矩阵 R 满足公式 (1)。可以直观想象不共面的两个点云(当然也需要点云本身不满足镜像性质)是不可能通过镜像变换映射成自身的。

2)点云 \(\mathcal{Q}\) 共面不共线

在这种情况下,满足公式 (1) 的解不唯一。一定同时存在旋转矩阵 R 和反射矩阵 R’ 都可以满足。说明如下:

对于 3×3 协方差矩阵 \(S=X W Y^{\top}\) 由于点 X 和 Y 是共面的,也就意味着存在平面法向 \(n = {n_x, n_y, n_z}\) 满足 \(nX = d\),那么 X 和 Y 均只有两个自由变量,S 的秩为 2。那么也就是:

\(S=U \Sigma V^{\top}=U \begin{pmatrix}\sigma_1 & 0 & 0\\ 0 & \sigma_2 & 0\\ 0 & 0 & \sigma_2\end{pmatrix} V^{\top}\tag{24}\)
由于 SVD 分解特征值是从大到小排序,则一定有:

\(S=\sigma_{1} u_{1} v_{1}^{\top}+\sigma_{2} u_{2} v_{2}^{\top}+\sigma_{3} u_{3} v_{3}^{\top}=\sigma_{1} u_{1} v_{1}^{\top}+\sigma_{2} u_{2} v_{2}^{\top}+0 \cdot u_{3} v_{3}^{\top},\\ \sigma_{1}>\sigma_{2}>\sigma_{3}=0\tag{25}\)
如果存在一个解 \(R=V U^{\top}=\left [ v_1, v_2, v_3 \right ] U^{\top}\) 满足以上 S 取得极大值,则一定有镜像变换:

\(R’=V’ U^{\top}=\left [ v_1, v_2, -v_3 \right ] U^{\top}=\left [ v_1, v_2, v_3 \right ] \begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\end{pmatrix}U^{\top}\tag{26}\)
满足上式 S 取得极大值。

所以当 \(\operatorname{det}\left(V U^{\top}\right)=-1\) 时,我们想要求得的 R 只需要去除中间乘的反射变换矩阵 \(\begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\end{pmatrix}\) 即可。也就是 \(R=V U^{\top}=V’ \begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1\end{pmatrix} T^{\top}\)。

如果是 \(\operatorname{det}\left(V U^{\top}\right)=1\) 时,那么本身所求的就是旋转矩阵 R 我们也可以写成 \(R=V U^{\top}=V \begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1\end{pmatrix} T^{\top}\)。

整理上述两种情况就可以统一成以下表达式:

\(R=V U^{\top}=V \begin{pmatrix}1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & \operatorname{det}\left(V U^{\top}\right)\end{pmatrix} T^{\top}\tag{27}\)

3)点云 \(\mathcal{Q}\) 共线

这种情况下存在不止一种解,无法用 SVD 求得旋转矩阵。

5 总结

经过上面的推导和镜像修正,我们可以总结出一套完整的使用 SVD 求解 ICP 问题的流程:

我们的问题是求解 R, t 使得下面的误差函数最小:

\(\sum_{i=1}^{n} w_{i}\left\|\left(R \mathbf{p}_{i}+\mathbf{t}\right)-\mathbf{q}_{i}\right\|^{2}\tag{27}\)
则进行如下步骤:

1)计算点云中心:

\(\overline{\mathbf{p}}=\frac{\sum_{i=1}^{n} w_{i} \mathbf{p}_{i}}{\sum_{i=1}^{n} w_{i}}, \quad \overline{\mathbf{q}}=\frac{\sum_{i=1}^{n} w_{i} \mathbf{q}_{i}}{\sum_{i=1}^{n} w_{i}}\tag{29}\)
2)对所有点做归一化:

\(\mathbf{x}_{i} :=\mathbf{p}_{i}-\overline{\mathbf{p}}, \quad \mathbf{y}_{i} :=\mathbf{q}_{i}-\overline{\mathbf{q}}, \quad i=1,2, \ldots, n\tag{30}\)
3)计算 \(d \times d\) 协方差矩阵(d代表数据维度):

\(S=X W Y^{\top}\tag{31}\)

其中:\(X=\begin{bmatrix}| & | & & |\\x_1 & x_2 & \cdots & x_n\\ | & | & & |\end{bmatrix}\),\(Y=\begin{bmatrix}| & | & & |\\y_1 & y_2 & \cdots & y_n\\ | & | & & |\end{bmatrix}\) 均是由 d 维列向量组成的 \(d \times n\) 维的矩阵;\(W=\left[\begin{array}{ccccc}{w_{1}} & {} & {} & {} & {} \\ {} & {w_{2}} & {} & {} & {} \\ {} & {} & {} & {\ddots} & {} \\ {} & {} & {} & {} & {w_{n}}\end{array}\right]\) 是权重矩阵。

4)计算 S 的 SVD 分解 \(S=U \Sigma V^{\top}\),则得到想要求的旋转矩阵 R 如下:

\(R=V\left(\begin{array}{ccccc}{1} \\ {} & {1} \\ {} & {} & {\ddots} & {} \\ {} & {} & {} & {\operatorname{det}\left(V U^{T}\right)}\end{array}\right) U^{\top}\tag{32}\)
5)计算平移向量 t:

\(\mathbf{t}=\overline{\mathbf{q}}-R \overline{\mathbf{p}}\tag{31}\)

参考文献

  1. https://blog.csdn.net/icameling/article/details/84103891
  2. https://blog.csdn.net/dfdfdsfdfdfdf/article/details/53213240
  3. https://blog.csdn.net/kfqcome/article/details/9356819https://blog.csdn.net/kfqcome/article/details/9356819
  4. https://math.stackexchange.com/questions/68119/why-does-ata-i-det-a-1-mean-a-is-a-rotation-matrix
  5. https://wangliangster.github.io/#/math/algra/rotatemat
  6. https://zhuanlan.zhihu.com/p/51362089

相关材料

原文 PDF:


Ubuntu 下监控并自动重启网卡

很多时候网站服务器挂掉也可能是因为网卡挂掉了,如果你网站不能访问时 SSH 也无效了一般都是这个问题。这时可以通过一个定时脚本监控网络并进行自动重启。相关文章:Linux 下如何监控并自动重启 apache

1 创建脚本 auto_restart_network.sh

内容如下:


上传此文件到服务器。

注释:
1)该脚本会 ping www.baidu.com 网站,当然你也可以改成其他的网站,比如你自己的页面。
2)根据 ping 的结果,如果不成功的话,则使用命令重启网络连接服务。

2 创建定时任务

创建定时任务有两种方法一种是直接编辑 /etc/crontab,运行:


在后面添加:


其中 */30 表示每 30 分钟执行一次脚本。
关于 crontab 时间说明如下:


另一种是使用命令:


在后面添加:


我的理解,这两种方式前者添加的是一个整个系统所有用户的定时任务,所以在时间后面需要指定用户 root。后者是当前用户的定任务。采用哪一种都可以。

重启定时任务:


如果顺利的话这一脚本就会每隔 30 分钟执行一次,并且可以在 /alidata/tools/network_restart_logs.txt 看到脚本执行的日志以及 /alidata/tools/logs.txt 看到每次执行的输出。因为频率低以及只输出不多的几个字符,这两个 log 文件通常不会增加很多,不过也应注意下及时清除。

附录:网络服务操作命令

Ubuntu 下网络服务操作命令如下:
开启网络服务


关闭网络服务


重新启动网络服务

参考材料

https://www.cyberciti.biz/faq/linux-restart-network-interface/

论文笔记:Dynamic Graph CNN for Learning on Point Clouds

DGCNN 是对 PointNet 的改进,PointNet 网络每个点单独提取特征缺乏局部关联。DGCNN 提出了 EdgeConv 就是对它的改进。

1 网络结构

DGCNN 网络结构如下图所示,可以看出其整体架构和 PointNet 是基本一致的,主要区别就是将其中的 MLP 替换成了 EdgeConv。

2 EdgeConv

2.1 EdgeConv 结构

上图是 EdgeConv 的示意图。假设一个F维点云 \(\mathbf{X}=\left\{\mathbf{x}_{1}, \dots, \mathbf{x}_{n}\right\} \subseteq \mathbb{R}^{F}\) 其中 F 表示每个点的维度,最简单的可能是 x, y, z 三维,另外还可能引入每个点颜色、法线等信息。给定一个有向图 \(\mathcal{G}=(\mathcal{V}, \mathcal{E})\) 用来表示点云的局部结构,其中顶点为 \(\mathcal{V}=\{1, \ldots, n\}\),边为 \(\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}\)。在一个简单例子中,我们用 \(\mathcal{G}\) 表示 k-nearest neighbor (k-NN) 图。针对上图有如下定义:

点和边缘:定义 \(x_i\) 周围最近 k 个点 \(x_{j_{i 1}}, \ldots, x_{j_{i k}}\) ,其有向边为 \(\left(i, j_{i 1}\right), \ldots,\left(i, j_{i k}\right)\)
边缘特征函数:定义边缘特征为 \(e_{i j}=h_{\Theta}\left(x_{i}, x_{j}\right)\),其中 \(h_{\Theta} : \mathbb{R}^{F} \times \mathbb{R}^{F} \rightarrow \mathbb{R}^{F^{\prime}}\) 是一些使用一些可学习的参数 \(\Theta\) 构成的非线性函数。
特征聚合函数:对于边缘特征进行聚合的函数定义为□,其定义是:

\(\begin{align*}\mathbf{x}_{i}^{\prime}=\square_{j :(i, j) \in \mathcal{E}} h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)\end{align*}\tag{1}\)

上述KNN图中的K值是一个超参,分类网络中K=20,而在分割网络中K=30。

关于公式中的 \(h\) 和 □ 有四种可能的选择:

第一种:认为 \(x_i\) 的特征是周围所有点的加权求和,这一点类似于图像的卷积操作,其中每个卷积核为 \(\theta_{m}\),\(\Theta=\left(\theta_{1}, \ldots, \theta_{M}\right)\) 表示所有卷积核的集合。其公式如下:

\(\begin{align*}x_{i m}^{\prime}=\sum_{j :(i, j) \in \mathcal{E}} \boldsymbol{\theta}_{m} \cdot \mathbf{x}_{j}\end{align*}\tag{2}\)
第二种:描述了一种仅考虑全局特征的形式,即只有关注点 \(x_i\)。其公式如下:

\(\begin{align*}h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)=h_{\Theta}\left(\mathbf{x}_{i}\right)\end{align*}\tag{3}\)
第三种:只关注周围特征形式。其公式如下:

\(\begin{align*}h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)=h_{\Theta}\left(\mathbf{x}_{j}\right)\end{align*}\tag{4}\)
其中:

\(\begin{align*}x_{i m}^{\prime}=\sum_{j \in \mathcal{V}}\left(h_{\boldsymbol{\theta}\left(\mathbf{x}_{j}\right)}\right) g\left(u\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)\right)\end{align*}\tag{5}\)
其中 \(g\) 表示高斯核,\(u\) 表示两个点的欧氏距离。当然在输入层是坐标之间的欧氏距离,后面的层级中,这个函数表示的是两个点特征向量的欧氏距离。

第四种:只关注局部特征形式,输入仅为点和周围点的差。其公式如下:

\(\begin{align*}h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)=h_{\Theta}\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right)\end{align*}\tag{6}\)
第五种:同时关注全局特征和局部特征,这也是本文中的主要形式:

\(\begin{align*}h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)=\overline{h}_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}-\mathbf{x}_{i}\right)\end{align*}\tag{7}\)
在本文中我们使用:

\(\begin{align*}e_{i j m}^{\prime}=\operatorname{ReLU}\left(\boldsymbol{\theta}_{m} \cdot\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right)+\boldsymbol{\phi}_{m} \cdot \mathbf{x}_{i}\right)\end{align*}\tag{8}\)
该项可以用 shared MLP 表示,最后使用 max 进行特征聚合:

\(\begin{align*}x_{i m}^{\prime}=\max _{j :(i, j) \in \mathcal{E}} e_{i j m}^{\prime}\end{align*}\tag{9}\)
全部超参数集合为: \(\Theta=\left(\theta_{1}, \ldots, \theta_{M}, \phi_{1}, \ldots, \phi_{M}\right)\)

4.2 EdgeConv 实现

输入为NxF的张量,其中包括原数据Xi和与其附近的K个点。首先原数据的每个点和附件的K个点,通过MLP生成K个NxM特征,即NxKxM,最后通过一个pooling操作输出NxM特征。

看它的代码更清晰:

其中的:

就是上面所说的 EdgeConv 部分。

3 动态图更新

实验表明,每次更新后重新计算每层点和周围的最近邻会取得更好的性能,这一动态图称为 Dynamic Graph CNN (DGCNN)。它的思想如下:

每一层都会得到一个不同的动态图,其中第 l 层动态图表示为:\(\mathcal{G}^{(l)}=\left(\mathcal{V}^{(l)}, \mathcal{E}^{(l)}\right)\)
每一层边缘 \(\left(i, j_{i 1}\right), \ldots,\left(i, j_{i k_{l}}\right)\) 表示 \(\mathbf{x}_{i}^{(l)}\) 周围 \(k_l\) 个最近邻点表示为 \(\mathbf{x}_{j_{i 1}}^{(l)}, \ldots, x_{j_{i k_{l}}}^{(l)}\)。每次更新后重新计算每层点和周围的最近邻重新得到边缘特征。
更新公式如下:

\(\begin{align*}x^{l+1}=\square_{j :(i, j) \in \epsilon^{t}} h_{\theta}^{l}\left(x_{i}^{l}, x_{j}^{l}\right)\end{align*}\tag{10}\)
相比一个普通的 CNN 操作,我们的动态图同时学习如何构建图结构本身,而不是仅仅把网络变成一个固定的参数

4 特性

排列不变性

对于点云来说,点的排列顺序不应该影响最终的结果,考虑到我们的每一层特征聚合公式:

\(\begin{align*}\mathbf{x}_{i}^{\prime}=\max _{j :(i, j) \in \mathcal{E}} h_{\Theta}\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)\end{align*}\tag{11}\)
在我们的网络中,由于每一层 global max pooling 的存在,使得排列不变性可以获得(这一点类似于 PointNet)。

平移不变性(部分)

假设我们将之前的公式中每个点云位置 x 加上平移 T ,那么会得到如下结果:

\(\begin{align*} e_{i j m}^{\prime} &=\theta_{m} \cdot\left(\mathbf{x}_{j}+T-\left(\mathbf{x}_{i}+T\right)\right)+\boldsymbol{\phi}_{m} \cdot\left(\mathbf{x}_{i}+T\right) \\ &=\boldsymbol{\theta}_{m} \cdot\left(\mathbf{x}_{j}-\mathbf{x}_{i}\right)+\boldsymbol{\phi}_{m} \cdot\left(\mathbf{x}_{i}+T\right) \end{align*}\tag{12}\)
假如我们只考虑局部特征 \(\mathbf{x}_{j}-\mathbf{x}_{i}\),也就是令 \(\phi_{m}=0\),则上市满足平移不变性。

5 与其他方法对比

在作者的分析中,将几种常见的点云网络进行的对比(分为两类:PointNet 系列:PointNet、PointNet++;Graph CNN 系列:MoNet、PCNN),并且说明其实他们可以看作是 EdgeConv 模型的特例(只是边缘函数 、聚合函数 □ 定义的不同),作者给出了这样一个表格:

6 实验结果

这部分建议看原文,作者做了非常丰富和详尽的实验。

个人小结

这篇文章的从建模到设计到数学分析非常到位,阅读起来很有价值,非常像传统方法中的一篇 Generalized-ICP ,将其他的方法建模到自己的框架里面,这样从理论上解释了为什么自己的方法会比其他方法更好。我觉得这是这篇文章最大的亮点。

参考文献

  1. Dynamic Graph CNN for Learning on Point Clouds, Wang, Yue and Sun, Yongbin and Liu, Ziwei and Sarma, Sanjay E. and Bronstein, Michael M. and Solomon, Justin M., 2019
  2. https://blog.csdn.net/legend_hua/article/details/79175315
  3. https://blog.csdn.net/hongbin_xu/article/details/85258278

相关材料

原文 PDF

代码
https://github.com/WangYueFt/dgcnn

论文笔记:GCNv2: Efficient Correspondence Prediction for Real-Time SLAM

GCNv2 是一个专门针对几何匹配的描述子网络,是对 GCN 的改进版主要工作如下:1)与常见深度学习特征匹配的性能并且显著减少了前向运算的时间;2)加入了二值化层,生成二值特征。

1 GCNv2 网络结构

GCNv2 网络结构如下图所示:

这一结构不必赘述,其实跟 SuperPoint 很像,不过 Descriptors 部分是把 Keypoints 直接拿来取了相应的部分。PixelShuffle (有人也叫作 Sub-pixel Convolution 或者亚像素卷积)在文章(《Real-Time Single Image and Video Super-Resolution Using an Efficient Sub-Pixel Convolutional Neural Network》)中讲过,它的大概思想如下图所示:

它的效果是将一个 \(H \times W\) 的低分辨率输入图像(Low Resolution),通过Sub-pixel操作将其变为 \(rH\times rW\) 的高分辨率图像(High Resolution)。

但是其实现过程不是直接通过插值等方式产生这个高分辨率图像,而是通过卷积先得到  个通道的特征图(特征图大小和输入低分辨率图像一致),然后通过周期筛选(periodic shuffing)的方法得到这个高分辨率的图像,其中 \(r\) 为上采样因子(upscaling factor),也就是图像的扩大倍率。

图中双层的表示连续两个卷积。具体可以看他的实现代码其实非常明晰:


为了在移动端上更快,作者还提出了 GCNv2-tiny,进一步压缩了参数(主要是把 conv2 和之后的卷积通道数减半),达到了 TX2 上实时(40Hz)的效果。GCNv2-tiny 的代码如下:

2 特征提取

2.1 定义

我们令 \(f\) 表示 convD 中提取的压缩尺度的特征图,\(o\) 表示 PixelShuffle 提取的原始大小的关键点热点图。当我们给出一个图像坐标系的坐标 \(x = (u, v)\),显然我们能够得到其在 \(f\left( \cdot \right)\) 和 \(\\o\left( \cdot \right)\) 的对应位置坐标,需要注意的是,由于 \(f\) 是缩略尺寸特征,其描述是通过双线性插值得到的。上角标表示输入的是 current 还是 target 帧(例如:\(f^{cur}\) 或 \(f^{tar}\))。另外,\(L_{feat}\) 表示 descriptor 的损失函数,\(L_{det}\) 表示 detector 的损失函数。

 2.2 二值特征

为了获取描述子的二值特征,我们定义了一个 Binary activation layer,它的定义如下(其实只是一个符号函数):

Forward: \(\boldsymbol{b}(\boldsymbol{x})=\operatorname{sign}(\boldsymbol{f}(\boldsymbol{x}))=\left\{\begin{array}{ll}{+1} & {\boldsymbol{f}(\boldsymbol{x}) \geq 0} \\ {-1} & {\text { otherwise }}\end{array}\right.\)

Backward: \( \frac{\partial b}{\partial f}=\mathbf{1}_{|f| \leq 1}\)

作者发现加入这样一个 Binary activation layer 比直接让网络输出 -1~+1 之间的数值效果更好。

这一二值网络来自于论文《Binarized neural networks: Training neural networks with weights and activations constrained to+ 1 or-1》,对其思想简述如下:

2.2.1 预测过程

前向运算

假设 \(W^{k}\) 是第 k 层的权重矩阵,\(x^{k}\) 是第 k 层的输入,则 Binary Layer 的计算过程如下:

\(x^{k}=\sigma\left(W^{k} \cdot x^{k-1}\right)\)

其中 \(\sigma( \cdot )\) 是二值化的激活函数其定义如下:

\(\sigma(x)=\operatorname{sign}(x)=\left\{\begin{array}{ll}{+1} & {\text { if } x \geq 0} \\ {-1} & {\text { otherwise }}\end{array}\right.\)

2.2.2 训练过程

由于网络参数已经二值化,直接正向或者负向调整会导致模型震荡 ,因此不能直接进行更新,为了解决这个问题,权值和梯度在训练过程中保持全精度。

前馈运算(神经网络后一层向前一层的传播)

作者给出了一个例子,以BN层的二值化为例,先将实数型权值参数二值化得到二值型权值参数,即 \(W_{k}^{b}=\operatorname{sign}\left(W_{k}\right)\)。然后利用二值化后的参数计算得到实数型的中间向量,该向量再通过 \(BN(⋅)\) 即 Batch Normalization 操作,得到实数型的隐藏层激活向量。如果不是输出层的话,就将该向量二值化。

\(x^{k}=\sigma\left(B N\left(W_{b}^{k} \cdot x^{k-1}\right)=\operatorname{sign}\left(B N\left(W_{b}^{k} \cdot x^{k-1}\right)\right)\right.\)

反向传播

由于 \(σ(⋅)\) 基本上是一个梯度处处为零的函数,其不适合用于反向传播,因此 Binarized neural networks 提出采用借鉴 Yoshua Benjio 在 2013年 的时候通过 “straight-through estimator” 随机离散化神经元来进行梯度计算和反向传播的工作。也就是将 \(σ(x)=sign(x)\) 近似为下面的函数:

\(\operatorname{Htanh}(x)=\operatorname{Clip}(x,-1,1)=\max (-1, \min (1, x))\)

本质上就是在 -1 与 +1 之间加上一个线性变换过程,导数为1,这样就变成了可导。

如图所示:

那么解决了 \(Binarize()\) 函数的可导问题之后,后面的就是关于 BN 层二值化的前向和反向传播过程:

2.3 描述子学习

误差函数

二值特征学习部分采用 triplet loss ,其定义如下:

\(L_{f e a t}=\sum_{i} \max \left(0, d\left(\boldsymbol{x}_{i}^{c u r}, \boldsymbol{x}_{i,+}^{\operatorname{tar}}\right)-d\left(\boldsymbol{x}_{i}^{c u r}, \boldsymbol{x}_{i,-}^{\operatorname{tar}}\right)+m\right)\)

其中 \(x^{cur}\)表示 anchor, 表示 \(x^{tar}_+\) positive,\(x^{tar}_-\) 表示 negative。距离函数 d 是汉明距离。其公式如下:

\(d\left(\boldsymbol{x}^{c u r}, \boldsymbol{x}^{t a r}\right)=\left\|\boldsymbol{b}^{c u r}\left(\boldsymbol{x}^{c u r}\right)-\boldsymbol{b}^{t a r}\left(\boldsymbol{x}^{t a r}\right)\right\|_{2}\)

样本生成

为了找到 target 的正样本,我们采用的是用 Pose 真值获取的方式,也就是将 current 帧的 feature 投影到 target 帧,然后通过寻找最近临的方式获得正样本,投影方式如下:

\(\boldsymbol{x}_{i,+}^{t a r}=\boldsymbol{\pi}^{-1}\left(\mathbf{R}_{g t} \boldsymbol{\pi}\left(\boldsymbol{x}_{i}^{c u r}, z_{i}\right)+\boldsymbol{t}_{g t}\right)\)

负样本生成(Exhaustive Negative Sample Mining)如下,该函数会根据阈值 c 返回负样本,也就是大于某个距离的认为是负样本(由于 k nearest 排序,返回的可以认为是最接近正样本的负样本):

其中的  \(b\) 表示量化成二值特征(256bit)的向量。

2.4 检测子学习

2.4.1 误差函数

GCN 没有像其他一些 local feature 网络一样学习一张热力图之类,而是直接生成了一张表示是否为关键点的二值图(0 表示不是关键点,1 表示是关键点),那么整个检测子就变成了一个分类问题。所以使用的就是分类问题的 Cross-entropy loss,由于正负样本数不均衡,所以使用 weighted cross-entropy loss。检测子 loss 如下:

\(L_{d e t}=L_{c e}\left(\boldsymbol{o}^{c u r}\left(\boldsymbol{x}^{c u r}\right)\right)+L_{c e}\left(\boldsymbol{o}^{tar}\left(\boldsymbol{x}_{+}^{tar}\right)\right)\)

其中 \(l_{ce}\) 就是 weighted cross-entropy loss 定义如下,其中 \(\alpha_1\) 和 \(\alpha_2\) 是正负样本权重系数:

\(\begin{aligned} L_{c e}(\boldsymbol{o}(\boldsymbol{x}))=&-\alpha_{1} \sum_{i} c_{\boldsymbol{x}_{i}} \log \left(\boldsymbol{o}\left(\boldsymbol{x}_{i}\right)\right) \\ &-\alpha_{2} \sum_{i}\left(1-c_{x_{i}}\right) \log \left(1-\boldsymbol{o}\left(\boldsymbol{x}_{i}\right)\right) \end{aligned}\)

2.4.2 样本生成

文章中使用 Shi-Tomasi corners (16×16)生成检测子,并对不同帧进行了 warp 操作获得另一帧上检测子位置。为什么不用每帧检测 Shi-Tomasi corners 而用投影方式呢?作者的解释是:This leads to better distribution of keypoints and the objective function directly reflects the ability to track the keypoints based on texture. 也就是说这样生成样本能够反映关键点的跟踪能力。也就是在不同帧关键点位置的鲁棒性。但另一方面,我没有特别明白的是,如果 warp 过去的位置,实际上被遮挡或者什么原因并不能匹配,这种方式会不会选择了错误的匹配对呢?

3 训练细节

总的 Loss 是 \(L_{feat}\) 和 \(L_{det}\) 的和,不过赋予了不同权重。在 GCNv2 中分别赋予权重 100 和 1。其他上述公式中提到的一些参数如下:

参数 赋值
\(m\)
1
\(\alpha_1\)
0.1
\(\alpha_2\)
1.0
Adam

初始值 \(10^{-4}\),每40 epoch 下降一半

4 实验结果

4.1 整体流程

GCNv2 整体流程如下图所示,其实就是整个代替了 ORB 的检测、描述与匹配过程。

在实验结果方面,GCNv2 取得了比 ORB-SLAM2 更高的精度,同时在 GPU 可以保持实时性能。

4.2 VO 跟踪结果

作者做了比较详细的比较,可以看到基本上精度秒杀 ORB ,即使和 SIFT 之类比也不逊色。比较有意思的是 SuperPoint 这种同样是学习的特征,其实换到 ORB SLAM 里面并没有明显优势,甚至很多数据集上还不如 ORB。这和我们之前实验的结果也比较相近。不过这里面有很多 trick,比如特征匹配时间影响 Local Mapping 优化速度,亚像素精度等等都是有影响的,并不能简单从这里说明描述子好效果就一定好。

4.3 加入闭环 SLAM 结果

可以看出加入闭环后,ORB SLAM2 效果还是不错的并没有 VO 那么明显,在部分测试集中表现甚至高于 GCN SLAM。不过整体而言 ORB SLAM2 有两个失败,而 GCN SLAM 表现稳定,显示出在鲁棒性方面还是有一定的优势。

4.4 匹配效率

作者比较了 GCNv2-tiny 和 ORB 的匹配效率:

可以看到 GCNv2-tiny 提取点个数更少而内点个数更多,表明其匹配效率较高。

点评

使用 Deep Learning 方式,但是实际上学习的是传统的检测子,感觉没有 MagicPoint 和 SuperPoint 定义更加贴近 ”语义特征点“。不过相对来说生成数据非常容易,仅依赖于 Pose 真值(往往是把 Vicon 或者激光这种当做真值)。可以理解为学习这些描述子是为了 ”fit“ pose(迁移性能有待验证,比如有些描述子遮挡了根本就不能匹配)。整体效果也还不错,而且二值描述子在存储空间和速度上有优势。

参考文献

  1. Binarized neural networks: Training neural networks with weights and activations constrained to+ 1 or-1笔记:https://zhuanlan.zhihu.com/p/24679842
  2. Real-Time Single Image and Video Super-Resolution Using an Efficient Sub-Pixel Convolutional Neural Network 笔记:https://blog.csdn.net/g11d111/article/details/82855946
  3. https://zhuanlan.zhihu.com/p/26600748
  4. https://zhuanlan.zhihu.com/p/60049140

EM估计(Expectation Maximization Algorithm)

这篇文章实际上是Ng在CS229上一个课程讲座。因为EM估计也很常用,特别是语义SLAM里面,这里根据一些粗浅的理解做一个总结备忘。

EM算法1是分为E(期望)和M(极大)的两步迭代优化算法,主要用于解决数据缺失情况下的概率估计问题。这个数据缺失其实不太好理解,有一个男女生身高的经典例子比较形象2,可以从ML到EM估计分别理解下。

1、极大似然估计(ML)

我们需要调查我们学校的男生和女生的身高分布。 假设你在校园里随便找了100个男生和100个女生。他们共200个人。将他们按照性别划分为两组,然后先统计抽样得到的100个男生的身高。假设他们的身高是服从正态分布的。但是这个分布的均值 \(\mu\) 和方差 \(\sigma^2\) 我们不知道,这两个参数就是我们要估计的。记作 \(\theta=[\mu,\sigma]^T\)。 这是一个典型的极大似然估计问题,我们知道100个男生,100个女生和他们的身高,我们想要估计男生的身高分布均值方差,女生的身高分布均值方差,都可以用极大似然估计来解决。

我们已知的有两个:样本服从的分布模型、随机抽取的样本;我们未知的有一个:模型的参数。根据已知条件,通过极大似然估计,求出未知参数。总的来说:极大似然估计就是用来估计模型参数的统计学方法。

因为我们知道哪个样本是男生哪个样本是女生,所以这一问题就分解成分别估计男生和女生的\(\theta=[\mu,\sigma]^T\)即可。那么以男生为例,设样本集 \(X={x_{1},x_{2},…,x_{N}}\),其中 \(N=100\),概率密度函数 \(p(x_{i}|\theta)\),首先写出似然函数也就是所有男生概率的乘积:

\(l(\theta)=L\left(x_{1}, \cdots, x_{n} ; \theta\right)=\prod_{i=1}^{n} p\left(x_{i} ; \theta\right), \theta \in \Theta \tag{1}\)

通常转为求解对数似然函数:

\(H(\theta)=\ln L(\theta)=\ln \prod_{i=1}^{n} p\left(x_{i} ; \theta\right)=\sum_{i=1}^{n} \ln p\left(x_{i} ; \theta\right) \tag{2}\)

极大似然估计就是求似然函数最大值(或者对数似然函数最小值):

\(\hat{\theta}=\arg \max _{\theta} l(\theta)=\arg \max _{\theta} \prod_{i=1}^{N} p\left(x_{i} | \theta\right) \tag{3}\)

这个极值自然就是通过求导得到:

\(\frac{d l(\theta)}{d \theta}=0 或者 \frac{d H(\theta)}{d \theta}=\frac{d \ln l(\theta)}{d \theta}=0 \tag{4}\)

特别的,对于正态分布样本 \(N\left(\mu, \sigma^{2}\right)\),其似然函数为:

\(L\left(\mu, \sigma^{2}\right)=\prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi} \sigma} e^{-\frac{\left(x_{i}-\mu\right)^{2}}{2 \sigma^{2}}}=\left(2 \pi \sigma^{2}\right)^{-\frac{n}{2}} e^{-\frac{1}{2 \sigma^{2}} \sum_{l=1}^{n}\left(x_{i}-\mu\right)^{2}} \tag{5}\)

对数为:

\(\ln L\left(\mu, \sigma^{2}\right)=-\frac{n}{2} \ln (2 \pi)-\frac{n}{2} \ln \left(\sigma^{2}\right)-\frac{1}{2 \sigma^{2}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2} \tag{6}\)

求导得方程组:

\(\left\{\begin{array}{l}{\frac{\partial \ln L\left(\mu, \sigma^{2}\right)}{\partial \mu}=\frac{1}{\sigma^{2}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)} \\ {\frac{\partial \ln L\left(\mu, \sigma^{2}\right)}{\partial \sigma^{2}}=-\frac{n}{2 \sigma^{2}}+\frac{1}{2 \sigma^{4}} \sum_{i=1}^{n}\left(x_{i}-\mu\right)^{2}=0}\end{array}\right. \tag{7}\)

则似然方程的唯一解为:

\(\left\{\begin{array}{l}{\mu^{*}=\overline{x}=\frac{1}{n} \sum_{i=1}^{n} x_{i}} \\ {\sigma^{* 2}=\frac{1}{n} \sum_{i=1}^{n}\left(x_{i}-\overline{x}\right)^{2}}\end{array}\right. \tag{8}\)

这个就是它的最大值点。也就是通过极大似然估计给出的符合当前样本的参数估计。

2、Jensen不等式

对于一个凸函数 \(f(X)\)(凸函数的定义为:对于X为实变量,\(f^{\prime \prime}(x) \geq 0 \text { (for all } x \in \mathbb{R} )\),对于X为向量,),其中X是随机变量则有如下不等式成立:

\(E[f(X)] \geq f(E X) \tag{9}\)

其中等号成立的条件当且仅当\(X=EX\),也就是随机变量为一个常数c:\(X=c\)。

Jensen不等式直观理解可以参见下图:

图中,实线f是凸函数,X是随机变量,有0.5的概率是a,有0.5的概率是b。(就像掷硬币一样)。X的期望值就是a和b的中值了,图中可以看到 \(E[f(X)] \geq f(E X)\) 成立。

对于f是凹函数的情况,则正相反:\(E[f(X)] \leq f(E X)\)。

3、最大期望估计(EM)

对比极大似然估计,EM估计理解起来非常简单,例如对于上述男女生的例子,我们现在有200个男女生的样本,但是我们不知道哪个样本属于男生哪个样本属于女生(信息缺失)。

3.1 EM估计公式推导

样本集 \(X={x_{1},x_{2},…,x_{m}}\),包含 m 个独立的样本;每个样本 \(x_i\) 对应的类别 \(z_i\)  是未知的(即上文中每个样本属于哪个分布是未知的);我们需要估计概率模型 \(p(x,z)\) 的参数 \(\theta\),即需要找到适合的 \(\theta\) 和 \(z\) 让 \(L(\theta)\) 最大。根据上文极大似然估计中的似然函数取对数所得 \(logL(\theta)\),可以得到如下式:

\(\begin{align*} \ell(\theta) &=\sum_{i=1}^{m} \log p(x ; \theta) \tag{10}\\ &=\sum_{i=1}^{m} \log \sum_{z} p(x, z ; \theta) \tag{11}\end{align*}\)

其中第二项 \(z\) 的求和在样本缺失情况下不容易做到。但是如果 \(z\) 的概率确定后,这个求和就很容易了,EM估计就是分成两个步骤来迭代实现了这样一个估计过程。

定义 \(Q_i(z)\) 表示 \(z\) 的概率,显然 \(\sum_{z} Q_{i}(z)=1, Q_{i}(z) \geq 0\)。将上述公式展开:

\(\begin{align*} \sum_{i} \log p\left(x^{(i)} ; \theta\right) &=\sum_{i} \log \sum_{z^{(i)}} p\left(x^{(i)}, z^{(i)} ; \theta\right) \tag{12} \\ &=\sum_{i} \log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{13} \\ & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{14}\end{align*}\)

式 (13) 就是分子分母同乘概率 \(Q_i(z)\),但是显然和的对数求导并不好求;式 (14) 根据 Jensen不等式而来,其中 \(f(x)=log(x)\) 是一个凹函数,同时设 \(X=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\),则 \(EX=\sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right)\left[\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}\right]\),根据 Jenson 不等式 \(f(EX) \geq E[f(X)]\),则有:

\(\log \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \geq \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{15}\)

再对i求和就得到了我们想要的公式。

根据之前的结论,仅当随机变量为常数时等号成立则有:

\(\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)}=c \tag{16}\)

显然 c 与 \(z_i\) 无关,则说明 \(Q_{i}\) 只和 \(p\left(x^{(i)}, z^{(i)} ; \theta\right)\) 有关,我们可以用后者作为前者的概率表达:

\(Q_{i}\left(z^{(i)}\right) \propto p\left(x^{(i)}, z^{(i)} ; \theta\right) \tag{17}\)

则显然可以推导:

\(\begin{align*} Q_{i}\left(z^{(i)}\right) &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{\sum_{z} p\left(x^{(i)}, z ; \theta\right)} \\ &=\frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{p\left(x^{(i)} ; \theta\right)} \\ &=p\left(z^{(i)} | x^{(i)} ; \theta\right) \tag{20} \end{align*} \)

至此,我们推出了在固定其他参数 \(\theta\) 后,\(Q_{i}\left(z^{(i)}\right)\) 的计算公式就是后验概率,解决了 \(Q_{i}\left(z^{(i)}\right)\) 如何选择的问题。这一步就是E步,建立 \(logL(\theta)\) 的下界。接下来的M步,就是在给定 \(Q_{i}\left(z^{(i)}\right)\) 后,调整\(\theta\),去极大化 \(logL(\theta)\) 的下界,是一个典型的极大似然估计过程(在固定 \(Q_{i}\left(z^{(i)}\right)\) 后,下界还可以调整的更大)。

如果继续用上面形象一点的例子就是:

E步:估计出男生和女生的概率分布;

M步:在上述概率分布的前提下,用极大似然估计得到当前最优的参数。

反复迭代直到收敛。

3.2 EM估计的基本步骤

基于以上推导,EM估计的典型迭代步骤如下:

循环重复直到收敛 {

(E-步骤) 对于每一个i,计算概率: \(Q_{\mathrm{i}}\left(z^{(\mathrm{i})}\right) :=p\left(z^{(i)} | x^{(\mathrm{i})} ; \theta\right)\)

(M-步骤) 使用极大似然估计最优参数:\( \theta :=\arg \max {\theta} \sum{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \)

}

3.3 EM估计有效性的证明

我们需要证明上述的步骤一定可以得到收敛的结果,也就是 \(\ell\left(\theta^{(t)}\right) \leq \ell\left(\theta^{(t+1)}\right)\) 我们整体的极大似然估计在每次迭代后一定单调增加,那么我们一定可以在一定次数的迭代后极大似然估计不再增加也就得到了整个系统的概率最大值的参数估计,这一参数就同时包含了极大似然估计 \(\theta\) 以及缺失的参数 \(z\)。

证明如下:

对于给定 \(\theta^{(t)}\) 我们在 E 步骤选定概率函数如下:

\(Q_{i}^{(t)}\left(z^{(i)}\right) :=p\left(z^{(i)} | x^{(i)} ; \theta^{(t)}\right) \tag{21}\)

根据之前证明,这一概率函数 \(Q_{i}^{(t)}\left(z^{(i)}\right)\) 使得前述公式中的 Jensen 不等式等号成立也就是:

\(\ell\left(\theta^{(t)}\right)=\sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{22}\)

则有如下推导:

\(\begin{align*} \ell\left(\theta^{(t+1)}\right) & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t+1)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{23} \\ & \geq \sum_{i} \sum_{z^{(i)}} Q_{i}^{(t)}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta^{(t)}\right)}{Q_{i}^{(t)}\left(z^{(i)}\right)} \tag{24} \\ &=\ell\left(\theta^{(t)}\right) \tag{25} \end{align*}\)

其中式 (23) 与前述 3.1 中的公式推导一样,只是 \(\theta^{(t)}\) 换成了 \(\theta^{(t+1)}\)。式 (24) 是因为 \(\theta^{(t+1)}\) 是对公式 (22) 进行极大似然估计的结果,也就是 \(\theta^{(t+1)}\) 是使得公式中 \(\theta^{(t)}\) 最大的那个值:

\(\theta^{t+1} = \arg \max _{\theta} \sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{26}\)

显然用 \(\theta^{(t+1)}\) 代替 \(\theta^{(t)}\) 是递增的。

根据上面的证明,可以知道 \(\ell\left(\theta^{(t)}\right) \leq \ell\left(\theta^{(t+1)}\right)\),那么我们一定能够通过 EM估计得到最优的结果。

如果我们定义:

\(J(Q, \theta)=\sum_{i} \sum_{z^{(i)}} Q_{i}\left(z^{(i)}\right) \log \frac{p\left(x^{(i)}, z^{(i)} ; \theta\right)}{Q_{i}\left(z^{(i)}\right)} \tag{27}\)

则EM估计可以看过是J的坐标上升法:,E步固定 \(\theta\),优化 \(Q\),M步固定 \(Q\) 优化 \(\theta\)。

参考资料

1.
Arthur D, Nan L, Donald R. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B. 1977;39(1):1–38.
2.
漫谈 Clustering (番外篇): Expectation Maximization « Free Mind. 漫谈 Clustering (番外篇): Expectation Maximization. http://blog.pluskid.org/?p=81. Accessed May 29, 2019.
3.
EM算法(Expectation Maximization Algorithm)详解. EM算法(Expectation Maximization Algorithm)详解. https://blog.csdn.net/zhihua_oba/article/details/73776553. Accessed May 29, 2019.
4.
Expectation–maximization algorithm. Expectation–maximization algorithm. https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm. Accessed May 29, 2019.
5.
极大似然估计详解. 极大似然估计详解. https://blog.csdn.net/zengxiantao1994/article/details/72787849. Accessed May 29, 2019.
6.
(EM算法)The EM Algorithm – JerryLead – 博客园. (EM算法)The EM Algorithm. https://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html. Accessed May 29, 2019.

论文材料

 

论文笔记:NetVLAD: CNN architecture for weakly supervised place recognition

NetVLAD1是一个较早的使用 CNN 来进行图像检索或者视频检索的工作,后续在此工作的基础上陆续出了很多例如 NetRVLAD、NetFV、NetDBoW 等等的论文,思想都是大同小异。

一、图像检索

VLAD 和 BoW、Fisher Vector 等都是图像检索领域的经典方法,这里仅简介下图像检索和 VLAD 的基本思想。

图像检索(实例搜索)是这样的一个经典问题:

1、我们有一个图像数据库 \(I_i\) 通过函数可以得到每一个图像的特征 \(f(I_i)\);

2、我们有一个待查询图像 \(q\) 通过函数得到它的特征 \(f(q)\);

3、则我们获得的欧氏距离 \(d(q, I) = \parallel  f(q) – f(I)\parallel\) 应该满足越相近的图像 \(d(q, I)\) 越小。

二、VLAD (Vector of Locally Aggregated Descriptors)

而 VLAD (Vector of Locally Aggregated Descriptors)则是图像检索中一个经典方法,也就可以理解为怎么获得上述的 \(f\) 的方法,记为 \(f_{VLAD}\)。通常在传统方法中我们会获得一系列的局部特征(SIFT、SURF、ORB)之类,假设为 N 个 D 维的局部特征(通常 N 可能比较大,而且每幅图特征多少不一,N 可能数量也不一定),我们希望通过这 N*D 维特征获得一个可以表示图像全局 K*D 维特征的方法(通常K是我们指定的数目,例如128维)。VLAD 的主要流程如下:

1、对全部 N*D 维局部特征进行 K-Means 聚类获得 K 个聚类中心,记为 \(C_k\)

2、通过以下公式将 N*D 维局部特征编写为一个全局特征 V,特征向量维数为 K*D,其中 \(k\in K\) \(j\in D\),公式如下:

\(\begin{equation}V(j, k) = \sum_{i=1}^N a_k(x_i) (x_i(j) – c_k(j))\end{equation}\)

其中 \(x_i\) 为第i个局部图像特征,\(c_k\) 为第k个聚类中心,\(x_i\) 和 \(c_k\) 都是 D 维向量。\(a_k(x_i)\) 是一个符号函数,当且仅当 \(x_i\) 属于聚类中心 \(c_k\) 时,\(a_k(x_i)=1\) ,否则 \(a_k(x_i)=0\)。

经过上面对于经典 VLAD 方法的解释,我们可以看出来 \(f_{VLAD}\) 是一个将若干局部特征压缩为一个特定大小全局特征的方法,通过聚类,实现了将特征降维,同时用特征与聚类中心的差值作为新的特征值,在实践中 VLAD 方法具有较好的检索效果。

三、NetVLAD

1、经典 VLAD 公式的可微化

经典 VLAD 方法显然是一个不可导的函数,其中主要不可导的地方在于 \(a_k(x_i)\) 这样一个符号函数,因此为了将 VLAD 变可训练的函数,必须将其变成可微计算。作者将 \(a_k(x_i)\) 平滑化:

\(\bar{a}_k(x_i)=\frac{e^{-\alpha\left \| x_i – c_k \right \|^2}}{\sum_{{k}’}e^{-\alpha\left \| x_i – c_{{k}’} \right \|^2}}\in (0,1)\)

这里面 \(\alpha\) 是一个正值参数,显然当 \(\alpha\rightarrow \infty \) 时,\(\bar{a}_k(x_i)\) 更趋近于 0 和 1 的两级。

再将 \(-\alpha \left \| x_i – c_k \right \|^2\) 展开,显然 \(e^{-\alpha\left \| x_i \right \|^2}\) 可以被约掉,因此得到:

\(\bar{\alpha}_k(x_i)=\frac{e^{w^T_kx_i+b_k}}{\sum_{{k}’}e^{w^T_{{k}’}+b_{{k}’}}}\)

其中 \(w_k=2\alpha c_k\) ,\(b_k=-\alpha\left \| c_k \right \|^2\),其实仔细观察这个公式,就是一个 softmax。

这样 VLAD 公式就被改写为:

\(\begin{equation}V(j, k) = \sum_{i=1}^N\frac{e^{w^T_kx_i+b_k}}{\sum_{{k}’}e^{w^T_{{k}’}+b_{{k}’}}}(x_i(j) – c_k(j))\end{equation}\)

显然这里面 \(w_k\) 、\(b_k\) 和 \(c_k\) 都是 NetVLAD 需要学习的参数。

2、NetVLAD 通过监督学习获得聚类中心的好处

我们这样将 \(c_k\) 作为学习参数有什么好处呢,论文中用了一幅图进行了直观解释:

传统 VLAD 的中心是聚类出来的,没有监督的标签数据 \(c_k^{VLAD}\) ,在聚类时我们使用了很多图像这些图像的描述符之间没有关系,那么也就很可能把本来不是一个物体的描述符聚为一类,使得我们原本期望的类内描述符都是一个物体的feature不太容易达到。而在使用监督数据进行训练时,我们可以已知图像数据属于同一物体,那么训练时就可以只把属于同一物体的特征聚在一起而把不是的划分到其他类别,这样就可能学习出一个更好的 \(c_k^{NetVLAD}\) 聚类中心,使得最终的特征更有区分度。

3、网络实现

首先,由于是 NN 的方法,我们这里使用 CNN Feature 代替了传统 VLAD 中的 N 个局部描述子,CNN 是一个全局的特征,它的 Feature Map 是 W*H*D 大小,那么类比于我们之前的传统方法 N*D,我们这里 NetVLAD 目标就是将 W*H*D (N=W*H)的特征转换为 K*D 的特征;

其次,我们将整个 NetVLAD 看做一个 pooling layer,它的作用是实现降最终和 VLAD 一样获得我们想要的 K*D 维描述子。

具体实现上,我们就分步骤来做,这部分作者的 PPT 里面有张图非常清晰:

1、实现 \(z_k=w_k^Tx_i+b_k\),也就是公式中的蓝色部分,论文里直接通过一个1×1的卷积来做。这也是1×1卷积的一个应用;

2、实现 \(\sigma_k (z)=\frac{e^{z_k}}{\sum_{{k}’} e^{z_{{k}’}} }\),也就是公式中的黄色部分,如之前所述这实际上就是一个 softmax 公式,论文里直接通过一个 softmax 来做;

3、实现 \(x_i (j)- c_k(j)\),也就是公式中的绿色部分,这部分就是一个减法,直接用一个 VLAD core来做;

4、1~3已经实现了 \(V(j,k)\) 的计算了,后面按照 All about VLAD 论文还要对 \(V(j,k)\) 做两步简单的归一化(这篇论文证明了归一化可以提升检索性能),包括:

1)intra-normalization

这个主要意思是将每一个 VLAD block(也就是每个聚类中心的所有残差)分别作 l2 normalization。

2) l2 normalization

这个是将所有 VLAD 特征一起做 l2 normalization。

四、弱监督训练

1、弱监督数据的获取

在图像检索中,训练通常需要明确的图像到底是哪个物体的标注。而论文针对的领域是地点识别,因此可以利用图片的位置信息,将 Google Street View 的数据作为训练图片,将同一地点、不同视角、不同时间的数据作为同一物体标签进行训练。相对并没有明确的标注图上面都是那个物体,只是说这几幅图中可能看到同样的一个物体,数据比较容易获取。这在机器人领域同样也是可行的。

2、弱监督 triplet ranking loss

我们将整个网络看做一个特征提取函数 \(f_{\theta}\) ,那我们的目标自然就是:对于一个需要检索的图像 \(q\) ,我们有一个数据库 \(I\),我们希望位置最近的图片\(I_{i*}\) 的特征欧氏距离,比其他所有图片\(I_{i}\) 的距离都要小:

\(d_{\theta}(q, I_{i*}) < d_{\theta}(q, I_i)\)

然而理想情况是这样,如刚才小节所述我们实际上不能获取最近的图片 \(I_{i*}\),但我们可以获取相近的图片集合(作为正样本) \(p_{i*}^q\),这些正样本集合一定包含了我们想要的 \(I_{i*}\),但是我们不知道是哪一个;同样,我们也可以获取绝对不可能是同一物体的负样本集合 \(n_j^q\) ,比如远离所选地点的图片。那么我们至少可以期望获得这样的结果:所有正样本集合的特征距离,应该要比负样本集合的特征距离要小:

\(d_{\theta}(q, p_{i}^q) < d_{\theta}(q, n_j^q)\)

并且正样本中距离最近的应该就是我们想要的最近图片:

\(p_{i*}^q=\underset{p_i^q}{argmin}d_{\theta}(q, p_i^q)\)

据此,我们构造如下 loss:

\(L_{\theta}=\underset{j}{\sum}l\left ( \underset{i}{\min} d_{\theta}^2 (q,p_i^q)+m-d_{\theta}^2(q,n_j^q)\right )\)

其中 \(l(x)=max(x,0)\) 也就是所说的 hinge loss,m是一个常量表示我们给负样本设置的 margin。这个 loss 与 triplet loss 是很像的。

3、训练&实现

作者在附录中给出了详细的实现细节,网络提取部分作者使用 VGG-16 但做了一些修改,对于 VLAD 的聚类中心使用 K=64,m=0.1。

训练 lr = 0.001 或 0.0001,优化器 SGD,momentun 0.9,weight decay 0.001,batchsize=4 turples。这里不一一例举。具体还是看作者的代码比较好。

论文下载

1、PDF:NetVLAD: CNN architecture for weakly supervised place recognition

2、PPT:cvpr16_NetVLAD_presentation.pptx

3、Code:

https://github.com/Relja/netvlad (作者的,使用 MatConv 框架)

https://github.com/sitzikbs/netVLAD/blob/master/netVLAD.py (第三方实现,不过只有 VLAD layer 部分)

https://github.com/shamangary/LOUPE_Keras (第三方实现,不过只有 VLAD layer 部分)

参考文献

  1. 1.
    Arandjelovic R, Gronát P, Torii A, Tomás Pajdla, Sivic J. NetVLAD: CNN architecture for weakly supervised place recognition. CoRR. 2015;abs/1511.07247.

Tensorflow C API 从训练到部署:使用 C API 进行预测和部署

前述博文 Tensorflow C++ 从训练到部署(2):简单图的保存、读取与 CMake 编译Tensorflow C++ 从训练到部署(3):使用 Keras 训练和部署 CNN 使用 Tensorflow/Keras 的 Python API 进行训练,并使用 C++ API 进行了预测。由于 C++ API 需要编译 Tensorflow 源码,还是比较麻烦的。而 Tensorflow 官方提供了 C API 编译好的库文件,相对部署上比较容易(直接复制库文件到自己的工程即可),本文将介绍使用 C API 进行预测的方法。对于 Python 训练部分,与前述文章相同不做赘述。

0、系统环境
Ubuntu 16.04
Tensorflow 1.12.0

1、安装依赖
1、GPU 支持安装(可选)
CUDA 9.0
cnDNN 7.x

2、Tensorflow 1.12.0
下载地址:
https://www.tensorflow.org/install/lang_c

其中 1.12.0 的下载地址如下(我这里提供了包含TX2 aarch64在内的几个版本):

TensorFlow C libraryCUDAcuDNNURL
Linux x86_64 CPUxxhttps://pan.baidu.com/s/1FDdXCgtJJlDJP8ziDs6dow
Linux x86_64 GPU9.07.xhttps://pan.baidu.com/s/1qxDntkQ-rcgvp1xxrSKW0w
macOS CPUxxhttps://pan.baidu.com/s/1F6NdNtCxg11P_EpEdsqttA
Linux aarch64 GPU (TX2)9.07.0.5https://pan.baidu.com/s/1mI76203wY9Nd5US4sH5-pg

将库解压到 third_party/libtensorflow 目录。

如果上面的版本都不符合你的需求,你可以参照这篇文章编译你需要的版本。

2、TFUtils 工具类
为了简便起见,我们首先将常用的 C API 封装为

1)文件 utils/TFUtils.hpp:

2)文件 utils/TFUtils.cpp:

3、简单图的读取与预测
在前述文章 Tensorflow C++ 从训练到部署(2):简单图的保存、读取与 CMake 编译 中我们已经介绍了一个 c=a*b 的简单“网络”是如何计算的。

其中 Python 构建网络和预测部分就不重复了,详见该文所述。这里直接给出 C API 的代码:

文件名:simple/load_simple_net_c_api.cc

简单解释一下:

这一行是加载 pb 文件。

这一段是创建两个输入 tensor 以及输入的 ops。注意这里的 CreateTensor 在后面都需要调用 DeleteTensors 进行内存释放。输出的 tensors 还没创建先定义为 nullptr。

这一行是运行网络。

这两行是从输出的 output_tensors 读取数据到一个二维vector const std::vector<:vector>>,我们这里输出只有 “c” 一个名字,而且只有一个索引 0,因此直接取出 data[0] 就是我们原本想要的输出。

编译运行这一文件,如果没有问题则会得到如下输出:

4、CNN的读取与预测
与刚才小节3相似,CNN网络也是一样的流程,还是以最基本的 fashion_mnist 为例,该网络的训练和保存流程请参考之前的文章。这里我们仅介绍 C API 进行预测的部分。由于我们这里需要读取一幅图并转化成 Tensor 输入网络,我们构造一个简单的函数 Mat2Tensor 实现这一转换:

1)Met2Tensor 部分文件:fashion_mnist/utils/mat2tensor_c_cpi.h

2)网络读取与预测,这部分与刚才的小节3基本一样,就不做解释了:

编译运行这一文件,如果没有问题,则会得到如下输出:

到此,我们就完成了使用 C API 运行 Tensorflow Model 的流程。

本文中的全部代码均已开源:
https://github.com/skylook/tensorflow_cpp

参考文献
[1] https://github.com/Neargye/hello_tf_c_api
[2] https://www.tensorflow.org/api_docs/python/tf/contrib/saved_model/save_keras_model
[3] https://www.tensorflow.org/install/lang_c