如何用线性代数优化矩阵快速幂

· · 个人记录

如何用线性代数优化矩阵快速幂

如果我们需要求维度为 c 的矩阵 A 的 m 次幂，那么一个显然的方法是用快速幂，复杂度为 c^3 \log m。不过当需要用到高精度的时候，该算法速度会大大下降。

考虑这么一个事情，如果可以找到一个矩阵 C 和一个对角矩阵 P（只有主对角线上的值不为 0），使得 A=C\times P \times C^{-1}，那么 A^m=C\times P \times C^{-1} \times C \times P \times C^{-1}...，不断的使用结合率，就可以得到 C \times P^m \times C^{-1}。因为 P 是对角矩阵，那么 P^m 本质上只要对 P 中的每个值跑一个实数的快速幂。这样就可以把复杂度优化到 c^3+c\log n，即实数的快速幂和两次矩阵乘法。

现在的问题就是怎么得到一个合法的 C。

引入两个概念：特征值和特征向量。

一个数 \lambda 是矩阵 A 的特征值，当且仅当存在一个非 0 向量 B，使得 A\times B=\lambda B，我们把这个向量 B 称作特征向量。需要注意的是，对于同一个 \lambda，B 可能不止一个。

如果我们求出了一组特征值 \lambda_1,\lambda_2...\lambda_c，以及对应的特征向量 B_1,B_2...B_c，那么就有 A\times [B_1,B_2...B_c]=[\lambda_1\cdot B_1,\lambda_2\cdot B_2...\lambda_c\cdot B_c]，设 C=[B_1,B_2...B_c]，右乘一个对角矩阵相当于每一列 i 乘上一个对角阵中的元素 a_{i,i}，所以可以构造对角矩阵 P=\begin{bmatrix} \lambda_1 & 0 & \cdots&0\\ 0 & \lambda_2 & \cdots&0\\ \vdots & \vdots&\ddots&\vdots \\ 0&\cdots&0&\lambda_c\end{bmatrix}，就可以得到 A\times C=C\times P\Leftrightarrow A=C\times P\times C^{-1} 。

考虑怎么求特征值，因为 A\times B=\lambda B 那么就有 \lambda B-A\times B=0，所以 (\lambda I-A)\times B=0，因为 B 向量非 0，所以必须有 \det(\lambda I-A)=0，那么可以通过行列式列方程求得 \lambda。

考虑使用反证法证明 \det(\lambda I-A)=0。设 G=\lambda I-A，那么如果 \det(G)\not=0，那么如果把 G 看作是一个方程组的系数矩阵，B 为未知数，那么 B 将会是唯一的。而因为等式右边等于 0 ，所以无论如何用高斯消元，最后总会有 B_c=0，那么一直往回迭代就会得到一个全 0 向量，与 B 向量非 0 冲突。而如果 det(G)=0，那么 G 至少有一个自由元，那么就可以通过自由元，得到至少一个非 0 向量 B。

由此我们就可以求出所有的特征值，对于某个特征值，把它带回矩阵，通过初等矩阵的行变换把它变成一个最简形式，然后把自由元带入任意一个值（一般为 1）就可以得到一个特征向量 B（当然对于任意的 k\cdot B,k\not=0 都是符合要求的）。

至此就可以得到所有的特征值及其对应的特征向量。那么就可以得到矩阵 C，求个逆就可以得到 C^{-1}（消元可以用辗转相除法，就可以避免求逆元），就可以用上述做法优化了。

下面就举个栗子，用此方法得到斐波那契数列的通项公式。

众所周知，斐波那契数列的递推矩阵是 A=\begin{bmatrix} 1 & 1 \\ 1 & 0 \\ \end{bmatrix} 。

先求出该矩阵的特征值，设其为 \lambda，根据定义，有 \det \left ( \begin{bmatrix} \lambda - 1 & -1 \\ -1 & \lambda \\ \end{bmatrix} \right) = 0 ，所以 (\lambda-1)\times \lambda-1=0，解出 \lambda=\frac{1\pm\sqrt{5}}{2} 。

再根据定理，列出方程 \begin{bmatrix} \lambda - 1 & -1 \\ -1 & \lambda \\ \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ \end{bmatrix}=\begin{bmatrix} 0 \\ 0 \\ \end{bmatrix} 可以解出两个特征向量 \lambda_1=\frac{1-\sqrt{5}}{2} 对应 B_1=\begin{bmatrix} 1 \\ \frac{-1-\sqrt{5}}{2} \\ \end{bmatrix}，\lambda_2=\frac{1+\sqrt{5}}{2} 对应 B_2=\begin{bmatrix} 1 \\ \frac{-1+\sqrt{5}}{2} \\ \end{bmatrix} 。那么根据定义，设 C=[B_1,B_2]，构造对角阵 P=\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \\ \end{bmatrix} ，就可以得到 A\times C=C\times P，所以有 A=C\times P \times C^{-1}。故 A^m=C\times P^m \times C^{-1}。

通过手算可以求得 C^{-1}= \begin{bmatrix} 1-\frac{1+\sqrt{5}}{2}\times \frac{\sqrt{5}}{5} & -\frac{\sqrt{5}}{5} \\ \frac{1+\sqrt{5}}{2}\times \frac{\sqrt{5}}{5} & \frac{\sqrt{5}}{5} \\ \end{bmatrix}。

对于斐波那契数列第 m 项，即是 A^m 的 a_{1,1}。显然 P^m=\begin{bmatrix} (\frac{1-\sqrt{5}}{2})^m & 0 \\ 0 & (\frac{1+\sqrt{5}}{2})^m \ \end{bmatrix} ，那么 C\times P^m \times C^{-1} 的 a_{1,1} 为 (\frac{1-\sqrt{5}}{2})^m\times(1-\frac{1+\sqrt{5}}{2}\times \frac{\sqrt{5}}{5})+(\frac{1+\sqrt{5}}{2})^m\times \frac{1+\sqrt{5}}{2}\times \frac{\sqrt{5}}{5}，整理之后就为 \frac{1}{\sqrt{5}}\left((\frac{1+\sqrt{5}}{2})^{m+1}-(\frac{1-\sqrt{5}}{2})^{m+1}\right) ，即为斐波那契数列的通项公式。

因为出现了根号，所以整个过程看似比较麻烦，不过在递推矩阵满足一些性质的时候有较快的方法求特征值和特征向量（比如在整个矩阵是一个上/下三角矩阵时，特征值就是主对角线上的值，自由元就是项数变少了某行，设它为 1，就可以在矩阵大小的复杂度内求得一个特征向量）。

整理自【线性代数】。

整理自【李超+线性代数在OI中的应用与题目讲解】。