『Python数値計算ノート』ではアフィリエイトプログラムを利用して商品を紹介しています。

グラム・シュミットの直交化法

≪【前の記事】完全正規直交系

グラム・シュミットの直交化法

前回記事では完全正規直交系の定義と性質について解説しました。今回は Rn の任意の基底 {a1, a2, , an} から完全正規直交系 {q1, q2, , qn} をつくるグラム・シュミットの直交化法を学びます。

直交化の手順

R3 を張る任意の基底 {a1, a2, a3} が与えられたとします。
この基底を使って完全正規直交系 {q1, q2, q3} をつくってみます。
最初に基底から適当なベクトルを 1 つ選んで u1 とします。
選び方は自由ですが、とりあえず a1 を選んで u1=a1 としましょう。
Gram–Schmidt orthonormalization
自身の大きさで割って ノルム1 にしたベクトルを q1 とします。
Gram–Schmidt の直交化2 u_1の正規化
次に a2 を選んで、q1 への 射影ベクトル p1 を求めます。
 (1)p1=q1Ta2q1Tq1q1
q1 のノルムは 1 なので、q1Tq1=1 となります。したがって、
 (2)p1=(q1Ta2)q1
a2 から p1 を差し引くと、q1 に垂直なベクトルが得られます。これを u2 とします。
 (3)u2=a2p1=a2(q1Ta2)q1
Gram–Schmidt の直交化3 a_2の射影ベクトル
これも u2 で割って、ノルム 1 のベクトル q2 をつくります。
Gram–Schmidt の直交化4 u2を単位ベクトルq2に変換
次は q1q2 の張る平面に対して a3 の射影ベクトル
 (4)p2=(q1Ta3)q1+(q2Ta3)q2
をつくって a3 から差し引くと、q1q2 に垂直なベクトル u3 を得ます。
 (5)u3=a3(q1Ta3)q1(q2Ta3)q2
Gram–Schmidt の直交化5 a3の正射影をつくる
u3 を自身の大きさで割って q3 を得ます。
グラム・シュミットの直交化6 完全正規直交系の完成
さらに次元の大きな空間の場合も同様にして、
 (6)uk=aki=1k1(qiTak)qi,qk=ukuk
によって順次 qk を求めていくことができます。新しい qk を計算するときには、それまで求めた {q1, q2, } をすべて使います。

Pythonで直交化法を実装する

グラム・シュミットの直交化法を使って、基底を直交化するコードです。

# python_gram_schmidt_orthonormalization

# In[1]

import numpy as np
from scipy import linalg

# グラム・シュミット直交化関数
def schmidt(arr):

    # 配列のデータ型を64バイト浮動小数点数型に変換
    arr = np.array(arr, dtype=np.float64)

    # 渡した配列の列数(基底に含まれるベクトル数)
    k = arr.shape[1]

    # 1列目のベクトルを選択
    u = arr[:,[0]]

    # uを正規化
    q = u / linalg.norm(u)

    # シュミットの直交化
    for j in range(1, k):
        u = arr[:,[j]]
        for i in range(j):
            u -= np.dot(q[:,i], arr[:,j]) * q[:,[i]]
        qi = u / linalg.norm(u)
        q = np.append(q, qi, axis=1)
    return q

schmidt() に適当な基底を渡して完全正規直交系をつくってみます。

# In[2]

A = np.array([[3, 1, 1],
              [1, 0, 1],
              [2, 5, 0]])

Q = schmidt(A)

print(Q)
# [[ 0.80178373 -0.47847438 -0.35805744]
#  [ 0.26726124 -0.24880668  0.93094934]
#  [ 0.53452248  0.8421149   0.07161149]]

出力結果が完全正規直交系の性質を満たしているか確認しておきましょう。各々のベクトルは互いに直交します。たとえば、1 列目のベクトルと 2 列目のベクトルの内積は 0 となっているはずです。

# In[3]

print(Q[:, 0] @ Q[:, 1])
# -5.551115123125783e-17

前回学んだように、直交行列 Q を適当なベクトルに作用させるとノルムが保存されます。

# In[4]

# ベクトルvを定義
v = np.array([[1],
              [2],
              [3]])

# Qvを計算
Qv = np.dot(Q, v)

# vとQvのノルム(大きさ)
v_norm = linalg.norm(v)
Qv_norm = linalg.norm(Qv)

print(v_norm)
print(Qv_norm)
# 3.7416573867739413
# 3.7416573867739418

≫【次の記事】QR分解
≫ Pythonで学ぶ線形代数トップページ

コメント

  1. HNaito より:

    下記は誤植と思われますので、ご確認ください。
    (3) 式の下の文章で、||u1|| で割って → ||u2||で割って
    In[3]のプログラムで、print(x[:, 0] @ x[:, 1]) → print(Q[:, 0] @ Q[:, 1])