线性代数 --- 计算斐波那契数列第n项的快速算法(矩阵的n次幂)-编程知识

计算斐波那契数列第n项的快速算法(矩阵的n次幂)

The n-th term of Fibonacci Numbers：

斐波那契数列的是一个古老而又经典的数学数列，距今已经有800多年了。关于斐波那契数列的计算方法不难，只是当我们希望快速求出其数列中的第100，乃至第1000项时，有没有又准又快的方法，一直是一个值得探讨和研究的问题。笔者(松下J27)在这篇文章中，就介绍了一种基于线性代数的快速算法，即，基于矩阵的n次幂找到斐波那契数列的第n项。这是我在MIT线性代数的公开课中看到的，并以此文记录下来。

（意大利数学家斐波那契，图片来源于参考文献【1】）

已知斐波那契数列的数学表示方式如下：

$Fibonacci:0,1,1,2,3,5,8,13,21,34......$

他始于数列的前两个初始值0，1，后面所有的项都是前两项的和，依此类推得到：

$0+1=1$ ， $1+1=2$ ， $1+2=3$ ， $2+3=5$ ， $3+5=8$ 。。。

基于上面的计算规律，假设我用 $F_{i}$ 来表示该数列的第i个数的话，那斐波那契数列用代数的方式可写成两个初值和一个递归公式的组合，即：

$F_{0}=0,F_{1}=1$

$F_{i+2}=F_{i}+F_{i+1},\: i=0,1,2...$

现在把上面的公式用矩阵和向量的方式来表示，以便于用线性代数来分析：

1，先把初值用向量u0来表示

$u_{0}=\begin{bmatrix} F_{0}\\ F_{1} \end{bmatrix}=\begin{bmatrix} 0\\ 1\end{bmatrix}$

2，用向量 $u_{i}$ 表示数列中的第n项。

$u_{i}=\begin{bmatrix} F_{i}\\ F_{i+1} \end{bmatrix}$

3，如此一来根据递归公式就能写出第i+1项 $u_{i+1}$

$u_{i+1}=\begin{bmatrix} F_{i+1}\\ F_{i+2} \end{bmatrix}=\begin{bmatrix} F_{i+1}\\ F_{i}+F_{i+1} \end{bmatrix}$

4，用矩阵来表示从第i项 $u_{i}$ 到第i+1项 $u_{i+1}$ 的过程

$u_{i+1}=\begin{bmatrix} F_{i+1}\\ F_{i+2} \end{bmatrix}=\begin{bmatrix} F_{i+1}\\ F_{i}+F_{i+1} \end{bmatrix}=\begin{bmatrix} 0 &1 \\ 1& 1\end{bmatrix}\begin{bmatrix} F_{i}\\ F_{i+1} \end{bmatrix}=\begin{bmatrix} 0 &1 \\ 1& 1\end{bmatrix}u_{i}$

其中，参与计算的矩阵A为：

$A=\begin{bmatrix} 0& 1\\ 1&1 \end{bmatrix}$

现在，我们知道了初始向量 $u_{0}$ ，也知道如果计算 $u_{i+1}$ ，我们就能通过线性代数中矩阵与向量的计算来计算斐波那契数列中的任意一项：

可以看到，如果你要计算第100项，不出意外，只要把上述步骤继续下去一定能够找到。另外，如果把上述计算中的三次计算合成一部，就是三个矩阵A连续相乘后，再乘以u0。

这就是说，斐波那契数列中第n项的计算公式应该是：

$u_{n}=A\cdot A\cdot A\cdot A...\cdot Au0=A^{n}u0$

$u_{n}=A^{n}u0$ ， (式1)

这里要注意的是，根据这种方法算出来的是一个向量，而我们需要的结果，即，第n项是该向量中的第一个元素。

引入矩阵的对角化：

熟悉线性代数的朋友看到A的n次幂时，马上就会想到是否可以通过特征向量矩阵X和特征值矩阵 $\Lambda$ 对A进行对角化。

线性代数 --- 矩阵的对角化以及矩阵的n次幂-CSDN博客文章浏览阅读1k次，点赞15次，收藏9次。本文从矩阵A的对角化开始，一直聊到了对角化的应用，并以一个A的n次幂为例子结尾。https://blog.csdn.net/daduzimama/article/details/138088128

因为，如果方阵A可以被对角化为：

$A=X\Lambda X^{-1}$

那么上面推导出来的斐波那契数列的第n项的计算公式就能简化为：

$u_{n}=A\cdot A\cdot A\cdot A...\cdot Au0=A^{n}u0\Rightarrow$

$u_{n}=(X\Lambda X^{-1})\cdot (X\Lambda X^{-1})\cdot (X\Lambda X^{-1})...\cdot (X\Lambda X^{-1})u0=X\Lambda ^{n}X^{-1}u0$

$u_{n}=X\Lambda ^{n}X^{-1}u0$ ， (式2)

现在我们试着用(式2)去计算斐波那契数列的第100项：

第一步，计算矩阵A的特征向量与特征值：

$\lambda _{1}=-0.61803399,x_{1}=\begin{bmatrix}-0.85065081 \\ 0.52573111\end{bmatrix}$

