はじめに

　『スタンフォード　ベクトル・行列からはじめる最適化数学』の学習ノートです。
　「数式の行間埋め」や「Pythonを使っての再現」によって理解を目指します。本と一緒に読んでください。

　この記事は5.4節「グラム・シュミット法」の内容です。
　2ベクトルの正規直交ベクトルを導出します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
グラム・シュミット法の導出
- 式の導出
参考書籍
おわりに

グラム・シュミット法の導出

　グラム・シュミット法(Gram-Schmidt algorithm)によるベクトルの正規直交化(orthogonalization)の計算を数式とグラフで確認します。
　正射影については「5.3：正射影ベクトルの導出【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照してください。

　次のライブラリを利用してグラフを作成します。不要であれば省略してください。

# 利用ライブラリ
import numpy as np
import matplotlib.pyplot as plt

式の導出

　グラム・シュミット法による正規直交化の計算式を導出します。2次元の場合を例として図示しますが、次元に関わらず成立します。

　ベクトル $\mathbf{a}_1 = (a_{1,1}, \cdots, a_{1,n})^{\top}$ と $\mathbf{a}_2 = (a_{2,1}, \cdots, a_{2,n})^{\top}$ を正規直交化したベクトル $\mathbf{q}_1 = (q_{1,1}, \cdots, q_{1,n})^{\top}$ と $\mathbf{q}_2 = (q_{2,1}, \cdots, q_{2,n})^{\top}$ を考えます。

　2つのベクトルとなす角の関係をグラフで確認します。

・作図コード(クリックで展開)

　2次元ベクトルを指定して、内積を計算します。

# ベクトルを指定
a1 = np.array([4.0, 1.0])
a2 = np.array([2.0, 3.0])

# ベクトルa1,a2のなす角を計算
dot_a1a2 = np.dot(a1, a2)
norm_a1 = np.linalg.norm(a1)
norm_a2 = np.linalg.norm(a2)
theta = np.arccos(dot_a1a2 / norm_a1 / norm_a2)
print(theta) # 0からπの値
print(theta / np.pi) # 0から1の値

0.7378150601204649
0.2348538278116319

　2つのベクトル $\mathbf{a}_1 = (a_{1,1}, a_{1,2})^{\top}, \mathbf{a}_2 = (a_{2,1}, a_{2,2})^{\top}$ をa1, a2として値を指定します。Pythonでは0からインデックスが割り当てられるので、 $a_{1,1}$ の値はa1[0]に対応します。

　ベクトル $\mathbf{a}_1, \mathbf{a}_2$ のなす角 $\theta = \arccos (\frac{\mathbf{a}_1^{\top} \mathbf{a}_2}{\|\mathbf{a}_1\| \|\mathbf{a}_2\|})$ を計算します。
　 $\arccos x$ は逆コサイン関数(コサイン関数 $\cos x$ の逆関数)(3.4節)、 $\|\mathbf{x}\|$ はユークリッドノルム(3.1節)でベクトルの大きさを表します。

　なす角を示す角マークと角ラベルを描画するための座標を作成します。

# ベクトルa1,a2とx軸のなす角を計算
sgn_a1 = 1.0 if a1[1] >= 0.0 else -1.0
theta_a1 = sgn_a1 * np.arccos(a1[0] / norm_a1)
sgn_a2 = 1.0 if a2[1] >= 0.0 else -1.0
theta_a2 = sgn_a2 * np.arccos(a2[0] / norm_a2)

# 差が180°を超える場合は、負の角度を優角に変換
if abs(theta_a1 - theta_a2) > np.pi:
    theta_a1 = theta_a1 if a1[1] >= 0.0 else theta_a1 + 2.0*np.pi
    theta_a2 = theta_a2 if a2[1] >= 0.0 else theta_a2 + 2.0*np.pi

# 角マーク用の値(ラジアン)を作成
rad_vals = np.linspace(start=theta_a1, stop=theta_a2, num=100)

# 角マークの座標を計算
r = 0.3
angle_arr = np.array(
    [r * np.cos(rad_vals), 
     r * np.sin(rad_vals)]
)
print(angle_arr[:, :5])

