本文是“图匹配入门学习笔记——以《Factorized Graph Matching》为例”系列的第二篇文章，主要介绍了FGM算法的实现过程，包括图匹配的因式分解形式，一种统一的图匹配框架和基于路径跟随策略的求解算法。

因式分解的图匹配（FGM）

关联矩阵KKK是图匹配问题中的关键，他编码了两幅图之间的一阶和二阶关系。但KKK的两个特性：满秩性（Full-rankness）和不定性（Indefinitenes）给许多传统方法（如谱方法和基于梯度的方法）带来了困难。上图中的第一行分别展示了三个数据集中矩阵KKK的秩(rank(K)rank(K)rank(K))和行数(n1n2n_1n_2n1n2)之间的比值，可以看到比值全为1，说明KKK通常是满秩的，这使得将KKK近似为秩为1的矩阵的谱方法不可避免地失去准确性。第二行展示了KKK的最大特征值eigmax(K)eig_{max}(K)eigmax(K)和最小特征值eigmin(K)eig_{min}(K)eigmin(K)之间的比值，可以看到比值均小于-0.4，即特征值异号，说明KKK是不定矩阵，这说明最大化vec(X)TKvec(X)vec(X)^TKvec(X)vec(X)TKvec(X)的过程是非凸的，使得基于梯度的方法容易陷入局部最优。了解了图匹配问题的定义和现有方法存在的问题之后，我们来介绍FGM算法是如何解决图匹配问题的。

由前文可知，关联矩阵KKK是由对角线元素（节点关联矩阵KpK_pKp）和非对角线元素（边关联矩阵KqK_qKq）构成的，如上图所示。除掉对角线元素，KKK是一个稀疏的块矩阵且具备以下三个结构特点：

KKK由n2×n2n_2\times n_2n2×n2个更小的块Kij∈Rn1×n1K_{ij}\in \mathbb{R}^{n_1\times n_1}Kij∈Rn1×n1构成；
有些块KijK_ijKij是零矩阵，这说明图G2\mathcal{G}_2G2中第iii个和第jjj个节点之间没有边连接，这些零矩阵可以通过G2H2TG_2H_2^TG2H2T来检索，即如果[G2H2T]ij=0[G_2H_2^T]_{ij}=0[G2H2T]ij=0，则对应的块Kij=0K_{ij}=0Kij=0
对于非零的块KijK_{ij}Kij有封闭形式的解G1diag(kcq)H1TG_1diag(k^q_c)H_1^TG1diag(kcq)H1T，其中ccc是图G2\mathcal{G}_2G2中第iii个和第jjj个节点之间边的下标。

基于上述结构特点，关联矩阵KKK可以被分解为：
K=diag(vec(Kp))+(G2⊗G1)diag(vec(Kq))(H1⊗H1)T(6)K=diag(vec(K_p))+(G_2\otimes G_1)diag(vec(K_q))(H_1\otimes H_1)^T \quad (6)K=diag(vec(Kp))+(G2⊗G1)diag(vec(Kq))(H1⊗H1)T(6)
其中⊗\otimes⊗表示矩阵的Kronecker积。利用上式大矩阵KKK被分解为6个小矩阵G1,H1,G2,H2,Kp,KqG_1,H_1,G_2,H_2,K_p,K_qG1,H1,G2,H2,Kp,Kq的组合，这使得算法复杂度由O(n12n22)O(n_1^2n_2^2)O(n12n22)下降到O(n1m1+n2m2+n1n2+m1m2)O(n_1m_1+n_2m_2+n_1n_2+m_1m_2)O(n1m1+n2m2+n1n2+m1m2)。将公式(6)带入公式(3)（见前文）中可得等价的图匹配目标函数
Jgm(X)=tr(KpTX)+tr(KqTY)(7)J_{gm}(X)=tr(K_p^TX)+tr(K_q^TY) \quad (7)Jgm(X)=tr(KpTX)+tr(KqTY)(7)其中Y=G1TXG2∘H1TXH2(8)Y=G_1^TXG_2\circ H_1^TXH_2\quad (8)Y=G1TXG2∘H1TXH2(8)
∘\circ∘表示矩阵的Hadamard积，Y∈{0,1}m1×m2Y\in\{0,1\}^{m_1\times m_2}Y∈{0,1}m1×m2表示边和边之间的对应矩阵，yc1c2=1y_{c_1c_2}=1yc1c2=1表示图G1\mathcal{G}_1G1中第c1c_1c1个边与第c2c_2c2个边匹配。

