Sparse Bundle Adjustment

概要

　Sparse Bundle Adjustment(SBA) [1] はStructure from Motion(SfM)の手法の一種である．

　SfMは一般的に再投影誤差の最小化問題として記述できるが，推定の対象となるランドマークやカメラ姿勢の数が膨大になり，計算コストが非常に大きくなってしまう．SBAは再投影誤差のヤコビ行列がスパースになることに着目し，計算量を大幅に削減した手法であり，次のような特徴がある．

利点

• 大規模な復元を高速に行うことができる
• ある視点からは見えないランドマークがあったとしても復元を行うことができる
• 内部行列を指定できるため，正投影以外にも多様なカメラモデルを用いることができる

欠点

• 問題に応じてハイパーパラメータを調整しなければならない
• 誤差関数のヤコビ行列を計算する際に， so(3) あるいは四元数などに関する微分が現れるため，手法が複雑である
• 誤差関数が減る方向にしかパラメータ探索を行わないため，容易に局所解に落ちてしまう

問題設定

入力

　SBAは，各視点から観測されたランドマークの像の集合 X を入力とする．

X=[x⊤11,…,x⊤1m,x⊤21,…,x⊤2m,…,x⊤n1,…,x⊤nm]

出力

• 3次元空間におけるランドマーク座標 bi,i=1,…,n
• カメラ姿勢 aj=[tj,ωj],j=1,…,m

　ただし tj∈R3 は並進を表すベクトルであり，ωj∈so(3) はカメラの回転を表す行列 R∈SO(3) に対応するリー代数の元である．

目的

　投影モデルを Q(ai,bj) とし，以下の誤差関数を最小化するような P=[a,b]=[a⊤1,…,a⊤m,b⊤1,…,b⊤n] を求める．

(11)E(P)=n∑i=1m∑j=1dΣxij(Q(aj,bi),xij)2

ここで dΣx は x に対応する分散共分散行列を Σx として

dΣx(x1,x2)=√(x1−x2)⊤Σ−1x(x1−x2)

で定義される距離関数である．

(12)ˆX=[ˆx⊤11,…,ˆx⊤1m,ˆx⊤21,…,ˆx⊤2m,…,ˆx⊤n1,…,ˆx⊤nm]⊤

(13)ˆxij=Q(aj,bi)

(14)ΣX=diag(Σx11,…,Σx1m,Σx21,…,Σx2m,…,Σxn1,…,Σxnm)

とおけば，誤差を次のように表現することができる．

E(P)=(X−ˆX)⊤Σ−1X(X−ˆX)

解法の概要

　SBAでは，誤差関数を最小化するような P を見つけるため， P(t) を逐次的に更新し，誤差関数を探索する．すなわち，時刻 t における P の更新量を δ(t)P=[δ⊤a1,…,δ⊤am,δ⊤b1,…,δ⊤bn] として，

(15)P(t+1)←P(t)+δ(t)P

というふうに P(t) を更新することで誤差関数を最小化するような P を見つける．

　更新量 δ(t)P の計算にはLM法を用いる．LM法の更新式は次のように表される．

(16)[J⊤Σ−1J+λI]δ(t)P=J⊤Σ−1[X−ˆX]

J は ˆX のヤコビ行列 J=∂ˆX∂P|P=P(t) であり， λ∈R,λ≥0 は damping parameter である．

　SBAでは，J の構造に着目し， (16) をより小さい複数の線型方程式に分解する．さらに，分解によって得られた方程式がスパースな行列によって構成されていることに着目し，計算を高速化している．

解法

線型方程式の分解

　まず J を分解する． P の定義より， A=∂ˆX∂a,B=∂ˆX∂b とおけば， J は

(17)J=∂ˆX∂P=∂ˆX∂(a,b)=[A,B]

と書ける．

　次に (16) の右辺を分解する． (17) を用いると， (16) の右辺は

ϵa=A⊤Σ−1(X−ˆX)ϵb=B⊤Σ−1(X−ˆX)

