northclimbの備忘録

徒然なるままに...です。

いっそのこと JavaScript の Array で行列を計算できるようにしてしまおう

初めに

稿は配列を配列のまま行列計算できるようにしたいという謎の考えの下で超絶適当な実装をしようとした筆者を供養するためのものです。ライブラリ使えよっていう

最新のコードはページ最下部のまとめに置いてます。

JavaScript を学び始めて1週間も経ってないのでエラーの出し方とかが気になる方もいるとは思いますが、そこは適宜変更を加えてください。

少し考えればLU分解とかも書けるのでしょうが疲れてしまったので気が向いたときに追記します。(たぶんそっちの方が重要なんだろうけど...)何か問題点を見つけた人はお手数ですがコメントまでお願いします。

関数群

配列(など)をJSONで厳密に比較

function JSON_equal(x, y) {
    return JSON.stringify(x) === JSON.stringify(y);
}

例えば、
JSON_equal([[1,2],[1,2]], [[1,2],[1,2]])true
JSON_equal([[1,2],[1,2]], [[1,2],[1,3]])false
を返します。

今回は配列を比較するために使いますが、それ以外にも適応可能で結構万能です。たぶん、null や undifined もきちんと区別します。

行列として扱えるか確認する

function is_matrix(array) {
    let row = array.length;
    let col = array[0].length;
    if (col === undefined) { return false }
    for (i in array) {
        if (array[i].length !== col) { return false }
    }
    return [row, col];
}

各行の要素数を比較し、一致しなければfalseを返します。要は、[[1,2], [3,4,5]]などを除外するために使います。ちなみに、[1,2,3]も除外します。行ベクトルを表現する場合は[[1,2,3]]です。

返り値は[行数, 列数]の配列です。

追記(今気づいた): 無駄に1回forを回しているので、気になる人は直しといてください。

m*nの零行列

function zeros(row, col) {
    if ((row < 1) || (col < 1)) { throw new Error('IndexError'); }
    let zeros = [];
    for (let i = 0; i < row; i++) {
        let tmp = [];
        for (let j = 0; j < col; j++) {
            tmp.push(0);
        }
        zeros.push(tmp);
    }
    return zeros;
}

その名の通りzeros(m,n)で m*n の零行列を返します。何のひねりもない素朴な実装です。

m*n の単位行列

function eye(row, col) {
    if ((row < 1) || (col < 1)) { throw new Error('IndexError'); }
    let eye = [];
    for (let i = 0; i < row; i++) {
        let tmp = [];
        for (let j = 0; j < col; j++) {
            tmp.push(Number(i === j));
        }
        eye.push(tmp);
    }
    return eye;
}

こちらも、その名の通りeye(m,n)で m*n の単位行列を返します。何のひねりもない(ry

対角行列

function diag(...c) {
    let diag = [];
    for (let i = 0; i < c.length; i++) {
        let tmp = [];
        for (let j = 0; j < c.length; j++) {
            if (i !== j) {
                tmp.push(0);
            } else {
                tmp.push(c[i]);
            }
        }
        diag.push(tmp);
    }
    return diag;
}

対角行列を求めます。例えば、diag(1,2)[[1,0],[0,2]]を返します。diag(1)もいけます。もちろん、[[1]]が返ってきます。

行列の和

function add_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('row or column does not match.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] + B[i][j];
        }
    }
    return tmp;
}

本当は X + Y みたいに書けるとベストなんだろうけど、そうすると Matrix オブジェクトを作って(たぶん)、演算子オーバーロードをして...とそれをやりたくないから Array のまま計算できるようにしたかったのだ!と本末転倒なので関数での実装に留めます。

行列の差

function sub_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('row or column does not match.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] - B[i][j];
        }
    }
    return tmp;
}

行列の差です。(見出しにそう書いてあるので、そうであることは当たり前)

和のときと同じです。特に言うことはありません。

行列のスカラー

function scalar_mul(A, x){
    let mn = is_matrix(A);
    if (mn) { throw new Error('Array is not matrixlike.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] * x;
        }
    }
    return tmp;
}

scalar_mul(A, x) x\cdot A を求めます。逆なのはあまり気にしないことにしましょう。スカラー倍なので。

行列の転置

function T(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    let tmp = zeros(mn[1], mn[0]);
    for (i in A) {
        for (j in A[i]) {
            tmp[j][i] = A[i][j];
        }
    }
    return tmp;
}

T(A)で行列 A の転置を求めます。

行列の乗算

