代数的ブロック化多色順序付け法のオープンソース実装を公開しました！

概要

今回公開した”AlgebraicBlockMulticolorOrdering”は、岩下・高橋・中島(2009)『代数ブロック化多色順序付け法による並列化ICCGソルバの性能評価』、情報処理学会研究報告、Vol.2009-HPC-121 No.11による「代数的ブロック化多色順序付け法」を実装したものになります。
詳しくは後述しますが簡単に言うと、この手法を使うと、ガウスザイデル法を並列化する時に元の格子の幾何形状を知らなくても収束性劣化を抑えることができるようになります。

下の図は、代数的ブロック化多色順序付け法を導入した時の効果を表したものです。
横軸が反復回数、縦軸がその時の残差（0だと完全に理論解と一致）を表しています。ですので、反復するほど残差が減少していきます。
ここで、緑色のSymGSというのは「対称ガウスザイデル法」というガウスザイデル法の1種であり、理想状態だと思ってください。
これを、通常の多色順序付け法(MulticolorOrdering)によって並列化した場合が茶色の線になります。ご覧のように同じ反復回数でも残差が大きく、収束性を劣化させている事がわかります。
一方で、ブロック化多色順序付け法(BlockMulticolorOrdering)を使った暗緑色の線では、収束性の劣化が抑えられている事が分かると思います。

このような効果をもたらすものが、代数的ブロック化多色順序付け法で、今回はこの手法を実装しオープンソースとして公開しました。
以下、この手法について「ガウスザイデル法をどのように並列化するか」という観点から、簡単ですが解説したいと思います。

ガウスザイデル法

ガウスザイデル法そのものについては、理学部・工学部であれば大学学部で習うと思いますが、ざっとおさらいしておきます。
「もちろん常識だよ！」という方は次の章まで飛ばしてもらって構いません。

ガウスザイデル法とは、線形方程式系の求解法の一種です。線形方程式系とは、連立一次方程式とも呼ばれている、中学校でも習う
\[ \begin{align*}
3x + 2y & = 1 \\
5x + 8y & = 9
\end{align*} \]
みたいなものです。ここでは一般化して、
\[ \sum_j a_{ij} x_j = b_i \]
と書くことにします。ここで、\(a_{ij}\)は↑の例だと3,2,5,8の係数のことで、\(x_j\)は\(x,y\)の未知変数、\(b_i\)は1,9の右辺定数に該当します。
\(\sum_j\)は\(j\)について総和するものです。

今、未知変数がn個あるとすると、方程式もn個あるはずな（というか無いと解けない）ので、iとjは0からn-1までの範囲の整数をとります。n=2の場合に書き下すと
\[ \begin{align*}
a_{00} x_0 + a_{01} x_1 & = b_0 \\
a_{10} x_0 + a_{11} x_1 & = b_1
\end{align*} \]
となります。式が2つあり、それぞれ上から順に式0と式1と呼ぶことにします。

さて、この線形方程式系を解く方法は世の中にたくさんありますが、ここでは何らかの反復式を繰り返すことで解を収束させることを目的にした反復法を考えます。
そのためには、↑の\(a_{ij} x_j = b_i\)を\(x’ = f(x)\)みたいな形にしなければなりません。

簡単にぱっと思いつくのは、式iを変数\(x_i\)の反復式に変形することです。つまり、
\[ a_{ii} x_i + \sum_{j \ne i} a_{ij} x_j = b_i \]
と、\(x_i\)の部分だけ取り出してやって、あとは変形して
\[ x_i = \frac{1}{a_{ii}} \left( b_i – \sum_{j \ne i} a_{ij} x_j \right) \]
という式を出します。
これを反復式にするためには、右辺を既知の値とすれば良いので、k番目の反復時の\(x_i\)を\(x_{i, k}\)と表記すると
\[ x_{i, k} = \frac{1}{a_{ii}} \left( b_i – \sum_{j \ne i} a_{ij} x_{j, k-1} \right) \]
という簡単な反復式ができます。これを繰り返せばなんとなく収束しそうですね。コードで書くと以下のようになります。

