SkyWT

10/2/2018

矩阵乘法在动态规划中的应用

This blog post is only available in Simplified Chinse.
[x11x12x13x21x22x23x31x32x33]\begin{bmatrix} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{bmatrix}

每增加一个维度，世界便会增加无限的美感。

从 Fibonacci 数列开始

矩阵乘法应用的入门题。

如果把 FnF_n 放在矩阵里，构造出一个 1×2 的矩阵 [Fn−1Fn−2]\begin{bmatrix} F_{n-1} & F_{n-2}\end{bmatrix}，就可以构造一个转移矩阵：

[Fn−1Fn−2]∗[1110]=[Fn−1+Fn−2Fn−1]=[FnFn−1]\begin{bmatrix} F_{n-1} & F_{n-2} \end{bmatrix} \ast \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} F_{n-1}+F_{n-2} & F_{n-1} \end{bmatrix} = \begin{bmatrix} F_{n} & F_{n-1} \end{bmatrix}

这里 [1110]\begin{bmatrix} 1 & 1 \\ 1 & 0\end{bmatrix} 就是转移矩阵。也就是说只需要将构造的 [Fn−1Fn−2]\begin{bmatrix} F_{n-1} & F_{n-2}\end{bmatrix} 乘以转移矩阵就可以得到 [FnFn−1]\begin{bmatrix} F_{n} & F_{n-1}\end{bmatrix}。显然可以用快速幂优化了。

其实也可以写成 4×4 的矩阵进行推导：

[Fn+1FnFnFn−1]=[1110]n=[1110][1110]⋯[1110]⏟n times\displaystyle \begin{bmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}^n = \underbrace{ \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\cdots \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix} }_ {\text{n times}}

拓展一下？

由上面的例子，我们可以扩展一下：对于这样的递推式：

Fn=Fn−1+Fn−2+⋯+Fn−kF_{n}=F_{n-1}+F_{n-2}+\dots+F_{n-k}

我们其实都可以通过构造矩阵的方法来优化，快速计算 FnF_n。 （不过很多简单的递推式也可以求通项公式……）

POJ 3233 Matrix Power Series

Description

Link： POJ 3233 Matrix Power Series

Given a n×nn × n matrix AA and a positive integer kk, find the sum S=A+A2+A3+⋯+AkS = A + A^2 + A^3 + \dots + A^k.

n⩽30n \leqslant 30, k⩽109k \leqslant 10^9 and m<104m < 10^4.

Hint

一道比较典型的矩阵乘法优化递推的题目。有很多方法。

Analysis #1

先把矩阵 AA 看成一个数字。由题可知：

F(n)=A+A2+A3+⋯+AnF(n)=A + A^2 + A^3 + \dots + A^n

那么

F(n)=A∗F(n−1)+AF(n)=A \ast F(n-1) + A

显然我们可以构造一个转移矩阵：

[An−1Fn−101]∗[AA01]=[AnFn01]\begin{bmatrix} A^{n-1} & F_{n-1} \\ 0 & 1 \end{bmatrix} \ast \begin{bmatrix} A & A \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} A^{n} & F_{n} \\ 0 & 1 \end{bmatrix}

也就是说我们发现：

[AA01]∗[AA01]=[A2A2+A01]\begin{bmatrix} A & A \\ 0 & 1 \end{bmatrix} \ast \begin{bmatrix} A & A \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} A^2 & A^2+A \\ 0 & 1 \end{bmatrix}

则：

S=[AA01]kS= \begin{bmatrix} A & A \\ 0 & 1 \end{bmatrix} ^k

写一个一百多行的矩阵套矩阵快速幂就可以了。

Analysis #2

这是最简单的方法。

F(n)=A+A2+A3+⋯+AnF(n)=A + A^2 + A^3 + \dots + A^n

  1. nn 为偶数：

F(n)=A+A2+⋯+An/2+An/2∗F(n/2)F(n)=A+A^2+\dots+A^{n/2} +A^{n/2}\ast F(n/2)

  1. nn 为奇数：