とおくことによって，

J⊤Σ−1(X−ˆX)=[ϵaϵb]

と書ける．

　さらに (16) の左辺を分解する．左辺の J⊤Σ−1J という項は大きく4つの行列に分解することができる．

(18)J⊤Σ−1J=[A⊤B⊤]Σ−1[AB]=[A⊤Σ−1AA⊤Σ−1BB⊤Σ−1AB⊤Σ−1B]=[UWW⊤V]

　以上の結果を用いると， (16) は

[[UWW⊤V]+[λI00λI]][δaδb]=[ϵaϵb]

という形にすることができる．さらに，

U∗=U+λIV∗=V+λI

とおけば，

[U∗WW⊤V∗][δaδb]=[ϵaϵb]

となる．この両辺に

[I−WV∗−10I]

という行列を左から作用させると，

(19)[I−WV∗−10I][U∗WW⊤V∗][δaδb]=[I−WV∗−10I][ϵaϵb]

(20)[U∗−WV∗−1W⊤0W⊤V∗][δaδb]=[ϵa−WV∗−1ϵbϵb]

という形にすることができる．ここから2つの方程式を取り出す．すると， (20) において左辺の行列の右上が 0 になったことから， δb を含まない δa についての式 (21) を得ることができる．

(21)(U∗−WV∗−1W⊤)δa=ϵa−WV∗−1ϵb

(22)V∗δb=ϵb−W⊤δa

したがって，(21) を先に解き，得られた δa を (22) に代入すれば δb を得ることができる．

具体的な計算

　前節では，LM法を分解し，より少ない計算量で更新量 δP を求める方法を述べた．ここでは，実際にヤコビ行列 J を計算し，その具体的なかたちを求める．

　まず，ヤコビ行列 J はスパースな行列になる．これは，∀j≠k について

∂Q(aj,bi)∂ak=0

∀i≠k について

∂Q(aj,bi)∂bk=0

が成り立つためである．

　例えば，n=4 ，m=3 のとき， Aij=∂Q(aj,bi)∂aj ， Bij=∂Q(aj,bi)∂bi とおけば，J は

(23)J=[A1100B110000A120B1200000A13B13000A21000B21000A2200B220000A230B2300A310000B3100A32000B32000A3300B330A4100000B410A420000B4200A43000B43]

となる．

　では Aij や Bij の具体的なかたちを求めてみよう．姿勢パラメータ aj=[tj,ωj] に関する微分 Aij=∂Q(aj,bi)∂aj は次のようになる．

∂ˆxij∂tj=∂π(p)∂p|p=R(ωj)bi+tj⋅∂(R(ωj)bi+v)∂v|v=tj=∂π(p)∂p|p=R(ωj)bi+tj

∂ˆxij∂ωj=∂π(p)∂p|p=R(ωj)bi+tj⋅∂(R(v)bi+tj)∂v|v=ωj=∂π(p)∂p|p=R(ωj)bi+tj⋅∂(R(v)bi)∂v|v=ωj

　ここで， ∂(R(v)bi)∂v はGallegoら [2] による計算結果を用いることができる．

∂(R(v)bi)∂v=−R(v)[bi]×vv⊤+(R(v)⊤−I)[v]×||v||2

　3次元点の座標 bi に関する微分 Bij=∂Q(aj,bi)∂bi は次のようになる．

∂ˆxij∂bi=∂π(p)∂p|p=R(ωj)bi+tj⋅∂(R(ωj)v+tj)∂v|v=bi=∂π(p)∂p|p=R(ωj)bi+tj⋅R(ωj)

　以上より， Aij と Bij の具体的なかたちを求めることができた．あとは，

1. 上記で得られた Aij と Bij (23) に代入して J を求める
2. (18) にしたがって U,V,W を求める
3. (21) と (22) によって姿勢パラメータ a と3次元点の座標 b それぞれについての更新量 δa と δb を求める

という3つのステップによって更新量を求めることができる．

