1.経緯

いま、数学の学び直しのため、世界標準MIT教科書ストラング：線形代数イントロダクション（邦訳第4版）（※以下、同著と略）を読んでいまして、つい先ごろ第4章まで読み終わったところで、4.4節「直交基底とグラム-シュミット法」に出てきたQR 分解のアルゴリズムをTypeScriptで実装したので共有します。
与えられた行列をプログラムでQR分解するだけなら、Pythonを使って Numpy のnumpy.linalg.qr を使えばよいし、JavaScript であれば、math.js の math.qr を使えばよいのですが、同著に掲載されているアルゴリズムの理解を深めるため自作してみた、という趣旨です。
使用した言語としては、目下の主業務であるフロントエンド開発で使い慣れているTypeScriptを使うことにしました。

2.設計メモ

説明のためのサンプルとして、同著P.248〜249に掲載の以下の例を使用します。

　　　　独立ではあるが直交しないベクトル $ a $ , $ b $ ,$ c $ が次のものであったとする。

a = 
\begin{pmatrix}
1 \\
-1 \\
0
\end{pmatrix}
 \hspace{15pt}
b = 
\begin{pmatrix}
2 \\
0 \\
-2
\end{pmatrix}
 \hspace{15pt}
c = 
\begin{pmatrix}
3 \\
-3 \\
3
\end{pmatrix}

　　　　上記の $ a $ , $ b $ ,$ c $ を列ベクトルとする行列 $ A $ のQR分解（の一例）は、以下のようになる。

A = 
\begin{pmatrix}
1 & 2  & 3 \\
-1 & 0 & -3 \\
0 & -2 & 3
\end{pmatrix}
= 
\begin{pmatrix}
1/\sqrt{2} & 1/\sqrt{6}  & 1/\sqrt{3} \\
-1/\sqrt{2} & 1/\sqrt{6} & 1/\sqrt{3} \\
0 & -2/\sqrt{6} & 1/\sqrt{3}
\end{pmatrix}

\begin{pmatrix}
\sqrt{2} & \sqrt{2}  & \sqrt{18} \\
0 & \sqrt{6} & -\sqrt{6} \\
0 & 0 & \sqrt{3}
\end{pmatrix}
= QR

　　　　上記の分解において、 $ Q $ は直交行列であり、$ R $は上三角行列である。

同著では、様々な概念の説明のときに、行列を（行を縦に並べたものではなく）列を横に並べたものとして見ることをしばしば求められるので、

A = 
\begin{pmatrix}
1 & 2  & 3 \\
-1 & 0 & -3 \\
0 & -2 & 3
\end{pmatrix}

という行列をコードで書くときに、行を表す配列の配列によって、

const A: Array<Array<number>> = [
  [1, 2, 3],
  [-1, 0, -3],
  [0, -2, 3]
];

とは書かないで、以下のように、行列 $A$ を構成する3つの列ベクトルを表す、 Vectorクラスのインスタンスを並べたものとして、行列 $A$ を表現する。

const a: Vector = new Vector([1, -1, 0]);
const b: Vector = new Vector([2,  0, -2]);
const c: Vector = new Vector([3, -3, 3]);

const A: Vector[] = [a, b, c];

上記のような Vector クラスをTypeScriptで実装し、これを使って QR分解する関数も実装する。

3.実装

3.1 Vector クラス

とり急ぎ、こんなの書きました。

// Vector class
class Vector {
    // properties
    [index: number]: number;

    // constructor
    constructor(elements: number[]) {
        for (let i=0; i < elements.length; ++ i) {
            this[i] = elements[i];
        }
    }

    // iterator
    [Symbol.iterator](): Iterator<number> {
        const that = this;
        let index = 0;
        return {
            next(): IteratorResult<number> {
                const value = that[index];
                if (Number.isFinite(value)) {
                    index ++;
                    return { value, done: false };
                } else {
                    return { value: NaN, done: true };
                }
            },
        };
    }

    // the norm of this vector
    get norm(): number {
        const sumOfSquares = [...this].reduce((s: number, v: number) => s + Math.pow(v, 2) , 0);
        return Math.sqrt(sumOfSquares);
    }

    // return the dot product of this vector and the other vector
    dot(other: Vector): number {
        return [...this].reduce((s, v, i) => s + v * other[i], 0);
    }

    // return the addition vector of this vector and the other vector
    add(other: Vector): Vector {
        return new Vector([...this].map((v, i) => v + other[i]));
    }

    // return the difference vector to subtract the other vector from this vector
    subtract(other: Vector): Vector {
        return new Vector([...this].map((v, i) => v - other[i]));
    }

    // return the vector of specified scalar 'k' multiplication of this vector
    mul(k: number) {
        return new Vector([...this].map(v => k * v));
    }

    // stringify this vector
    toString(): string {
        return JSON.stringify([...this]);
    }
}

以下のような感じで使います。