F(n)=A∗F(n−1)+AF(n)=A\ast F(n-1)+A

这样其实就是分治，可以 log⁡k\log k 求解。

Analysis #3

看看 F(n)=A∗F(n−1)+AF(n)=A\ast F(n-1)+A 这个递推式，考虑用它推通项公式。

F(n)+λ=A∗F(n−1)+A+λF(n)+\lambda = A\ast F(n-1)+A +\lambda

（我们用 II 表示单位矩阵）

F(n)+λ=A∗(F(n−1)+I+λA)F(n)+\lambda = A \ast (F(n-1)+I+\frac {\lambda} {A})

得 λ=I+λA\displaystyle \lambda = I+\frac {\lambda} A ，解得 λ=A∗(A−I)−1\displaystyle \lambda =A\ast (A-I)^{-1}

(A−I)−1(A-I)^{-1} 就是 A−IA-I 的逆矩阵。因为这个给你的矩阵是 n∗nn\ast n 的方阵，所以其实可以直接用高斯消元求其逆矩阵。 接下来令 G(n)=F(n)+A∗(A−I)−1G(n)=F(n)+A\ast (A-I)^{-1}，那么

G(n)=A∗G(n−1)G(n)=A\ast G(n-1)

即

G(n)=An−1∗G(1)G(n)=A^{n-1}\ast G(1)

矩阵快速幂即可求解。

Code

我只写了第一种思路的代码……

这题代码有个坑点，如果全都开 long long 会过不了，必须开 int 让其自然溢出……可能数据有点问题[1]。

（不过谁能想到有人居然所有代码前面都加 #define int long long 的……）

/*
 * Vjudge CONTEST244508 矩阵乘法专题训练
 * POJ 3233
 * C - Matrix Power Series
 * 180929 By SkyWT
 */

#include<cstdio>
#include<iostream>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#include<vector>
#include<queue>
#include<stack>
#include<map>

using namespace std;

#define memset_clear(_) memset(_,0,sizeof(_))
#define memset_clear_tre(_) memset(_,1,sizeof(_))
#define memset_clear_reg(_) memset(_,-1,sizeof(_))
#define memset_clear_max(_) memset(_,0x3f,sizeof(_))
#define memset_clear_min(_) memset(_,0x80,sizeof(_))
#define sqr(_) ((_)*(_))

// #define int long long

const int maxn=35;
int n,k,tt;

inline int read(){
	int ret=0,f=1;char ch=getchar();
	while (ch<'0'||ch>'9'){if (ch=='-') f=-1;ch=getchar();}
	while (ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar();
	return ret*f;
}

struct matrix{
	int a[maxn][maxn];
	void init(){ // 构造一个单位矩阵
		memset(a,0,sizeof(a));
		for (int i=0;i<n;i++) a[i][i]=1;
	}
	void clear(){
		memset(a,0,sizeof(a));
	}
	matrix operator +(matrix b){
		matrix c;c.clear();
		for (int i=0;i<n;i++)
		for (int j=0;j<n;j++)
			c.a[i][j]=a[i][j]+b.a[i][j];
		return c;
	}
	matrix operator *(matrix b){
		matrix c; c.clear();
		for (int i=0;i<n;i++)
		for (int j=0;j<n;j++)
		for (int k=0;k<n;k++)
			c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%tt)%tt;
		return c;
	}
	matrix operator ^(int b){
		matrix ret; ret.init();
		matrix w;   w.clear();
		for (int i=0;i<2;i++)
		for (int j=0;j<2;j++)
			w.a[i][j]=a[i][j];
		while (b){
			if (b&1) ret=ret*w;
			w=w*w;b>>=1;
		}
		return ret;
	}
}fst;