void Jacobi(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(n)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x_prev[j];
            }
        }
        x_next[i] = x_i / a[i][i];
    }
}

この方法はヤコビ法と呼ばれていて、最も簡単で、かつ各式を並列で計算できるという長所があります。
ただしこのヤコビ法には欠点があって、なにしろ収束がとても遅いのです。下手するとどこかで飽和して収束しないということにもなります。

これを改良するための方法が、ガウスザイデル法です。
方法は簡単です。今、i番目の変数\(x_i\)の計算時には、右辺は全てひとつ前の反復の値を使っています。
しかし、逐次実行していれば、i番目の式を計算している時には既に0からi-1の今の反復した後の値が分かっているはずですので、それを使うのです。つまり
\[ x_i = \frac{1}{a_{ii}} \left( b_i – \sum_{j < i} a_{ij} x_j – \sum_{j > i} a_{ij} x_j\right) \]
とします。コードで書くと以下のようになります。

void GaussSeidel(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < i; j++)
        {
            x_i -= a[i][j] * x_next[j];
        }
        for(auto j = i + 1; j < n; j++)
        {
            x_i -= a[i][j] * x_prev[j];
        }
        x_next[i] = x_i / a[i][i];
    }
}

あるいは単純にxを1つにまとめて以下のように書くほうが簡潔で分かりやすいかもしれません。

void GaussSeidel(const double** a, double* x, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
}

この結果、下の図のように、収束性が格段に上がります。

このように、ガウスザイデル法の方が収束性が良い一方、欠点もあります。
それは、先に述べた通り、この方法があくまで逐次実行することを前提としているという問題です。つまり、このままだと並列化できません。

ということで、太古の昔は、並列性を取るならヤコビ法、収束性を取るならガウスザイデル法、という時期もあったようです。
しかしなんとかしてガウスザイデル法も高速化したい・・・そのための方法が、次に説明する多色順序付け法です。

多色順序付け法

多色順序付け法は、Multi-color OrderingやMCと呼ばれるものです。

多色順序付け法の説明に入る前に。まず、先の線形方程式系
\[ \sum_j a_{ij} x_j = b_i \]
は、
\[A \boldsymbol{x} = \boldsymbol{b} \]
と、係数行列\(A\)、未知変数ベクトル\(\boldsymbol{x}\)、右辺定数ベクトル\(\boldsymbol{b}\)で表せることは簡単に分かると思います。

ここで、線形方程式系が現実世界でどんな状況で出てくるか考えてみると、大体が\(A\)は疎行列になることがほとんどです。疎行列とは、\(A\)のうち殆どの要素が0であるような行列です。
どういうことかというと、例えば、熱伝導方程式（あるいはラプラス方程式）
\[ \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} = 0 \]
を解きたい時、2階偏微分が
\[\begin{align*}
\left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i,j,k} & = \frac{\phi_{i-1, j, k} – 2 \phi_{i, j, k} + \phi_{i+1, j, k}}{{\Delta x}^2} \\
\left. \frac{\partial^2 \phi}{\partial y^2} \right|_{i,j,k} & = \frac{\phi_{i, j-1, k} – 2 \phi_{i, j, k} + \phi_{i, j+1, k}}{{\Delta y}^2} \\
\left. \frac{\partial^2 \phi}{\partial z^2} \right|_{i,j,k} & = \frac{\phi_{i, j, k-1} – 2 \phi_{i, j, k} + \phi_{i, j, k+1}}{{\Delta z}^2} \\
\end{align*}\]
という形で中心差分で離散化できるので、3次元空間を正方形で区切った時に各頂点の値を求めるのには、上下左右奥手前の6点とだけつながっているモデルを作ることができます。
そのような時にそこから生成される線形方程式系は、1つの式が1つの頂点に対応しており、つまり行列で言えば、1行が1点に対応します。
とすると、1点はたかだか6点としかつながっていないことから、1行は6列だけ非ゼロで、あとは全てゼロとなります。

熱伝導方程式とはちょっと違いますが、平面で周囲8点（上下左右斜め）から影響を受けるという、もうちょっとまともに離散化したモデルだと行列は以下のような感じに見えます。

