基于嵌套剖分的位姿圖分層優化算法

摘 要：位姿圖優化 （pose graph optimization，PGO）是一種在同時定位與地圖構建（simultaneous localization and mapping，SLAM）后端優化中常用的高維非凸優化算法，通常建模成極大似然估計。由于目前的PGO算法優化大規模大噪聲數據集時很難在保證精度的同時提升速度，所以提出了一種基于嵌套剖分的位姿圖分層優化算法。該算法首先建立不同距離度量的χ2檢驗模型，進而剔除異常值點。然后利用嵌套剖分算法將位姿圖分割成一組子圖，再從這些子圖中提取出一個表示原SLAM問題的抽象拓撲的骨架圖，從而優化該骨架圖，完成初始化。最后在模擬和真實的位姿圖數據集上進行實驗評估，結果表明該算法在不影響精度的情況下，可以提高算法的計算速度，具有可伸縮性。

關鍵詞：位姿圖優化； 嵌套剖分； 噪聲； χ2檢驗; 最大似然估計

中圖分類號：TP391 文獻標志碼：A

文章編號：1001-3695（2024）06-046-1916-05

doi：10.19734/j.issn.1001-3695.2023.08.0386

Hierarchical pose graph optimization algorithm based on nested dissection

Abstract：PGO is a high-dimensional non-convex optimization algorithm commonly used in the back-end optimization of SLAM， which is usually modeled as maximum likelihood estimation. Since the current PGO algorithm faces difficulties in improving speed while ensuring accuracy when optimizing large-scale noise datasets， this paper proposed a hierarchical pose graph optimization algorithm based on nested dissection algorithm for large noise datasets. The algorithm firstly established a χ2 test model with different distance measures， and then removed outlier points. Secondly， it used nested dissection method to split the original pose graph into a set of subgraphs and extracted a skeleton graph from these subgraphs. The skeleton represented the abstract topology of the original SLAM problem. Then， the algorithm optimized the skeleton graph and completed the initialization. Finally， experimental evaluation on simulated and real pose graph datasets show that the proposed algorithm can improve the calculation speed and scalability of the algorithm without affecting the accuracy.

Key words：pose graph optimization; nested dissection; noise; χ2 test; maximum-likelihood estimation

0 引言

在缺乏環境先驗信息的情況下，SLAM使機器人能依據自身搭載的傳感器實現自主定位和導航［1］，其本質是從相對測量中估計機器人位姿。位姿圖是該問題的圖形表示，其中節點表示位姿，通常用特殊歐氏群SE（3）上的元素表示；邊表示節點對之間的空間約束，即相對測量［2］，通常用概率的方式來表示由外部觀測引起的不確定性。由于位姿圖優化 （pose graph optimization，PGO） 具有非凸性，所以已有研究［3～5］將其建模成一個非凸極大似然估計問題。傳統的求解方法是將該問題轉為非線性最小二乘問題，通過高斯牛頓［6］、列文伯格-馬夸爾特［7，8］等迭代優化算法來獲得與概率約束最大一致的節點估計［9～11］。在實際應用中，隨著問題規模的擴大和復雜性的提升，PGO算法的研究開始向分層策略遷移。Bosse等人［12］首次將分治法［13］應用于SLAM領域，并解決了城市尺度區域的空中三角測量和測繪問題。文獻［14］提出一種用于大規模環境的精確局部建圖方法，通過緩存分解的子圖來提升運行的速度。之后，Ni等人［15］將嵌套剖分算法運用到原始的SLAM圖上，將其遞歸地劃分為多個層次的子圖，極大地減少了兩棵子樹之間的依賴，并沿相應的聚類樹采用自底向上的推理對原始的SLAM圖優化。同時，Grisetti等人［16］提出了一種層級式的局部圖，在得到一個新的觀測值時，位于高層級的局部圖會被立刻修改，而在低層級的局部圖中只有受影響的區域才會被更改。Tang等人［17］提出了一種高效的分層位姿圖優化算法，該方法利用一種修改后的標準切割算法自動將位姿圖分割成在關鍵節點條件下相互獨立的子圖。然后通過位姿圖稀疏算法將關鍵節點生成一個上層圖，并在優化上層圖后再獨立解決每個子圖的優化問題。Tian等人［18］提出一種基于段的分層優化方法，該方法提出了一種可靠的軌跡分割方法以提高優化的效率，并首次提出一種緩沖機制來提高分割的魯棒性。 Korovko等人［19］將部分優化與分層位姿相結合，提出一種構建層次結構的方法，允許快速添加新節點進入層次結構，無須重建，極大地提升了運行速度。Guadagnino等人［20］提出 HiPE 算法，該算法通過結合來自輸入局部區域的最大似然估計，將問題進行幾何圖形的抽象表示，再將其編碼為粗粒度圖。相較于Cauchy 算法［21］，該算法通過利用這種表示的稀疏性，能更加穩健地初始化姿態圖。

