はじめに

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

　この記事は3.4節「角度」の内容です。
　2つのベクトルのなす角のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2つのベクトルのなす角の可視化
- 2次元の場合
- 3次元の場合
参考書籍
おわりに

2つのベクトルのなす角の可視化

　ベクトルのなす角をグラフで確認します。
　「【Python】3.4：ベクトルとx軸のなす角の可視化【『スタンフォード線形代数入門』のノート】 - からっぽのしょこ」では、1つのベクトルとx軸のなす角のグラフを作成しました。
　この記事では、任意の2つのベクトル $\mathbf{a}, \mathbf{b}$ のなす角のグラフを作成します。重複する内容は省略するので、前回の記事も参照してください。

　利用するライブラリを読み込みます。

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

2次元の場合

　まずは、2次元空間(平面)上でベクトルのなす角を可視化します。

　2次元ベクトルと係数を指定します。

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

# 係数を指定
alpha = 1.0
beta = 1.0

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

　2つのベクトルのなす角を計算します。

# ベクトルa,bのなす角を計算
dot_ab = np.dot(alpha*a, beta*b)
norm_a = np.linalg.norm(alpha * a)
norm_b = np.linalg.norm(beta * b)
theta = np.arccos(dot_ab / norm_a / norm_b)
print(theta) # ラジアン
print(theta / np.pi) # ラジアンの目安
print(theta * 180.0 / np.pi) # 角度

2.8198420991931505
0.8975836176504332
161.56505117707798

　 $\alpha \mathbf{a}$ と $\beta \mathbf{b}$ のなす角 $\theta$ は、次の式で定義されます。

$\displaystyle \theta = \arccos \left( \frac{ \alpha \beta \mathbf{a}^{\top} \mathbf{x} }{ \|\alpha \mathbf{a}\| \|\beta \mathbf{x}\| } \right)$

　コサイン関数 $\cos x$ の逆関数(逆コサイン関数) $\arccos x$ はnp.arccos()、ユークリッドノルム $|\mathbf{x}|$ はnp.linalg.norm()で計算できます。
　ただし、 $\alpha \mathbf{a}$ または $\beta \mathbf{x}$ が0ベクトルだと0除算になるため定義できません(計算結果がnanになります)。

　2つのベクトル(線分)によって劣角 $0 \leq \theta \leq \pi$ と優角 $2 \pi - \theta$ の2つの角ができます。なす角 $\theta$ は劣角の角度に対応します。
　 $\theta$ はラジアン(弧度法における角度)であり、円周率 $\pi$ はnp.piで扱えます。

　なす角を示す角マークを描画するための座標を作成します。なす角の計算自体には不要な処理です。

# ベクトルaのなす角を計算
norm_a = np.linalg.norm(alpha * a)
sgn_a = 1.0 if alpha*a[1] >= 0.0 else -1.0
theta_a = sgn_a * np.arccos(alpha*a[0] / norm_a)
print(theta_a / np.pi)

# ベクトルbのなす角を計算
norm_b = np.linalg.norm(beta * b)
sgn_b = 1.0 if beta*b[1] >= 0.0 else -1.0
theta_b = sgn_b * np.arccos(beta*b[0] / norm_b)
print(theta_b / np.pi)

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

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

# 角マーク用の座標を計算
r = 0.2
angle_x_vals = np.cos(rad_vals) * r
angle_y_vals = np.sin(rad_vals) * r

 # 角ラベルの座標を計算
r = 0.3
angle_x = np.cos(np.median(rad_vals)) * r
angle_y = np.sin(np.median(rad_vals)) * r

0.25000000000000006
-0.8524163823495667
0.25000000000000006
1.1475836176504333

　 $\alpha \mathbf{a}, \beta \mathbf{b}$ それぞれのベクトルとx軸(標準単位ベクトル $\mathbf{e}_1 = (1, 0)^{\top}$ )とのなす角 $\theta_a, \theta_b$ を計算します。ただし、第1・2象限側にできる角が優角の場合(第3象限にまで跨る場合)は、第3・4象限側にできる劣角を負の角度として扱います。詳しくは前回の記事を参照してください。