この行列はサイズが64×64の正方行列で、対角項が8（紫色）、非対角項には8点だけ-1が入っており（赤色）、他の要素は0（白色）になります。

さて、ということは、ガウスザイデル法の計算方法を再掲しますが

void GaussSeidel(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < i; j++)
        {
            x_i -= a[i][j] * x_next[j];
        }
        for(auto j = i + 1; j < n; j++)
        {
            x_i -= a[i][j] * x_prev[j];
        }
        x_next[i] = x_i / a[i][i];
    }
}

で、ある行x[i]を計算するにあたっては、実はたかだか8個のx[j]からしか影響を受けません。
例えば、0点（行）目の計算では、1,8,9点（列）目からしか影響を受けないことが分かります。
また、2点（行）目は、1,3,9,10,11点（列）目からしか影響を受けません。
ということは、0行目と2行目はお互いに影響しておらず、実は並列で計算しても問題ないということが分かります。

このようにお互いに影響しない点同士を集めてグループにして、そのグループの中では並列で実行する方法を多色順序付け法と呼びます。
順序付け(ordering)と呼んでいるのは、計算順を上から順ではなくて、お互いに影響しない点のグループ順にする、という意味です。

多色(multicolor)の方は、実は歴史的経緯で、元々は、平面で上下左右の4点からしか影響を受けないような問題について考えた時に、
以下の図のように、黒色の点を計算する時には赤色からのみ影響を受け、赤色は黒色からのみ影響を受ける、と点を赤と黒の2グループに分けたことに起因します（この場合だけを赤黒法、Red-Black GaussSeidelとも呼ばれます）

この多色順序付けですが、単に計算順序を変えるだけだと、飛び飛びの行を計算することになり、キャッシュ効率等が悪くなってしまいます。
そこで、実際には、計算順序に沿って行列・ベクトルそのものを並び変えることが多いです。

下の図は、先ほどの行列を多色順序付けして並び替えたものです。

若干分かりづらいのですが、この行列は4色（グループ）に分けることができました。元々のサイズが64だったので、16行（点）ずつに分かれたことになります。
実際、先頭の色である0-15点（行）目は、同じ色の0-15点（列）目との係数はすべて0になっていることが分かります。次の16-31行目も同様です。
なので、この行列は、16行ずつなら並列化することが可能、ということになります。
16という絶対数だけ聞くと少なく聞こえるかもしれません。
しかし、問題サイズが大きくなっても、点と点の接続状況が変わらないかぎり色の数も同じなので、例えば32万点だった場合、8万点ずつ並列で計算できることになるので、計算時間がつらくなる問題サイズが大きいものでは十分に効果を発揮します。

多色順序付け法の欠点

ということで、ガウスザイデル法を多色順序付けを使うことで並列化することができました。めでたしめでたし。
と言いたいところなのですが、実は多色順序付けすると、元の逐次実行に比べて欠点があったりします。
それはベクトルデータのランダムアクセス化と、収束性の悪化です。

ベクトルへのランダムアクセス

元の行列と多色順序付けした後の行列の様子を並べて再掲します。

ここで、ガウスザイデル法は、基本的には疎行列・ベクトル積と同じく、係数行列\(a_{ij}\)に入出力ベクトル\(x_j\)をかける演算でした。
つまり、各行は、その行の非ゼロ（図でいうと色の付いている）部分に該当するベクトルxの値をアクセスします。

ということは、元の行列だと、例えば、0行目は、自分自身(0)と1,8,9列目のベクトルxの値が必要で、1行目は0,1,2,8,9,10列目が必要です。
という感じで見ていくと、元の行列のままだと、i行目の計算をする時には、基本的には中央のi-1,i,i+1列目と、そこから6つずつ左右に離れた3列ずつ、という3つの塊でアクセスしていることが分かります。
このように固まっている上、上から少しずつずれていっているので、隣の行を計算するときには既にキャッシュに乗っている可能性が高いです。
つまり、キャッシュには割と当たりやすい形をしています。

