からっぽのしょこ

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

5.3:正射影ベクトルの導出【『スタンフォード線形代数入門』のノート】

はじめに

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

 この記事は5.3節「正規直交ベクトル」の内容です。
 2ベクトル間の正射影ベクトルを導出します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

正射影ベクトルの導出

 正射影ベクトル(orthographic projection of a vector・vector projection)を数式とグラフで確認します。
 なす角については「【Python】3.4:ベクトルとx軸のなす角の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」、作図については「【Python】5.3:正射影ベクトルの可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」を参照してください。

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

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


式の導出

 正射影ベクトルの計算式を導出します。2次元の場合を例として図示しますが、次元に関わらず成立します。

 ベクトル \mathbf{a} = (a_1, \cdots, a_n)^{\top}に対してベクトル \mathbf{b} = (b_1, \cdots, b_n)^{\top}を正射影したベクトル \mathbf{p} = (p_1, \cdots, p_n)^{\top}を考えます。

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

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

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

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

# ベクトルa,bのなす角を計算
dot_ab = np.dot(a, b)
norm_a = np.linalg.norm(a)
norm_b = np.linalg.norm(b)
theta = np.arccos(dot_ab / norm_a / norm_b)
print(theta) # 0からπの値
print(theta / np.pi) # 0から1の値
0.7378150601204649
0.2348538278116319

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

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

# ベクトルa,bとx軸のなす角を計算
sgn_a = 1.0 if a[1] >= 0.0 else -1.0
theta_a = sgn_a * np.arccos(a[0] / norm_a)
sgn_b = 1.0 if b[1] >= 0.0 else -1.0
theta_b = sgn_b * np.arccos(b[0] / norm_b)

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

# 角マーク用の値(ラジアン)を作成
rad_vals = np.linspace(start=theta_a, stop=theta_b, 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, a[0], b[0]])) - 1.0
x_max =  np.ceil(np.max([0.0, a[0], b[0]])) + 1.0
y_min = np.floor(np.min([0.0, a[1], b[1]])) - 1.0
y_max =  np.ceil(np.max([0.0, a[1], b[1]])) + 1.0

# 2Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.quiver(0, 0, *a, 
          color='red', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