　 $\theta_a$ から $\theta_b$ までのラジアン $t$ を作成して、角マークの座標を計算します。
　ただし、 $\theta_a - \theta_b$ が優角になる場合( $\theta_a, \theta_b$ の一方が正の角度でもう一方が負の角度で、差の絶対値が $\pi$ を超える場合)は、負の角度の方に $2 \pi$ を足して優角に変更します。
　theta_aからtheta_bまでの値を作成してrad_valsとします。

　rad_valsの各要素を円形の座標に変換します。
　座標 $x, y$ は、半径を $r$ として、それぞれ次の式で計算できます。

$\displaystyle \left\{ \begin{aligned} x &= r \cos t \\ y &= r \sin t \end{aligned} \right.$

　半径rで角マークのサイズを調整します。

　角ラベルは、rad_valsの中点の座標に表示することにします。

　2次元空間上にベクトルのなす角を作図します。

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

# 2Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(6, 6), facecolor='white')
ax.quiver(0, 0, *alpha*a, color='red', 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\alpha a = \\alpha ('+', '.join(map(str, a.round(2)))+')$') # ベクトルαa
ax.quiver(0, 0, *beta*b, color='blue', 
          angles='xy', scale_units='xy', scale=1, 
          label='$\\beta b = \\beta ('+', '.join(map(str, b.round(2)))+')$') # ベクトルβb
ax.plot(angle_x_vals, angle_y_vals, color='black', linewidth=1) # 角マーク
ax.text(angle_x, angle_y, s='$\\theta$', ha='center', va='center') # 角ラベル
ax.set_xticks(ticks=np.arange(x_min, x_max+1))
ax.set_yticks(ticks=np.arange(y_min, y_max+1))
ax.grid()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('$\\alpha='+str(alpha)+', ' + 
             '\\beta='+str(beta)+', ' + 
             '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
             '\\theta^\circ='+str((theta * 180.0 / np.pi).round(1))+'$', loc='left')
fig.suptitle('angle a, b', fontsize=20)
ax.legend(loc='upper left')
ax.set_aspect('equal')
plt.show()

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

　係数の値を変えてみます。

　係数の大小はなす角に影響せず、係数の正負はベクトルの向きに影響するためなす角に影響するのが分かります。

　ベクトルの値を変化させたアニメーションを作成します。

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

# フレーム数を設定
frame_num = 600

# ベクトルとして利用する値を指定
theta_a_vals = np.linspace(start=0.5*np.pi, stop=24.5*np.pi, num=frame_num+1)[:frame_num]
a_vals = np.array(
    [-np.cos(theta_a_vals), 
     np.sin(theta_a_vals)]
).T
theta_b_vals = np.linspace(start=0.5*np.pi, stop=2.5*np.pi, num=frame_num+1)[:frame_num]
b_vals = np.array(
    [-np.cos(theta_b_vals), 
     np.sin(theta_b_vals)]
).T

# 係数を指定
alpha = 1.0
beta = 1.0

# 作図用の値を設定
x_min = np.floor(np.min([0.0, *alpha*a_vals[:, 0], *beta*b_vals[:, 0]])) - 1
x_max = np.ceil(np.max([0.0, *alpha*a_vals[:, 0], *beta*b_vals[:, 0]])) + 1
y_min = np.floor(np.min([0.0, *alpha*a_vals[:, 1], *beta*b_vals[:, 1]])) - 1
y_max = np.ceil(np.max([0.0, *alpha*a_vals[:, 1], *beta*b_vals[:, 1]])) + 1

# 作図用のオブジェクトを初期化
fig, ax = plt.subplots(figsize=(6, 6), facecolor='white')
fig.suptitle('angle a, b', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # i番目のベクトルを取得
    a = a_vals[i]
    b = b_vals[i]
    
    # ベクトルa,bのなす角を計算
    dot_ab = np.dot(alpha*a, beta*b)
    norm_a = np.linalg.norm(alpha * a)
    norm_b = np.linalg.norm(beta * b)
    theta = np.arccos(dot_ab / norm_a / norm_b)
    
    # ベクトルaのなす角を計算
    norm_a = np.linalg.norm(alpha * a)
    sgn_a = 1.0 if alpha*a[1] >= 0.0 else -1.0
    theta_a = sgn_a * np.arccos(alpha*a[0] / norm_a)
    
    # ベクトルbのなす角を計算
    norm_b = np.linalg.norm(beta * b)
    sgn_b = 1.0 if beta*b[1] >= 0.0 else -1.0
    theta_b = sgn_b * np.arccos(beta*b[0] / norm_b)
    
    # 差が180°を超える場合は、負の角度を優角に変換
    if abs(theta_a - theta_b) > np.pi:
        theta_a = theta_a if alpha*a[1] >= 0.0 else theta_a + 2.0*np.pi
        theta_b = theta_b if beta*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.2
    angle_x_vals = np.cos(rad_vals) * r
    angle_y_vals = np.sin(rad_vals) * r

     # 角ラベルの座標を計算
    r = 0.3
    angle_x = np.cos(np.median(rad_vals)) * r
    angle_y = np.sin(np.median(rad_vals)) * r
    
    # 2Dベクトルのなす角を作図
    ax.quiver(0, 0, *alpha*a, color='red', 
              angles='xy', scale_units='xy', scale=1, 
              label='$\\alpha a = \\alpha ('+', '.join(map(str, a.round(2)))+')$') # ベクトルαa
    ax.quiver(0, 0, *beta*b, color='blue', 
              angles='xy', scale_units='xy', scale=1, 
              label='$\\beta b = \\beta ('+', '.join(map(str, b.round(2)))+')$') # ベクトルβb
    ax.plot(angle_x_vals, angle_y_vals, color='black', linewidth=1) # 角マーク
    ax.text(angle_x, angle_y, s='$\\theta$', ha='center', va='center') # 角ラベル
    ax.set_xticks(ticks=np.arange(x_min, x_max+1))
    ax.set_yticks(ticks=np.arange(y_min, y_max+1))
    ax.grid()
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title('$\\alpha='+str(alpha)+', ' + 
                 '\\beta='+str(beta)+', ' + 
                 '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
                 '\\theta^\circ='+str((theta * 180.0 / np.pi).round(2))+'$', loc='left')
    ax.legend(loc='upper left')
    ax.set_aspect('equal')

# gif画像を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100)

# gif画像を保存
ani.save('angle_ab_2d.gif')

　作図処理をupdate()として定義して、FuncAnimation()でgif画像を作成します。

　この例では、角マークの座標の処理と同様にして、単位円上を通るように $\mathbf{a}, \mathbf{b}$ の値を指定しています。
　円周は $2 \pi$ で1周期なので、 $2 \pi$ の $n$ 倍の範囲(サイズ)をフレーム数+1個の等間隔の値に切り分けて最大値を除くと、最後のフレームと最初のフレームの繋がりが良くなります。

3次元の場合

　続いて、3次元空間上でベクトルのなす角を可視化します。

　3次元ベクトルと係数を指定します。

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

# 係数を指定
alpha = 1.0
beta = 1.0

　 $\mathbf{a} = (a_1, a_2, a_3)^{\top}$ をa、 $\mathbf{b} = (b_1, b_2, b_3)^{\top}$ をbとして値を指定します。

　2つのベクトルのなす角を計算します。

# ベクトルa,bのなす角を計算
dot_ab = np.dot(alpha*a, beta*b)
norm_a = np.linalg.norm(alpha * a)
norm_b = np.linalg.norm(beta * b)
theta = np.arccos(dot_ab / norm_a / norm_b)
print(theta) # ラジアン
print(theta / np.pi) # ラジアンの目安
print(theta * 180.0 / np.pi) # 角度

1.9106332362490186
0.6081734479693928
109.47122063449069

　「2次元の場合」のコードで計算できます。

　参考として、xy面(グラフの底面)にx軸とy軸におけるなす角を示す角マークを描画するための座標を作成します。

# xy面上のベクトルaのなす角を計算
sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
theta_a_xy = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
print(theta_a_xy / np.pi)

# xy面上のベクトルbのなす角を計算
sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
theta_b_xy = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
print(theta_b_xy / np.pi)

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

# xy面上の角マーク用の値(ラジアン)を作成
rad_xy_vals = np.linspace(start=theta_a_xy, stop=theta_b_xy, num=100)

# xy面上の角マーク用の座標を計算
r = 0.2
angle_xy_x_vals = np.cos(rad_xy_vals) * r
angle_xy_y_vals = np.sin(rad_xy_vals) * r

0.25000000000000006
0.75
0.25000000000000006
0.75

　「2次元の場合」のときと同様に処理します。

　同様に、xz面(デフォルトのグラフの奥面)にx軸とz軸におけるなす角についても作成します。

# xz面上のベクトルaのなす角を計算
sgn_z = 1.0 if alpha*a[2] >= 0.0 else -1.0
theta_a_xz = sgn_z * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 2]]))
print(theta_a_xz / np.pi)