一方で、多色順序付けをして並び替えた方は、例えば0行目は0,16,32,48列目となり飛び飛びの場所のベクトルの値を必要とします。
1行目も同様に1,16,17,33,48,49と飛び飛びです。
これは当然といえば当然で、なぜなら多色順序付けは並列化のために格子の上で隣り合った（接続している）点を別のグループに分割するという手法だからです。
なので、今回の例のように4色に分かれていれば、基本的に必ず他の色のどこかにはアクセスする必要があるので、離れた4箇所にアクセスする必要が出てきます。
一般的に、n色に分けられた場合はベクトルはn箇所に分散します。
そのせいでキャッシュ効率が下がり、メモリアクセス速度がとても遅いなど、場合によっては並列化すると逆に遅くなる、ということも起きたりします。

というわけで、多色順序付けは並列化できる一方で、原理的にキャッシュ効率の悪化を避けられない、という問題があります。

収束性の悪化

キャッシュ効率が落ちるのは問題ですが、どれほど落ちるのかはやってみないと分からない上に、その対象とする問題や実効環境によっては別にあまり気にならなかったりします。

一方、こちらの収束性の悪化は結構悲惨です。 収束性の悪化も、最終的には、やはり元の問題（行列）の性質によるのですが、基本的にはどんな問題でも少しは悪化するはずです。

収束性が顕著に悪化する例として、ここでは対称ガウスザイデル法を例に挙げます。
通常のガウスザイデル法を見て気づいた方もいると思いますが、基本的には後に計算した値の方が前に計算した値より改善された（より解に近い）値を得ることができます。
ということは、0行目の結果とn-1行目の結果では、明らかに後者だけがとても解に近く、前者は離れたところにいるという非対称性を生むことになります。

非対称になると困る場合はいくつかあり、これを解説するとそれだけで記事1つできてしまうので詳しくは述べませんが。
ただ、多くの場合は、ガウスザイデル法をなんらかの他の反復法の前処理として使う場合に、前処理後に対称性が崩れてほしくないということに起因して対称性を残したいことが多いです
（※たまに、対称ガウスザイデル法はガウスザイデル法より収束性が高いと勘違いしている人も多いですが、それは嘘です。
対称ガウスザイデル法は↑のように2回ループが回るので1回のガウスザイデル法より収束性が高いように見えがちですが、実は素直に2回ガウスザイデル法を回すほうが収束性が高いです）。

このような非対称性を改善するものが以下の対称ガウスザイデル法です。

void SymmetricGaussSeidel(const double** a, double* x, const double* b, const std::size_t n)
{
    // forward sweep
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
    // backward sweep
    for(auto i = static_cast<std::make_signed_t<decltype(n)>>(n)-1; i > 0; i--)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
}

この対称ガウスザイデル法では、0行目からn-1行目まで順に計算した後、逆にn-1行目から0行目に向かって逆順にも計算する手法です。
このように逆回しも実行することで、どの値もちょうど同じぐらい改善されており、結果が対称になります。

と、対称ガウスザイデル法の話で少し脱線してしまいましたが、では実際に多色順序付けをするとどれほどガウスザイデル法の収束性が劣化するのかを示したのが下の図です。

まず、なにもしないガウスザイデル法（GaussSeidel、赤色）と多色順序付け（MulicoloOrdering-GS、黄緑色）はそれほど収束性に差がありません。
一方で、なにもしない対称ガウスザイデル法（SymGS、緑色）に対して多色順序付けした場合（MulicoloOrdering-SymGS、暗赤色）は大きく収束性が劣化して線が随分上に行ってしまっているのが分かります。
なお、ついでなのでガウスザイデル法を素直に2回まわした結果も乗せました（GaussSeidel x2、黄色）。
ご覧のとおり、SymGSよりこちらのほうが収束性が高いことが分かります。

今回の例ではこういう結果になりましたが、先述の通り実際には、多色順序付けすると普通のガウスザイデル法だけでも収束性がもっと劣化する場合もあります。
逆に、対称ガウスザイデル法でもあまり劣化しないこともあります。
が、いずれにせよ多色順序付けすると多かれ少なかれ収束性が劣化するという性質があることは間違いありません。

