三維空間中如何尋找和一組給定直線垂直程度最高的直線

這是個挺有意思的小問題，給定一組直線(至少兩條不平行)，希望能找到和這組直線盡可能垂直的直線。打個比方，比如在三維空間中，如下圖(forked from wiki)

a和b分別是在一個平面上不平行的兩條直線上，那么顯而易見與a和b所在直線垂直程度最高的就是與a和b倆倆垂直的豎線，也就是叉積axb方向平行的直線。兩條直線可以用叉積，那么多於兩條的情況呢？想象如果又有了一條直線，其一端的方向可以用向量c表示，如果c和a、b在一個平面上的話，那么好辦，還是axb，或者bxc也行；但是如果c和a、b不在一個平面上的話，那么叉積的辦法就不管用了。回到問題的最開始，其實是如何定義和其他直線最大程度垂直。還是考慮只有a和b的情況，axb其實就是在a和b上的投影均為零，也就是點積為零，所以不妨把在其他所有向量上投影最少的向量定義為最大程度垂直的方向。

那么考慮有n條直線，對應n個單位向量代表每條直線的一個方向，其中第i個向量可以表示為\({{v}_{i}}=\left( {{x}_{i}},{{y}_{i}},{{z}_{i}} \right)\)，那么對於一個向量\(v=\left( x,y,z \right)\)，在\({{v}_{i}}\)上投影的長度為：

\[\begin{align}
 & \left| \left\langle v,{{v}_{i}} \right\rangle  \right|=\left| x{{x}_{i}}+y{{y}_{i}}+z{{z}_{i}} \right| \\
\end{align}\]

那么一個很直接的辦法就是用優化算法來求出使投影絕對值之和最小的v：

\(\text{min: }\sum\limits_{i=1}^{n}{\left| \left\langle v,{{v}_{i}} \right\rangle  \right|}\)
\(\begin{align}
 & \text{subject to: }\left| {{v}_{i}} \right|=1\text{ and }\left| v \right|=1 \\
\end{align}\)

\(\left| {{v}_{i}} \right|=1\)是為了是每一個向量上的投影能夠相加比較，\(\left| v \right|=1\)是為了優化，否則最后的結果就是\(\left( 0,0,0 \right)\)。這個方法固然直觀，但是考慮到1)要利用絕對值相加，2)要利用優化算法還是覺得有些麻煩。考慮一下用平方和作為優化目標，當然這會和絕對值的辦法有些許不一樣，主要是每一個投影大小對結果影響的sensitivity變了，不過總體而言是差不多的，總之，要求的目標變成了：

\( \sum\limits_{i=1}^{n}{{{\left\langle v,{{v}_{i}} \right\rangle }^{2}}} \)
\( =\sum\limits_{i=1}^{n}{{{\left( x{{x}_{i}}+y{{y}_{i}}+z{{z}_{i}} \right)}^{2}}} \)
\( =\sum\limits_{i=1}^{n}{{{x}^{2}}{{x}_{i}}^{2}+{{y}^{2}}{{y}_{i}}^{2}+{{z}^{2}}{{z}_{i}}^{2}+2xy{{x}_{i}}{{y}_{i}}+2yz{{y}_{i}}{{z}_{i}}+2zx{{z}_{i}}{{x}_{i}}} \)
\( ={{x}^{2}}\sum\limits_{i=1}^{n}{{{x}_{i}}^{2}}+{{y}^{2}}\sum\limits_{i=1}^{n}{{{y}_{i}}^{2}}+{{z}^{2}}\sum\limits_{i=1}^{n}{{{z}_{i}}^{2}}+xy\sum\limits_{i=1}^{n}{2{{x}_{i}}{{y}_{i}}}+yz\sum\limits_{i=1}^{n}{2{{y}_{i}}{{z}_{i}}}+zx\sum\limits_{i=1}^{n}{2{{z}_{i}}{{x}_{i}}} \)
\(\begin{align}
&\\
\end{align}\)

注意觀察最后的展開，如果右邊加個等於常數項的話，就是個橢球嘛，而且沒有一次項，所以是個中心在原點的橢球。這樣一來就很簡單了，考慮之前提出的限制條件\(\left| v \right|=1\)，相當於有一個半徑為1，球心在原點的球面，和一個中心也在原點，大小未定，但是三個極軸比例確定的橢球面。我們要尋找這么一種情況：球面和橢球面有交點(v存在)，並且橢球的軸長盡可能小(\(\sum\limits_{i=1}^{n}{{{\left\langle v,{{v}_{i}} \right\rangle }^{2}}}\)最小)。很顯然這種情況就是橢球內切於單位球面，而此時v對應的就是長軸的方向，所以問題就更簡單了，我們都不需要求出橢球的具體表達式，而是找出長軸的方向向量就可以了。回顧一下橢球的知識：任意一個橢球可以表達如下：

