LoginSignup
7
1

More than 5 years have passed since last update.

WebGLAdvent Calendar 2017Day 3
@nyamadandan

WebGLでガウスの消去法を解いてみた

Last updated at Posted at 2017-12-03

はじめに

ガウスの消去法は掃き出し法などとも呼ばれ、連立一次方程式を解くための多項式アルゴリズムです。
ガウスの消去法を使用すれば逆行列などを求めることもできるのですが、こちらをWebGL 1.0を用いたGPGPUを使用して解いてみます。

ガウスの消去法を使用した逆行列の求め方

ガウスの消去法を使用した逆行列ですが逆行列を求めたい行列mat1と単位行列mat2を並べて書きます。

mat1 = \left(
\begin{array}{cccc}
4.000 & 0.000 & -1.000 & 2.000 \\
1.000 & 0.000 & 0.000 & 0.000 \\
2.000 & 0.000 & 0.000 & 1.000 \\
-7.000 & 1.000 & 1.000 & -3.000 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
1.000 & 0.000 & 0.000 & 0.000 \\
0.000 & 1.000 & 0.000 & 0.000 \\
0.000 & 0.000 & 1.000 & 0.000 \\
0.000 & 0.000 & 0.000 & 1.000 
\end{array}
\right)

上記２行列に対して以下の操作を加えてmat1を単位行列へと変換していきます。

  • ある行と行を入れ替える
  • ある行をn倍する
  • ある行をn倍した結果を別の行に加える

例えば１行目と４行目を入れ替える操作を行った場合はmat1及びmat2は次のように変換されます。

mat1 = \left(
\begin{array}{cccc}
-7.000 & 1.000 & 1.000 & -3.000 \\
1.000 & 0.000 & 0.000 & 0.000 \\
2.000 & 0.000 & 0.000 & 1.000 \\
4.000 & 0.000 & -1.000 & 2.000 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & 0.000 & 0.000 & 1.000 \\
0.000 & 1.000 & 0.000 & 0.000 \\
0.000 & 0.000 & 1.000 & 0.000 \\
1.000 & 0.000 & 0.000 & 0.000 
\end{array}
\right)

実は機械的にmat1を単位行列に変換していく方法があります。
まずは対角成分である１行１列に着目し、１行１列の成分の絶対値が最も大きくなるように入れ替えます。
ちょうど先程の１行目と４行目の交換がそれに当たります。
次に１列目の対角成分より下方向を０にするように１行目を加算していきます。
例えば２行目であれば１行目の1/7倍を加算し、３行目であれば2/7倍を加算すればよいわけです。
この操作を繰り返すと１列目は次のように整理され、合わせてmat2も変形されていきます。


mat1 = \left(
\begin{array}{cccc}
-7.000 & 1.000 & 1.000 & -3.000 \\
0.000 & 0.143 & 0.143 & -0.429 \\
0.000 & 0.286 & 0.286 & 0.143 \\
0.000 & 0.571 & -0.429 & 0.286 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & 0.000 & 0.000 & 1.000 \\
0.000 & 1.000 & 0.000 & 0.143 \\
0.000 & 0.000 & 1.000 & 0.286 \\
1.000 & 0.000 & 0.000 & 0.571 
\end{array}
\right)

まとめると以下の操作を繰り返すことになります

  • 対角成分より下方向に最も絶対値が大きくなる成分を検索する
  • 対角成分が最も大きくなるように行を入れ替える
  • 注目している対角成分を含む行を、対角成分の下方向が0になるように加算する

処理は１行１列、２行２列、３行３列、４行４列のように対角成分に注目して進めていきます。
３行３列まで処理を繰り返すと次のようになります。

mat1 = \left(
\begin{array}{cccc}
-7.000 & 1.000 & 1.000 & -3.000 \\
0.000 & 0.571 & -0.429 & 0.286 \\
0.000 & 0.000 & 0.500 & 0.000 \\
0.000 & 0.000 & 0.000 & -0.500 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & 0.000 & 0.000 & 1.000 \\
1.000 & 0.000 & 0.000 & 0.571 \\
-0.500 & 0.000 & 1.000 & 0.000 \\
0.000 & 1.000 & -0.500 & 0.000 
\end{array}
\right)