# xz面上のベクトルbのなす角を計算
sgn_z = 1.0 if beta*b[2] >= 0.0 else -1.0
theta_b_xz = sgn_z * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 2]]))
print(theta_b_xz / np.pi)

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

# xz面上の角マーク用の値(ラジアン)を作成
rad_xz_vals = np.linspace(start=theta_a_xz, stop=theta_b_xz, num=100)

# xz面上の角マーク用の座標を計算
r = 0.2
angle_xz_x_vals = np.cos(rad_xz_vals) * r
angle_xz_z_vals = np.sin(rad_xz_vals) * r

-0.25000000000000006
0.75
-0.25000000000000006
0.75

　xy面のときと同様に処理します。

　空間上に、なす角を示す角マークを描画するための座標を作成します。この処理も、なす角の計算自体には不要です。

# 角マーク用の角aの値(球面座標)を計算
psi_a = np.arccos(alpha*a[2] / np.linalg.norm(alpha * a))
sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
phi_a = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
print(psi_a / np.pi)
print(phi_a / np.pi)

# 角マーク用の角bの値(球面座標)を計算
psi_b = np.arccos(beta*b[2] / np.linalg.norm(beta * b))
sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
phi_b = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
print(psi_b / np.pi)
print(phi_b / np.pi)

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