const a: Vector = new Vector([1, -1, 0]);

console.log(`${a}`);  // => [1,-1,0]
console.log([...a]);  // => [ 1, -1, 0 ]
console.log(a.norm);  // => 1.4142135623730951

for (let x of a) {
    console.log(x);  // 一行ずつ、1, -1, 0 と表示される。
}

3.2 QR分解

同著P.250に掲載されている、下記のアルゴリズム

を TypeScript で書きました。

const qr = (A: Vector[]): [Vector[], Vector[]] => {
    const n: number = A.length;
    const Q: Vector[] = Array(n).fill(null);
    const R: Vector[] = [...Array(n)].map(_ => new Vector([...A[0]].fill(0)));

    for (let j = 0; j < n; ++ j) {
        let v = A[j];
        for (let i =0; i < j; ++ i) {
            R[j][i] = A[j].dot(Q[i]);
            v = v.subtract(Q[i].mul(R[j][i]));
        }
        R[j][j] = v.norm;
        Q[j] = v.mul(1 / R[j][j]);
    }

    return [Q, R];
}

なお同著掲載のコードでは、上三角行列 $R$ の i行j列成分を自然に R(i, j) と書けていますが、上記の自作コードだと、ゼロ始まりのインデクスjとi を使って $R$ のi行j列成分を(iとjが逆の)R[j][i] と書かなければならないところが、やや（いや、かなり?）気持ちワルい。けれども、同著のポリシーに従って、行列を（行の配列ではなく）列の配列として表しているので、そこは（今のところ）妥協。

4.テスト実行

先の例に挙げた、

A = 
\begin{pmatrix}
1 & 2  & 3 \\
-1 & 0 & -3 \\
0 & -2 & 3
\end{pmatrix}
= 
\begin{pmatrix}
1/\sqrt{2} & 1/\sqrt{6}  & 1/\sqrt{3} \\
-1/\sqrt{2} & 1/\sqrt{6} & 1/\sqrt{3} \\
0 & -2/\sqrt{6} & 1/\sqrt{3}
\end{pmatrix}

\begin{pmatrix}
\sqrt{2} & \sqrt{2}  & \sqrt{18} \\
0 & \sqrt{6} & -\sqrt{6} \\
0 & 0 & \sqrt{3}
\end{pmatrix}
= QR

を確認します。

const a: Vector = new Vector([1, -1, 0]);
const b: Vector = new Vector([2,  0, -2]);
const c: Vector = new Vector([3, -3, 3]);

const A: Vector[] = [a, b, c];

const [Q, R] = qr(A);

console.log(`A: ${A}`);
console.log(`Q: ${Q}`);
console.log(`R: ${R}`);

上記によって、以下のように出力されます。

A: [1,-1,0],[2,0,-2],[3,-3,3]
Q: [0.7071067811865475,-0.7071067811865475,0],[0.4082482904638632,0.40824829046386296,-0.8164965809277261],[0.5773502691896262,0.5773502691896255,0.5773502691896258]
R: [1.4142135623730951,0,0],[1.414213562373095,2.449489742783178,0],[4.242640687119285,-2.449489742783178,1.7320508075688772]

上記の Q と R が含む小数点以下を持つ各値は、先の $Q$ および $R$ の成分の近似値となっていることを確認できました。（※ 上記のコードのテスト結果では、A だけではなく、Q と R も列ベクトルの並びとして表示されます。）

5.備考

5.1 Numpyを使ってみると

Numpy のnumpy.linalg.qr を使ってみると、以下のように Q の第1列と第2列の符号（±）が異なり、それにあわせて、R の第1行と第2行の符号（±）が異なる結果になった。

$ python
Python 3.8.3 (default, Jul  2 2020, 11:26:31) 
[Clang 10.0.0 ] :: Anaconda, Inc. on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> A = np.array([
...      [1, 2, 3],
...      [-1, 0, -3],
...      [0, -2, 3]
... ])
>>> Q, R = np.linalg.qr(A)
>>> Q
array([[-0.70710678, -0.40824829,  0.57735027],
       [ 0.70710678, -0.40824829,  0.57735027],
       [-0.        ,  0.81649658,  0.57735027]])
>>> R
array([[-1.41421356, -1.41421356, -4.24264069],
       [ 0.        , -2.44948974,  2.44948974],
       [ 0.        ,  0.        ,  1.73205081]])
>>>

同著の P.250 の一番下に

告白　　MATLABやOctave, Python のような良質のシステムで使われるLAPACKのような良質のソフトウェアではこのグラムーシュミット法のプログラムは使われない。今ではより良い方法がある。

と書かれているから、アルゴリズムが異なるのだろう。

5.2 GitHubレポジトリ

今回作成したものは、github.com/jun68ykt/linalg-ts に上げておいた。同著を先に進めていくにあたって、また別のナントカ分解を学ぶ際に、理解の確認のためにコードを追加するかもしれない。

以上