次は逆方向に４行４列成分から上方向に加算を繰り返し、対角成分よりも上方向を0にするように対角成分を含む行を加算します。
たとえば、２行目に着目した時、４行目の$0.286/0.500$倍を加算すると２行４列成分を０とすることができます。同様に４行目を$-0.300/0.500$倍したものを１行目に加算します。


mat1 = \left(
\begin{array}{cccc}
-7.000 & 1.000 & 1.000 & 0.000 \\
0.000 & 0.571 & -0.429 & 0.000 \\
0.000 & 0.000 & 0.500 & 0.000 \\
0.000 & 0.000 & 0.000 & -0.500 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & -6.000 & 3.000 & 1.000 \\
1.000 & 0.571 & -0.286 & 0.571 \\
-0.500 & 0.000 & 1.000 & 0.000 \\
0.000 & 1.000 & -0.500 & 0.000 
\end{array}
\right)

その後、４行目を対角成分の-0.500で割り、対角成分を１に正規化します。

mat1 = \left(
\begin{array}{cccc}
-7.000 & 1.000 & 1.000 & 0.000 \\
0.000 & 0.571 & -0.429 & 0.000 \\
0.000 & 0.000 & 0.500 & 0.000 \\
0.000 & 0.000 & 0.000 & 1.000 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & -6.000 & 3.000 & 1.000 \\
1.000 & 0.571 & -0.286 & 0.571 \\
-0.500 & 0.000 & 1.000 & 0.000 \\
0.000 & -2.000 & 1.000 & 0.000 
\end{array}
\right)

まとめると

  • 対角成分の上の行に対角成分の行を加算し、対角成分よりも上方向の成分を０にする
  • 対角成分を１に正規化する

これを３行３列、２行２列、１行１列の順で処理します。するとmat1は単位行列となりmat2に逆行列が求まります。

mat1 = \left(
\begin{array}{cccc}
1.000 & 0.000 & 0.000 & 0.000 \\
0.000 & 1.000 & 0.000 & 0.000 \\
0.000 & 0.000 & 1.000 & 0.000 \\
0.000 & 0.000 & 0.000 & 1.000 
\end{array}
\right)
,mat2 = \left(
\begin{array}{cccc}
0.000 & 1.000 & 0.000 & 0.000 \\
1.000 & 1.000 & 1.000 & 1.000 \\
-1.000 & 0.000 & 2.000 & 0.000 \\
0.000 & -2.000 & 1.000 & 0.000 
\end{array}
\right)

逆行列が求まる過程のデモ

WebGLで実装するには

まずは浮動小数点テクスチャに行列を書き込みます。
その際にr成分には行列の成分を格納し、gb成分にはUV座標を格納します。

絶対値が最大となる対角成分の検索ですが、列の長さの半分のテクスチャに書き込みを繰り返すことにより検索します。
４行４列の場合はその半分のサイズである1x2サイズのテクスチャを作成し、それに書き込む事になります。
例えば先程の行列だと１列目は$(4, 1, 2, -7)$ですがこれを半分のサイズのテクスチャに書き込む場合、隣接する4と1を比較して絶対値の大きい4を書き込み、同様に隣接する2と-7を比較して-7を書き込みます。
その結果$(4, 1, 2, -7)$は1x2テクスチャに$(4, -7)$が書き込まれることとなります。
同じ処理をもう一度繰り返します。
書き込むテクスチャは1x2テクスチャの更に半分の1x1テクスチャとなり、4と-7の比較が行われた結果、1x1テクスチャには-7とそのuv座標が書き込まれることになります。
この際、上方向に検索をしないようにoffsetよりも上方向の成分は選択されないようにしています。

precision highp float;
uniform sampler2D dataTexture;
uniform vec2 resolution;
uniform vec2 offset;