それは、多色順序付けで収束性が劣化する原因が、主にはincompatible nodeの増加によるものだからです。
incompatible nodeとは、「周囲の点/行がひとつも改善されていない状態から始める点/行」のことです。
詳細な説明等は中島先生の資料等に譲りますが、直感的には、なるべく周囲に改善された（解に近い）値があったほうが自分自身の解はより良いものになるので、その数が少ないと収束性が劣化する、ということです。

なぜincompatible nodeが増えるのかというのは、これも多色順序付けが何をしているのか考えれば当然で、同じグループの中の行は互いに独立して計算ができるようになっているので、
最初に計算する色（グループ）に属する行は全部周囲の点はまだ一度も計算されていないことになり、incompatible nodeに該当してしまうからです。
対して、何も並び替えず逐次で実行する場合は、基本的には最初の1行/点だけで、後続の処理は基本的には全部、周囲のどこかは既に改善された点になっているはずです。
つまり、n次元正方行列がc色に多色順序付けされたとすると、incompatible nodeの数は必ず1点からn/c点に増加します。

ということで、多色順序付けをすると、原理的に収束性を劣化させることにつながる（かもしれない）のです。
せっかく多色順序付けをして並列化できるようにしたのに、収束までに必要な反復回数が増えてしまってはあまり効果がありません。

ブロック化多色順序付け法

さて、通常の多色順序付け法ではこのように原理的に避けようがない問題が2点ありました。さてどうしましょうか。これを避ける方法がブロック化です。

ブロック化の説明のために、先に紹介した2色で塗り分けられる5点ステンシルの以下の例を挙げて説明します。

この場合、黒→赤の順で計算が実行されるとすると、この状態だと黒色が全てincompatible nodeになり合計で8個あります。
また、黒色が必要とする周囲の赤色の点は全て遠いところにあるのでキャッシュ効率も下がっているという先ほどの欠点があります。

しかし、これを、2×2の点でブロック化して見てください。つまり行番号のままで言うと

• 左下のブロック#0：0,8,10,2
• 右下のブロック#1：1,9,11,3
• 左上のブロック#2：4,12,14,16
• 右上のブロック#3：5,13,15,7

の4つのブロックに分けます（※ブロック≠色のグループです、気をつけてください）。
そうすると、ブロック単位でみると、左下のブロック#0は右上のブロック#3とは繋がっていないので独立して計算できることが分かります。また、右下のブロック#1と左上のブロック#2は独立しています。

と、いうことは、ブロックに分けた状態でも独立しているブロック同士は並列で計算できます（※ブロックの中の点同士は接続されているので、逐次実行する必要があります）。
もうちょっと言うと、ブロックに分けた後に、ブロック間の接続状況で多色順序付けすれば、ブロック同士は並列計算できるようになります。これが、ブロック化多色順序付けです。

ブロック化の効果

ブロック化すると、先に上げた多色順序付けの問題が緩和されます。

先の例だと、まずincompatible nodeはブロック化することによって8から4に半減します。
なぜなら、ブロックの中は逐次実行されるので、基本的にはincompatible nodeは「最初に計算される色に属するブロックの中の最初の点」だけになるからです。
例えばどのブロックも左下から計算するのであれば、incompatible nodeは、以下の4点だけになります。

• ブロック#0の点0
• ブロック#1の点1
• ブロック#2の点4
• ブロック#3の点5

実際に収束性が改善されている様子は、先に通常の多色順序付け法の収束性を論じていた時と同じ対象にブロック化を適用して実験してみた以下の結果を見ればよく分かると思います。

ブロック化していない普通の多色順序付け（MulticolorOrdering-SymGS、暗赤色）に比べて、ブロック化多色順序付け（BlockMulticolorOrdering-SymGS、水色・暗緑色）は、
理想的な逐次板対称ガウスザイデル法（SymGS、明緑色）にぐっと近くなっている様子がわかります。
また、2×2の4点ブロック（水色）にするより、ブロックのサイズ（1ブロックに含まれる点の数）を4×4に増やした方（暗緑色）が、より効果的であることが分かります。

ブロック化すると行列はどうなっているかというとこんな感じ↓になっています。