在實際應用中，光照、天氣、溫度變化等環境因素會對傳感器性能造成影響，使得傳感器獲取的測量數據具有較大的噪聲，從而導致HiPE算法迭代次數顯著增加，魯棒性較低，極大地增加了陷入局部極小值無法收斂到全局最優解的可能。因此，本文針對大噪聲數據集提出了一種基于嵌套解剖的位姿圖分層優化方法ND-PGO。其貢獻如下：a）提出了一種基于不同距離度量的χ2檢驗模型，對大噪聲數據集進行處理，檢測剔除大噪聲數據集中的異常值點；b）利用嵌套剖分算法將數據處理后的位姿圖分割成一組子圖，再提取一個骨架圖來抽象表示待解決的SLAM問題，然后通過對該圖骨架非線性優化來初始化位姿變量，最后，將優化后的圖骨架位姿變量放入原始位姿圖中固定，對剩余變量進行優化；c）在PGO數據集Sphere5000、Parking-garage、torus3D和grid3D上進行對比實驗。結果表明，相較于其他先進的位姿圖算法，ND-PGO算法在不犧牲精度的情況下，顯著提高了優化運行速度。

1 PGO算法介紹

1.1 PGO模型

位姿圖通常用一個有向圖G={V，E}表示（圖1）。其中V是非空集合，稱為節點集；E是V中元素構成的無序二元組集合，稱為邊集。相鄰節點之間的邊為里程約束（圖1實線），其余邊為閉環約束［22］ （圖1虛線）。位姿節點i，i∈V可用SE（3）中的元素表示：

PGO問題的目標是在位姿圖所有相對測量的約束下找到絕對位姿最小二乘誤差的節點構型。采用極大似然估計（maximum likelihood estimation，MLE）方法，即求解：

假設測量值是獨立同分布的，則：

1.1.1 相對測量分布模型

在MLE中最常見的選擇是假設測量值服從高斯分布。然而，通過高斯分布獲得的測量值在求代價函數的過程中引入了對數映射，從而導致該優化問題具有很高的非凸性。為了克服這一問題，Rosen等人［23］提出了SE（3）上相對測量分布的解耦：

1.1.2 極大似然估計

式（6）在負對數下有較好的數學形式，將其取負對數得：

求解位姿的MLE等價于最小化負對數似然，則式（4）變為

其中：‖·‖F為Frobenius范數。將式（7）代入式（8）中，得到PGO問題的解為

1.2 PGO初始化

式（9）通常采用迭代的方法求解，但是迭代的方法非常依賴初值。一個良好的初始化結果不僅可以大大降低收斂到次優配置的風險，還可以加快迭代的收斂速度，減少優化的計算時間。因此求解式（9）非常關鍵的一步是尋找盡可能接近最優值的初始值。然而PGO初始值是通過前端里程計測量獲得的，精度較低，迭代求解時容易陷入局部最小值，甚至難以收斂。因此研究者們提出了一些經典的PGO初始化方法，其中最廣泛使用的算法就是弦初始化算法［24］。該方法首先通過弦松弛估計旋轉矩陣，然后利用這個旋轉估計來計算平移向量的估計量。最后將這個位姿的粗略估計用作迭代優化算法的初始值。弦初始化算法將式（9）的旋轉部分替換為［25］