图匹配的统一框架

前文中我们提到图匹配有两种形式化方式：Lawler’s QAP和Koopmans-Beckmann’s QAP，两者之间的区别与联系我们也做了简单的讨论。FGM基于因式分解的方法提出一种统一的图匹配框架。首先，边关联矩阵KqK_qKq可以通过SVD的方式分解为Kq=UVTK_q=UV^TKq=UVT，利用KqK_qKq本身的低秩结构特点，公式(7)可以重写为：

前文说过关联矩阵KKK具有满秩性，为啥这里又说KqK_qKq具备低秩结构特点呢。我理解是因为KKK的满秩性主要是来自于对角线元素KpK_pKp，KpK_pKp表示的是节点之间的相似性。而边的关联矩阵KqK_qKq本身是比较稀疏的，因为不是所有的节点之间都有边连接，导致它的秩也比较低。

Jgm(X)=tr(KPTX)+tr((UVT)T(G1TXG2∘H1TXH2))=tr(KPTX)+tr((∑i=1cuiviT)T(G1TXG2∘H1TXH2))=tr(KPTX)+∑i=1ctr(Ai1XAi2XT)(9)J_{gm}(X)=tr(K_P^TX)+tr((UV^T)^T(G_1^TXG_2\circ H_1^TXH_2))\\ \quad\quad\quad\quad\quad=tr(K_P^TX)+tr((\sum_{i=1}^cu_iv_i^T)^T(G_1^TXG_2\circ H_1^TXH_2))\\ \quad\quad\quad\quad\quad=tr(K_P^TX)+\sum_{i=1}^ctr(A_i^1XA_i^2X^T)\quad\quad\quad\quad\quad\quad(9)Jgm(X)=tr(KPTX)+tr((UVT)T(G1TXG2∘H1TXH2))=tr(KPTX)+tr((i=1∑cuiviT)T(G1TXG2∘H1TXH2))=tr(KPTX)+i=1∑ctr(Ai1XAi2XT)(9)
其中Ai1=G1diag(ui)H1T∈Rn1×n1,Ai2=G2diag(vi)H2T∈Rn2×n2A_i^1=G_1diag(u_i)H_1^T\in\mathbb{R}^{n_1\times n_1},A_i^2=G_2diag(v_i)H_2^T\in\mathbb{R}^{n_2\times n_2}Ai1=G1diag(ui)H1T∈Rn1×n1,Ai2=G2diag(vi)H2T∈Rn2×n2可以理解为两个图对应的邻接矩阵。我们可以看到公式(7)原本是由公式(3)（Lawler’s QAP）推导而来，经过上面的转化却可以写成类似Koopmans-Beckmann’s QAP的形式。当c=1c=1c=1时，公式(9)就退化为公式(4)也就是Koopmans-Beckmann’s QAP。而当c>1c>1c>1时，求解Lawler’s QAP可以看作是ccc个Koopmans-Beckmann’s QAP的组合优化问题。

路径跟随算法（Path-Following Algorithm）求解FGM

上文利用因式分解的方式得到一种新的统一的图匹配形式化方法，下面要介绍的就是如何求解上面定义的优化问题了。传统的求解图匹配问题通常分成两个步骤：1. 将GM进行连续性松弛，求解松弛后的问题得到一个近似的解；2.将近似解取整得到一个二元解。

解释一下，前文我们介绍过的双随机矩阵松弛就是一种连续性松弛的方式，他将XXX由一个只有{0,1}\{0,1\}{0,1}元素的离散矩阵，松弛为一个取值范围为[0,1][0,1][0,1]的连续矩阵，这样就可以用一些基于梯度的优化方法来求解了。但是得到的结果XXX也是连续的近似解，而我们需要的是节点和节点之间明确的对应关系，也就是说我们必须要将连续近似解再还原为一个只有{0,1}\{0,1\}{0,1}的二元解，且解要是一个置换矩阵（每行每列有且只有一个元素为1，其余为0）。

但这类方法有两个限制：第一，虽然经过了连续性松弛，但由于问题是非凸的因此容易陷入局部最优；第二，圆整的步骤不可避免地会导致精度损失。为了解决上述问题，作者提出一种基于路径跟随算法地迭代优化策略，这一策略有三个理论优势：1. 优化的效果与初始化条件无关；2.最终解确保可以收敛到一个整数解，因此无需额外的圆整步骤；3.迭代更新的步骤可以利用一些数值型连续求解方法的思想，这些方法已经广泛用于求解非线性系统的问题。