# 角ラベルの座標を計算
r = 0.45
angle_vec = np.array(
    [r * np.cos(np.median(rad_vals)), 
     r * np.sin(np.median(rad_vals))]
)
print(angle_vec)

[[0.29104275 0.29049241 0.28992594 0.28934336 0.28874471]
 [0.07276069 0.07492769 0.07709054 0.0792491  0.08140327]]
[0.367837  0.2592218]

　詳しくは「円周上の2点を結ぶ曲線を作図したい(今書いてる)」を参照してください。

　なす角のグラフを作成します。

# グラフサイズ用の値を設定
x_min = np.floor(np.min([0.0, a1[0], a2[0]])) - 1.0
x_max =  np.ceil(np.max([0.0, a1[0], a2[0]])) + 1.0
y_min = np.floor(np.min([0.0, a1[1], a2[1]])) - 1.0
y_max =  np.ceil(np.max([0.0, a1[1], a2[1]])) + 1.0

# 2Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.quiver(0, 0, *a1, 
          color='red', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a_1=('+', '.join(map(str, a1))+')$') # ベクトルa1
ax.quiver(0, 0, *a2, 
          color='blue', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a_2=('+', '.join(map(str, a2))+')$') # ベクトルa2
ax.plot(*angle_arr, 
        color='black', linewidth=1) # なす角マーク
ax.text(*angle_vec, 
        s='$\\theta$', size=15, ha='center', va='center') # なす角ラベル
ax.text(*0.5*a1, s='$a_1$', size=15, ha='left', va='top') # ベクトルa1ラベル
ax.text(*0.5*a2, s='$a_2$', size=15, ha='left', va='top') # ベクトルa2ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
ax.grid()
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('$\\theta=' + str((theta / np.pi).round(2)) + '\pi, ' + 
             '\|a_1\|=' + str(norm_a1.round(2)) + ', ' + 
             '\|a_2\|=' + str(norm_a2.round(2)) + '$', loc='left')
fig.suptitle('angle between two vectors', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

　axes.quiver()でベクトルを描画します。第1・2引数に始点の座標、第3・4引数にベクトルのサイズ(移動量)を指定します。その他に、指定した値の通りに(調整せずに)矢印を描画するための設定をしています。
　配列a1, a2の前に*を付けてアンパック(展開)して座標を指定しています。始点は原点とします。

　左図はなす角が鋭角のとき、右図は鈍角のときのグラフです。(左右の図の作成時に引数の設定を少し変更しています。)

　 $\mathbf{a}_1, \mathbf{a}_2$ のなす角を $\theta$ とします。

　2つのベクトルと正射影ベクトル・垂線ベクトルの関係をグラフで確認します。

・作図コード(クリックで展開)

　正射影ベクトルとその直交ベクトルを計算します。

# ベクトルを正規化
q1 = a1 / norm_a1
print(q1)

# 正射影ベクトルを計算
p = dot_a1a2 / norm_a1**2 * a1
print(p)

# 垂線ベクトルを計算
tilde_q2 = a2 - p
print(tilde_q2)

[0.9701425  0.24253563]
[2.58823529 0.64705882]
[-0.58823529  2.35294118]

　直角マークを描画するための座標を作成します。

# ベクトルを正規化
q2 = tilde_q2 / np.linalg.norm(tilde_q2)

# 直角マークの座標を計算
r = 0.2
rightangle_pq2_arr = np.array(
    [[p[0]+r*q2[0], p[0]+r*(q2[0]+q1[0]), p[0]+r*q1[0]], 
     [p[1]+r*q2[1], p[1]+r*(q2[1]+q1[1]), p[1]+r*q1[1]]]
)
print(rightangle_pq2_arr)

[[2.53972817 2.73375667 2.78226379]
 [0.84108732 0.88959445 0.69556595]]

　正射影ベクトルと垂線ベクトルの交点(垂線の足) $\mathbf{p}$ の座標から、ベクトル $\tilde{\mathbf{q}}_2$ とベクトル $\mathbf{q}_1$ 方向に少しずつ移動した点の座標を配列に格納します。

　原点と点 $\mathbf{a}_1$ を通る(ベクトル $\mathbf{a}_1$ と平行な)直線を描画するための座標を作成します。

# 原点を通るベクトルaと平行な直線の座標を計算
line_x_vals = np.linspace(start=x_min, stop=x_max, num=101)
line_arr = np.array(
    [line_x_vals, 
     a1[1]/a1[0] * line_x_vals]
)
print(line_arr[:, :5])

[[-1.    -0.94  -0.88  -0.82  -0.76 ]
 [-0.25  -0.235 -0.22  -0.205 -0.19 ]]

　ベクトル $\mathbf{a}_1$ の傾き $\frac{a_{1,2}}{a_{1,1}}$ を用いて、直線の座標 $y = \frac{a_{1,2}}{a_{1,1}} x$ を計算します。ただし、a1[0]( $a_{1,1}$ )が0だと0除算になるため計算できません。

　正射影のグラフを作成します。

# 2Dの正射影ベクトルを作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.quiver(0, 0, *q1, 
          color='hotpink', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$q_1=('+', '.join(map(str, q1.round(2)))+')$') # 正規化ベクトルq1
ax.quiver(*p, *tilde_q2, 
          color='chocolate', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\tilde{q}_2=('+', '.join(map(str, tilde_q2.round(2)))+')$') # 未正規化垂線ベクトルq2
ax.quiver(0, 0, *p, 
          color='orange', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # 正射影ベクトルp
ax.plot(*line_arr, 
        color='red', linestyle=':', zorder=0) # 直線a1
ax.quiver(0, 0, *a1, 
          fc='none', ec='red', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a_1=('+', '.join(map(str, a1.round(2)))+')$') # ベクトルa1
ax.quiver(0, 0, *a2, 
          fc='none', ec='blue', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a_2=('+', '.join(map(str, a2))+')$') # ベクトルa2
ax.plot(*rightangle_pq2_arr, 
        color='black', linewidth=1) # 直角マーク
ax.plot(*angle_arr, 
        color='black', linewidth=1) # なす角マーク
ax.text(*angle_vec, 
        s='$\\theta$', size=15, ha='center', va='center') # なす角ラベル
ax.text(*0.5*q1, s='$q_1$', 
        size=15, ha='center', va='top') # 正規化ベクトルq1ラベル
ax.text(*0.5*p, s='$p = \\frac{a_1 \cdot a_2}{\|a_1\|^2} a_1$', 
        size=15, ha='left', va='top') # 正射影ベクトルpラベル
ax.text(*0.5*(a2+p), s='$\\tilde{q}_2 = a_2-p$', 
        size=15, ha='left', va='center') # 未正規化垂線ベクトルq2ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
ax.grid()
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('$\|q_1\|=' + str(np.linalg.norm(q1).round(2)) + ', ' + 
             '\|p\|=' + str(np.linalg.norm(p).round(2)) + ', ' + 
             '\|\\tilde{q}_2\|=' + str(np.linalg.norm(tilde_q2).round(2)) + '$', loc='left')
fig.suptitle('vector projection', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

　正射影ベクトルは、次の式で定義されました(5.3節)。

$\displaystyle \mathbf{p} = \frac{ \mathbf{a}_1^{\top} \mathbf{a}_2 }{ \|\mathbf{a}_1\|^2 } \mathbf{a}_1$

　ベクトル $\mathbf{a}_1$ を正規化したベクトルを

$\displaystyle \mathbf{q}_1 = \frac{\mathbf{a}_1}{\|\mathbf{a}_1\|}$

と置いて、正射影ベクトルの式を変形します。

$\displaystyle \begin{aligned} \mathbf{p} &= \left( \Bigl( \frac{\mathbf{a}_1}{\|\mathbf{a}_1\|} \Bigr)^{\top} \mathbf{a}_2 \right) \frac{\mathbf{a}_1}{\|\mathbf{a}_1\|} \\ &= (\mathbf{q}_1^{\top} \mathbf{a}_2) \mathbf{q}_1 \end{aligned}$

　ベクトル $\mathbf{q}_1$ や $\mathbf{a}_1$ と直交するベクトル $\tilde{\mathbf{q}}_2$ として、次の式で計算できます。

$\displaystyle \tilde{\mathbf{q}}_2 = \mathbf{a}_2 - \mathbf{p}$

　2つの正規直交ベクトルの関係をグラフで確認します。

・作図コード(クリックで展開)

　正規直交ベクトルを計算します。

# ベクトルを正規化
q2 = tilde_q2 / np.linalg.norm(tilde_q2)
print(q2)

[-0.24253563  0.9701425 ]

　直角マークを描画するための座標を作成します。

# 直角マークの座標を計算
r = 0.2
rightangle_q1q2_arr = np.array(
    [[r*q2[0], r*(q2[0]+q1[0]), r*q1[0]], 
     [r*q2[1], r*(q2[1]+q1[1]), r*q1[1]]]
)
print(rightangle_q1q2_arr)

[[-0.04850713  0.14552138  0.1940285 ]
 [ 0.1940285   0.24253563  0.04850713]]

　原点の座標から、1つ目の正規直交ベクトル $\mathbf{q}_1$ 方向と2つ目の正規直交ベクトル $\mathbf{q}_2$ 方向に少しずつ移動した点の座標を配列に格納します。

　正規直交ベクトルのグラフを作成します。

# 2Dの正規直交ベクトルを作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.plot(*line_arr, 
        color='red', linestyle=':', zorder=0) # 直線a1
ax.quiver(0, 0, *q1, 
          color='hotpink', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$q_1=('+', '.join(map(str, q1.round(2)))+')$') # 正規化ベクトルq1
ax.quiver(0, 0, *q2, 
          color='chocolate', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$q_2=('+', '.join(map(str, q2.round(2)))+')$') # 正規化垂線ベクトルq2
ax.quiver(0, 0, *tilde_q2, 
          fc='none', ec='chocolate', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\tilde{q}_2=('+', '.join(map(str, tilde_q2.round(2)))+')$') # 未正規化垂線ベクトルq2
ax.quiver(*p, *tilde_q2, 
          fc='none', ec='chocolate', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1) # 未正規化垂線ベクトルq2
ax.quiver(0, 0, *p, 
          fc='none', ec='orange', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # 正射影ベクトルp
ax.quiver(*tilde_q2, *p, 
          fc='none', ec='orange', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1) # 正射影ベクトルp
ax.quiver(0, 0, *a2, 
          fc='none', ec='blue', ls='--', lw=1.0, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a_2=('+', '.join(map(str, a2))+')$') # ベクトルa2
ax.plot(*rightangle_q1q2_arr, 
        color='black', linewidth=1) # 直角マーク
ax.plot(*rightangle_pq2_arr, 
        color='black', linewidth=1) # 直角マーク
ax.text(*0.5*q1, s='$q_1$', 
        size=15, ha='center', va='top') # 正規化ベクトルq1ラベル
ax.text(*0.5*q2, s='$q_2$', 
        size=15, ha='right', va='center') # 正規化ベクトルq2ラベル
ax.text(*0.5*p, s='$p = (q_1 \cdot a_2) q_1$', 
        size=15, ha='left', va='top') # 正射影ベクトルpラベル
ax.text(*0.5*(a2+p), s='$\\tilde{q}_2 = a_2-p$', 
        size=15, ha='left', va='center') # 未正規化垂線ベクトルq2ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+0.1))
ax.set_yticks(ticks=np.arange(y_min, y_max+0.1))
ax.grid()
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_title('$\|q_1\|=' + str(np.linalg.norm(q1).round(2)) + ', ' + 
             '\|q_2\|=' + str(np.linalg.norm(q2).round(2)) + ', ' + 
             'q_1 \cdot q_2=' + str(np.dot(q1, q2).round(2)) + '$', loc='left')
fig.suptitle('orthonormal vectors', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()