定義rki為RTi（k=1，2，3）的第k列，則式（10）為

放松SO（3）約束，則式（11）轉換為一個無約束問題：

求解式（12）可得到n個矩陣Ai=［r1i，r2i，r3i］。這些矩陣不一定是旋轉矩陣，但可以通過Ai找到滿足旋轉約束的矩陣Ri：

式（13）可通過SVD分解算法［26］得到一個封閉的解。對矩陣Ai進行奇異值分解Ai=USVT，則有

R*i=U diag（［1 1 det（UVT）］）VT（14）

已知旋轉矩陣，可以通過線性最小二乘估計平移向量，即求解：

該算法在實踐中表現良好。此外，由于初始化是由兩個線性問題組成的，所以計算會更高效。

2 本文方法

2.1 基于不同距離度量的χ2檢驗模型

a）測地線距離（又稱角距離）：表示相對旋轉RTaRb對應的旋轉角度，公式如下：

dang （Ra，Rb）=‖Log（RTaRb）‖=‖Log（RTbRa）‖（17）

其中：Log（R） 表示SO（3）在恒等式處的對數映射。

b）弦距離：Ra-Rb的Frobenius范數：

dchord（Ra，Rb）=‖Ra-Rb‖F=‖RTaRb-I3‖F（18）

c）四元數距離：

dquat（Ra，Rb）=min（‖qa－qb‖，‖qa+qb‖）（19）

其中：qa、qb是旋轉矩陣Ra、Rb的四元數表示。因為qa與-qa表示同一旋轉，所以式（19）中的最小運算符用于解決符號歧義。

分別將式（17）～（19）引入式（16）可得基于不同距離度量的χ2檢驗模型如下：

假設里程計中無異常值，則不含異常值的回環邊須滿足以下條件［27］：

rij<χ2α，dij（21）

其中：dij是i、j之間的約束自由度；α為χ2檢驗顯著性水平。

2.2 基于嵌套剖分的骨架圖提取算法

針對式（9），一般求解其泰勒展開的線性估計的最小值。在線性代數［28］和圖論［29］文獻中，已經證明了線性系統求解的效率在很大程度上取決于消去順序。將同樣的思想應用于SLAM問題［11］，良好的消去排序在圖三角化過程中產生小團，并在矩陣分解過程中引入較少的非零填充。雖然找到一個最優排序是NP完全問題［30］，但是找到一個有效的排序有最小度［31，32］和嵌套剖分［33，34］兩種方法，其中最小度方法又包括多重最小度［31］（multiple minimum degree，MMD）、近似最小度［32］ （approximate minimum degree， AMD）等變種。針對基于子映射的方法，嵌套剖分比最小度剖分的性質更好。雖然這兩種方法生成的消除順序質量相當，但嵌套剖分傾向于在大型圖上執行得更好，并且它的消除樹通常更低、更平衡，因此更適合本文的分層劃分。

本文利用嵌套剖分算法遞歸地將輸入劃分為代表問題空間連通區域的分區。劃分過程旨在最小化分區之間的連接，并使得每個分區具有相對均衡的負載。嵌套剖分算法首先初始化圖結構和變量索引，再根據變量索引構造近似的Hessian矩陣模式。通常情況下，近似的Hessian矩陣模式是一個稀疏矩陣，用于表示變量之間的相關性。使用近似的Hessian矩陣模式，對問題的圖進行嵌套分解。嵌套分解是一種遞歸分割圖的過程，類似于樹的分解。在每一次遞歸中，選擇一個分割點，將圖分成兩個子圖。繼續對每個子圖進行遞歸分解，直到子圖足夠小，不再進行分解。對于這些分區，使用局部圖中可用的測量值提取一個由變量子集組成的骨架圖，該骨架圖必須在保持問題空間結構的同時保持其稀疏性。圖2展示了骨架圖的優化效果，其中（a）為原始圖，（b）為利用嵌套剖分算法分區后提取出的骨架圖，（c）為利用弦初始化算法優化后的骨架圖。

