三重対角行列アルゴリズム

数値線形代数において、三重対角行列アルゴリズム（トーマスアルゴリズムとも呼ばれる。ルウェリン・トーマスにちなんで命名）は、ガウス消去法の簡略化された形式であり、三重対角方程式系を解くために使用できる。n個の未知数を持つ三重対角系は次のように記述できる。

${\displaystyle a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},}$

どこ${\displaystyle a_{1}=0}$そして${\displaystyle c_{n}=0}$。

${\displaystyle {\begin{bmatrix}b_{1}&c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}}={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}.}$

このようなシステムの場合、解は次のように得られます。${\displaystyle O(n)}$代わりに操作${\displaystyle O(n^{3})}$ガウス消去法で必要とされる。最初のスイープで除去されるのは${\displaystyle a_{i}}$'s と、(簡略化された) 後方代入によって解が得られます。このような行列の例は、 1 次元ポアソン方程式の離散化や自然三次スプライン補間からよく発生します。

トーマスのアルゴリズムは一般には安定していませんが、行列が対角優位（行または列のいずれか）である場合や対称正定値である場合など、いくつかの特殊なケースでは安定しています。^{[ 1 ] [ 2 ]}トーマスのアルゴリズムの安定性のより正確な特徴付けについては、ハイアムの定理9.12を参照してください。^{[ 3 ]}一般的なケースで安定性が必要な場合は、代わりに部分ピボット付きガウス消去法（GEPP）が推奨されます。^{[ 2 ]}

順方向掃引は、以下のように新しい係数を計算することから成り、新しい係数にはプライム記号が付されます。

${\displaystyle c'_{i}={\begin{cases}{\cfrac {c_{i}}{b_{i}}},&i=1,\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n-1\end{cases}}}$

そして

${\displaystyle d'_{i}={\begin{cases}{\cfrac {d_{i}}{b_{i}}},&i=1,\\{\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n.\end{cases}}}$

解は、逆代入によって得られる。

${\displaystyle x_{n}=d'_{n},}$

${\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1},\quad i=n-1,n-2,\ldots ,1.}$

上記の方法では元の係数ベクトルは変更されませんが、新しい係数も記録しておく必要があります。係数ベクトルが変更される可能性がある場合、記録管理の手間が少ないアルゴリズムは次のようになります。

のために${\displaystyle i=2,3,\dots ,n,}$する

${\displaystyle w:={\cfrac {a_{i}}{b_{i-1}}},}$

${\displaystyle b_{i}:=b_{i}-wc_{i-1},}$

${\displaystyle d_{i}:=d_{i}-wd_{i-1},}$

続いてバック交代

${\displaystyle x_{n}={\cfrac {d_{n}}{b_{n}}},}$

${\displaystyle x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}}\quad {\text{for }}i=n-1,n-2,\dots ,1.}$

acの入力値を変更しないようにスクラッチ領域を使用するC関数としての実装であり、入力値の再利用を可能にする。

void thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double scratch [ restrict X ]) { /*  Ax = d を解きます。ここで、A はベクトル a、b、c からなる三重対角行列です。X  = 方程式の数 。x[] = 最初は入力 d を含み、x を返します。 a[] = サブダイアゴナル、[1, ..., X - 1] からインデックス付け b[] = メインダイアゴナル、[0, ..., X - 1] からインデックス付け c[] = スーパーダイアゴナル、[0, ..., X - 2] からインデックス付け scratch[] = 呼び出し元によって提供される長さ X のスクラッチ領域。これにより、a、b、c を const にすることができます。 この例では実行されません: 手動でコストのかかる共通部分式の削除 */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; x [ 0 ] = x [ 0 ] / b [ 0 ];/* 1からX-1までループ */ for ( int ix = 1 ; ix < X ; ix ++ ) { if ( ix < X -1 ){ scratch [ ix ] = c [ ix ] / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); } x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); }/* X - 2 から 0 までループ */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) x [ ix ] -= scratch [ ix ] * x [ ix + 1 ]; }

三重対角行列アルゴリズムの導出は、ガウス消去法の特殊なケースである。

未知数が${\displaystyle x_{1},\ldots ,x_{n}}$そして、解くべき方程式は以下のとおりです。

${\displaystyle {\begin{alignedat}{4}&&&b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&&&=d_{n}\,.\end{alignedat}}}$

2番目を修正することを検討してください（${\displaystyle i=2}$最初の式との方程式は次のとおりです。

