はじめに
『スタンフォード ベクトル・行列からはじめる最適化数学』の学習ノートです。
「数式の行間埋め」や「Pythonを使っての再現」によって理解を目指します。本と一緒に読んでください。
この記事は5.4節「グラム・シュミット法」の内容です。
2ベクトルの正規直交ベクトルを導出します。
【前の内容】
【他の内容】
【今回の内容】
グラム・シュミット法の導出
グラム・シュミット法(Gram-Schmidt algorithm)によるベクトルの正規直交化(orthogonalization)の計算を数式とグラフで確認します。
正射影については「5.3:正射影ベクトルの導出【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照してください。
次のライブラリを利用してグラフを作成します。不要であれば省略してください。
# 利用ライブラリ import numpy as np import matplotlib.pyplot as plt
式の導出
グラム・シュミット法による正規直交化の計算式を導出します。2次元の場合を例として図示しますが、次元に関わらず成立します。
ベクトルとを正規直交化したベクトルとを考えます。
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つのベクトルをa1, a2
として値を指定します。Pythonでは0からインデックスが割り当てられるので、の値はa1[0]
に対応します。
ベクトルのなす角を計算します。
は逆コサイン関数(コサイン関数の逆関数)(3.4節)、はユークリッドノルム(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
の前に*
を付けてアンパック(展開)して座標を指定しています。始点は原点とします。
左図はなす角が鋭角のとき、右図は鈍角のときのグラフです。(左右の図の作成時に引数の設定を少し変更しています。)
のなす角をとします。
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]]
正射影ベクトルと垂線ベクトルの交点(垂線の足)の座標から、ベクトルとベクトル方向に少しずつ移動した点の座標を配列に格納します。
原点と点を通る(ベクトルと平行な)直線を描画するための座標を作成します。
# 原点を通るベクトル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 ]]
ベクトルの傾きを用いて、直線の座標を計算します。ただし、a1[0]
()が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節)。
ベクトルを正規化したベクトルを
と置いて、正射影ベクトルの式を変形します。
ベクトルやと直交するベクトルとして、次の式で計算できます。
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つ目の正規直交ベクトル方向と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()
ベクトルを正規化したベクトルをとします。
は、それぞれノルムがで、直交します(内積がになります)。
3次元の場合をこれから書きます…
この記事では、グラム・シュミット法の計算式を導出しました。次の記事では、グラム・シュミット法による正規直交化を可視化します。
参考書籍
- Stephen Boyd・Lieven Vandenberghe(著),玉木 徹(訳)『スタンフォード ベクトル・行列からはじめる最適化数学』講談社サイエンティク,2021年.
おわりに
5章の中では最初に書いた(内容的には次の)記事から1か月以上かかってしまいましたが、この記事で5章完了です。
2023年5月5日は、私立恵比寿中学のメジャーデビュー11周年&ココユノノカ加入2周年の日です!!!
エビ中の魅力に気付いてまだ1年未満の新参で絶賛ハマってる途中ですが、これからも楽しみです🦐
【次の内容】