矩阵乘法优化线性递推

矩阵乘法是线性代数中一块很重要的内容.矩阵乘法的定义很奇怪^[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$,$nle 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$,$nle 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[ egin{array} {lcr}
0 & 1 \
1 & 1
end{array} ight]$

初始值为:

$A = initial = left[ egin{array} {lcr}
1 \
1
end{array} ight]$ //非习惯性将向量写作竖格.因为这是一个$2 imes 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 imes 2$矩阵,可以认为时间为$ ext{Time}(TMultiply imes 8+TModulo imes 4)$(求余是必需的)

而二分快速幂复杂度为$ ext{O}left( log_2{n} ight)$

让我们干一件无聊的事: 估算可以达到的上界.

首先,假定CPU每秒运算20亿次$TGeneral$,$TMultiply=TGeneral imes 8, TModulo=TGeneral imes 20$

而二分的位移和位与时间均为$TGeneral$,四次赋值和判断是否应该退出为$TGeneral imes 5$,则一次迭代时间大约$TGeneral imes 151$

算上各种损耗和一些指令集的并发优化, 记一次迭代为$TFibItr$大约$TGeneral imes 200$

可执行迭代次数为$frac{2 imes 10^9}{200}=10^7$,即最大可计算到$2^{10^7}$,这是一个大约$3010300$位的大数!但是输入都达不到这么大.

这就是快速幂的威力.