矩阵乘法优化线性递推
矩阵乘法是线性代数中一块很重要的内容.矩阵乘法的定义很奇怪[1],但正是这种奇怪的性质,让矩阵乘法成为在除了线性代数和其衍生学科(还有诸如矩阵力学之类)外最广泛使用的关于矩阵变换的应用.(什么?FFT不属于矩阵变换吧...)
注:
[1]: 矩阵乘法有另外的很多定义,如未说明,指的是中间不带符号的矩阵乘法,即一般矩阵乘积.另有 标量乘积(即所有数乘上一个固定的数),阿达马乘积等,没有那么诡异,但是在大多数问题的用途上也不大.
你不会矩阵乘法?没关系,下一篇会写到的.快速幂也是.
矩阵乘法的本质
矩阵乘法的本质是线性组合.
一个线性组合是像这样的:
$c=a_1b_1 + a_2b_2 + ... + a_nb_n$
其中$a_n$这个数列是一个常量,则乘c是数列b的一个线性组合.同理对于a.
So, 如何应用到线性递推中?
让我们来看看最经典的题目,求斐波那契数列.
当你们刚学递归的时候,Fib的题目是这样的:
求$Fib(n)$,其中$Fib(n)=Fib(n-1)+Fib(n-2)$,$Fib(0)=Fib(1)=1$,$n\le 30$
你的程序大约会是这样的:
int fib(int n){
if(n==0||n==1) return 1;
return fib(n-1)+fib(n-2);
}//你大概注意到了,这样写的效率很低,在n不大时就崩溃了.(n=30还是可以承受的.)事实上,这个时间是以指数级别增长的.
然后你们学习了记忆化搜索(其实这是一个非常好的思想,几乎是递推和动规生命的早期).这时你们的题目会是这样的:
求$Fib(n) \mod 10^9+7$,其中$Fib(n)=Fib(n-1)+Fib(n-2)$,$Fib(0)=Fib(1)=1$,$n\le 100000$
你的程序会是这样的:
int fibb[100010];
int fib(int n){//不严密.可能(几乎一定)会爆栈
if(fibb[n]!=0){
return fibb[n];
}
if(n==0||n==1) return (fibb[n]=1);
return (fibb[n]=(fib(n-1)+fib(n-2))%1000000007);
}
相比起纯粹递归来说,这是非常好的算法.
当你们学习了递推后,程序变得快而简洁:
#include <cstdio>
int n,i,f[10000000];//大概可以处理一千万的数据.
int main(){//本文给出的第一个完整程序.
scanf("%d",&n);
f[0]=f[1]=1;
for(i=2;i<=n;++i){
f[i]=(f[i-1]+f[i-2])%1000000007;
}
printf("%d",f[n]);
return 0;
}
你又知道了滚动数组.空间已经不是问题.
#include <cstdio>
int n,i,a,b,c;
int main(){
scanf("%d",&n);
a=b=1;
for(i=2;i<=n;++i){
c=(a+b)%1000000007;
a=b;
b=c;
}
printf("%d",b);
return 0;
}
这是最优算法了么?
看起来是的,但其实这只是线性递推的初步.
线性递推是什么呢?
线性递推是一类线性逐步推进的问题.
何为线性?
注意到上面提到的线性组合.线性递推指用线性方程递推的问题.
很容易写出Fibonacci数列的线性递推方程:
$a_1=0a_1+a_2\\
a_2=a_1+a_2$ //注意这是一个阶段,当$a_2$求解时$a_1$并未更改.
而矩阵乘法的本质是线性组合.于是,一个简单的矩阵可以代替这个公式:
$T = tr = \left[ \begin{array} {lcr}
0 & 1 \\
1 & 1
\end{array} \right]$
初始值为:
$A = initial = \left[ \begin{array} {lcr}
1 \\
1
\end{array} \right]$ //非习惯性将向量写作竖格.因为这是一个$2 \times 1$矩阵.向量就是竖的.虽然可能不符合直觉,但是矩阵乘法也不符合直觉.
那么一次递推即是:
$TA$ //注意不可调换位置.切记.
两次递推就是:
$T(TA)$
三次则是:
$T(T(TA))$
注意矩阵乘法满足结合率,那么:
$T(T(TA))=(T^3)A$
因为结合率,于是我们可以使用二分快速幂,以十次递归为例:
$T(T(T(T(T(T(T(T(T(TA)))))))))=(T^{10})A=(T^2T^8)A$
这里$T^2$和$T^8$可以翻倍生成,即$(T)^2=T^2,(T^2)^2=T^4,(T^4)^2=T^8...$
矩阵运算时间复杂度为常数,因为固定$2\times 2$矩阵,可以认为时间为$\text{Time}(TMultiply\times 8+TModulo\times 4)$(求余是必需的)
而二分快速幂复杂度为$\text{O}\left( \log_2{n}\right)$
让我们干一件无聊的事: 估算可以达到的上界.
首先,假定CPU每秒运算20亿次$TGeneral$,$TMultiply=TGeneral \times 8, TModulo=TGeneral \times 20$
而二分的位移和位与时间均为$TGeneral$,四次赋值和判断是否应该退出为$TGeneral \times 5$,则一次迭代时间大约$TGeneral\times 151$
算上各种损耗和一些指令集的并发优化, 记一次迭代为$TFibItr$大约$TGeneral \times 200$
可执行迭代次数为$\frac{2 \times 10^9}{200}=10^7$,即最大可计算到$2^{10^7}$,这是一个大约$3010300$位的大数!但是输入都达不到这么大.
这就是快速幂的威力.