struct matrix_matrix{
	matrix a[3][3];
	void init(matrix number,bool flag){
		if (flag){
			a[0][0].init();a[1][1].init();
			a[0][1].clear();a[1][0].clear();
		} else {
			a[0][0]=a[0][1]=number;
			a[1][0].clear();a[1][1].init();
		}
	}
	void clear(){
		a[0][0].clear();a[0][1].clear();
		a[1][0].clear();a[1][1].clear();
	}
	matrix_matrix operator *(matrix_matrix b){
		matrix_matrix c; c.clear();
		for (int i=0;i<2;i++)
		for (int j=0;j<2;j++)
		for (int k=0;k<2;k++)
			c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]);
		return c;
	}
	matrix_matrix operator ^(int b){
		matrix tmp;
		matrix_matrix ret;ret.init(tmp,true);
		matrix_matrix w;w.init(fst,false);
		while (b){
			if (b&1) ret=ret*w;
			w=w*w;b>>=1;
		}
		return ret;
	}
	void write(){
		for (int i=0;i<n;i++){
			for (int j=0;j<n;j++) printf("%d ",a[0][1].a[i][j]%tt);
			printf("\n");
		}
	}juzhegn
};

signed main(){
	n=read();k=read();tt=read();
	for (int i=0;i<n;i++)
	for (int j=0;j<n;j++)
		fst.a[i][j]=read()%tt;
	matrix_matrix now;
	now.init(fst,false);
	now=now^k;
	now.write();
	return 0;
}

HDU 2065 "红色病毒"问题

Description

Link：HDU 2065 "红色病毒"问题

医学界发现的新病毒因其蔓延速度和 Internet 上传播的“红色病毒”不相上下，被称为“红色病毒”。经研究发现，该病毒及其变种的 DNA 的一条单链中，胞嘧啶、腺嘧啶均是成对出现的。
现在有一长度为N的字符串，满足以下条件：

  1. 字符串仅由 A,B,C,D 四个字母组成；
  2. A 出现偶数次（也可以不出现）；
  3. C 出现偶数次（也可以不出现）。

计算满足条件的字符串个数。 当 N=2 时，所有满足条件的字符串有如下 6 个：BB,BD,DB,DD,AA,CC。 由于这个数据肯能非常庞大，你只要给出最后两位数字即可。

Hint

矩阵优化 DP 的典型。

对于前面一题递推的题目，我们尚可直接用求通项的数学方法，但是这题似乎不行了……

Analysis #1

首先可以想到 F[i] 表示长度为 i 的字符串，偶数个 A 并且偶数个 C 的数量。但是直接这样定义显然是没法直接状态转移的…… 我们自然想到分别在 DP 状态里记录下 A 和 C 了。 为了方便放进矩阵里，我们尽量定义成一个二维的：F[i][0/1/2/3]，i 表示字符串长度为 i，第二维讨论一下：

  1. A 和 C 均出现了偶数次；
  2. A 出现了偶数次，C 出现了奇数次；
  3. A 出现了奇数次，C 出现了偶数次；
  4. A 和 C 均出现奇数次。

那么就有如下转移方程：

F(i,0)=2∗F(i−1,0)+F(i−1,1)+F(i−1,2)F(i,1)=2∗F(i−1,1)+F(i−1,0)+F(i−1,3)F(i,2)=2∗F(i−1,2)+F(i−1,0)+F(i−1,3)F(i,3)=2∗F[i−1,3]+F(i−1,1)+F(i−1,2)F(i,0)=2\ast F(i-1,0) + F(i-1,1) + F(i-1,2) \\ F(i,1)=2\ast F(i-1,1) + F(i-1,0) + F(i-1,3) \\ F(i,2)=2\ast F(i-1,2) + F(i-1,0) + F(i-1,3) \\ F(i,3)=2\ast F[i-1,3] + F(i-1,1) + F(i-1,2)

让我们对齐一下，以便于后续的推导：

