牛顿方法在高阶马尔可夫随机场推理中的应用

Newton-type Methods for Inference in Higher-Order Markov Random FieldsHariprasad KannanCentraleSupélec-INRIA SaclayUniversité Paris-Saclayhkannan@gmail.comNikos KomodakisEcole des Ponts ParisTechUniversité Paris Estnikos.komodakis@enpc.frNikos ParagiosCentraleSupélec-INRIA SaclayUniversité Paris-Saclaynikos.paragios@ecp.fr173000摘要0线性规划松弛在离散马尔可夫随机场中的 MAP推断中起着核心作用。正确求解拉格朗日对偶是这种方法的关键组成部分。在本文中，我们研究了使用牛顿类型方法来求解问题的平滑版本的拉格朗日对偶的好处。我们研究了它们实现优越的收敛行为和更好地处理问题的病态性质的能力，与一阶方法相比。我们表明，对于广泛范围的 MAP推断问题，确实可以高效地应用信任区域牛顿方法。在本文中，我们提出了一个可证明的全局高效框架，其中包括（i）在 Hessian矩阵构造方面在计算复杂性和精度之间的卓越折中，（ii）有助于高效优化的阻尼策略，（iii）与共轭梯度的通用预处理器相结合的截断策略，（iv）用于稀疏团势的高效求和计算。高阶马尔可夫随机场的结果展示了这种方法的潜力。01. 引言0许多计算机视觉问题可以使用马尔可夫随机场（MRF）进行建模。在 MRF中，最大后验（MAP）估计为节点分配标签，以最大化它们的联合概率分布。然而，对于一般图形，MAP 推断是NP-hard的，因此存在多种方法来获得近似解-图割、信念传播和基于线性规划松弛的方法[17]。近年来，高阶 MRF在各种应用中取得了出色的结果，因为它们模拟了节点之间的广泛相互作用。虽然对于成对 MRF的推理问题已经有了很好的研究，但是对于高阶模型的可扩展和高效技术的开发仍在不断发展[15]，[21]，[18]，[10]，[20]，[40]。线性规划松弛方法已经导致了最先进的算法，并且为 MAP 推断提供了理论基础。0MAP 推断的话题。这种方法的一个吸引人的特性是它很容易适用于高阶 MRF[21]，[39]的推理。在原始问题中求解线性规划松弛是不可扩展的，更好的方法是通过利用问题的图结构来求解对偶问题[37]，[41]。求解对偶的各种方法可以分为坐标优化和基于梯度的方法。块坐标方法收敛快，但在优化非平滑对偶时可能陷入次优解的角落[39]，[20]。由于在接近非平滑对偶的最优解时可以恢复更好的整数原始标签，这可能导致较差的解决方案。另一方面，超梯度方法在理论上可以达到全局最优解[22]。然而，它们的迭代次数为 O ( 10� 2  ) 收敛速度趋向于一个 �-准确解的速度。这些缺点可以通过用一个平滑版本来近似对偶来解决。平滑对偶和应用加速梯度技术是由[16]引入的，并且由[33]，[34]进一步研究。这些算法在 O ( 10�)。块坐标方法也可以处理平滑目标[26]，[14]，[28]，这使得这些算法可以避免次优解。其中一些方法[14]在平滑度减小时仍然可能陷入困境，以便获得更准确的结果。与基于梯度的方法相比，块坐标优化方法在并行化方面的空间更有限。乘子法的交替方向方法（ADMM）灵感方法[25]，[27]使用增广拉格朗日函数，这也是一个平滑目标。然而，这些方法的收敛速度分析尚未得到解决，并且在实践中也观察到这一点。例如，AD3[25]在许多中等规模的问题上效果很好，但在某些情况下可能无法收敛。此外，迄今为止，将这些方法扩展到具有高阶团的大型视觉问题尚未成功。近年来，牛顿类型方法在各种机器学习问题中取得了最先进的结果[36]，[24]，[23]。这些方法能够通过考虑曲率信息来选择更好的方向，并且具有θc(xc) +n�i=1θi(xi)(1)�xcθc(xc)φc(xc) +�i�xiθi(xi)φi(xi)φc(xc) = φi(xi),∀c, i ∈ c, xi�xiφi(xi) = 1,∀i;�xcφc(xc) = 1,∀cmax.δ�cmin.xc�θc(xc) −�i:i∈cδci(xi)�+�imin.xi�θi(xi) +�c:i∈cδci(xi)�(3)Here the dual variables are δci(xi), for each label xi of eachnode i within each clique c. If all the cliques have k nodeseach, with the nodes taking l labels, then the total numberof dual variables is |C|.k.l. We will denote this number asN and δ ∈ RN , is the vector of dual variables. We notethat the size of the dual is much lesser than the primal (2),hence, there has been considerable research effort towardssolving the MAP problem through the dual.In MAP inference, there is an interplay between the sizesof the graph and the sub-graphs within dual decomposi-tion [21] [42].For medium sized graphs, decomposinginto individual cliques, leads to excellent practical conver-gence. However, for large graphs, only bigger sub-graphslike chains of cliques lead to practical convergence [21],[29].We would like to emphasize that the formulationshown in (3) decomposes the graph according to the cliquesbut in general, the sub-graphs can be bigger regions.2.1. Smooth dualThe optimization of the non-smooth dual (3) has slowconvergence and can stall at a point away from the opti-73010当接近最优解时，具有二次收敛速度。在解决任何优化问题时，其中一个挑战是目标的条件。直观地说，如果目标值的变化由于变量值的扰动在方向上变化很大，则问题是病态的。人们发现，当条件数较小或中等时，一阶方法更快，但是二阶牛顿类型的方法对于病态问题的性能更好[11]。在MAP推断中，随着接近最优解，平滑性必须减小，导致一个相当病态的对偶目标。牛顿类型的方法的主要缺点是需要计算Hessian矩阵，这可能非常昂贵。然而，值得研究的是，对于手头的问题，是否确实如此。此外，使用 Hessian近似的拟牛顿方法也可以得到最先进的结果[36]，[5]。我们的贡献总结如下：(i)我们展示了对于基于平滑性的方法，可以高效地计算广泛类别问题的 Hessian 矩阵和 Hessian-向量乘积。(ii)我们研究了如何将牛顿类型的方法适应高阶 MRF 中的 MAP推断。(iii)我们展示了在稀疏模式的高阶团链中如何高效地执行求和-乘积推断。这极大地提高了基于平滑性的方法在高阶 MRF推断中的适用性。(iv)我们展示了牛顿类型的方法如何在高阶数据集上击败一阶方法。本文的其余部分组织如下：第2节概述了 MAP推断的概念和基于平滑性的方法。本文的主要贡献在第3节和第4节中介绍，实验验证在第5节中进行，最后进行讨论。02. 通过优化平滑对偶问题进行 MAP 推断0考虑一个图 G，其中 V 是节点集，表示随机变量，C是团集，强制它们包含的节点之间存在某种关系（例如平滑性）。每个节点从离散集合 l中取一个标签。例如，可以使用 MRF来进行目标检测，其中每个节点对应于对象的一部分，团表示部分之间的几何约束，标签表示图像中的位置。MAP推断是找到最大化这些随机变量的联合概率分布的标签的任务。它可以表示为等价的能量最小化问题，如下所示，0最小化 �0其中 θ i ( x i ) 是节点 i 的标签 x i 的势函数，θ c ( x c )是团 c 的标签 x c的势函数。具有两个以上节点的团称为高阶团。此外，这个能量最小化问题可以表示为整数线性规划（ILP），如下所示，0最小化 �0满足 �0φ i ( x i )  ∈ { 0 ,  1 } ,  � i, x i ; φ c ( x c )  ∈ { 0 ,  1 }� c,  x c. (2) 这里，φ i ( x i ) 和 φ c ( x c ) 是给定节点 i (或者团 c)的指示向量，对于标签 x i (或者 x c)，它们是离散优化变量。约束条件表示一个多面体，称为边际多面体。放松整数约束（最后一行），得到线性规划松弛问题，该问题定义在局部多面体上。通过对偶分解[22]，[38]可以得到该线性规划松弛问题的拉格朗日对偶。其基本思想是将图分解为易处理的子图，以推导出拉格朗日函数。这导致了一个凹、无约束和非光滑的优化问题。特别地，如果将每个团视为易处理的子图，则[38]中给出的方法会导致以下对偶问题。mum. In order to reach closer to the optimum of the non-smooth dual, optimization over smoother versions are car-ried out. The smooth dual can be obtained by adding tothe primal objective (2), entropies corresponding to all thenodes and cliques. The dual of this modified problem canbe derived based on duality theory (refer [28], [16] for fur-ther details). The smooth dual looks very similar to thenon-smooth one, where each minimum operation in (3) isreplaced by the negative soft-max function. The smoothoptimization problem takes the following form,max.δ�csminxc�θc(xc) −�i:i∈cδci(xi); τ�+�isminxi�θi(xi) +�c:i∈cδci(xi); τ�(4)We will denote the smooth dual function as g(δ). Thenegative soft-max function is defined as sminx (f(x); τ) =−1τ log �x exp(−τf(x)). As τ increases, (4) is a closerapproximation of (3) and we get closer to the optimum ofthe non-smooth dual. However, a smooth approximationbecomes increasingly ill-conditioned as it approaches theshape of the non-smooth problem (§7, [31]). Hence, it iscostlier to optimize for larger values of τ. It is better to startwith a very smooth version and switch to less smoothness,with warm start from the previous smoother version [34].3. Trust-Region Newton for MAP InferenceThe central concept in Newton-type methods is the useof curvature information to compute the search direction.In each iteration, a local quadratic approximation is con-structed and a step towards the minimum of this quadraticis taken. For the smooth dual g(δ), the equations to obtainthe Newton direction are shown below. In these equations,B(δ) ∈ RN XN carries the curvature information and canbe either the exact Hessian (∇2g(δ)) or a modified Hessianor a Hessian approximation (e.g., quasi-Newton).g(δ + p) ≈ g(δ) + ∇g(δ)T p + pT B(δ)p(5)q(p) = ∇g(δ)T p + pT B(δ)p(6)p∗ = argminp∇g(δ)T p + pT B(δ)p(7)∴ B(δ)p∗ = −∇g(δ)(8)Here, the minimum p∗ of the quadratic approximation q(p)will be a descent direction for positive definite B(δ) and aline search along this direction determines the step size. Forunconstrained problems, the Newton direction is found bysolving the linear system (8). Since, N is large for com-puter vision problems, we use Conjugate Gradients (CG) toobtain the Newton direction in our work.73020图1.Hessian的两个组件。在每个组件中，相同颜色的块具有相同的值。在组件一中，有与团相同的唯一块。在组件二中，每个块行/列都有相同的块。0对于牛顿类型的方法，如果B(δ)是精确的Hessian(�2g(δ))，并且如果回溯线搜索首先测试全步长，那么当接近最优解时，可以实现二次收敛[4]。与一阶方法相比，这是可取的。如第2.1节所述，随着算法的进行，平滑度降低，导致病态。牛顿类型的方法对病态具有一定的鲁棒性，因为它们具有仿射不变性，即对于某些函数f(x)和¯f(¯x)=f(Ax)，其中A是可逆方阵，牛顿类型的方法的迭代将相关联，即¯pk=Apk。因此，从理论上讲，这些方法不受病态的影响，并且在有限精度下，牛顿类型的方法表现出鲁棒行为的条件数范围比一阶方法更大。值得研究牛顿类型的方法如何处理MAP推断。03.1. Hessian相关计算0尽管牛顿类型的方法具有优势，但填充和解线性系统(8)可能是昂贵的。共轭梯度(CG)算法的计算量较大的步骤是：计算Hessian-向量乘积和定期构建和应用预处理矩阵。通常情况下，如果Hessian难以填充，复杂步长微分(CSD)[1]可以提供非常准确的Hessian-向量乘积，每个乘积的成本大约是一个梯度计算，这在CG中使用可能是昂贵的。因此，我们研究了高效填充Hessian和获得快速Hessian-向量乘积的可能性。如第2节所述，对于中等规模的问题，根据团进行分解并实现实际收敛是足够的。我们现在将展示，对于根据像单个团这样的小子图进行分解，Hessian可以非常高效地计算。实际上，计算梯度和Hessian只需要大约两倍的梯度计算时间。如果我们仔细观察Hessian，我们会发现它可以被写成两个组件的和，如下所示µci(xi) =�xc:xc(i)=xiexp�τ.(θc(xc) − �n:n∈cδcn(xi))��xcexp�τ.(θc(xc) − �n:n∈cδcn(xn))�(9)µi(xi) =exp�τ.(θi(xi) + �k:i∈kδki(xi))��xlexp�τ.(θi(xl) + �k:i∈kδki(xl))�(10)(12)73030图1. 两个组件都具有块结构。考虑以下数量，0组件一的元素可以写为Hc,ij(x i,x j)=τ�µcij(x i,x j)−µci(xi)µcj(x j)�，其中µcij(x i,x j)可以像µci(xi)一样通过固定两个节点的标签来计算。这导致了一个块对角矩阵，有与团一样多的唯一对称块。组件二的元素可以写为Hst,i(x i,x i)=τµi(x i)(1−µi(x i))和Hst,i(x i,x j)=−τµi(xi)µi(xj)。它有与节点一样多的唯一对称块，并且给定的行或列有相同块的重复副本。这里c、s、t是团，i、j是成员节点，xi、xj是节点标签，xc是团标签。而且，正是组件二包含了由于重叠团而引起的非对角块。因此，Hessian的元素可以通过仅迭代一次遍历所有团和节点来计算。避免了与重叠团对应的迭代。对于梯度，需要通过一次遍历所有团和节点来计算值的数组。对于Hessian，我们需要在每个团和节点上计算（对称）块的值。实际上，由于缓存重用，这只需要一个梯度计算的开销时间。这些块可以方便地用于并行化Hessian-向量乘法，我们使用简单的OpenMP代码来实现。因此，高效的Hessian-向量乘积导致了我们的情况下高效的CG例程。对于一些问题，Hessian是稀疏的，因为一个团只与少数其他团重叠。然而，由于特殊的结构，没有必要利用这种稀疏性来计算Hessian的唯一块。此外，如果有许多重叠的团，计算Hessian-向量乘积也不会受到不利影响，因为每行组件二矩阵中只会有更多相同块的副本，对应于共享节点。因此，计算机内存中的数据移动受到限制。此外，这些数据结构可以方便地用于构建预处理器。03.2. 阻尼矩阵方法0光滑的对偶问题（4）只是严格凸的，而不是强凸的，即在最优点附近，Hessian矩阵只是半正定的。在这种情况下，需要采用信任区域方法来进行有意义的牛顿步骤。在这里，0使用相同的二次近似（6）来最小化约束∥p∥ ≤∆，其中∆是信任半径。在开发用于MRF推断的信任区域牛顿方法（TRN-MRF）时，我们解决了几个问题：如何强制信任区域，如何通过处理病态问题来提高收敛速度，以及如何设计合适的预处理器。正是这些选择的组合导致了一个可用的算法，并且我们注意到，对于MAP推断有效的方法可能对其他任务无效，反之亦然。Steihaug方法似乎是尝试解决大型问题的第一种方法。它具有将信任区域根据景观塑造成椭球的良好特性。这是通过最小化（6）并在约束∥p∥M ≤∆的情况下实现的，其中M是一个预处理器，∥p∥M =pTMp[7]。然而，在缺乏良好的预处理器的情况下，在给定的外部牛顿迭代中，算法在计算出一个好的方向之前很快达到信任半径，导致多次外部迭代。Levenberg-Marquardt算法通常用于最小二乘问题，它提供了另一种实施信任区域的方法。我们借用的思想是阻尼矩阵，它是一个正则化器，用于解决严格凸性和病态性。该矩阵被添加到Hessian矩阵中，以获得修改后的Hessian矩阵[8]。阻尼矩阵将CG返回的牛顿方向限制在信任区域内。Levenberg-Marquardt算法使用缩放后的Hessian对角线作为阻尼矩阵。这个选择在MAP推断中效果很差。相反，我们对Hessian进行如下修改：B(δ)= �2g(δ) + λI，其中λ >0。虽然Steihaug明确地使用信任半径∆工作，但我们通过λ隐式地施加它。这可以通过以下方程看出，这些方程将阻尼参数λ与信任半径∆联系起来。0(�2g(δ) + λI)p� = -�g(δ)0λ(∆ - ∥p�∥) = 0 (11)0在一个足够正定的区域（接近最优点），λ接近于零，牛顿方向p�是用真实的Hessian计算的。如果λ > 0，则∥p�∥ =∆，即方向受到信任区域的限制。在病态区域，λ将会很大，强制施加的信任半径将会很小。因此，每次迭代后，λ都会适应，并且可以按照以下方式进行：0ρ < 0.25: λ ← 2λ; 0.25 < ρ < 0.5: λ ← λ; 0.5 < ρ< 0.9: λ ← 0.5λ; 0.9 < ρ: λ ← 0.25λ0其中，ρ = g(δ + p) - g(δ)0q(p) - q(0)0ρ表示方程（6）的二次近似对双对偶问题（g(δ)）的逼近程度。随着算法接近最优点，λ变得非常小，算法在没有任何扰动的情况下达到真正的最优点。Having found out the suitability of damping matrix basedapproach, it is still necessary to properly address the ill-conditioning caused by annealing.Trust-region Newtonproceeds by taking approximate Newton steps, where ineach outer iteration the run of CG iterations is truncated bya suitable criterion. However, as the algorithm approachesoptimality, it is critical to solve the linear systems to greateraccuracy and get better Newton steps [35]. Otherwise, thealgorithm will take too long to converge or will not con-verge at all. Generally, at an outer iteration k, CG can betruncated at iteration j according to the following condition,∥rj∥ ≤ ηk∥∇g(δk)∥. Here, rj = B(δk)δjk + ∇g(δk), isthe residual of equation 8, at iteration j of CG and ηk isreferred to as the forcing sequence. Through ηk we reducethe residual and achieve more accurate Newton direction.For MAP inference, we found textbook choices of the forc-ing sequence leading to Newton iterations not convergingbecause the residual never becomes low enough to achievethe more accurate directions required for further progress.Hence, for TRN-MRF, the following criterion has been de-730403.3. CG截断的强制序列0∥�g(δk)∥）。这个条件自然地对后续迭代的残差施加了更强的条件，而且，随着退火的进行，项（�τ）的值变小，确保获得足够准确的牛顿步长，因为平滑减少了。03.4. 基于团的预处理器和回溯搜索0为了进一步提高TRN-MRF的效率，我们将解决一个影响信任区域牛顿方法的重要方面。这涉及共轭梯度例程的计算效率。CG迭代的收敛性取决于线性系统的不同特征值簇的数量，而良好预处理的系统（M-1Ax =M-1b）将具有较少的簇。矩阵M应尽可能与A相似，并且应构造和求逆高效。标准的预处理方法，如不完全Cholesky分解、拟牛顿法和多重网格法，无法处理通用的MAP推断问题。基于团的块对角预处理器在对角线上具有特定团（Hc）的双导数项，即Hc i,xi,j,xj =∂2g(δ)0∂δ ci ( x i ) ∂δ cj ( x j)，其中i，j是属于团的节点，xi，xj是它们的标签。由于其结构与MAP推理问题密切相关，因此它的性能相当好。计算成本是分别计算和求逆这些块。这是在计算梯度和Hessian矩阵的同时进行的。应用这个预处理器对应于涉及块逆的矩阵-向量乘法。对于足够大的λ值，修改后的Hessian矩阵为0自动良好条件和CG将快速收敛。当λ减小到接近于零的值时，CG运行的迭代次数会增加。我们观察到使用这个预处理器的CG性能有了显著的改进，本文中的结果都是基于此的。此外，接近最优解时，会有CG运行达到最大允许迭代次数的情况。我们将其设置为250，适用于所有实验。最后一个方面涉及到在CG例程计算出牛顿方向后的回溯搜索。如果方程（12）中的ρ小于一个小值（比如ερ），这意味着函数值要么减少得很少，要么增加。因此，我们不能直接沿着这个方向迈出一步。然而，在TRN-MRF的每个外部迭代中，希望能够迈出足够减少的一步，因此在这些情况下我们进行回溯搜索[30]。回溯可以沿着一条直线或沿着一个曲线路径（CG迭代的一个子集）进行。尽管在其他问题中，沿着曲线路径进行回溯可以得到良好的结果[24]，但我们观察到在MAP推理中效果不佳：最终方向非常接近最陡下降方向。另一方面，沿着CG计算出的方向进行回溯线搜索，可以大大加快TRN-MRF的速度。我们实现了一个基于三次插值的搜索，非常高效。03.5. 退火计划和停止条件0[34]建议通过定期计算平滑问题的原始对偶（PD）间隙来退火τ。对于中间对偶变量，他们通过为每个团解决一个小型线性规划（称为运输问题）来恢复可行的原始变量。由于他们的方法适用于成对图，他们在少于一千次的预测调用（预测调用是迭代或计算PD间隙）中获得结果。然而，对于高阶MRF，计算可行的原始变量和原始目标是昂贵的。因此，每隔几次迭代计算平滑问题的PD间隙会极大地影响计算效率。我们使用了一个简单但直观的判断标准来判断何时进行退火。由于我们正在处理一个凹函数，梯度欧几里得范数较小的区域保证比梯度欧几里得范数较大的区域更接近最优解。因此，如果�g(δk)表示在给定τ下运行k次迭代后的梯度，我们更新τ←ατ，其中1<α，如果∥�g(δk)∥2<γτ。如果我们设置一个足够强的阈值γτ，我们可以保证对于特定的τ获得足够的改进，并且可以执行退火。我们定义γτ=β∥�g(δ)∥，就在τ退火之后。与[34]中的证明类似，这种退火方法适用于任何对于固定τ值收敛到全局最优解的优化算法。我们测试的所有算法，在平滑、凹问题上都有这个保证。此外，我们确保τ已经达到了一个large enough value τmax in order to obtain accurate results.In order to exit the least smooth problem, [34] use thenon-smooth PD gap. We have observed that TRN-MRF, dueto its quadratic convergence rate, can exit based on classi-cal gradient based condition itself. More precisely, with thel∞ norm and a threshold of ζ = 10−3 (§8, [12]), TRN-MRFachieves good exit behaviour. However, first order meth-ods take too long a time to achieve this gradient based exitcondition and many times never do so. Hence, we have im-plemented a PD gap based approach, so that all algorithmscan exit gracefully. The approach in [34], is available onlyfor pairwise graphs in the openGM library and implement-ing small LP solvers for all the higher order cliques lookschallenging. [28] proposed a method involving only closedform calculations and we have implemented their approach.Our complete trust-region Newton method, within an an-nealing framework, is described in Algorithm 1.In ourexperiments, we set, λ0 = 1, α = 2, β =16, ǫρ =10−4, ζ = 10−3, τmax = 213 and ǫτ = 0.1 if τ <τmax4, 0.01 if τmax4< τ < τmax2, 0.001 if τ7:set ηk = min( ǫτk ,�∥∇g(δk)∥)8:Run CG while ∥rj∥ > ηk∥∇g(δk)∥9:obtain Newton direction p and calculate ρ10:update λ according to equation (12)11:if ρ < ǫρ then backtracking line search along p toobtain Newton step pk12:δk+1 = δk + pk13: end while4. Scalable Smoothing based approachEven though smoothing based approach has been scaledfor large pairwise graphs [34], higher order MRFs with largelabel spaces and/or large graphs are less studied. In thissection we show an efficient way to compute the smoothingoperation with large label spaces for a very useful class ofclique potentials. Next we demonstrate the use of quasi-Newton methods for problems with large graphs.4.1. Pattern-based SmoothingIn the smooth dual (4), we denote the first term as g1(δ).It corresponds to cliques and involves log-sum-exp calcula-tions over all possible labellings for each clique. This scalesexponentially (lk) with clique size k, where each node takesl labels.Hence, computing the gradient of this dual iscomputationally heavy, especially for higher order cliques.[16] observed that these terms in the gradient correspond tomarginal probabilities in a suitably defined graphical model.Hence, they can be computed using the sum-product algo-rithm [19]. Still, O(lk) complexity remains.Sparse, pattern-based clique potentials have been veryuseful in computer vision [21], [32]. In these potentials, abig majority of the labellings take a constant (usually high)value and a small subset take other significant values. Oncloser observation many clique potentials fit this descrip-tion. [21] showed an efficient way to perform max-productcomputation in chains of such cliques. However, it is notclear how to extend their work to sum-product computation.[8] showed a sum-product approach in pairwise MRFs withsparse, pattern-based potentials. They achieved a complex-ity of O(2l + s), instead of O(l2). They mention that theiridea can be used for higher order cliques with a factor graphrepresentation but don’t go into details.We have foundthat with a factor graph representation, their approach hasa complexity of O(lk−1 + l + s). Instead of a factor graph,if a cliq