計算量の削減

　前節までで更新量の計算 (16) を2つの計算 (21) (22) に分解する過程を見た．(16) ， (21) ， (22) はいずれも線型方程式とみなすことができる．

　線型方程式 y=Ax,x∈Rn,y∈Rm,A∈Rn×m の解は

x=(A⊤A)−1A⊤y=K−1A⊤yK=A⊤A,K∈Rn×n

によって得られるが，行列 K のサイズが大きくなると解を求めるための計算量が急激に増加する．これは， n×n 行列の逆行列を計算するアルゴリズムが O(n2.3) 〜 O(n3) 程度の計算量をもつことに起因する [4] ．したがって，線型方程式を高速に解くには，問題の構造を見極め， K の逆行列を直接計算することを避けて計算量を減らす必要がある．

　SBAでは， (16) を直接解くのではなく，それを分割して得た (21) と (22) をそれぞれ解くことで δP を得ている．さらに， V∗ がスパースであるという性質に基づいて計算量を大幅に削減している．(23) で定義された J を用いて V∗ を計算すると次のようになる．

V∗=[V∗10000V∗20000V∗30000V∗4]

ただし V∗i は

Vi=m∑j=1B⊤ijΣ−1ijBijV∗i=Vi+λI.

である．

　(21) には V∗ の逆行列が両辺に含まれている．また， (22) を解いて δb を得る際にも両辺に左から V∗ の逆行列をかける必要がある．V∗ のサイズが大きいとその逆行列を求めるのに多大なコストがかかってしまう．しかし， V∗ がスパースな行列であることに着目すると， V∗ の逆行列は

(24)V∗−1=[V∗−110000V∗−120000V∗−130000V∗−14]

となるため， V∗i,i=1,…,n のそれぞれについて逆行列を求めればよいことがわかる．結果として V∗ の逆行列の計算量はランドマーク数 n に対して線型に増加することになり， V∗ の逆行列を直接求めるのと比較すると計算量を一気に削減できる．

　δa を求める際には， S=U∗−WV∗−1W⊤ の逆行列を (21) の両辺に左からかける必要がある．しかし，一般的にランドマーク数 n よりもカメラの視点数 m の方が圧倒的に小さい (m≪n) ため， S のサイズは V∗ と比べると圧倒的に小さい．したがって， S の逆行列を求める処理は全体の計算量にはほとんど影響しない．

　問題のサイズ(視点数や復元対象となるランドマークの数)が大きいときは， (16) を直接解いて δP を得るよりも， (21) (22) (24) によって δa と δb をそれぞれ計算し結合することで δP を得るほうが圧倒的に高速である．

改良

　Agarwalらは inexact Newton method とPCG(Preconditioned Conjugate Gradients)法を組み合わせることでより高速に更新量を求める手法を提案している [5] ．

　SBAでは，誤差関数の更新則 (16) を変形し， (21) (22) という2つの線型方程式を解く問題に落とし込んでいる．このうち (22) は V∗ のスパース性を利用して高速に解くことができたが， (21) は S の逆行列を直接計算する必要があった．SBAでは (21) と (22) を解くことで各ステップにおける”厳密な”更新量 δP を求めている．これに対してAgarwalら [5] は必ずしも (21) (22) の厳密な解を求める必要はなく，より高速な近似的計算によって厳密解を代替できることを主張している． すなわち，最終的な目的は誤差関数 (11) を十分小さくするような解を見つけることであり，もしそれが達成できるのであれば，必ずしも各ステップにおいて厳密な更新量を見つける必要はないのである．各ステップにおいてより少ない計算量で近似的に更新量を求められれば，最適解に達するまでのステップ数が増えたとしても，解に到達するまでの計算量の総和を小さくすることができる可能性がある．

LM法

概要

　勾配の2次微分の情報を利用する最適化手法の一種Gauss-Newton法は収束性が保証されていない．LM法 [3] はGauss-Newton法と最急降下法を組み合わせることで収束性を保証したアルゴリズムである [6] ．