void main()
{
    vec2 texelSize = 0.5 / resolution;
    vec2 uv = vec2(offset.x, (gl_FragCoord.y - 0.25) / resolution.y);

    vec2 u0 = uv;
    vec2 u1 = uv + vec2(0.0, 1.0) * texelSize;
    vec4 v0 = texture2D(dataTexture, u0);
    vec4 v1 = texture2D(dataTexture, u1);

    vec4 v;

    // 対角成分よりも上なのでv1を選択
    if(offset.y > u0.y) {
        v = v1;
    } else {
        // 絶対値が大きい方を選択
        if(abs(v1.r) < abs(v0.r)) {
            v = v0;
        } else {
            v = v1;
        }
    }

    gl_FragColor = v;
}

こうして求められた絶対値が最大となる成分を行の交換シェーダに渡して交換を行います。

precision highp float;
uniform sampler2D dataTexture;
uniform sampler2D pivotTexture;
uniform vec2 resolution;
uniform float v0;

void main()
{
    vec2 uv = gl_FragCoord.xy / resolution;
    float v1 = texture2D(pivotTexture, vec2(0.5, 0.5)).b; // 先程検索した行のv座標
    vec4 color = texture2D(dataTexture, uv);

    // 対象の行ならば値の交換を行う
    if(abs(v0 - uv.y) < 0.0001) {
        color.r = texture2D(dataTexture, vec2(uv.x, v1)).r;
    } else if(abs(v1 - uv.y) < 0.0001) {
        color.r = texture2D(dataTexture, vec2(uv.x, v0)).r;
    }

    gl_FragColor = color;
}

交換が終わったら次に対角成分の下方向の成分を0とします。
GLSLを使用すれば1パスで１列目の対角成分以下を0にできます。

precision highp float;
uniform sampler2D dataTexture;
uniform vec2 resolution;
uniform vec2 pivot;

void main()
{
    vec2 uv = gl_FragCoord.xy / resolution;
    float p = texture2D(dataTexture, pivot).r;
    vec4 color = texture2D(dataTexture, uv);

    if(pivot.y + 0.5 / resolution.y < uv.y) {
        float r0 = texture2D(dataTexture, vec2(pivot.x, uv.y)).r;
        float r1 = texture2D(dataTexture, vec2(uv.x, pivot.y)).r;
        // 対角成分の行をr0 / p倍したものを引く
        color.r -= r0 * r1 / p;
    }

    gl_FragColor = color;
}

そして最後に対角成分よりを上方向を0にするとともに対角成分を1で正規化します。

precision highp float;
uniform sampler2D dataTexture;
uniform vec2 resolution;
uniform vec2 pivot;

void main()
{
    vec2 uv = gl_FragCoord.xy / resolution;
    float p = texture2D(dataTexture, pivot).r;
    vec4 color = texture2D(dataTexture, uv);

    if(abs(pivot.y - uv.y) < 0.0001) {
        // 対角成分を1にする
        color.r = color.r / p;
    } else if(pivot.y - 0.5 / resolution.y > uv.y) {
        float r0 = texture2D(dataTexture, vec2(pivot.x, uv.y)).r;
        float r1 = texture2D(dataTexture, vec2(uv.x, pivot.y)).r;
        // 対角成分の行をr0 / p倍したものを引く
        color.r -= r0 * r1 / p;
    }

    gl_FragColor = color;
}

これでWebGLによって逆行列を求めることができました。
最後にこの数値をJavaScript側から読み出したいのですがglReadPixelsは浮動小数点テクスチャには対応していないのでそのままでは読み出すことができません。
そのため、こちらのシェーダを使用して浮動小数点テクスチャからRGBA24の整数テクスチャに変換します。

実行結果のデモ

おわりに

今回はWebGLを使用して逆行列を求めてみました。気になるパフォーマンスですが今回は実装にThree.jsを使用したのもあり十分な最適化ができていませんのであまり期待はできません。
行列のサイズも今回のような小さなものではあまり効果は期待はできないと思います。

ですが絶対値の検索や浮動小数点テクスチャの読み出しはWebGLでGPGPUをする上では非常に便利なテクニックだと思うので今後活用できればいいなと思います。

7
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
7
1