F(i,0)=2∗F(i−1,0)+1∗F(i−1,1)+1∗F(i−1,2)+0∗F(i−1,3)F(i,1)=1∗F(i−1,0)+2∗F(i−1,1)+0∗F(i−1,2)+1∗F(i−1,3)F(i,2)=1∗F(i−1,0)+0∗F(i−1,0)+2∗F(i−1,2)+1∗F(i−1,3)F(i,3)=0∗F(i−1,0)+1∗F(i−1,1)+1∗F(i−1,2)+2∗F(i−1,3)F(i,0)=2\ast F(i-1,0) + 1\ast F(i-1,1) + 1\ast F(i-1,2) + 0\ast F(i-1,3) \\ F(i,1)=1\ast F(i-1,0) + 2\ast F(i-1,1) + 0\ast F(i-1,2) + 1\ast F(i-1,3) \\ F(i,2)=1\ast F(i-1,0) + 0\ast F(i-1,0) + 2\ast F(i-1,2) + 1\ast F(i-1,3) \\ F(i,3)=0\ast F(i-1,0) + 1\ast F(i-1,1) + 1\ast F(i-1,2) + 2\ast F(i-1,3)

现在考虑通过矩阵转移。假设我们把 F(i,j)F(i,j) 看成如下矩阵：

[F(i,0)F(i,1)F(i,2)F(i,3)]\begin{bmatrix} F(i,0) \\ F(i,1) \\ F(i,2) \\ F(i,3) \end{bmatrix}

现在我们需要转移成这样子：

[2∗F(i−1,0)+1∗F(i−1,1)+1∗F(i−1,2)+0∗F(i−1,3)1∗F(i−1,0)+2∗F(i−1,1)+0∗F(i−1,2)+1∗F(i−1,3)1∗F(i−1,0)+0∗F(i−1,0)+2∗F(i−1,2)+1∗F(i−1,3)0∗F(i−1,0)+1∗F(i−1,1)+1∗F(i−1,2)+2∗F(i−1,3)]\begin{bmatrix} 2\ast F(i-1,0) + 1\ast F(i-1,1) + 1\ast F(i-1,2) + 0\ast F(i-1,3) \\ 1\ast F(i-1,0) + 2\ast F(i-1,1) + 0\ast F(i-1,2) + 1\ast F(i-1,3) \\ 1\ast F(i-1,0) + 0\ast F(i-1,0) + 2\ast F(i-1,2) + 1\ast F(i-1,3) \\ 0\ast F(i-1,0) + 1\ast F(i-1,1) + 1\ast F(i-1,2) + 2\ast F(i-1,3) \end{bmatrix}

其实转移矩阵已经很显然了，就是那几个系数构成的矩阵：