function mul_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (Amn[1] !== Bmn[0]) { throw new Error('multiplication is not defined.') }
    let tmp = zeros(Amn[0], Bmn[1]);
    for (i in tmp) {
        for (j in tmp[i]) {
            for (let k = 0; k < Amn[1]; k++) {
                tmp[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return tmp;
}

ここまで来てやっと一般的には  AB \neq BA であることを確認できますね。

対角成分の抽出

function extract_diag(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = []
    for (i in A) {
        tmp.push(A[i][i]);
    }
    return diag(...tmp);
}

正方行列 A の対角成分を抽出した正方行列を返します。
例えば、extract_diag([[1,2],[3,4]])[[1,0],[0,4]]になります。

下三角行列

function lower(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = zeros(...mn);
    for (i in A) {
        for (let j = 0; j <= i; j++) {
            tmp[i][j] = A[i][j];
        }
    }
    return tmp;
}

正方行列 A の下三角行列を返します。
例えば、lower([[1,2],[3,4]])[[1,0],[3,4]]になります。

上三角行列

function upper(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = zeros(...mn);
    for (i in A) {
        for (let j = i; j < mn[0]; j++) {
            tmp[i][j] = A[i][j];
        }
    }
    return tmp;
}

正方行列 A の上三角行列を返します。
例えば、upper([[1,2],[3,4]])[[1,2],[0,4]]になります。
だんだん冗長になってきましたね?

行列の対角和(トレース)

function tr(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = 0;
    for (i in A) {
        tmp += A[i][i];
    }
    return tmp;
}

行列 A の対角和(トレース)を返します。
例えば、tr([[1,2],[3,4]])なら1+4で5になります。

行列の内積

function inner_prod(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('inner product is not defined.') }
    let tmp = 0;
    for (i in A) {
        for (j in A) {
            tmp += A[i][j] * B[i][j]
        }
    }
    return tmp;
}

読んで字のごとく内積〈A, B〉を求めます。

まとめ

上記のコードをまとめておきます。見てもよくわからない人は取り敢えず以下をコピペすれば上で紹介した関数をすべて使えるようになります。使用者がこのページへのリンクを布教すると、筆者は泣いて喜びます。

//配列(など)をJSONで厳密に比較
function JSON_equal(x, y) {
    return JSON.stringify(x) === JSON.stringify(y);
}

//行列として扱えるか?扱えないならfalse, 扱えるなら[row, col]を返す。
function is_matrix(array) {
    let row = array.length;
    let col = array[0].length;
    if (col === undefined) { return false }
    for (i in array) {
        if (array[i].length !== col) { return false }
    }
    return [row, col];
}
//零行列
function zeros(row, col) {
    if ((row < 1) || (col < 1)) { throw new Error('IndexError'); }
    let zeros = [];
    for (let i = 0; i < row; i++) {
        let tmp = [];
        for (let j = 0; j < col; j++) {
            tmp.push(0);
        }
        zeros.push(tmp);
    }
    return zeros;
}
//単位行列
function eye(row, col) {
    if ((row < 1) || (col < 1)) { throw new Error('IndexError'); }
    let eye = [];
    for (let i = 0; i < row; i++) {
        let tmp = [];
        for (let j = 0; j < col; j++) {
            tmp.push(Number(i === j));
        }
        eye.push(tmp);
    }
    return eye;
}
//対角行列
function diag(...c) {
    let diag = [];
    for (let i = 0; i < c.length; i++) {
        let tmp = [];
        for (let j = 0; j < c.length; j++) {
            if (i !== j) {
                tmp.push(0);
            } else {
                tmp.push(c[i]);
            }
        }
        diag.push(tmp);
    }
    return diag;
}
//行列の和 A+B
function add_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('row or column does not match.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] + B[i][j];
        }
    }
    return tmp;
}
//行列の差 A-B
function sub_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('row or column does not match.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] - B[i][j];
        }
    }
    return tmp;
}
//行列のスカラー倍 xA
function scalar_mul(A, x) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.'); }
    let tmp = zeros(...mn);
    for (i in A) {
        for (j in A[i]) {
            tmp[i][j] = A[i][j] * x;
        }
    }
    return tmp;
}
//行列の転置 A^t
function T(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    let tmp = zeros(mn[1], mn[0]);
    for (i in A) {
        for (j in A[i]) {
            tmp[j][i] = A[i][j];
        }
    }
    return tmp;
}
//行列の乗算 AB
function mul_asmatrix(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (Amn[1] !== Bmn[0]) { throw new Error('multiplication is not defined.') }
    let tmp = zeros(Amn[0], Bmn[1]);
    for (i in tmp) {
        for (j in tmp[i]) {
            for (let k = 0; k < Amn[1]; k++) {
                tmp[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return tmp;
}
//行列の対角成分を抽出する
function extract_diag(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = []
    for (i in A) {
        tmp.push(A[i][i]);
    }
    return diag(...tmp);
}
//行列Aの下三角行列
function lower(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = zeros(...mn);
    for (i in A) {
        for (let j = 0; j <= i; j++) {
            tmp[i][j] = A[i][j];
        }
    }
    return tmp;
}
//行列Aの上三角行列
function upper(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = zeros(...mn);
    for (i in A) {
        for (let j = i; j < mn[0]; j++) {
            tmp[i][j] = A[i][j];
        }
    }
    return tmp;
}
//行列のトレース(対角和)
function tr(A) {
    let mn = is_matrix(A);
    if (mn === false) { throw new Error('Array is not matrixlike.') }
    if (mn[0] !== mn[1]) { throw new Error('square matrix only')}
    let tmp = 0;
    for (i in A) {
        tmp += A[i][i];
    }
    return tmp;
}
//行列の内積 <A,B>
function inner_prod(A, B) {
    let Amn = is_matrix(A); Bmn = is_matrix(B);
    if (Amn === false || Bmn === false) { throw new Error('Array is not matrixlike.'); }
    if (!JSON_equal(Amn, Bmn)) { throw new Error('inner product is not defined.') }
    let tmp = 0;
    for (i in A) {
        for (j in A) {
            tmp += A[i][j] * B[i][j]
        }
    }
    return tmp;
}