凸松弛和凹松弛

为了应用路径跟随算法，我们首先要找到Jgm(X)J_{gm}(X)Jgm(X)的一个凸松弛和凹松弛，为了实现这个要求，我们引入了一个辅助函数
Jcon(X)=∑i=1ctr(Ai1TXXTAi1)+tr(Ai2XTXAi2T)(10)J_{con}(X)=\sum_{i=1}^ctr({A_i^1}^TXX^TA_i^1)+tr({A_i^2}X^TX{A_i^2}^T)\quad (10)Jcon(X)=i=1∑ctr(Ai1TXXTAi1)+tr(Ai2XTXAi2T)(10)
由于置换矩阵XXX也是一个正交矩阵，因此XTX=IX^TX=IXTX=I所以Jcon(X)J_{con}(X)Jcon(X)是一个常数
Jcon(X)=γ=∑i=1ctr(Ai1TAi1)+tr(Ai2Ai2T)(11)J_{con}(X)=\gamma=\sum_{i=1}^ctr({A_i^1}^TA_i^1)+tr({A_i^2}{A_i^2}^T)\quad (11)Jcon(X)=γ=i=1∑ctr(Ai1TAi1)+tr(Ai2Ai2T)(11)
将12Jcon(X)\frac{1}{2}J_{con}(X)21Jcon(X)带入到Jgm(X)J_{gm}(X)Jgm(X)中可以构造得到两个等价的目标函数：
Jvex(X)=Jgm(X)−12Jcon(X)=tr(KPTX)−12∑i=1c∥XTAi1−Ai2XT∥F2(12)Jcav(X)=Jgm(X)+12Jcon(X)=tr(KPTX)+12∑i=1c∥XTAi1+Ai2XT∥F2(13)J_{vex}(X)=J_{gm}(X)-\frac{1}{2}J_{con}(X)\quad\quad\quad\quad\quad\quad\\\quad\quad\quad\quad\quad\quad\quad\quad=tr(K_P^TX)-\frac{1}{2}\sum_{i=1}^c\|X^TA_i^1-A_i^2X^T\|_F^2\quad (12)\\J_{cav}(X)=J_{gm}(X)+\frac{1}{2}J_{con}(X)\quad\quad\quad\quad\quad\quad\\\quad\quad\quad\quad\quad\quad\quad\quad=tr(K_P^TX)+\frac{1}{2}\sum_{i=1}^c\|X^TA_i^1+A_i^2X^T\|_F^2\quad (13)Jvex(X)=Jgm(X)−21Jcon(X)=tr(KPTX)−21i=1∑c∥XTAi1−Ai2XT∥F2(12)Jcav(X)=Jgm(X)+21Jcon(X)=tr(KPTX)+21i=1∑c∥XTAi1+Ai2XT∥F2(13)
对于包含任意置换矩阵和正交矩阵的Jgm(X)J_{gm}(X)Jgm(X)，优化过程中上述两个函数都将达到相同的结果，且差值为一个常数γ2\frac{\gamma}{2}2γ。由于Jvex(X)J_{vex}(X)Jvex(X)和Jcav(X)J_{cav}(X)Jcav(X)的Hessian矩阵（二阶导数矩阵）分别为：
▽X2Jvex(X)=−∑i=1c(I⊗Ai1−Ai2⊗I)T(I⊗Ai1−Ai2⊗I)(14)▽X2Jcav(X)=∑i=1c(I⊗Ai1+Ai2⊗I)T(I⊗Ai1+Ai2⊗I)(15)\bigtriangledown_X^2J_{vex}(X)=-\sum_{i=1}^c(I\otimes A_i^1 - A_i^2\otimes I)^T(I\otimes A_i^1 - A_i^2\otimes I)\quad (14)\\\bigtriangledown_X^2J_{cav}(X)=\sum_{i=1}^c(I\otimes A_i^1 + A_i^2\otimes I)^T(I\otimes A_i^1 + A_i^2\otimes I)\quad (15)▽X2Jvex(X)=−i=1∑c(I⊗Ai1−Ai2⊗I)T(I⊗Ai1−Ai2⊗I)(14)▽X2Jcav(X)=i=1∑c(I⊗Ai1+Ai2⊗I)T(I⊗Ai1+Ai2⊗I)(15)
分别是半负定和半正定的，因此最大化Jvex(X)J_{vex}(X)Jvex(X)和Jcav(X)J_{cav}(X)Jcav(X)分别是凸优化和凹优化问题。