[2110120110210112]\begin{bmatrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{bmatrix}

现在还有个微小的问题，这个矩阵是 4×4 的，但是原来是 4×1 的……4 行 1 列的矩阵可不能与 4 行 4 列的矩阵相乘。我们只能把 4×1 的矩阵写成 1×4 的：

[F(i,0)F(i,1)F(i,2)F(i,3)]\begin{bmatrix} F(i,0) & F(i,1) & F(i,2) & F(i,3) \end{bmatrix}

（可以发现这是一个沿”主对角线“对称的矩阵） 接下来考虑下初始矩阵，不难得出答案：

[2110]∗[2110120110210112]n−1=[F(n,0)F(n,1)F(n,2)F(n,3)]\begin{bmatrix} 2 & 1 & 1 & 0 \end{bmatrix} \ast \begin{bmatrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{bmatrix} ^{n-1} = \begin{bmatrix} F(n,0) & F(n,1) & F(n,2) & F(n,3) \end{bmatrix}

矩阵快速幂优化即可。

Analysis #2

网上题解搜得：泰勒级数？不会。

Code

需要注意题目的输入格式：

每组输入的第一行是一个整数 T，表示测试实例的个数。下面是 T 行数据，每行一个整数 N(1⩽N<264)N (1\leqslant N<2^{64})，当 T=0 时结束。

N 达到了 2642^{64}，所以我们需要开 unsigned long long。 （正常的操作当然是#define int unsigned long long 啦～）

解法 #1 的代码：

/*
 * Vjudge CONTEST 244508
 * HDU 2065
 * B - "红色病毒"问题
 * 181002 By SkyWT
 */

#include<cstdio>
#include<iostream>
#include<cstring>
#include<string>
#include<cmath>
#include<algorithm>
#include<vector>
#include<queue>
#include<stack>
#include<map>

using namespace std;

#define memset_clear(_) memset(_,0,sizeof(_))
#define memset_clear_tre(_) memset(_,1,sizeof(_))
#define memset_clear_reg(_) memset(_,-1,sizeof(_))
#define memset_clear_max(_) memset(_,0x3f,sizeof(_))
#define memset_clear_min(_) memset(_,0x80,sizeof(_))
#define sqr(_) ((_)*(_))

#define int unsigned long long

inline int read(){
	int ret=0,f=1;char ch=getchar();
	while (ch<'0'||ch>'9') {if (ch=='-') f=-1;ch=getchar();}
	while (ch>='0'&&ch<='9') ret=ret*10+ch-'0',ch=getchar();
	return ret*f;
}

const int maxn=5,tt=100;
const int transfer[4][4]={{2,1,1,0},
						  {1,2,0,1},
						  {1,0,2,1},
						  {0,1,1,2}};

struct matrix{
	int n,m,a[maxn][maxn];
	void set_identity(){
		memset_clear(a);
		for (int i=0;i<min(n,m);i++) a[i][i]=1;
	}
	void init(){
		memset_clear(a);
	}
	void give(matrix &b){
		b.n=n;b.m=m;
		for (int i=0;i<n;i++)
		for (int j=0;j<m;j++)
			b.a[i][j]=a[i][j];
	}
	matrix operator +(matrix b){
		if (b.n!=n||b.m!=m) printf("ERROR: TWO DEFFERENT MATRIX MAKE +\n");
		matrix c; c.init();
		c.n=n;c.m=m;
		for (int i=0;i<n;i++)
		for (int j=0;j<m;j++)
			c.a[i][j]=(a[i][j]+b.a[i][j])%tt;
		return c;
	}
	matrix operator -(matrix b){
		if (b.n!=n||b.m!=m) printf("ERROR: TWO DEFFERENT MATRIX MAKE -\n");
		matrix c; c.init();
		c.n=n;c.m=m;
		for (int i=0;i<n;i++)
		for (int j=0;j<m;j++)
			c.a[i][j]=(a[i][j]-b.a[i][j]+tt)%tt;
		return c;
	}
	matrix operator *(matrix b){
		if (m!=b.n) printf("ERROR: TWO DEFFERENT MATRIX MAKE *\n");
		matrix c; c.init();
		c.n=n;c.m=b.m;
		for (int i=0;i<  n;i++)
		for (int j=0;j<  m;j++)
		for (int k=0;k<b.m;k++)
			c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%tt)%tt;
		return c;
	}
	matrix operator ^(int b){
		matrix ret;ret.n=n;ret.m=m;ret.set_identity();
		matrix w;w.n=n;w.m=m;give(w);
		while (b){
			if (b&1) ret=ret*w;
			w=w*w;b>>=1;
		}
		return ret;
	}
};

signed main(){
	int T=0,n=0;
	scanf("%llu",&T);
	while (T){
		for (int k=0;k<T;k++){
			scanf("%llu",&n);

			matrix now;now.init();
			now.n=1;now.m=4;
			now.a[0][0]=2;now.a[0][1]=now.a[0][2]=1;

			matrix tomul;
			tomul.n=4;tomul.m=4;
			for (int i=0;i<4;i++)
			for (int j=0;j<4;j++)
				tomul.a[i][j]=transfer[i][j];

			tomul=tomul^(n-1);
			now=now*tomul;
			printf("Case %llu: %llu\n",k+1,now.a[0][0]);
		}
		printf("\n");
		scanf("%llu",&T);
	}
	return 0;
}
10/2/2018 10:54:00 AM

Post a New Comment

Please login to leave a comment.