2.3 ND-PGO算法流程

對于含有大噪聲的數據，ND-PGO方法使用2.1節中提出的χ2檢驗模型篩選和剔除相對觀測值中的異常值。然后利用嵌套剖分算法遞歸地劃分分區，形成分區集后，使用局部估計來形成骨架圖，然后對輸入分兩個階段進行初始化。

算法1 ND-PGO算法

形成骨架圖首先需要對每個分區進行優化，獲得其變量的局部估計（算法1第8行）。由于形成原問題抽象表示的骨架圖僅可使用相對度量，所以需要從每個分區的變量中選擇一個錨點變量作為后續優化過程中固定變量的值。本文選擇度數最高（連接的因子數量最多）的變量作為錨點變量（算法1第6行）。此外，本算法在錨點變量和每個邊界變量之間添加虛擬測量，從而得到局部圖的稀疏近似（算法1第10行）。整體骨架將由一組錨點變量和邊界變量χS以及相應的虛擬測量值Euclid Math OneZApS組成。

分階段進行初始化時，首先通過極大似然估計計算骨架圖中變量的初始猜測值，合并塊中的局部信息。由于邊界變量在不同分區中的估計不一致，所以需要對骨架的一致估計對齊塊（算法1第19行）。然后再固定得到的骨架變量的初始值，優化得到剩余變量χR=χ＼χS的初始猜測（算法1第21行）。其中χS和χR均是通過1.2節中的弦初始化方法來引導收斂的（算法第18、20行）。

3 實驗

本文采用PGO公開數據集，首先在低噪聲環境下的公共數據集上驗證算法的性能，接著向這些數據集中添加噪聲。本文選取g2o［10］和HiPE［20］算法作為實驗對比方案。其中g2o算法是最廣泛使用的SLAM框架后端優化算法，具有很高的魯棒性；HiPE算法是最新的分層算法，在運行速度方面有很大提升。因此將提出的ND-PGO與g2o、HiPE算法在大噪聲數據集上進行實驗對比，評估ND-PGO算法在大噪聲數據集下的性能。本文使用歸一化χ2度量作為實驗的評估標準。

其中：m為相對測量邊數；n為位姿點數。一般將χ2度量數值作為評價標準，χ2度量數值越低，機器人軌跡優化結果越精確，算法越魯棒。實驗的硬件條件為Intel Xeon CPU E5-1620 v3 @ 3.50 GHz×8、NVIDIA Corporation GP104，軟件環境為Ubuntu 20.04，編程語言為C++。

3.1 數據集

在3DSLAM中，PGO公共數據集主要有Sphere5000、Parking-garage、torus3D和grid3D。其中Sphere5000、torus3D和grid3D數據集均為模擬生成的數據集，分別模擬了機器人在球面、環面表面和三維網格世界上運動時的軌跡。Parking-garage是Vertigo收集的用于研究自動停車的斯坦福停車場的3D地圖，為真實場景數據集［35］。本文將含有大量錯誤和異常的數據稱之為大噪聲數據集，包括相對觀測值中含有較大的隨機誤差以及偏離正常數據較遠的異常值點。

表1記錄了每個數據集的具體參數，θnoise表示相對旋轉測量的噪聲（單位：弧度（rad））。mean、std分別為其均值與標準差，n為位姿節點數量，m為相對測量邊的數量。

