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

This blog post is only available in Simplified Chinse.

\begin{bmatrix} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{bmatrix}

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

从 Fibonacci 数列开始

矩阵乘法应用的入门题。

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

\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}

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

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

\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}}

拓展一下？

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

$F_{n}=F_{n-1}+F_{n-2}+\dots+F_{n-k}$

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

POJ 3233 Matrix Power Series

Description

Link： POJ 3233 Matrix Power Series

Given a $n × n$ matrix $A$ and a positive integer $k$ , find the sum $S = A + A^2 + A^3 + \dots + A^k$ .

$n \leqslant 30$ , $k \leqslant 10^9$ and $m < 10^4$ .

Hint

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

Analysis #1

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

$F(n)=A + A^2 + A^3 + \dots + A^n$

那么

$F(n)=A \ast F(n-1) + A$

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

\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}

也就是说我们发现：

\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= \begin{bmatrix} A & A \\ 0 & 1 \end{bmatrix} ^k

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

Analysis #2

这是最简单的方法。

$F(n)=A + A^2 + A^3 + \dots + A^n$

$n$ 为偶数：

$F(n)=A+A^2+\dots+A^{n/2} +A^{n/2}\ast F(n/2)$

$n$ 为奇数：

$F(n)=A\ast F(n-1)+A$

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

Analysis #3

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

$F(n)+\lambda = A\ast F(n-1)+A +\lambda$

（我们用 $I$ 表示单位矩阵）

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

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

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

$G(n)=A\ast G(n-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的字符串，满足以下条件：

字符串仅由 A,B,C,D 四个字母组成；

A 出现偶数次（也可以不出现）；

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，第二维讨论一下：

A 和 C 均出现了偶数次；
A 出现了偶数次，C 出现了奇数次；
A 出现了奇数次，C 出现了偶数次；
A 和 C 均出现奇数次。

那么就有如下转移方程：

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\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)$ 看成如下矩阵：

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

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

\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}

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

\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 的：

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

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

\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\leqslant N<2^{64})$ ，当 T=0 时结束。

N 达到了 $2^{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

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

从 Fibonacci 数列开始

拓展一下？

POJ 3233 Matrix Power Series

Description

Hint

Analysis #1

Analysis #2

Analysis #3

Code

HDU 2065 "红色病毒"问题

Description

Hint

Analysis #1

Analysis #2

Code

Post a New Comment