# 角マーク用の値(ラジアン)を作成
rad_psi_vals = np.linspace(start=psi_a, stop=psi_b, num=100)
rad_phi_vals = np.linspace(start=phi_a, stop=phi_b, num=100)

# 角マーク用の座標を計算
r = 0.2
angle_x_vals = np.sin(rad_psi_vals) * np.cos(rad_phi_vals) * r
angle_y_vals = np.sin(rad_psi_vals) * np.sin(rad_phi_vals) * r
angle_z_vals = np.cos(rad_psi_vals) * r

# 角ラベル用の座標を計算
r = 0.4
angle_x = np.sin(np.median(rad_psi_vals)) * np.cos(np.median(rad_phi_vals)) * r
angle_y = np.sin(np.median(rad_psi_vals)) * np.sin(np.median(rad_phi_vals)) * r
angle_z = np.cos(np.median(rad_psi_vals)) * r

0.6959132760153038
0.25000000000000006
0.3040867239846964
0.75
0.6959132760153038
0.25000000000000006

　 $\alpha \mathbf{a}, \beta \mathbf{b}$ それぞれのベクトルとx軸(標準単位ベクトル $\mathbf{e}_1$ )との球面座標系における2つの角度を $\psi_a, \phi_a$ 、 $\psi_b, \phi_b$ とします。詳しくは前回の記事を参照してください。
　角度 $\psi, \phi$ は、それぞれ次の式で計算できます。

