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

2024-04-26 23:09
32人已看

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

The n-th term of Fibonacci Numbers：

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

By the way：似乎这也是唯一一个不需要迭代就能直接找到斐波那契数列的第n项的算法！

（意大利数学家斐波那契，图片来源于参考文献【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}$ 表示数列中的第i项和第i+1项。

$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 Au_{0}=A^{n}u_{0}$

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

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

这样一来，计算斐波那契数列第n项的问题就变成了求解矩阵A的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项的计算公式 (式1) 就能简化为：

$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)

其中，矩阵X是特征向量矩阵，他的每一列都是一个特征向量。矩阵 $\Lambda$ 是特征值矩阵，这是一个对角阵，只有主对角线上才有元素，且主对角线上的元素都是特征值 $\lambda$ 。

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

第一步，计算矩阵A的特征向量与特征值，根据特征值的个数判断矩阵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}$

通过计算得出了两个不同的特征值，等于矩阵的维度2，说明矩阵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)

把 (式3) 代入斐波那契的第n项的计算公式(式1)就变成了：

$u_{n}=A^{n}u_{0}=A^{n}Xc=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$

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

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

第三步，求和。把第二步的结果加在一起，求出最终的结果。在本例中，因为 $c_{1}\lambda _{1}^{100}x_{1}$ 中 $\lambda _{1}$ 的100次方是一个很小的数，因此这一项的计算结果可以忽略不计，所以他们二者的和几乎约等于 $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$

因此，斐波那契数列的第n项最简化的数值计算公式为：

$u_{n}=c_{2}\lambda _{2}^{n}x_{2}=-0.85065081*1.618033988749895^{n}*\begin{bmatrix} -0.52573111\\ -0.85065081\end{bmatrix}=\begin{bmatrix} F_{n}\\ F_{n+1} \end{bmatrix}$