$\lambda _{2}=1.61803399,x_{2}=\begin{bmatrix}-0.52573111\\ -0.85065081\end{bmatrix}$

通过计算得出了两个不同的特征值，说明矩阵A可以对角化。

第二步，构建特征向量矩阵X，特征值矩阵 $\Lambda$ 和特征向量矩阵的逆：

$X=\begin{bmatrix} -0.85065081&-0.52573111 \\ 0.52573111&-0.85065081 \end{bmatrix}$

$\Lambda =\begin{bmatrix} -0.61803399 &0 \\ 0 & 1.61803399 \end{bmatrix}$

$X^{-1}=\begin{bmatrix} -0.85065081&0.52573111 \\ -0.52573111&-0.85065081 \end{bmatrix}$

第三步，计算特征值矩阵 $\Lambda$ 的100次幂：

$\Lambda ^{100}=\begin{bmatrix} -0.61803399^{100} &0 \\ 0 & 1.61803399^{100} \end{bmatrix}$

$=\begin{bmatrix}1.26251334E-21 &0 \\ 0 & 7.92070840E+20 \end{bmatrix}$

第四步，基于(式2)去计算第100项的值：

$u_{100}=\begin{bmatrix} 3.54224848E+20\\ 5.73147844E+20 \end{bmatrix}$

其中，第100项的值为：

$F_{100}=3.54224848E+20$

进一步简化：

如果我能把u0表示成所有特征向量的线性组合的话，就能更进一步简化计算：

$u_{0}=c_{1}x_{1}+c_{2}x_{2}...+c_{n}x_{n}=Xc$ ， (式3)

其中，权重系数向量c等于：

$c=X^{-1}u_{0}$ ， (式4)

如此一来，斐波那契的第n项的计算公式(式1)就变成了：

$u_{n}=A^{n}(c_{1}x_{1}+c_{2}x_{2}...+c_{n}x_{n})=c_{1}A^{n}x_{1}+c_{2}A^{n}x_{2}...+c_{n}A^{n}x_{n}$ ， (式5)

又因为这里的x都是特征向量，根据特征向量的性质，我们有Key Equation:

$Ax=\lambda x$

对于上式中的第一项 $c_{1}A^{n}x_{1}$ 有：

$c_{1}A^{n}x_{1}=c_{1}A\cdot A...\cdot A\cdot Ax_{1}=c_{1}A\cdot A...\cdot A\cdot \lambda _{1}x_{1}$

因为 $\lambda _{1}$ 是一个系数，所以我们把他提到前面去：

$c_{1}A^{n}x_{1}=c_{1}\lambda _{1}A\cdot A...\cdot Ax_{1}$

这样一来又出现了一个 $\lambda _{1}$ ，继续：

$c_{1}A^{n}x_{1}=c_{1}\lambda _{1}A\cdot A...\cdot \lambda _{1}x_{1}$

如此反复，一直到第n个乘法：

$c_{1}A^{n}x_{1}=c_{1}\lambda _{1}^{n}A$

依此类推，(式5)可改写成：

$u_{n}=c_{1}\lambda _{1}^{n}x_{1}+c_{2}\lambda _{2}^{n}x_{2}...+c_{n}\lambda _{n}^{n}x_{n}$ ， (式6)

(式6) 也可以写成：

$u_{n}=c_{1}\lambda ^{n}x_{1}+c_{2}\lambda ^{n}x_{2}...+c_{n}\lambda ^{n}x_{n}=X\Lambda ^{n}c$ ， (式7)

这一点也可以通过把 (式3)代入(式2)得到(式7)：

$u_{n}=X\Lambda ^{n}X^{-1}u0=X\Lambda ^{n}X^{-1}(Xc)=X\Lambda ^{n}c$

这里面最重要的是(式6)，因为他把计算分离开了，分成了一个个常数与向量的乘积之和。

现在针对这一简化算法做一个小结，并以求斐波那契的第100项为例：

第一步，优先使用(式4)算出权重向量c

$c=X^{-1}u_{0}=\begin{bmatrix} 0.52573111\\ -0.85065081 \end{bmatrix}$

$\Rightarrow c_{1}=0.52573111, c_{2}=-0.85065081$

第二步，分别计算 $c_{1}\lambda _{1}^{100}x_{1}$ 和 $c_{2}\lambda_{2} ^{100}x_{2}$ ，其中c和 $\lambda$ 都是常数。

第三步，求和。把第二步的结果加在一起，求出最终的结果。在本例中，因为 $c_{1}\lambda _{1}^{100}x_{1}$ 中的元素几乎为0，所以他们二者的和几乎约等于 $c_{2}\lambda_{2} ^{100}x_{2}$ 。

$u_{100}=c_{1}\lambda _{1}^{100}x_{1}+c_{2}\lambda _{2}^{100}x_{2}=\begin{bmatrix} F_{100}\\ F_{101} \end{bmatrix}$

最终，得到了同样正确的结果：

$F_{100}=3.54224848E+20$

（全文完）

--- 作者，松下J27

参考文献(鸣谢)：

1，https://zh.wikipedia.org/wiki/%E6%96%90%E6%B3%A2%E9%82%A3%E5%A5%91

2， Lec22_对角化和矩阵乘幂_哔哩哔哩_bilibili

（配图与本文无关）