$\displaystyle \begin{aligned} \psi &= \arccos \left( \frac{ \beta x_3 }{ \|\beta \mathbf{x}\| } \right) \\ \phi &= \mathrm{sgn}(x_2) \arccos \left( \frac{ \beta x_1 }{ \|\beta (x_1, x_2)\| } \right) \\ \mathrm{sgn}(x) &= \begin{cases} 1 & (x \geq 0) \\ -1 & (x < 0) \end{cases} \end{aligned}$

　ただし、 $\beta \mathbf{x}$ や $\beta (x_1, x_2)$ が0ベクトルだと0除算になるため定義できません(計算結果がnanになります)。

　 $\psi_a$ から $\psi_b$ までのラジアン $t$ と、 $\phi_a$ から $\phi_b$ までのラジアン $u$ を作成して、角マークの座標を計算します。
　ただし、 $\phi_a - \phi_b$ が優角になる場合( $\phi_a, \phi_b$ の一方が正の角度でもう一方が負の角度で、差の絶対値が $\pi$ を超える場合)は、負の角度の方に $2 \pi$ を足して優角に変更します。
　 $t$ の値rad_psi_vals、 $u$ の値rad_phi_valsを作成して、各要素を球形の座標に変換します。
　座標 $x, y, z$ は、半径を $r$ として、それぞれ次の式で計算できます。

$\displaystyle \left\{ \begin{aligned} x &= \sin t \cos u \\ y &= \sin t \sin u \\ z &= \cos t \end{aligned} \right.$

　半径rで角マークのサイズを調整します。

　角ラベルは、rad_*_valsの中点の座標に表示することにします。

　3次元空間上にベクトルのなす角を作図します。

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

# 3Dベクトルのなす角を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.quiver(0, 0, 0, *alpha*a, color='red', arrow_length_ratio=0.1, 
          label='$\\alpha a = \\alpha ('+', '.join(map(str, a.round(2)))+')$') # ベクトルαa
ax.quiver(0, 0, 0, *beta*b, color='blue', arrow_length_ratio=0.1, 
          label='$\\beta x = \\beta ('+', '.join(map(str, b.round(2)))+')$') # ベクトルβb
