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

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

The n-th term of Fibonacci Numbers：

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

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

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

他始于数列的前两个初始值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

（配图与本文无关）