In \(D\)-dimensions, \(p_1\)-plane and \(p_2\)-plane has intersection of \((D-p_1-p_2)\)-plane at least or with more dimensions. If \(D-p_1-p_2 <0\), they can have 0 to \(\mathrm{min}\{p_1,p_2\}\) dimension proximity object or intersections. (i.e. no intersection guaranteed). This is because parameter of each dimension of plane gives constraint to the \(D\)-dimensional coordinate equation.
The two objects may be written by
\[\vec{x}(a) = \vec{v}_0 + \sum_{k=1}^{p_1}a_k \vec{v}_k,\]
\[\vec{y}(b) = \vec{w}_0 + \sum_{k=1}^{p_2}b_k \vec{w}_k.\]
The proximity or intersection maybe given by
\[\mathrm{min}\left|\vec{x} - \vec{y}\right|,\]
which may be represented by this extrema equation:
\[\begin{cases}
0=\frac{\partial}{\partial a_k}\left|\vec{x} - \vec{y}\right|^2 = 2\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial}{\partial a_k}\left(\vec{x} - \vec{y}\right)= 2\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial \vec{x}}{\partial a_k}
\\
0\frac{\partial}{\partial b_k}\left|\vec{x} - \vec{y}\right|^2 = 2\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial}{\partial b_k}\left(\vec{x} - \vec{y}\right) = -2\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial \vec{y}}{\partial b_k}
\end{cases}\]
\[\Rightarrow\begin{cases}
0=\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial \vec{x}}{\partial a_k}=\vec{v}_k\cdot\left(\vec{x} - \vec{y}\right)
\\
0=\left(\vec{x} - \vec{y}\right)\cdot\frac{\partial \vec{y}}{\partial b_k}=\vec{w}_k\cdot\left(\vec{x} - \vec{y}\right)
\end{cases}\]
\[\Rightarrow\begin{cases}
0=\vec{v}_k\cdot\left(\vec{x} - \vec{y}\right) =\vec{v}_k \cdot\left( \vec{v}_0 - \vec{w}_0\right) + \sum_{l=1}^{p_1}\vec{v}_k\cdot\vec{v}_l a_l - \sum_{l=1}^{p_1}\vec{v}_k\cdot\vec{w}_l b_l
\\
0=-\vec{w}_k\cdot\left(\vec{x} - \vec{y}\right) =-\vec{w}_k \cdot\left( \vec{v}_0 - \vec{w}_0\right) - \sum_{l=1}^{p_2}\vec{w}_k\cdot\vec{v}_l a_l+ \sum_{l=1}^{p_2}\vec{w}_k\cdot\vec{w}_l b_l
\end{cases}\]
The equations may be re-written by
\[\begin{pmatrix}*&*&*\\\vec{v}_k\cdot\left( \vec{v}_0 - \vec{w}_0\right)&\vec{v}_k\cdot\vec{v}_l&-\vec{v}_k\cdot\vec{w}_l \\-\vec{w}_k\cdot\left( \vec{v}_0 - \vec{w}_0\right)&-\vec{w}_k\cdot\vec{v}_l&\vec{w}_k\cdot\vec{w}_l\end{pmatrix}\begin{pmatrix}1 \\ a_l \\ b_l \end{pmatrix}=\begin{pmatrix}* \\ 0 \\ 0 \end{pmatrix}\]
or
\[\begin{pmatrix}\vec{v}_k\cdot\vec{v}_l&-\vec{v}_k\cdot\vec{w}_l \\-\vec{w}_k\cdot\vec{v}_l&\vec{w}_k\cdot\vec{w}_l\end{pmatrix}\begin{pmatrix}a_l \\ b_l \end{pmatrix}=\begin{pmatrix} -\vec{v}_k\cdot\left( \vec{v}_0 - \vec{w}_0\right) \\ \vec{w}_k\cdot\left( \vec{v}_0 - \vec{w}_0\right) \end{pmatrix}\]
With convention of column vector \(\vec{v}\)'s, and
\[{\mathbf{m}_0}_{(D \times 1)} := \vec{v}_0-\vec{w}_0 ,\quad \mathbf{M}_{(D \times (p_1 + p_2))} := \begin{pmatrix} \vec{v}_k & -\vec{w}_k \end{pmatrix},\quad \mathbf{a}_{((p_1 + p_2 )\times 1)}=\begin{pmatrix}a_k \\ b_k \end{pmatrix},\]
the equations will be re-written by
\[ (\mathbf{M}^T \mathbf{M}) \mathbf{a}= -\mathbf{M}^T \mathbf{m}_0. \]
So, the solution of parameters is
\[\mathbf{a}= -(\mathbf{M}^T \mathbf{M})^{-1} \mathbf{M}^T \mathbf{m}_0\]
unless \(\mathbf{M}^T \mathbf{M}\) is singular.
If singular, to separate singular (freed) and non-singular (constraint) components, we need to diagonalize it. We may try SVD to avoid catastrophic result.
\[(\mathbf{M}^T \mathbf{M}) =: \mathbf{U}\mathbf{\Sigma} \mathbf{V}^T\]
where \(U\) and \(V\) are orthogonal (Hermitian) and \(\Sigma\) is a nonnegative diagonal matrix.
The parameters will be rotated by the orthogonal matrices to diagonalize the matrix to see the free variable.
\[\mathbf{U}\mathbf{\Sigma} \mathbf{V}^T\mathbf{a}= - \mathbf{M}^T \mathbf{m}_0 \;\Rightarrow \;\mathbf{\Sigma} \mathbf{V}^T\mathbf{a}= -\mathbf{U}^T \mathbf{M}^T \mathbf{m}_0=: \begin{pmatrix}m_r \\ 0\end{pmatrix}\]
where in the rows with the zero values in \(\Sigma\), the vector component is also zeros, so the solution for that row is indefinite. We know indefinite because it is zero divide by zero.
\[ \mathbf{V}^T\mathbf{a}=\mathbf{\Sigma}^{-1} \begin{pmatrix}m_r \\ z\end{pmatrix} = \begin{pmatrix}\Sigma_r^{-1} m_r \\ 0^{-1} \cdot 0\end{pmatrix}\].
Now, substitute the indefinite numbers to free parameters \(z_i\).
\[ \mathbf{V}^T\mathbf{a}= \begin{pmatrix}\Sigma_r^{-1} m_r \\ z_i\end{pmatrix}\;\Rightarrow \;\mathbf{a}= \mathbf{V}\begin{pmatrix}\Sigma_r^{-1} m_r \\ z_i\end{pmatrix}.\]
Also, we can find the shortest distance using the equation and the solution. The squared distance is
\[\left|\vec{x} - \vec{y}\right|^2 = \left| \mathbf{m}_0 + \mathbf{M} \mathbf{a}\right|^2 = \mathbf{m}_0^T \mathbf{m}_0 +2 \mathbf{a}^T \mathbf{M}^T \mathbf{m}_0 +\mathbf{a}^T \mathbf{M}^T\mathbf{M}\mathbf{a} \]
\[= \mathbf{m}_0^T \mathbf{m}_0 + \mathbf{a}^T \mathbf{M}^T \mathbf{m}_0 = \mathbf{m}_0^T \left(\mathbf{I} - \mathbf{M}(\mathbf{M}^T\mathbf{M})^{-1} \mathbf{M}^T\right)\mathbf{m}_0\]