ax.quiver(0, y_max, 0, alpha*a[0], 0, alpha*a[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルαa
ax.quiver(0, 0, z_min, alpha*a[0], alpha*a[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルαa
ax.quiver(0, y_max, 0, beta*b[0], 0, beta*b[2], color='blue', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルβb
ax.quiver(0, 0, z_min, beta*b[0], beta*b[1], 0, color='blue', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルβb
ax.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
        color='black', linewidth=1) # xy面上の角マーク
ax.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
        color='black', linewidth=1) # xz面上の角マーク
ax.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
ax.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
ax.quiver([0, alpha*a[0], beta*b[0], 0, alpha*a[0], beta*b[0]], 
          [0, alpha*a[1], beta*b[1], y_max, y_max, y_max], 
          [z_min, z_min, z_min, 0, alpha*a[2], beta*b[2]], 
          [0, 0, 0, 0, 0, 0], 
          [0, 0, 0, -y_max, alpha*a[1]-y_max, beta*b[1]-y_max], 
          [-z_min, alpha*a[2]-z_min, beta*b[2]-z_min, 0, 0, 0], 
          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_ylim(y_min, y_max)
ax.set_zlim(z_min, z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$\\alpha='+str(alpha)+', ' + 
             '\\beta='+str(beta)+', ' + 
             '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
             '\\theta^\circ='+str((theta * 180.0 / np.pi).round(1))+'$', loc='left')
fig.suptitle('angle a, b', fontsize=20)
ax.legend(loc='upper left')
ax.set_aspect('equal')
#ax.view_init(elev=90, azim=270) # xy面
#ax.view_init(elev=0, azim=270) # xz面
plt.show()

　axes.quiver()の第1・2・3引数に始点の座標、第4・5・6引数にベクトルのサイズを指定します。原点を始点とします。

　係数の値を変えてみます。

　ベクトルの値を変化させたアニメーションを作成します。

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

# フレーム数を設定
frame_num = 300

# ベクトルとして利用する値を指定
psi_a_vals = np.linspace(start=0.5*np.pi, stop=2.5*np.pi, num=frame_num+1)[:frame_num]
phi_a_vals = np.linspace(start=-2.0*np.pi, stop=4.0*np.pi, num=frame_num+1)[:frame_num]
a_vals = np.array(
    [np.sin(psi_a_vals) * np.cos(phi_a_vals), 
     np.sin(psi_a_vals) * np.sin(phi_a_vals), 
     np.cos(psi_a_vals)]
).T
psi_b_vals = np.linspace(start=-3.5*np.pi, stop=1.5*np.pi, num=frame_num+1)[:frame_num]
phi_b_vals = np.linspace(start=-2.5*np.pi, stop=2.5*np.pi, num=frame_num+1)[:frame_num]
b_vals = np.array(
    [-np.sin(psi_b_vals) * np.cos(phi_b_vals), 
     np.sin(psi_b_vals) * np.sin(phi_b_vals), 
     -np.cos(psi_b_vals)]
).T

# 係数を指定
alpha = 1.0
beta = 1.0

# 作図用の値を設定
x_min = np.floor(np.min([0.0, *alpha*a_vals[:, 0], *beta*b_vals[:, 0]])) - 1
x_max = np.ceil(np.max([0.0, *alpha*a_vals[:, 0], *beta*b_vals[:, 0]])) + 1
y_min = np.floor(np.min([0.0, *alpha*a_vals[:, 1], *beta*b_vals[:, 1]])) - 1
y_max = np.ceil(np.max([0.0, *alpha*a_vals[:, 1], *beta*b_vals[:, 1]])) + 1
z_min = np.floor(np.min([0.0, *alpha*a_vals[:, 2], *beta*b_vals[:, 2]])) - 1
z_max = np.ceil(np.max([0.0, *alpha*a_vals[:, 2], *beta*b_vals[:, 2]])) + 1

# 作図用のオブジェクトを初期化
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
fig.suptitle('angle a, b', fontsize=20)

# 作図処理を関数として定義
def update(i):
    
    # 前フレームのグラフを初期化
    plt.cla()
    
    # i番目のベクトルを取得
    a = a_vals[i]
    b = b_vals[i]
    
    # ベクトルa,bのなす角を計算
    dot_ab = np.dot(alpha*a, beta*b)
    norm_a = np.linalg.norm(alpha * a)
    norm_b = np.linalg.norm(beta * b)
    theta = np.arccos(dot_ab / norm_a / norm_b)
    
    # xy面・xz面上のベクトルaのなす角を計算
    sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
    theta_a_xy = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
    sgn_z = 1.0 if alpha*a[2] >= 0.0 else -1.0
    theta_a_xz = sgn_z * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 2]]))
    
    # xy面・xz面上のベクトルbのなす角を計算
    sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
    theta_b_xy = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
    sgn_z = 1.0 if beta*b[2] >= 0.0 else -1.0
    theta_b_xz = sgn_z * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 2]]))
    
    # 差が180°を超える場合は、負の角度を優角に変換
    if abs(theta_a_xy - theta_b_xy) > np.pi:
        theta_a_xy = theta_a_xy if alpha*a[1] >= 0.0 else theta_a_xy + 2.0*np.pi
        theta_b_xy = theta_b_xy if beta*b[1] >= 0.0 else theta_b_xy + 2.0*np.pi
    if abs(theta_a_xz - theta_b_xz) > np.pi:
        theta_a_xz = theta_a_xz if alpha*a[2] >= 0.0 else theta_a_xz + 2.0*np.pi
        theta_b_xz = theta_b_xz if beta*b[2] >= 0.0 else theta_b_xz + 2.0*np.pi
    
    # xy・xz面面上の角マーク用の値(ラジアン)を作成
    rad_xy_vals = np.linspace(start=theta_a_xy, stop=theta_b_xy, num=100)
    rad_xz_vals = np.linspace(start=theta_a_xz, stop=theta_b_xz, num=100)

    # xy面・xz面上の角マーク用の座標を計算
    r = 0.2
    angle_xy_x_vals = np.cos(rad_xy_vals) * r
    angle_xy_y_vals = np.sin(rad_xy_vals) * r
    angle_xz_x_vals = np.cos(rad_xz_vals) * r
    angle_xz_z_vals = np.sin(rad_xz_vals) * r
    
    # 角マーク用の角aの値(球面座標)を計算
    psi_a = np.arccos(alpha*a[2] / np.linalg.norm(alpha * a))
    sgn_y = 1.0 if alpha*a[1] >= 0.0 else -1.0
    phi_a = sgn_y * np.arccos(alpha*a[0] / np.linalg.norm(alpha * a[[0, 1]]))
    
    # 角マーク用の角bの値(球面座標)を計算
    psi_b = np.arccos(beta*b[2] / np.linalg.norm(beta * b))
    sgn_y = 1.0 if beta*b[1] >= 0.0 else -1.0
    phi_b = sgn_y * np.arccos(beta*b[0] / np.linalg.norm(beta * b[[0, 1]]))
    
    # 差が180°を超える場合は、負の角度を優角に変換
    if abs(phi_a - phi_b) > np.pi:
        phi_a = phi_a if alpha*a[1] >= 0.0 else phi_a + 2.0*np.pi
        phi_b = phi_b if beta*b[1] >= 0.0 else phi_b + 2.0*np.pi
    
    # 角マーク用の値(ラジアン)を作成
    rad_psi_vals = np.linspace(start=psi_a, stop=psi_b, num=100)
    rad_phi_vals = np.linspace(start=phi_a, stop=phi_b, num=100)

    # 角マーク用の座標を計算
    r = 0.2
    angle_x_vals = np.sin(rad_psi_vals) * np.cos(rad_phi_vals) * r
    angle_y_vals = np.sin(rad_psi_vals) * np.sin(rad_phi_vals) * r
    angle_z_vals = np.cos(rad_psi_vals) * r

    # 角ラベル用の座標を計算
    r = 0.4
    angle_x = np.sin(np.median(rad_psi_vals)) * np.cos(np.median(rad_phi_vals)) * r
    angle_y = np.sin(np.median(rad_psi_vals)) * np.sin(np.median(rad_phi_vals)) * r
    angle_z = np.cos(np.median(rad_psi_vals)) * r
    
    # 3Dベクトルのなす角を作図
    ax.quiver(0, 0, 0, *alpha*a, color='red', arrow_length_ratio=0.1, 
              label='$\\alpha a = \\alpha ('+', '.join(map(str, a.round(2)))+')$') # ベクトルαa
    ax.quiver(0, 0, 0, *beta*b, color='blue', arrow_length_ratio=0.1, 
              label='$\\beta x = \\beta ('+', '.join(map(str, b.round(2)))+')$') # ベクトルβb
    ax.quiver(0, y_max, 0, alpha*a[0], 0, alpha*a[2], color='red', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルαa
    ax.quiver(0, 0, z_min, alpha*a[0], alpha*a[1], 0, color='red', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルαa
    ax.quiver(0, y_max, 0, beta*b[0], 0, beta*b[2], color='blue', linestyle='--', arrow_length_ratio=0.0) # xy面上のベクトルβb
    ax.quiver(0, 0, z_min, beta*b[0], beta*b[1], 0, color='blue', linestyle='--', arrow_length_ratio=0.0) # xz面上のベクトルβb
    ax.plot(angle_xy_x_vals, angle_xy_y_vals, np.repeat(z_min, len(angle_xy_x_vals)), 
            color='black', linewidth=1) # xy面上の角マーク
    ax.plot(angle_xz_x_vals, np.repeat(y_max, len(angle_xz_x_vals)), angle_xz_z_vals, 
            color='black', linewidth=1) # xz面上の角マーク
    ax.plot(angle_x_vals, angle_y_vals, angle_z_vals, color='black', linewidth=1) # 角マーク
    ax.text(angle_x, angle_y, angle_z, s='$\\theta$', size=15, ha='center', va='center') # 角ラベル
    ax.quiver([0, alpha*a[0], beta*b[0], 0, alpha*a[0], beta*b[0]], 
              [0, alpha*a[1], beta*b[1], y_max, y_max, y_max], 
              [z_min, z_min, z_min, 0, alpha*a[2], beta*b[2]], 
              [0, 0, 0, 0, 0, 0], 
              [0, 0, 0, -y_max, alpha*a[1]-y_max, beta*b[1]-y_max], 
              [-z_min, alpha*a[2]-z_min, beta*b[2]-z_min, 0, 0, 0], 
              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_ylim(y_min, y_max)
    ax.set_zlim(z_min, z_max)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_title('$\\alpha='+str(alpha)+', ' + 
                 '\\beta='+str(beta)+', ' + 
                 '\\theta='+str((theta / np.pi).round(2))+'\pi, ' + 
                 '\\theta^\circ='+str((theta * 180.0 / np.pi).round(1))+'$'+'\n'+
                 '$\psi_a='+str((psi_a/np.pi).round(2))+'\pi, \psi_b='+str((psi_b/np.pi).round(2))+
                 '\pi, \phi_a='+str((phi_a/np.pi).round(2))+'\pi, \phi_b='+str((phi_b/np.pi).round(2))+'\pi$', loc='left')
    ax.legend(loc='upper left')
    ax.set_aspect('equal')
    #ax.view_init(elev=90, azim=270) # xy面
    #ax.view_init(elev=0, azim=270) # xz面

# gif画像を作成
ani = FuncAnimation(fig=fig, func=update, frames=frame_num, interval=100)

# gif画像を保存
ani.save('angle_ab_3d.gif')

　「2次元の場合」と同様にして、アニメーションを作成します。

　ちなみに、点 $\alpha \mathbf{a}, \beta \mathbf{b}$ (矢印の先端)は、次の軌道をしています。

# ベクトルの軌道を作図
fig, ax = plt.subplots(figsize=(8, 8), facecolor='white', 
                       subplot_kw={'projection': '3d'})
ax.plot(*alpha*a_vals.T, color='red', label='$\\alpha a$') # ベクトルαaの軌道
ax.plot(*beta*b_vals.T, color='blue', label='$\\beta b$') # ベクトルβbの軌道
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_ylim(y_min, y_max)
ax.set_zlim(z_min, z_max)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title('$\\beta='+str(beta)+'$', loc='left')
fig.suptitle('trace $\\alpha a, \\beta b$', fontsize=20)
ax.legend()
ax.set_aspect('equal')
plt.show()

　a_vals, b_valsを転置して、次元ごとにアンパックして描画します。

　この記事では、2つのベクトルのなす角を作図しました。次の記事では、相関係数を考えます。

参考書籍

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

おわりに

　前回にも増して装飾用の処理が氾濫しています。必要なところだけ読んでください。
　苦労した割に、意図通りに角度マーク(正しい呼び方が分からない)を描画できなくて残念です。真上や真横から見ると、変な方向にカーブしているのが分かると思います。
　2次元だと円で描けるのだから、3次元だと球で描けるはずなのですが上手くできませんでした。大円なるものを計算できればいいようなのですが、今回は撤退しました。行列の章まで進めば、回転行列なるものでも実現できるかもしれません。
　現在は一周回って2次元的に見れば今の方が分かりやすいのかもしれないと思ったり思わなかったりしてます。

　2023年3月16日は、BEYOOOOONDSの高瀬くるみさんの24歳のお誕生日です！