读Matlab7.7的subspace函数（一）——两向量的夹角

我觉得这是一个名不副实的函数，名为“子空间”，实际作用是求解两个子空间的夹角（Angle between  two subspace）。所谓两个子空间的夹角，对于不怎么接触矩阵理论的人，听起来可能有点费解。其实说简单一点，像空间几何中，两个向量的夹角，向量与平面的夹角，两个平面的夹角等都属于子空间夹角问题的范畴。只是这里的subspace函数能够求解任意维的两个子空间的夹角。不局限于三维空间。也就是说，subspace所使用的算法，是一种通用的算法。而像空间几何中使用的求两个向量夹角的方法，利用内积的特性，利用反余弦函数求解，并不能应用于向量和平面之间夹角的计算，所以不是通用的算法。根据subspace的文档，文献Numerical methods for computing angles between linear subspaces和On angles between subspaces of a finite dimensional inner product space指出了subspace所使用的算法的理论依据。我看过其中一篇，不过就看了几行。

     可以通过一些严格的数学方法说明subspace算法的原理，但是这里先不打算这样做，而是首先将几种特例拿来分析。这里不贴完整的subspace.m的代码，只要打开Matlab，输入edit subspace就可以看了。

     最简单的特例就是两个一维子空间的夹角，即两个向量之间的夹角。当subspace的两个参数都是向量时，代码一开始的
A = orth(A);
B = orth(B);
将会得到两个新的向量，这两个向量的长度为1，新的A与原来的A同向或反向，新的B也是这样。

http://xqnu7q.blu.livefilestore.com/y1pbvT4lMG8vigM7BnLij3qjMv5BhdEUhy1187VQ5EMxW1cGvVPQgq6biAX-sS0BzC3JVqEsDrA6XBMjazq75qqc50L5tlmyDu_/20100915_1.jpg

     如上图，我们用向量OC代表向量A，OD代表向量B。

     因为A和B都代表一维子空间（向量），所以代码中的循环退化成：
B = B - A * (A' * B);
这是比较关键的一个步骤。代码的精华就在这里。

     首先是A' * B，即计算A和B的内积。因为此时A和B的长度都为1，因此A' * B等于A到B上投影的长度，即OF的长度（CF为C到OD的垂线）。

     然后是A * (A' * B)，这样会得到向量OE，它的方向与A（OC）同，长度与OF同，那么现在看看两个三角形OFC和OED，它们两个是≌的（SAS）。因为CF垂直于OD，所以DE也垂直于OC，OED为直角三角形。

     最后，B - A * (A' * B)，即向量OD减去向量OE，结果为向量ED。

     往下看代码，最后一行可能有点费解，其实去掉基于代码稳定性的考虑，完全可以简化为：
theta = asin(norm(B));
即向量ED的长度的反正弦值。因为三角形为直角三角形，而OD的长度为1，所以ED长度的反正弦值就是∠COD的角度。COD就是两个向量的夹角。