Gram–Schmidt の正規直交化法

一年生の線形代数で学習したGram–Schmidt (GS) の正規直交化法のプログララムを書いていきましょう．

GS の正規直交化法とは，内積の定義されたベクトル空間 $V$ において，ベクトルの集合 $\{a_1,a_2,\dots,a_n \}$ が与えられたとき，これらの貼る部分空間の正規直交基底を具体的に構成する手続きです．この講義では，簡単のため $V =\mathbb{R}^m$ とし，その内積は $a,b\in\mathbb{R}^m$ に対して $(a,b) := a^\top b$ であるとします（$a^\top$ は $a$ の転置を表します：${}^t a$ と同じ意味です）．基本的なアルゴリズムは次の通りです（プログラムを書くことを意識し，$=$ は右辺で左辺の変数を定義する（あるいは再定義する）という意味です）

    $q_1 = a_1$

    $\displaystyle q_1 = \frac{q_1}{\|q_1\|}$

    $q_2= a_2 - (q_1,a_2)q_1$

    $\displaystyle q_2= \frac{q_2}{\|q_2\|}$

    $q_3= a_3 - (q_1,a_3)q_1 - (q_2,a_3)q_2$

    $\displaystyle q_3= \frac{q_3}{\|q_3\|}$

    $\vdots$

    $\displaystyle q_n= a_n - \sum_{j=1}^{n-1} (q_j,a_n)q_j$

    $\displaystyle q_n = \frac{q_n}{\|q_n\|}$

以上を擬似コード的に書けば次のようになるでしょう．

for i = 1:n
    q_i = a_i
    for j = 1:i-1
        q_i -= (q_j,a_i) * q_j
    end
    q_i = q_i / norm(q_i)
end

あとは，Julia の文法に沿って，ベクトルの集合を入力として受け取り，ベクトルの集合を出力する関数を定義すればよいわけです．しかし，ここでいくつか注意すべきことがあります．

まず，内積やノルムはどうすればよいのでしょうか？実は，LinearAlgebraというパッケージには，内積やノルムを計算する関数が既に用意されています．内積は dot で計算でき，ノルムは norm で計算できます．ノルムの項を読むと分かるように，ユークリッドノルム以外のノルムにも対応しています．

ではこれで準備万端かというと，実はまだ考えなければならないことがあります．入力したベクトルの列が線形従属だとすると何が起きるでしょうか？そのような場合，どこかの i で q_i が零ベクトルとなってしまい，このときノルムは $0$ ですから 零割り が発生してしまいます．このような場合の対応策はいくつか考えられますが（GS の正規直交化法を使って「何をしたいのか」に強く依存します），少なくとも 「q_i が零ベクトルかどうかの判定」は必要そうです．どのように判定すればよいでしょうか？ベクトルの全ての成分が零かどうかを判定してもよいですが，（零ベクトルでなければ）その後すぐに正規化することを考えると，norm(q_i) が 0 かどうか判定すればよさそうです．では，if norm(q_i) == 0 のような条件分岐を書けばよいのでしょうか？実は，このような条件分岐は「数学としては正しい」ですが，コンピュータで計算する際には丸め誤差が避けられないことを踏まえると（数学的には 0 でも，実際の計算では丸め誤差分のズレが生じます），倍精度で計算する場合には例えば tol = 1e-10 のように閾値を設定して，if norm(q_i) < tol のように条件分岐を書けばよいということになります（1e-10 は $10^{-10}$ のことです ）．

サンプルファイルはこちら