${\displaystyle ({\mbox{equation 2}})\cdot b_{1}-({\mbox{equation 1}})\cdot a_{2}}$

これにより、以下の結果が得られます。

${\displaystyle (b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}.}$

ご了承ください${\displaystyle x_{1}}$は2番目の式から消去されました。修正された2番目の式を3番目の式に同様の手法で適用すると、次のようになります。

${\displaystyle (b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,}$

この時${\displaystyle x_{2}}$が除去されました。この手順を繰り返すと、${\displaystyle n^{th}}$行; (修正済み)${\displaystyle n^{th}}$方程式には未知数が 1 つだけ含まれます。${\displaystyle x_{n}}$これは解くことができ、その後、解くために使用できます。${\displaystyle (n-1)^{th}}$方程式を立てていき、すべての未知数が解かれるまでこれを繰り返します。

明らかに、修正された方程式の係数を明示的に示すと、ますます複雑になる。手順を検討することで、修正された係数（チルダで表記）は、代わりに再帰的に定義することができる。

${\displaystyle {\tilde {a}}_{i}=0\,}$

${\displaystyle {\tilde {b}}_{1}=b_{1}\,}$

${\displaystyle {\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,}$

${\displaystyle {\tilde {c}}_{1}=c_{1}\,}$

${\displaystyle {\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,}$

${\displaystyle {\tilde {d}}_{1}=d_{1}\,}$

${\displaystyle {\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,}$

解決プロセスをさらに加速するために、${\displaystyle {\tilde {b}}_{i}}$（ゼロ除算のリスクがない場合）割り算が可能であれば、それぞれプライム記号で表記される新しい修正係数は次のようになります。

${\displaystyle a'_{i}=0\,}$

${\displaystyle b'_{i}=1\,}$

${\displaystyle c'_{1}={\frac {c_{1}}{b_{1}}}\,}$

${\displaystyle c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,}$

${\displaystyle d'_{1}={\frac {d_{1}}{b_{1}}}\,}$

${\displaystyle d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,}$

これにより、上記の元の未知数と係数で定義された同じ未知数と係数を持つ、以下のシステムが得られます。

${\displaystyle {\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,}$

最後の式には未知数が1つしか含まれていません。これを解くと、次の最後の式も未知数が1つに減算されるため、この逆代入法を用いてすべての未知数を求めることができます。

${\displaystyle x_{n}=d'_{n}\,}$

${\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad$ ; i=n-1,n-2,...,1.}

場合によっては、特に周期境界条件が関係するような状況では、三重対角系のわずかに摂動した形式を解く必要があるかもしれません。

$${\displaystyle {\begin{alignedat}{4}&a_{1}x_{n}&&+b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&+c_{n}x_{1}&&=d_{n}\,.\end{alignedat}}}$$

この場合、シャーマン・モリソン公式を利用することで、ガウス消去法の追加演算を回避しつつ、トーマスアルゴリズムを使用することができます。この方法では、入力と疎な補正ベクトルの両方について、システムの非巡回型修正版を解き、それらの解を組み合わせる必要があります。純粋な三重対角行列アルゴリズムの順方向部分を共有できるため、両方の解を同時に計算すれば効率的に処理できます。

次のように示す場合： $${\displaystyle A={\begin{bmatrix}b_{1}&c_{1}&&&a_{1}\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\c_{n}&&&a_{n}&b_{n}\end{bmatrix}},x={\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}},d={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}}$$

すると、解くべき連立方程式は次のようになります。$${\displaystyle Ax=d}$$

この場合、係数は${\displaystyle a_{1}}$そして${\displaystyle c_{n}}$一般的に、 はゼロではないため、それらが存在するとトーマスアルゴリズムを直接適用することはできません。したがって、 を考慮することができます。${\displaystyle B\in \mathbb {R} ^{n\times n}}$そして${\displaystyle u,v\in \mathbb {R} ^{n}}$以下のとおりです。 $${\displaystyle B={\begin{bmatrix}b_{1}-\gamma &c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}-{\frac {c_{n}a_{1}}{\gamma }}\end{bmatrix}},u={\begin{bmatrix}\gamma \\0\\0\\\vdots \\c_{n}\end{bmatrix}},v={\begin{bmatrix}1\\0\\0\\\vdots \\a_{1}/\gamma \end{bmatrix}}.}$$ どこ${\displaystyle \gamma \in \mathbb {R} }$は選択すべきパラメータである。行列Aは次のように再構成できる。${\displaystyle A=B+uv^{\mathsf {T}}}$解は次のようにして得られます。^{[ 4 ]}まず、トーマスアルゴリズムを適用して、2 つの三重対角方程式系を解きます。 $${\displaystyle By=d\qquad \qquad Bq=u}$$

次に、シャーマン・モリソン公式を用いて解xを再構築する。 $${\displaystyle {\begin{aligned}x&=A^{-1}d=(B+uv^{T})^{-1}d=B^{-1}d-{\frac {B^{-1}uv^{T}B^{-1}}{1+v^{T}B^{-1}u}}d=y-{\frac {qv^{T}y}{1+v^{T}q}}\end{aligned}}}$$

acの入力値を変更しないようにスクラッチ領域を使用するC関数としての実装であり、入力値の再利用を可能にする。

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double u [ restrict X ]) {/* Ax = v を解く。ここで A はベクトル a、b、c からなる巡回三重対角行列である。 X = 方程式の数 x[] = 最初は入力vが格納され、xを返します。インデックスは[0, ..., X - 1]です。 a[] = サブダイアゴナル、インデックスは [1, ..., X - 1] から規則的に取得され、a[0] は左下隅です。 b[] = 主対角線、インデックスは[0, ..., X - 1] c[] = 上対角線、インデックスは [0, ..., X - 2] から規則的に付けられ、c[X - 1] は右上隅です。 cmod[]、u[] = それぞれ長さXのスクラッチベクトル *//* それぞれ循環三対角システムの左下隅と右上隅 */const double alpha = a [ 0 ];const double beta = c [ X - 1 ];/* 任意だが、ゼロ除算を回避するように選択されている */const double gamma = - b [ 0 ];cmod [ 0 ] = c [ 0 ] / ( b [ 0 ] - gamma );u [ 0 ] =ガンマ/ ( b [ 0 ] -ガンマ);x [ 0 ] /= ( b [ 0 ] - gamma );/* 1からX-2までルー​​プする */for ( int ix = 1 ; ix + 1 < X ; ix ++ ) {const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]);cmod [ ix ] = c [ ix ] * m ;u [ ix ] = ( 0.0f - a [ ix ] * u [ ix - 1 ]) * m ;x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ;}/* X - 1 を処理する */const double m = 1.0 / ( b [ X - 1 ] - alpha * beta / gamma - a [ X - 1 ] * cmod [ X - 2 ]);u [ X - 1 ] = ( alpha - a [ X - 1 ] * u [ X - 2 ]) * m ;x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;/* X - 2 から 0 までループする */for ( int ix = X - 2 ; ix >= 0 ; ix -- ) {u [ ix ] -= cmod [ ix ] * u [ ix + 1 ];x [ ix ] -= cmod [ ix ] * x [ ix + 1 ];}const double fact = ( x [ 0 ] + x [ X - 1 ] * alpha / gamma ) / ( 1.0 + u [ 0 ] + u [ X - 1 ] * alpha / gamma );/* 0からX-1までループする */for ( int ix = 0 ; ix < X ; ix ++ )x [ ix ] -= fact * u [ ix ];}

上記で検討した三重対角系のわずかに摂動した形式を解く別の方法もある。^{[ 5 ]}次元の2つの補助線形システムを考えてみよう。${\displaystyle (n-1)\times (n-1)}$: $${\displaystyle {\begin{aligned}\qquad \ \ \ \ \ b_{2}u_{2}+c_{2}u_{3}&=d_{2}\\a_{3}u_{2}+b_{3}u_{3}+c_{3}u_{4}&=d_{3}\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=d_{i}\\\dots \\a_{n}u_{n-1}+b_{n}u_{n}\qquad &=d_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1\qquad \qquad {\begin{aligned}\qquad \ \ \ \ \ b_{2}v_{2}+c_{2}v_{3}&=-a_{2}\\a_{3}v_{2}+b_{3}v_{3}+c_{3}v_{4}&=0\\a_{i}v_{i-1}+b_{i}v_{i}+c_{i}v_{i+1}&=0\\\dots \\a_{n}v_{n-1}+b_{n}v_{n}\qquad &=-c_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1}$$

便宜上、さらに定義します${\displaystyle u_{1}=0}$そして${\displaystyle v_{1}=1}$解決策を見つけることができる${\displaystyle \{u_{2},u_{3}\dots ,u_{n}\}}$そして${\displaystyle \{v_{2},v_{3}\dots ,v_{n}\}}$トーマスアルゴリズムを2つの補助三重対角システムに適用する。

解決策${\displaystyle \{x_{1},x_{2}\dots ,x_{n}\}}$は次のように表すことができます。 $${\displaystyle x_{i}=u_{i}+x_{1}v_{i}\qquad i=1,2,\dots ,n}$$

実際、第2補助システムの各方程式に${\displaystyle x_{1}}$第一補助システムの対応する方程式を加え、表現を用いる${\displaystyle x_{i}=u_{i}+x_{1}v_{i}}$方程式の数がすぐにわかります元のシステムの2からnまでの条件は満たされています。あとは方程式番号を満たすだけです。1.そのためには、次の式を検討してください。${\displaystyle i=2}$そして${\displaystyle i=n}$そして代替する${\displaystyle x_{2}=u_{2}+x_{1}v_{2}}$そして${\displaystyle x_{n}=u_{n}+x_{1}v_{n}}$元のシステムの最初の式に代入します。これにより、次の1つのスカラー方程式が得られます。${\displaystyle x_{1}}$: $${\displaystyle b_{1}x_{1}+c_{1}(u_{2}+x_{1}v_{2})+a_{1}(u_{n}+x_{1}v_{n})=d_{1}}$$

したがって、以下のことが分かります。 $${\displaystyle x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}}$$

acの入力値を変更しないようにスクラッチ領域を使用するC関数としての実装であり、入力値の再利用を可能にする。

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double v [ restrict X ]) {/* まず、ix == 0 を無視して、長さ X - 1 の連立方程式を 2 つの右辺について解きます */cmod [ 1 ] = c [ 1 ] / b [ 1 ];v [ 1 ] = - a [ 1 ] / b [ 1 ];x [ 1 ] = x [ 1 ] / b [ 1 ];/* 2からX-1までループする */for ( int ix = 2 ; ix < X - 1 ; ix ++ ) {const double m = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ]);cmod [ ix ] = c [ ix ] * m ;v [ ix ] = ( 0.0f - a [ ix ] * v [ ix - 1 ]) * m ;x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ;}/* X - 1 を処理する */const double m = 1.0 / ( b [ X - 1 ] - a [ X - 1 ] * cmod [ X - 2 ]);cmod [ X - 1 ] = c [ X - 1 ] * m ;v [ X - 1 ] = ( - c [ 0 ] - a [ X - 1 ] * v [ X - 2 ]) * m ;x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;/* X - 2 から 1 までループする */for ( int ix = X - 2 ; ix >= 1 ; ix -- ) {v [ ix ] -= cmod [ ix ] * v [ ix + 1 ];x [ ix ] -= cmod [ ix ] * x [ ix + 1 ];}x [ 0 ] = ( x [ 0 ] - a [ 0 ] * x [ X - 1 ] - c [ 0 ] * x [ 1 ]) / ( b [ 0 ] + a [ 0 ] * v [ X - 1 ] + c [ 0 ] * v [ 1 ]);/* 1からX-1までループする */for ( int ix = 1 ; ix < X ; ix ++ )x [ ix ] += x [ 0 ] * v [ ix ];}

どちらの場合も、解くべき補助システムは真に三重対角行列であるため、システムを解く全体の計算複雑度は${\displaystyle Ax=d}$システムの次元nに関して線形性を維持する、つまり${\displaystyle O(n)}$算術演算。

その他の状況では、方程式系はブロック三重対角行列（ブロック行列を参照）になる場合があり、その場合、より小さな部分行列が上記の行列系の個々の要素として配置されます（例：2Dポアソン問題）。このような状況のために、ガウス消去法の簡略化された形式が開発されています。^{[ 6 ]}

アルフィオ・クアルテローニ、サッコ、サレリ共著の教科書『数値数学』には、一部の除算を回避し（代わりに乗算を使用する）、一部のコンピュータアーキテクチャで有利となる、アルゴリズムの改良版が掲載されている。

並列三重対角ソルバーは、GPU を含む多くのベクトルおよび並列アーキテクチャ向けに公開されています^{[ 7 ] [ 8 ] 。}

並列三重対角ソルバーとブロック三重対角ソルバーの詳細な解説については^{[ 9 ]を参照のこと。}

• Conte, SD; de Boor, C. (1972).初等数値解析. McGraw-Hill、ニューヨーク. ISBN 0070124469。
• この記事は、 GFDLライセンスの下で公開されているCFD-Wikiの記事「Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)」のテキストを組み込んでいます。

• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). 「セクション 2.4」 . Numerical Recipes: The Art of Scientific Computing (第 3 版). ニューヨーク: Cambridge University Press. ISBN 978-0-521-88068-8。