\(\begin{align}
 & {{\left( v-{{v}_{0}} \right)}^{T}}M\left( v-{{v}_{0}} \right)=1 \\
\end{align}\)

其中v0是中心坐標，當然對於本文的情況這是個零向量，M的特征向量就是三個極軸的方向，特征值就是三個軸半長的平方分之一，所以我們只要求出M，然后找到M對應特征值絕對值最小的那個特征向量，該向量就是我們最終要求的方向向量v。所以，將橢球表達式展開，為簡便，首先假設M如下：

\(\begin{align}
&M=\left[ \begin{matrix}
   A & D & F  \\
   D & B & E  \\
   F & E & C  \\
\end{matrix} \right]\\
\end{align}\)

於是有：

\( {{v}^{T}}Mv \)
\( ={{\left[ \begin{matrix}
   x  \\
   y  \\
   z  \\
\end{matrix} \right]}^{T}}\left[ \begin{matrix}
   A & D & F  \\
   D & B & E  \\
   F & E & C  \\
\end{matrix} \right]\left[ \begin{matrix}
   x  \\
   y  \\
   z  \\
\end{matrix} \right] \)
\( =A{{x}^{2}}+B{{y}^{2}}+C{{z}^{2}}+2Dxy+2Eyz+2Fzx \)
\[\begin{align}
&\\
\end{align}\]

把這個結果和前面的公式(3)對照，得到A、B、C、D、E和F的值，代入M，得到如下：

\(\begin{align}
&M=\left[ \begin{matrix}
   \sum\limits_{i=1}^{n}{{{x}_{i}}^{2}} & \sum\limits_{i=1}^{n}{{{x}_{i}}{{y}_{i}}} & \sum\limits_{i=1}^{n}{{{z}_{i}}{{x}_{i}}}  \\
   \sum\limits_{i=1}^{n}{{{x}_{i}}{{y}_{i}}} & \sum\limits_{i=1}^{n}{{{y}_{i}}^{2}} & \sum\limits_{i=1}^{n}{{{y}_{i}}{{z}_{i}}}  \\
   \sum\limits_{i=1}^{n}{{{z}_{i}}{{x}_{i}}} & \sum\limits_{i=1}^{n}{{{y}_{i}}{{z}_{i}}} & \sum\limits_{i=1}^{n}{{{z}_{i}}^{2}}  \\
\end{matrix} \right]\\
\end{align}\)

對這個矩陣求絕對特征值最小的特征向量，就得到了最大程度垂直於給定直線的向量了。需要注意的是當所有給定直線在一個平面時，M不滿秩，不過這種情況恰好對應絕對值最小的特征值就是0，所以應用到這個算法里還是沒有問題。這個辦法還可以類似地推廣到二維或者更高維，只不過對三維以上的情況矩陣不滿秩的處理會比較麻煩一些。

下面是三維情況下的一段Python驗證程序：

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Distance from x to the plane with v as the normal vector, passing through (0, 0, 0)
d = lambda x, v: abs(x[0]*v[0]+x[1]*v[1]+x[2]*v[2]) / np.linalg.norm(v)

# Random normal vector for test
norm_v_plane = np.random.randn(3)
N = 200

# Generate test unit vectors that within 0.2 to the test plane
v_test = [x for x in [x/np.linalg.norm(x) for x in np.random.randn(N, 3)] if d(x, norm_v_plane) < 0.2]

# Plot test vectors
fig = plt.figure('Minimum dot product vector')
ax = fig.add_subplot(111, projection='3d')
for v in v_test:
    ax.plot([0, v[0]], [0, v[1]], [0, v[2]], 'b')

# Calculate the coefficients matrix of the arbitrary ellipsoid
A, B, C, D, E, F = [0] * 6
for v in v_test:
    x, y, z = v
    A += x * x
    B += y * y
    C += z * z
    D += x * y
    E += y * z  
    F += x * z
M = np.array([[A, D, F],
              [D, B, E], 
              [F, E, C]])

# Pick the eigenvector with the minimum egienvalue
u, v = np.linalg.eig(M)
vp = v[:, np.argmin(np.abs(u))]
if vp[2] < 0:
    vp = -vp

ax.plot([0, vp[0]], [0, vp[1]], [0, vp[2]], 'r')

plt.show()

 

這段小程序會隨機選取一個過原點的平面，然后在距離平面0.2的距離范圍內產生一些隨機向量，然后找到最大程度垂直於這些向量的向量，比如下圖：