為評估算法的抗噪能力，在數據集的相對方向測量中分別加入均值μR=0，標準差σR為0.1 rad、0.2 rad、0.3 rad、0.4 rad、0.5 rad的噪聲，以比較算法在大噪聲數據集下的魯棒性。

3.2 實驗結果

在PGO數據集中加入不同水平的噪聲，重復實驗50次取平均值。此次實驗采用模型中的rij（χ）ang、rij（χ）chord、rij（χ）quat作為χ2檢驗函數的ND-PGO算法，分別記為ND-PGO+a、ND-PGO+c和ND-PGO+q。圖3記錄了Parking-garage、torus3D、Sphere5000和grid3D數據集在無噪聲和有噪聲的情況下，ND-PGO+a、ND-PGO+c、ND-PGO+q、HiPE和g2o算法對數據集優化所用的時間。從圖3可得，ND-PGO算法明顯快于g2o算法，相較于HiPE算法，運行速度也有顯著優勢。在基于ND-PGO的三種算法中，ND-PGO+c稍微快于其他兩種算法。

表2～5分別記錄了各個算法對四個數據集優化得到的χ2度量值和迭代次數（表中加黑數據說明對應的算法最優）。從表2～5中的數據可知，分層優化算法會比傳統算法的效率更高，其中在中小噪聲集中，HiPE算法在運行精度上會比本文算法更有優勢。但隨著噪聲增大，HiPE算法需花費時間考慮異常值的影響，計算速度受到影響，ND-PGO的運行時間明顯少于HiPE算法。此外，基于ND-PGO的三種算法的χ2度量值基本低于HiPE算法和g2o算法。其中ND-PGO+c算法在Parking-garage、torus3D和grid3D上的表現稍優于其他兩種基于ND-PGO的算法。ND-PGO+q在Sphere5000數據集上表現最佳。ND-PGO算法首先采用χ2檢驗來剔除異常值，這在一定程度上可以幫助算法在大噪聲情況下避免局部極小值，保證算法精度的同時提升后續優化速度。此外，本文算法優先優化骨架圖，再利用優化后的骨架圖優化剩余變量以提升運行速度。

圖4是各個算法對Parking-garage數據集優化的軌跡圖。HiPE算法和g2o算法可以優化出較為精確的軌跡結果，從圖4中可以看出ND-PGO算法優化的軌跡圖并不亞于HiPE算法與g2o算法。這表明ND-PGO在保證算法精度的同時，提升了運行速度，總體性能優于HiPE法和g2o法。

4 結束語

本文提出了一種基于嵌套剖分的分層位姿圖優化方法（ND-PGO）。首先利用χ2檢驗剔除外點后提取出一個骨架圖，該骨架圖對圖形幾何的高級表示進行編碼，以快速地初始化變量。為驗證本文算法的有效性，在多個不同的數據集上進行實驗，并將其與最先進的方法進行比較。實驗結果表明，與現有的優化方法相比，ND-PGO在不損失優化精度的情況下具有更高的優化效率。在以后的工作中可以將χ2檢驗改為自適應的，根據數據集自身的特點選擇最優的距離度量。需要指出的是，χ2檢驗剔除一些異常值后會損失一部分信息，未來可進一步研究適合該算法的魯棒核函數，通過降低異常值權重的方式來減少異常值對優化精度的影響。最后，希望能將本文算法整合到一個完整的SLAM系統中，以進一步驗證和增強算法的實際性能。

參考文獻：

［1］Durrant W H， Bailey T. Simultaneous localization and mapping： part I［J］. IEEE Robotics and Automation Magazine， 2006，13（2）： 99-108.

［2］王苗苗，魏國亮，蔡潔，等. 基于偏差矩陣的3D SLAM位姿圖優化算法［J］. 信息與控制， 2023，52（3）： 334-342. （Wang Miaomiao， Wei Guoliang， Cai Jie， et al. Deviation matrix based for 3D SLAM pose graph optimization［J］. Information and Control， 2023， 52（3）： 334-342.）