この例は2×2（つまり4点）でブロック化したもので、ブロックごと4色に分割されています。
この行列の様子を見ると分かる通り、上から順に実行したとして、周囲にまったく改善された点が含まれていない点(incompatible node)に該当するのは、0,4,8,12行目だけです。
なにもしない多色順序付けだと最初の色（0-31行目）が全てincompatible nodeで32点あるのに対して、4点まで減らすことに成功しています。
そしてまた、ブロックの中は接続されている上、他のブロックの点を必要としないor必要とする色の数が少ない点があるので、キャッシュ効率も上げられていることが分かります。

一般的に、n次元正方行列をブロックサイズbにブロック化してc色に分けられた場合、incompatible nodeはn/(bc)にできます。なので、可能な限りブロックサイズは大きくしたほうが得策ですね！

代数的ブロック化多色順序付け法

さて、ブロック化できましためでたしめでたし・・・と行きたいところなのですが、
↑の方法は、格子の状態を見て説明をした通り、基本的は行列を生成する元の問題の形状に依存してブロックに分類します。
そのため、幾何的ブロック化多色順序付け(Geometic Block Multicolor Ordering, GBMC)と呼ばれます。

でも、元の問題なんてわからないことも多いですし、元の問題によってブロック化の方法を変えるのって面倒ですよね。

そこで、行列の接続状況からだけでブロック化する方法があります。
それが代数的ブロック化多色順序付け(Algebraic Block Multicolor Ordering, ABMC)です。
今回公開したものは、このABMCを実装したものになります。

ABMCの方法については、先述の岩下ら(2009)の論文に詳しく書かれています。
英語でよければ、同じ著者による、IEEE Trans.の論文や
スライド資料もあります（発表年順からすると、IEEE Trans.の方が元のようです）。
詳しくはこれらの論文等を読んでもらうとして、おおざっぱな方法は以下のようになっています。

1. まだどのブロックにも入れていない行（点）iを探す
2. 点iを、新しいブロックbに入れる。
3. 点iに接続している列（点）のうち、まだどのブロックにも入れていない列（点）jを探す
4. 点jをブロックbに入れる
5. ブロックbが一杯になったら最初から繰り返す
6. まだ空きがあるなら、今度は点jに接続している点kを探す
7. 点kを新しいブロックbに入れる
8. ブロックbが一杯になったら最初から繰り返す
9. 以下、今度は点kに接続している列（点）を探してまたbに入れて・・・を繰り返す
10. 全点がブロックに入ったら、今度はブロック間の接続状況を確認して、普通の多色順序付けをすれば、完成

今回公開したコードでは、Abmc.hppの321行目からの関数がこの部分に該当します。
文章で説明されるよりコードを見たほうが分かりやすい方はぜひそちらをご覧ください。

まとめ

ということで、ガウスザイデル法の復習から、代数的ブロック化多色順序付け(ABMC)法までを簡単ですが解説しました。

この法に限らず、論文に手法自体は載っていても、実際の実装は公開されていることはほぼ稀で、仕方ないので自分で実装する時に苦労することがよくあります（特に数値計算分野）。
今回このコードを公開することで、後に続く方々の参考になれば、とても嬉しく思います。

なお、今回のコードは、これ自体は特に高速化等をあまり考えて作ってありません。
ひとつは先述のように「他の方の参考に」という意図を含めて分かりやすさを重視しているのと、またそもそも行列の前処理部分は実際の計算部分より小さく無視して良いことが多いから、という理由からです。
ただ、とても挑戦しがいのある問題だとは思いますので、またそのうち機会があればやってみた結果をご報告させていただくことができるかもしれません（アイディアはあるので）。
もし、並列化実装してみたとか、あるいはそれ以外でもバグ等見つけられた方は、遠慮無くissueやpull requestを投げてください。
お待ちしています。

あと最後に一応宣伝です。
今回のような「アルゴリズムそのものを見直す」という部分も、フィックスターズの得意とする「高速化」に含まれており、アルゴリズムからハードウェア特化したものまで、総合的なな解決方法を提供しています。
もし今現在何か「これを速くしてみたいんだけど」とお困りのプロジェクトがありましたら、ぜひお問い合わせください。
また、そういった仕事に興味がある方は、ぜひ一緒に働いてみませんか？