路径跟随策略

作者通过迭代优化如下一系列的凸-凹问题来求解非凸的图匹配问题，
max⁡X∈DJα(X)=(1−α)Jvex(X)+αJcav(X)(16)\max_{X\in \mathcal{D}}J_{\alpha}(X)=(1-\alpha)J_{vex}(X)+\alpha J_{cav}(X)\quad (16)X∈DmaxJα(X)=(1−α)Jvex(X)+αJcav(X)(16)α∈[0,1]\alpha\in [0,1]α∈[0,1]是用于平衡凸松弛和凹松弛之间的平衡系数，在迭代过程中从0逐步增长到1。路径跟随策略的优势体现在三个方面：

当α=0\alpha =0α=0时该问题是一个纯粹的凸优化问题，数值型的优化器可以得到与初始化无关的全局最优解；
当α=1\alpha =1α=1时（迭代的终点），该问题是一个凹优化问题，局部的最优点通常在约束集合的极值点上。根据Birkhoff-von Neumann定理，所有的形成Birkhoff多面体的双随机矩阵，其顶点都是置换矩阵。也就是说局部最优点会收敛为置换矩阵，这也就保证了收敛的解为二元解而无需进行圆整操作；
α\alphaα平滑的从0增长到1的过程，与基于梯度的方法相比，路径跟随算法更有可能找到更好的局部最优解。

如上图所示，当α=0\alpha =0α=0和α=1\alpha =1α=1时，JαJ_{\alpha}Jα分别等于JvexJ_{vex}Jvex和JcavJ_{cav}Jcav。当α\alphaα从0到1的增长过程中，JvexJ_{vex}Jvex和JcavJ_{cav}Jcav分别是单调递减和单调递增的，表明JcavJ_{cav}Jcav部分对于JαJ_{\alpha}Jα的影响不断增强。当α=12\alpha =\frac{1}{2}α=21时，JαJ_{\alpha}Jα和JgmJ_{gm}Jgm的曲线相交。当α=1\alpha =1α=1时，JgmJ_{gm}Jgm与JvexJ_{vex}Jvex和JcavJ_{cav}Jcav两个松弛项之间的差值为一个常数γ2\frac{\gamma}{2}2γ。

子问题优化

前文中我们将原本的图匹配问题分别松弛为一个凸优化问题和凹优化问题，然后采用了路径跟随策略，通过逐步增大平衡系数，使得整个优化过程能够兼顾凸优化和凹优化的优势。接下来，需要解决的就是在特定的α\alphaα下，如何优化求解JαJ_{\alpha}Jα的问题。
本文采用的是Frank-Wolfe’s algorithm (FW)算法，首先给定一个初始的X0X_0X0，FW算法不断更新X∗=X0+ηDX^*=X_0+\eta DX∗=X0+ηD直至收敛。在每次迭代过程中，需要计算两个内容：优化的方向D∈DD\in \mathcal{D}D∈D和优化步长η∈[0,1]\eta\in [0,1]η∈[0,1]。通过求解以下线性规划问题，来计算DDD，
max⁡D∈Dtr(▽Jα(X0)T(D−X0))(17)\max_{D\in\mathcal{D}}tr(\bigtriangledown J_{\alpha}(X_0)^T(D-X_0))\quad (17)D∈Dmaxtr(▽Jα(X0)T(D−X0))(17)
上述目标函数是Jα(X)J_{\alpha}(X)Jα(X)在X0X_0X0处的一阶近似，▽Jα(X0)\bigtriangledown J_{\alpha}(X_0)▽Jα(X0)可以采用矩阵操作计算：
▽Jα(X)=▽Jgm(X)+(α−12)▽Jcon(X)(18)\bigtriangledown J_{\alpha}(X)=\bigtriangledown J_{gm}(X)+(\alpha-\frac{1}{2})\bigtriangledown J_{con}(X)\quad (18)▽Jα(X)=▽Jgm(X)+(α−21)▽Jcon(X)(18)▽Jgm(X)=Kp+H1(G1TXG2∘Kq)H2T+G1(H1TXH2∘Kq)G2T(19)\bigtriangledown J_{gm}(X)=K_p+H_1(G_1^TXG_2\circ K_q)H_2^T+G_1(H_1^TXH_2\circ K_q)G_2^T\quad (19)\\ ▽Jgm(X)=Kp+H1(G1TXG2∘Kq)H2T+G1(H1TXH2∘Kq)G2T(19)▽Jcon(X)=2(G1(H1TH1∘UUT)G1T+G2(H2TH2∘VVT)G2T)X(20)\bigtriangledown J_{con}(X)=2(G_1(H_1^TH_1\circ UU^T)G_1^T+G_2(H_2^TH_2\circ VV^T)G_2^T)X\quad (20)▽Jcon(X)=2(G1(H1TH1∘UUT)G1T+G2(H2TH2∘VVT)G2T)X(20)
上述线性规划问题（公式(17)）可采用Hungarian算法求解。得到梯度后，最优的η\etaη是下述抛物线的最大值点：Jα(X0+ηD)=aη2+bη+const(21)J_{\alpha}(X_0+\eta D)=a\eta^2+b\eta+const\quad (21)Jα(X0+ηD)=aη2+bη+const(21)
根据抛物线的特点，最大值点可能为η∗∈{0,1,−b2a}\eta^*\in\{0,1,-\frac{b}{2a}\}η∗∈{0,1,−2ab}。