［3］Carlone L， Rosen D M， Calafiore G， et al. Lagrangian duality in 3D SLAM： verification techniques and optimal solutions［C］//Proc of IEEE International Conference on Intelligent Robots and Systems. Piscataway，NJ： IEEE Press， 2015： 125-132.

［4］Tron R， Rosen D M， Carlone L. On the inclusion of determinant constraints in Lagrangian duality for 3D SLAM［C］//Robotics： Science and Systems Workshop. 2015.

［5］Briales J， Gonzalez J J. Convex global 3D registration with Lagrangian duality［C］//Proc of IEEE Conference on Computer Vision and Pattern Recognition. Piscataway，NJ：IEEE Press， 2017： 4960-4969.

［6］郝亞東， 張奇志， 周亞麗. 基于高斯牛頓的局部優化 SLAM 系統［J］. 傳感器世界， 2018， 24（3）： 7-11. （Hao Yadong， Zhang Qizhi， Zhou Yali. Local optimization SLAM system based on Gauss-Newton method［J］. Sensor World， 2018， 24（3）： 7-11.）

［7］Dou Xinglei， Huo Yuzhou， Liu Yongchang， et al. An unmanned aerial vehicle pose estimation system based on SLAM［C］//Proc of IOP Conference Series： Materials Science and Engineering. London： IOP Press， 2019.

［8］曾俊飛， 楊海清， 吳浩. 面向三維重建的自適應列文伯格-馬夸爾特點云配準方法［J］. 計算機學， 2020， 47（3）： 137-142. （Zeng Junfei， Yang Haiqing， Wu Hao. Adaptive Levenberg-Marquardt cloud registration method for 3D reconstruction［J］. Computer Science， 2020， 47（3）： 137-142.）

［9］Grisetti G， Guadagnino T， Aloise I， et al. Least squares optimization： from theory to practice［J］. Robotics， 2020， 9（3）： 51.

［10］Kyummerle R， Grisetti G， Strasdat H， et al. g2o： a general framework for graph optimization［C］//Proc of IEEE International Confe-rence on Robotics and Automation. Piscataway，NJ： IEEE Press， 2011： 3607-3613.

［11］Dellaert F， Kaess M. Square root SAM： simultaneous localization and mapping via square root information smoothing［J］. The Internatio-nal Journal of Robotics Research， 2006， 25（12）： 1181-1203.

［12］Bosse M， Newman P， Leonard J， et al. Simultaneous localization and map building in large-scale cyclic environments using the Atlas framework［J］. The International Journal of Robotics Research， 2004， 23（12）： 1113-1139.

［13］Paz L M， Tardos J D， Neira J. Divide and conquer： EKF SLAM in o （n）［J］. IEEE Trans on Robotics， 2008， 24（5）： 1107-1120.

［14］Ni K， Steedly D， and Dellaert F. Tectonic SAM： exact， out-of-core， submap-based SLAM［C］//Proc of IEEE International Conference on Robotics and Automation. Piscataway，NJ： IEEE Press， 2007： 1678-1685.

［15］Ni K， Dellaert F. Multi-level submap based SLAM using nested dissection［C］//Proc of IEEE/RSJ International Conference on Intelligent Robots and Systems. Piscataway，NJ： IEEE Press， 2010： 2558-2565.

［16］Grisetti G， Kyummerle R， Stachniss C， et al. Hierarchical optimization on manifolds for online 2D and 3D mapping［C］//Proc of IEEE International Conference on Robotics and Automation. Piscataway，NJ： IEEE Press， 2010： 273-278.

［17］Tang Hengbo， Liu Yunhui， Li Luyang. Pose graph optimization with hierarchical conditionally independent graph partitioning［C］//Proc of IEEE/RSJ International Conference on Intelligent Robots and Systems. Piscataway，NJ： IEEE Press， 2016： 3255-3260.