ax.quiver(0, 0, *b, 
          color='blue', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$b=('+', '.join(map(str, b))+')$') # ベクトルb
ax.plot(*angle_arr, 
        color='black', linewidth=1) # なす角マーク
ax.text(*angle_vec, 
        s='$\\theta$', size=15, ha='center', va='center') # なす角ラベル
ax.text(0, 0, s='O', size=15, ha='right', va='top') # 原点ラベル
ax.text(*a, s='A', size=15, ha='left', va='top') # 点aラベル
ax.text(*b, s='B', size=15, ha='left', va='bottom') # 点bラベル
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$', loc='left')
fig.suptitle('angle between two vectors', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

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

ベクトルのなす角

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

 原点を O、原点からそれぞれ \mathbf{a}, \mathbf{b}移動した点を点 A, Bとします。

 2つのベクトル \mathbf{a}, \mathbf{b}のなす角 \thetaは、次の式で定義されます(3.4節)。ただし、 \mathbf{a}または \mathbf{b}が0ベクトルだと0除算になるため定義できません。

 \displaystyle
\theta
    = \arccos \left(
          \frac{
              \mathbf{a}^{\top} \mathbf{b}
          }{
              \|\mathbf{a}\| \|\mathbf{b}\|
          }
      \right)

  \arccos xは逆コサイン関数(コサイン関数 \cos xの逆関数)(3.4節)、 \|\mathbf{x}\|はユークリッドノルム(3.1節)でベクトルの大きさを表します。
  \thetaは、弧度法のおける角度(ラジアン)です。逆コサイン関数の定義より、 0 \leq \theta \leq \piの値(度数法による角度だと 0^{\circ} \leq \theta^{\circ} \leq 180^{\circ}に対応)をとります。

 また内積は、( \mathbf{x}^{\top} \mathbf{y} = \sum_{i=1}^n x_i y_iとは別に)、次の式で定義されます(3.4節)。

 \displaystyle
\mathbf{a}^{\top} \mathbf{b}
    = \|\mathbf{a}\| \|\mathbf{b}\|
      \cos \theta

  \mathbf{a}, \mathbf{b}が直交する(直角に交わる)とき、なす角が \theta = 0であり、 \cos \theta = 0なので、 \mathbf{a}^{\top} \mathbf{b} = 0です。

 「上の式の両辺のコサインをとる」または「下の式をコサインについて整理する」と、コサイン関数は、次の式で表わせます。

 \displaystyle
\cos \theta
    = \frac{
          \mathbf{a}^{\top} \mathbf{b}
      }{
          \|\mathbf{a}\| \|\mathbf{b}\|
      }
\tag{1}

 コサイン関数は、 0 \leq \cos \theta \leq 1の値をとります。

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

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

 (導出前ですが)正射影ベクトルを計算して、直角マークを描画するための座標を作成します。

# 正射影ベクトルを計算
p = dot_ab / norm_a**2 * a
print(p)

# ベクトルを正規化
e_a  = a / norm_a
e_bp = (b - p) / np.linalg.norm(b - p)

# 直角マークの座標を計算
r = 0.2
rightangle_arr = np.array(
    [[p[0]+r*e_bp[0], p[0]+r*(e_bp[0]-e_a[0]), p[0]-r*e_a[0]], 
     [p[1]+r*e_bp[1], p[1]+r*(e_bp[1]-e_a[1]), p[1]-r*e_a[1]]]
)
print(rightangle_arr)
[2.58823529 0.64705882]
[[2.53972817 2.34569967 2.39420679]
 [0.84108732 0.7925802  0.5985517 ]]

 点 Pの座標から、垂線方向とベクトル \mathbf{p}方向に少しずつ移動した点の座標を配列に格納します。

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

# 原点を通るベクトルaと平行な直線の座標を計算
line_x_vals = np.linspace(start=x_min, stop=x_max, num=101)
line_arr = np.array(
    [line_x_vals, 
     a[1]/a[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}の傾き \frac{a_2}{a_1}を用いて、直線の座標 y = \frac{a_2}{a_1} xを計算します。ただし、a[0]( a_1)が0だと0除算になるため計算できません。

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

# 2Dベクトルの正射影を作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.quiver(0, 0, *a, 
          color='red', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
ax.plot(*line_arr, 
        color='red', linestyle=':', zorder=0) # 直線a
ax.quiver(0, 0, *b, 
          color='blue', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$b=('+', '.join(map(str, b))+')$') # ベクトルb
ax.quiver(0, 0, *p, 
          fc='none', ec='black', ls='--', lw=1.5, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # 正射影ベクトルp
ax.plot([b[0], p[0]], [b[1], p[1]], 
        color='black', linestyle=':') # 点bから直線aへの垂線
ax.plot(*rightangle_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, 0, s='O', size=15, ha='right', va='top') # 原点ラベル
ax.text(*p, s='P', size=15, ha='left', va='top') # 垂線の足ラベル
ax.text(*a, s='A', size=15, ha='left', va='top') # 点aラベル
ax.text(*b, s='B', size=15, ha='left', va='bottom') # 点bラベル
ax.text(*0.5*a, s='$\|a\|$', 
        size=15, ha='center', va='bottom') # ベクトルaのノルムラベル
ax.text(*0.5*b, s='$\|b\|$', 
        size=15, ha='right', va='center') # ベクトルbのノルムラベル
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('$\|a\|=' + str(norm_a.round(2)) + ', ' + 
             '\|b\|=' + str(norm_b.round(2)) + '$', loc='left')
fig.suptitle('vector projection', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

ベクトルの正射影

 点 Bからベクトル \mathbf{a}またはその延長線(原点と点 Aを通る直線)上への垂線を考えます。その交点(垂線の足)を点 Pとします。ベクトル \overrightarrow{OP}を正射影ベクトルと言います。ここでは、正射影ベクトルを \mathbf{p}とします。
  \mathbf{p} \mathbf{a}と平行なベクトルなので、 \betaをスカラーとして \mathbf{p} = \beta \mathbf{a}で表せます(1.3節)。

 正射影ベクトルと平行なベクトルの大きさの関係をグラフで確認します。

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

 ベクトル \mathbf{a}を正規化します。

# ベクトルを正規化
tilde_a = a / norm_a
print(tilde_a)
[0.9701425  0.24253563]


 正射影ベクトルのグラフを作成します。

# ラベル用の文字列を設定
sgn_p = '-' if theta > 0.5*np.pi else ''

# 2Dの正射影ベクトルを作図
fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.quiver(0, 0, *tilde_a, 
          color='hotpink', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # 正規化ベクトルa
ax.plot(*line_arr, 
        color='red', linestyle=':', zorder=0) # 直線a
ax.quiver(0, 0, *b, 
          color='blue', width=0.01, headwidth=4, headlength=5, 
          angles='xy', scale_units='xy', scale=1, 
          label='$b=('+', '.join(map(str, b))+')$') # ベクトルb
ax.quiver(0, 0, *p, 
          fc='none', ec='black', ls='--', lw=1.5, width=0.001, headwidth=40, headlength=50, 
          angles='xy', scale_units='xy', scale=1, 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # 正射影ベクトルp
ax.plot([b[0], p[0]], [b[1], p[1]], 
        color='black', linestyle=':') # 点bから直線aへの垂線
ax.plot(*rightangle_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, 0, s='O', size=15, ha='right', va='top') # 原点ラベル
ax.text(*p, s='P', size=15, ha='left', va='top') # 垂線の足ラベル
ax.text(*b, s='B', size=15, ha='left', va='bottom') # 点bラベル
ax.text(*0.5*b, s='$\|b\|$', 
        size=15, ha='right', va='center') # 斜辺ラベル
ax.text(*0.5*p, s='$' + sgn_p + '\|p\| = \|b\| \cos \\theta$', 
        size=15, ha='left', va='top') # 底辺ラベル
ax.text(*0.5*(b+p), s='$\|b-p\| = \|b\| \sin \\theta$', 
        size=15, ha='left', va='center') # 垂線ラベル
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, ' + 
             '\|\\tilde{a}\|=' + str(np.linalg.norm(tilde_a).round(2)) + ', ' + 
             '\|\\tilde{p}\|=' + str(np.linalg.norm(p).round(2)) + '$', loc='left')
fig.suptitle('orthographic projection vector', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

正射影ベクトル

 ベクトル \mathbf{a}を正規化したベクトルを \tilde{\mathbf{a}} = \frac{\mathbf{a}}{\|\mathbf{a}\|}と置きます。正規化したベクトルは \|\tilde{\mathbf{a}}\| = 1になります(ノルムを1にすることを正規化すると言います)。
 ノルムはスカラーなので、 \mathbf{a} \tilde{\mathbf{a}}は平行なベクトルです。つまり、 \mathbf{p} \tilde{\mathbf{a}}も平行なベクトルです。
 よって、なす角が鋭角 0 \leq \theta \leq \frac{\theta}{2}のとき、 \mathbf{p}は、 \tilde{\mathbf{a}}と平行でノルムが \|\mathbf{p}\|のベクトルと言えるので、次の式で表わせます(正規化の反対の操作です)。

 \displaystyle
\begin{aligned}
\mathbf{p}
   &= \|\mathbf{p}\| \tilde{\mathbf{a}}
\\
   &= \|\mathbf{p}\|
      \frac{\mathbf{a}}{\|\mathbf{a}\|}
\end{aligned}

 また、なす角が鈍角 \frac{\theta}{2} \leq \theta \leq \piのとき、 \mathbf{p}は、 \tilde{\mathbf{a}}と平行で向きが逆でノルムが \|\mathbf{p}\|のベクトルと言えるので、次の式で表わせます。

 \displaystyle\begin{aligned}
\mathbf{p}
   &= - \|\mathbf{p}\| \tilde{\mathbf{a}}
\\
   &= - \|\mathbf{p}\|
      \frac{\mathbf{a}}{\|\mathbf{a}\|}
\end{aligned}

 なす角が直角 \theta = \frac{\theta}{2}のとき、 \|\mathbf{p}\|が0になるので、どちらの式でも(符号に関わらず)成り立ちます。
 まとめると、次の式で表わせます。

 \displaystyle
\mathbf{p}
    = \begin{cases}
          \|\mathbf{p}\|
          \frac{\mathbf{a}}{\|\mathbf{a}\|}
             &\quad
                (0 \leq \theta \leq \frac{\theta}{2}) \\
          - \|\mathbf{p}\|
          \frac{\mathbf{a}}{\|\mathbf{a}\|}
             &\quad
                (\frac{\theta}{2} < \theta \leq \pi)
      \end{cases}

  \pm \|\mathbf{p}\|を符号付きの長さと言います。

 ベクトル \mathbf{b}を斜辺とする直角三角形 BOPを考えます。
 底辺 OPの符号付きの長さは、ベクトル \mathbf{p}のノルムを使って \pm \|\mathbf{p}\|であり、またコサインの定義(底辺を斜辺で割ると \cos x)より \|\mathbf{b}\| \cos \thetaなので、場合分けした式を置き替えます( \theta \cos \thetaの符号の関係は上の式の条件と一致します)。

 \displaystyle
\mathbf{p}
    = \|\mathbf{b}\| \cos \theta
      \frac{\mathbf{a}}{\|\mathbf{a}\|}

 さらに、式(1)より \cos \theta = \frac{\mathbf{a}^{\top} \mathbf{b}}{\|\mathbf{a}\| \|\mathbf{b}\|}で置き換えて、式を整理します。

 \displaystyle
\begin{aligned}
\mathbf{p}
   &= \|\mathbf{b}\|
      \frac{
          \mathbf{a}^{\top} \mathbf{b}
      }{
          \|\mathbf{a}\| \|\mathbf{b}\|
      }
      \frac{\mathbf{a}}{\|\mathbf{a}\|}
\\
   &= \frac{
          \mathbf{a}^{\top} \mathbf{b}
      }{
          \|\mathbf{a}\|^2
      }
      \mathbf{a}
\end{aligned}

 またこの式は、ノルムの定義よりノルムの2乗は内積 \|\mathbf{x}\|^2 = \mathbf{x}^{\top} \mathbf{x}なので(3.1節)、次の式でも表せます。

 \displaystyle
\mathbf{p}
    = \frac{
          \mathbf{a}^{\top} \mathbf{b}
      }{
          \mathbf{a}^{\top} \mathbf{a}
      }
      \mathbf{a}

 正射影ベクトルの計算式が得られました。

 3次元の場合もグラフで確認します。

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

 3次元ベクトルを指定して、同様に処理します。

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

# ベクトルa,bのなす角を計算
dot_ab = np.dot(a, b)
norm_a = np.linalg.norm(a)
norm_b = np.linalg.norm(b)
theta = np.arccos(dot_ab / norm_a / norm_b)
print(theta) # 0からπの値
print(theta / np.pi) # 0から1の値

# ベクトルを正規化
tilde_a = a / norm_a
print(tilde_a)
2.2870084014997913
0.7279773839827717
[-0.94280904 -0.23570226 -0.23570226]

 2つのベクトル \mathbf{a} = (a_1, a_2, a_3)^{\top}, \mathbf{b} = (b_1, b_2, b_3)^{\top}a, bとして値を指定します。

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


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

# 角マーク用の値(ラジアン)を作成
rad_vals = np.linspace(start=0.0, stop=0.5*np.pi, num=100)
rad_arr  = np.stack([rad_vals, rad_vals, rad_vals])
print(rad_arr.shape)

# 混合率(線形結合の係数)を作成
beta1_arr = np.cos(rad_arr)
beta2_arr = np.sin(rad_arr)

# 角マークの座標を計算
r = 0.5
angle_arr = (beta1_arr.T * a/np.linalg.norm(a) + beta2_arr.T * b/np.linalg.norm(b)).T # 角マーク用の曲線を計算
angle_arr *= r / np.linalg.norm(angle_arr, axis=0) # 曲線のサイズを調整
print(angle_arr[:, :5])

# 角ラベル用の値(角マークの中点のラジアン)を作成
rad = 0.25 * np.pi

# 混合率(線形結合の係数)を作成
beta1 = np.cos(rad)
beta2 = np.sin(rad)

# 角ラベルの座標を計算
r = 0.8
angle_vec = beta1 * a/np.linalg.norm(a) + beta2 * b/np.linalg.norm(b) # なす角の中点の座標を計算
angle_vec *= r / np.linalg.norm(angle_vec) # 表示位置を調整
print(angle_vec)
(3, 100)
[[-0.47140452 -0.47335498 -0.47527588 -0.47716395 -0.47901562]
 [-0.11785113 -0.11461698 -0.11129601 -0.10788466 -0.10437923]
 [-0.11785113 -0.11312828 -0.10828683 -0.10332213 -0.09822936]]
[-0.5515516   0.31021006  0.48944924]

 詳しくは「球面上の2点を結ぶ弧を作図したい【Matplotlib】 - からっぽのしょこ」を参照してください。

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

# 正射影を計算
p = dot_ab / norm_a**2 * a

# ベクトルを正規化
e_a  = a / norm_a
e_bp = (b - p) / np.linalg.norm(b - p)

# 直角マーク用の座標を計算
r = 0.2
rightangle_arr = np.array(
    [[p[0]+r*e_bp[0], p[0]+r*(e_bp[0]-e_a[0]), p[0]-r*e_a[0]], 
     [p[1]+r*e_bp[1], p[1]+r*(e_bp[1]-e_a[1]), p[1]-r*e_a[1]], 
     [p[2]+r*e_bp[2], p[2]+r*(e_bp[2]-e_a[2]), p[2]-r*e_a[2]]]
)
print(rightangle_arr)
[[3.26768454 3.45624635 3.52189514]
 [0.94001263 0.98715308 0.88047379]
 [0.98924922 1.03638967 0.88047379]]

 点 Pの座標から、垂線方向とベクトル \mathbf{p}方向に少しずつ移動した点の座標を配列に格納します。

 直線を描画するための座標を作成します。

# 作図用の値を設定
x_min = np.floor(np.min([0.0, a[0], b[0], p[0]])) - 1.0
x_max =  np.ceil(np.max([0.0, a[0], b[0], p[0]])) + 1.0
y_min = np.floor(np.min([0.0, a[1], b[1], p[1]])) - 1.0
y_max =  np.ceil(np.max([0.0, a[1], b[1], p[1]])) + 1.0
z_min = np.floor(np.min([0.0, a[2], b[2], p[2]])) - 1.0
z_max =  np.ceil(np.max([0.0, a[2], b[2], p[2]])) + 1.0

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

# 描画範囲外の要素を除去
bool_idx = (line_arr[1] < y_min) | (line_arr[1] > y_max) | (line_arr[2] < z_min) | (line_arr[2] > z_max)
bool_idx = [not bl for bl in bool_idx]
line_arr = line_arr[:, bool_idx]
print(line_arr[:, :5])
[[-5.    -4.9   -4.8   -4.7   -4.6  ]
 [-1.25  -1.225 -1.2   -1.175 -1.15 ]
 [-1.25  -1.225 -1.2   -1.175 -1.15 ]]

 「2次元の場合」ときに加えて、z軸方向の傾き \frac{a_3}{a_1}を用いて、直線のz軸の座標 z = \frac{a_3}{a_1} xを計算します。ただし、a[0]( a_1)が0だと0除算になるため計算できません。
 また、描画範囲外の要素を取り除きます。

 正射影ベクトルのグラフを作成します。

# 矢のサイズを指定
l = 0.5

# ラベル用の文字列を設定
sgn_p = '-' if theta > 0.5*np.pi else ''

# 3D正射影を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.quiver(0, 0, 0, *a, 
          color='red', lw=1, arrow_length_ratio=l/np.linalg.norm(a), 
          label='$a=('+', '.join(map(str, a))+')$') # ベクトルa
ax.quiver(0, 0, 0, *tilde_a, 
          color='hotpink', lw=3, arrow_length_ratio=l/np.linalg.norm(tilde_a), 
          label='$\\tilde{a}=('+', '.join(map(str, tilde_a.round(2)))+')$') # 正規化ベクトルa
ax.quiver(0, 0, 0, *b, 
          color='blue', lw=2, arrow_length_ratio=l/np.linalg.norm(b), 
          label='$b=('+', '.join(map(str, b))+')$') # ベクトルb
ax.quiver(0, 0, 0, *p, 
          color='black', lw=2, ls='--', arrow_length_ratio=l/np.linalg.norm(p), 
          label='$p=('+', '.join(map(str, p.round(2)))+')$') # 正射影ベクトルp
ax.plot(*line_arr, 
        color='red', linestyle=':') # ベクトルaの延長線
ax.plot([b[0], p[0]], [b[1], p[1]], [b[2], p[2]], 
        color='black', linestyle=':') # 点bから直線aへの垂線
ax.plot(*rightangle_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', zorder=20) # なす角ラベル
ax.text(0, 0, 0, s='O', size=15, ha='right', va='top', zorder=20) # 原点ラベル
ax.text(*p, s='P', size=15, ha='left', va='top') # 垂線の足ラベル
ax.text(*a, s='A', size=15, ha='right', va='top', zorder=20) # 点aラベル
ax.text(*b, s='B', size=15, ha='left', va='bottom', zorder=20) # 点bラベル
ax.text(*0.5*b, s='$\|b\|$', 
        size=12, ha='right', va='center') # 斜辺ラベル
ax.text(*0.5*p, s='$' + sgn_p + '\|p\| = \|b\| \cos \\theta$', 
        size=12, ha='center', va='top') # 底辺ラベル
ax.text(*0.5*(b+p), s='$\|b-p\| = \|b\| \sin \\theta$', 
        size=12, ha='left', va='center') # 垂線ラベル
ax.quiver([0, a[0], b[0], p[0]], 
          [0, a[1], b[1], p[1]], 
          [0, a[2], b[2], p[2]], 
          [0, 0, 0, 0], 
          [0, 0, 0, 0], 
          [z_min, z_min-a[2], z_min-b[2], z_min-p[2]], 
          color='gray', arrow_length_ratio=0, linestyle=':') # 座標用の補助線
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.set_zticks(ticks=np.arange(z_min, z_max+1))
ax.set_zlim(bottom=z_min, top=z_max)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$x_3$')
ax.set_title('$\\theta=' + str((theta / np.pi).round(2)) + '\pi, ' + 
             '\|a\|=' + str(norm_a.round(2)) + ', ' + 
             '\|\\tilde{a}\|=' + str(np.linalg.norm(tilde_a).round(2)) + ', ' + 
             '\|b\|=' + str(norm_b.round(2)) + ', ' + 
             '\|p\|=' + str(np.linalg.norm(p).round(2)) + '$', loc='left')
fig.suptitle('vector projection', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

 axes.quiver()でベクトルを描画します。第1・2・3引数に始点の座標、第4・5・6引数にベクトルのサイズ(移動量)を指定します。
 配列a, bの前に*を付けてアンパック(展開)して座標を指定しています。始点は原点です。

3次元ベクトルの正射影

 同様に、 n次元ベクトルに対しても成り立ちます。

 この記事では、正射影ベクトルの計算式を導出しました。次の記事では、正射影ベクトルを可視化します。

参考書籍

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

おわりに

 導出の記事でも自然とPythonでお絵描きするようになってしまった。 先に書いた可視化の方の記事と内容が被ってしまう。そしてコードが長い。
 と言いつつ丁寧に図で確認しながらやれば微妙に分からなかった式変形を追えました。 \|\mathbf{b}\| \cos \thetaはどこから出てきた?と思ったらコサインの定義そのままだったり、逆に \frac{1}{\|\mathbf{a}\|^2}が幾何的意味というより式変形によるものだったり、最終的な計算式が2つあるなと思ったらノルムの2乗は内積で表現できるんだったとか、これ辺の長さがマイナスということ?と調べたら符号付き長さなる概念に辿り着いたとか。

 線形代数の勉強をしているのか三角関数の勉強をしているのか段々分からなくなってきました。

【次の内容】

www.anarchive-beta.com