からっぽのしょこ

読んだら書く!書いたら読む!同じ事は二度調べ(たく)ない

5.4:グラム・シュミット法の導出【『スタンフォード線形代数入門』のノート】

はじめに

 『スタンフォード ベクトル・行列からはじめる最適化数学』の学習ノートです。
 「数式の行間埋め」や「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()

正規直交ベクトル

 ベクトル \tilde{\mathbf{q}}_2を正規化したベクトルを \mathbf{q}_2とします。

 \displaystyle
\mathbf{q}_2
    = \frac{\tilde{\mathbf{q}}_2}{\|\tilde{\mathbf{q}}_2\|}

  \mathbf{q}_1, \mathbf{q}_2は、それぞれノルムが \|\mathbf{q}_1\| = \|\mathbf{q}_2\| = 1で、直交します(内積が \mathbf{q}_1^{\top} \mathbf{q}_2 = 0になります)。

 3次元の場合をこれから書きます…

 この記事では、グラム・シュミット法の計算式を導出しました。次の記事では、グラム・シュミット法による正規直交化を可視化します。

参考書籍

  • Stephen Boyd・Lieven Vandenberghe(著),玉木 徹(訳)『スタンフォード ベクトル・行列からはじめる最適化数学』講談社サイエンティク,2021年.

おわりに

 5章の中では最初に書いた(内容的には次の)記事から1か月以上かかってしまいましたが、この記事で5章完了です。

 2023年5月5日は、私立恵比寿中学のメジャーデビュー11周年&ココユノノカ加入2周年の日です!!!

 エビ中の魅力に気付いてまだ1年未満の新参で絶賛ハマってる途中ですが、これからも楽しみです🦐

【次の内容】

www.anarchive-beta.com