总结

本文从分析关联矩阵KKK的满秩性和不定性出发，分析了现有方法存在的不足和问题。然后引出了一种新的基于因式分解的图匹配形式化方法，并提出一种统一的图匹配框架。最后，设计了一种基于路径跟随策略的图匹配优化方法，通过对目标函数进行凸松弛和凹松弛，并不断调整平衡系数α\alphaα，将优化目标由凸优化问题转移到凹优化问题。在每个迭代步骤中，都采用FW算法进行优化求解。FGM的主要内容在本文中已经做了详细的阐述，下篇文章将分析一下与图匹配问题与点配准（Point Registration）问题的区别和联系，进而介绍如何在FGM中引入额外的几何约束。

图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（一）
图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（二）
图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（三）

图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（二）相关推荐

图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（一）
这篇文章本身是图匹配经典论文<Factorized Graph Matching>的阅读笔记,后来发现该文介绍并串联了许多图匹配相关的知识,甚至可以看作一个小小的综述性文章,因此就作为图匹 ...
激光SLAM入门学习笔记
激光SLAM入门学习笔记激光SLAM入门学习笔记一.推荐阅读书籍二.推荐公众号.知乎.博客 1.公众号 2.知乎 3.博客三.推荐阅读论文&代码(参考泡泡机器人) 2D激光SLAM 3 ...
原创 OpenCV3编程入门学习笔记（总）
版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明. 本文链接:https://blog.csdn.net/qq_36163358/article/ ...
OpenCV3编程入门学习笔记（总）
OpenCV3编程入门学习笔记 2018.12.12-2018.12.29 此博客为在看过毛星云版<OpenCV3编程入门>后所总结的一本笔记,可供复习使用. 文章目录 OpenCV3编 ...
Altium Designer入门学习笔记4：PCB设计中各层的含义
Altium Designer入门学习笔记4:PCB设计中各层的含义阻焊层:solder mask,是指板子上要上绿油的部分:因为它是负片输出,所以实际上有solder mask的部分实际效果并不上 ...
机器学习入门学习笔记：（3.2）ID3决策树程序实现
前言之前的博客中介绍了决策树算法的原理并进行了数学推导(机器学习入门学习笔记:(3.1)决策树算法).决策树的原理相对简单,决策树算法有:ID3,C4.5,CART等算法.接下来将对ID3决策树算法 ...
node入门-学习笔记
文章目录 node入门-学习笔记 node 启动node服务重构express-run node入门-学习笔记 node 为什么后端要用node, 因为它是js js运行时环境主要使用expres ...
【带着canvas去流浪（11）】Three.js入门学习笔记
[摘要] three.js 入门学习笔记示例代码托管在:http://www.github.com/dashnowords/blogs 一. 资料推荐及建议 1.官方文档很详细,但是API部分单独 ...
RT-Thread Nano入门学习笔记（2）
RT-Thread系列 Keil模拟器 STM32F103上手指南学习笔记 RT-Thread Studio快速上手 RT-Thread Nano入门学习笔记(1) 项目实战:快速打造一个桌面mini ...

图匹配（Graph Matching）入门学习笔记——以《Factorized Graph Matching》为例（二）

目录