［18］Tian Yuxin， Wang Yujie， Ouyang Ming， et al. Hierarchical segment-based optimization for SLAM［C］//Proc of IEEE/RSJ International Conference on Intelligent Robots and Systems. Piscataway，NJ： IEEE Press， 2021： 6573-6580.

［19］Korovko A， Robustov D. Partial hierarchical pose graph optimization for SLAM［EB/OL］. （2021）. https：//arxiv.org/abs/2110. 08639.

［20］Guadagnino T， Di Giammarino L， Grisetti G. HiPE： hierarchical initialization for pose graphs［J］. IEEE Robotics and Automation Letters， 2021， 7（1）： 287-294.

［21］Hu G， Khosoussi K， Huang Shoudong. Towards a reliable SLAM back-end［C］//Proc of IEEE/RSJ International Conference on Intelligent Robots and Systems. Piscataway，NJ ： IEEE Press， 2013： 37-43.

［22］Juric＇ A， Kende F， Markovic＇ I， et al. A comparison of graph optimization approaches for pose estimation in SLAM［C］//Proc of the 44th International Convention on Information， Communication and Electronic Technology. Piscataway，NJ： IEEE Press， 2021： 1113-1118.

［23］Rosen D M， Carlone L， Bandeira A S， et al. SE-Sync： a certifiably correct algorithm for synchronization over the special Euclidean group［J］. The International Journal of Robotics Research， 2019， 38（2-3）： 95-125.

［24］Martinec D，Pajdla T. Robust rotation and translation estimation in multiview reconstruction［C］//Proc of IEEE Conference on Computer Vision and Pattern Recognition. Piscataway，NJ： IEEE Press， 2007： 1-8.

［25］Carlone L， Tron R， Daniilidis K， et al. Initialization techniques for 3D SLAM： a survey on rotation estimation and its use in pose graph optimization［C］//Proc of IEEE International Conference on Robotics and Automation. Piscataway，NJ： IEEE Press， 2015： 4597-4604.

［26］Hartley R， Trumpf J， Dai Yuchao， et al. Rotation averaging［J］. International Journal of Computer Vision， 2013， 103： 2LMRmduukxA045i57pkxWZQ==67-305.

［27］Latif Y， Cadena C， Neira J. Robust loop closing over time［M］//Robotics： Science Systems Ⅷ. Cambridge，MA： MIT Press， 2013： 233-240.

［28］George A， Liu J W. Computer solution of large sparse positive definite［M］.[S.l.]：Prentice Hall， 1981.

［29］Blair J swS0yeEuWKBl8kXUqGktXg==R S， Peyton B. An introduction to chordal graphs and clique trees［M］//Graph Theory and Sparse Matrix Computation. New York： Springer， 1993： 1-29.

［30］Yannakakis M. Computing the minimum fill-in is NP-complete［J］. SIAM Journal on Algebraic Discrete Methods， 1981， 2（1）： 77-79.

［31］Liu J W H. Modification of the minimum-degree algorithm by multiple elimination［J］. ACM Trans on Mathematical Software， 1985， 11（2）： 141-153.

［32］Amestoy P R， Davis T A， Duff I S. An approximate minimum degree ordering algorithm［J］. SIAM Journal on Matrix Analysis and Applications， 1996，17（4）： 886-905.

［33］George A. Nested dissection of a regular finite element mesh［J］. SIAM Journal on Numerical Analysis， 1973，10（2）： 345-363.

［34］Lipton R J， Rose D J， Tarjan R E. Generalized nested dissection［J］. SIAM Journal on Numerical Analysis， 1979， 16（2）： 346-358.

［35］李雨潔， 魏國亮， 蔡潔，等. 基于特征分解的快速位姿圖優化算法［J］. 計算機應用研究， 2022， 39（10）： 3065-3070. （Li Yujie， Wei Guoliang， Cai Jie，et al. Fast pose optimization algorithm based on eigen decomposition［J］. Application Research of Computers， 2022， 39（10）： 3065-3070.）