　β をパラメータとするベクトル値関数 f(β) と，目標値ベクトル y について，次で定義される誤差 d2Σ(y,f(β)) を最小化するような β を見つける問題を考える．

(25)d2Σ(y,f(β))=(y−f(β))⊤Σ−1(y−f(β))

　LM法はGauss-Newton法と最急降下法を組み合わせた手法であると解釈することがすることができる．

　J を関数 f のヤコビ行列 ∂f∂β ， δ を β の更新量として，Gauss-Newton法，最急降下法，LM法それぞれによる δ を示す．

δGN=(J⊤Σ−1J)−1J⊤Σ−1[y−f(β)]δGD=J⊤Σ−1[y−f(β)]δLM=(J⊤Σ−1J+λI)−1J⊤Σ−1[y−f(β)]

I は単位行列であり， λ∈R,λ>0 は damping parameter と呼ばれる値である．それぞれの式を見比べると，

• LM法による更新量の計算方法はGauss-Newton法と最急降下法を組み合わせたものである
• Gauss-Newton法と最急降下法のどちらの性質を強くするかを damping parameter がコントロールしている

ということがわかる．Damping parameter を大きくすると最急降下法の性質が強くなり，小さくするとGauss-Newton法の性質が強くなる(誤差が発散する可能性が高くなる)．

　時刻 t におけるパラメータ β の値を β(t) とする．このとき，LM法は次に示す規則にしたがってパラメータ β を更新する．

• 誤差が減少する (f(β(t)+δ)<f(β(t))) ならばパラメータを β(t+1)←β(t)+δ と更新する．
• 誤差が減少しない (f(β(t)+δ)≥f(β(t))) ならば λ の値を大きくし，再度更新量 δ を計算し直す．誤差が減少するような δ が見つかるまでこれを繰り返す．

　LM法は，damping parameter を変化させながら誤差が必ず減少するような更新量 δ を探し出すことで収束を保証している．

導出

　Σ を分散共分散行列とし，誤差をmahalanobis距離によって次のように定義する．

(26)d2Σ(y,f(β+δ))=(y−f(β+δ))⊤Σ−1(y−f(β+δ))

　関数 f を f(β+δ)≈f(β)+Jδ と近似すると， (26) は

d2Σ(y,f(β+δ))≈(y−f(β)−Jδ)⊤Σ−1(y−f(β)−Jδ)=(y−f(β))⊤Σ−1(y−f(β))−2(y−f(β))⊤Σ−1Jδ+δ⊤J⊤Σ−1Jδ

となる．これを δ で微分して 0 とおくと，

J⊤Σ−1Jδ=J⊤Σ−1[y−f(β)]

が得られる．左辺に λI という項を組み込んでしまえば，即座にLM法が得られる．

(J⊤Σ−1J+λI)δ=J⊤Σ−1[y−f(β)]

[1]	Lourakis, Manolis IA, and Antonis A. Argyros. “SBA: A software package for generic sparse bundle adjustment.” ACM Transactions on Mathematical Software (TOMS) 36.1 (2009): 2.

[2]	Gallego, Guillermo, and Anthony Yezzi. “A compact formula for the derivative of a 3-D rotation in exponential coordinates.” Journal of Mathematical Imaging and Vision 51.3 (2015): 378-384.

[3]	Levenberg, Kenneth. “A method for the solution of certain non-linear problems in least squares.” Quarterly of applied mathematics 2.2 (1944): 164-168.

[4]	Coppersmith, Don, and Shmuel Winograd. “Matrix multiplication via arithmetic progressions.” Journal of symbolic computation 9.3 (1990): 251-280.

[5]	(1, 2) Agarwal, Sameer, et al. “Bundle adjustment in the large.” European conference on computer vision. Springer, Berlin, Heidelberg, 2010.

[6]	Wright, Stephen, and Jorge Nocedal. “Numerical optimization.” Springer Science 35.67-68 (1999): 7.