0 向量与矩阵基础
向量是一个 \(n\) 维有方向的量记为\(\alpha=(a_1,a_2,\dots,a_n)\),\(a_i\) 是其在第 \(i\) 维上的分量。
向量可以定义加法(两个 \(n\) 维向量 \(\alpha,\beta\)):\(\alpha+\beta=(a_1+b_1,a_2+b_2,...a_n+b_n)\)。
减法即是加法的逆运算,也可以定义。
向量的数乘:\(c \times \alpha=(ca_1,ca_2,...,ca_n)\)。
向量也分行向量与列向量,\(\alpha^T\) 记为 \(\alpha\) 的转置。
向量点积:\(\alpha \cdot \beta=a_1b_1+a_2b_2+...+a_nb_n\)。
向量的模:\(|\alpha|=\sqrt{a1^2+a2^2+...+an^2}\)
在二维平面若向量夹角为 \(\theta\) 。则 \(\alpha \times \beta=\mid\alpha\mid\mid\beta\mid\cos{\theta}\)
二维平面内向量叉积:\(\alpha \times \beta=a_1b_2-a_2b_1\)。方向是垂直于平面的 \(n\) 维向量,由右手定则判断。
矩阵是一个 \(n*m\) 的数组(第 \(i\) 行第 \(j\) 列记为 \(a+{i,j}\))。矩阵转置即把 \(j,i\) 互换。
这 \(n\times m\) 个数称为矩阵的元。\(a_{i,j}\)称为 \((i,j)\) 元。以数 \(a_{i,j}\) 为 \((i,j)\) 元的矩阵可记为 \(a(ij)\) 或 \(a(ij)_{n\times m}\)。
矩阵的乘法要求 \(A\) 是 \(n\times m\) 而 \(B\) 是 \(m\times r\) 的矩阵,得到 \(C\) 是 \(n\times r\) 的矩阵:\(C_{i,j}=\sum_{k=1}^{m}A_{i,k}B_{k,j}.\)
1 矩阵快速幂优化数列递推
我们来看最基础的斐波那契数列: \(f_n=f_{n-1}+f_{n-2}\)。
暴力做是 \(O(n)\) 的我们可以设计矩阵来优化——设计一个初始矩阵 \(A\) 表示数列的初值,然后再设计一个矩阵 \(B\) 使得每次初始矩阵乘上其之后都会变成下一个状态。
一般转移有几步就会设多大的矩阵。
比如一个斐波那契数列的转移,我们可以都早这样一个 \(A= \begin{bmatrix} 1&1\\ 0&0 \end{bmatrix} .\)
再构造这样一个 \(B=\begin{bmatrix}1&1\\1&0\end{bmatrix}.\)
我们发现 \(A \times B=\begin{bmatrix}2&1\\0&0\end{bmatrix}.\)
可是这样有什么用呢?
应为矩阵的乘法也满足交换律和结合律,所以也可以用快速幂来优化乘法。
比如我们要求一个递推数列的第 \(M\) 项,这个数列的递推式有 \(N\) 项,则我们可以在 \(O(N^3 \times \log{M})\) 的时间复杂度内求得。
例题:数列
Code:AC记录
2 矩阵快速幂优化dp
很多dp题都是递推数列,所以就可以用矩阵快速幂来做。
比如随便举个dp的例子:
dp方程:\(dp[i][j] += dp[i-1][k]\times [w(j,k)==1]\)。要求 \(dp[n][m]\),\(n\le 10^9,m\le 100\)。
很显然这是个递推数列(从状态 \(i-1\) 推到 \(i\)),所以我们可以用矩阵乘法来优化。
我们把所有排斥的部分设为 \(0\),其余设为 \(1\) 然后做矩阵乘法即可。
3 图上的矩阵快速幂
我们最开始学存图的时候肯定都是用邻接矩阵存的,矩阵,矩阵!!!
所以图上一个点能到的点就构成一个矩阵结构。
例题:(NOI online3# T2)魔法值(矩阵快速幂只能拿暴力的 \(40pts\),可以先尝试写一下)
4 一些很妙的优化
1 常数优化
如果将做矩阵乘法时的代码改成k在j前循环,则内存的访问是连续的,可以有一点优化效果(虽然大部分时候是用不到的)。
2 一些玄学优化、
照刚刚的思想,我们发现如果 \(A[i][k]=0\) 则最后一步j的循环时不需要的,可以直接 continue
。
这样做看似没有什么大作用,其实很有用。
因为原本的的 \(A\) 矩阵大部分时间其实是个向量,也就是用 \(A \times B\) 的时间只是 \(O(N^2)\)。
本质:向量乘矩阵是平方级别的。
3 倍增处理 \(B\) 矩阵
针对多组数据。
如果把 \(B\) 矩阵先用倍增处理好,然后对于每一个查询都直接用 \(A\) 矩阵乘上那些处理好的矩阵则复杂度变为 \(O(N^2 \log{M})\)。
可以过掉上面的例题
#include <iostream>
#include <cstdio>
using namespace std;
typedef long long ll;
const ll N = 110, M = 10010;
struct Matrix{
ll size;
ll a[N][N];
}ans1;
Matrix ans2[35];
struct node{
ll pre, to;
}edge[M << 1];
ll head[N], tot;
ll n, m, q;
ll p[N];
ll f[N];
void init(Matrix &x) {
for (ll i = 1; i <= x.size; i++) {
for (ll j = 1; j <= x.size; j++) {
x.a[i][j] = 0;
}
}
}
Matrix Mul(Matrix x, Matrix y) {
Matrix ret;
ret.size = x.size;
init(ret);
for (ll i = 1; i <= x.size; i++) {
for (ll k = 1; k <= x.size; k++) {
if (x.a[i][k] == 0) continue;
for (ll j = 1; j <= x.size; j++) {
ret.a[i][j] ^= x.a[i][k] * y.a[k][j];
}
}
}
return ret;
}
Matrix ksm(Matrix x, ll y) {
Matrix ret;
ret.size = x.size;
init(ret);
for (ll i = 1; i <= n; i++) {
ret.a[i][i] = 1;
}
while (y) {
if (y & 1) ret = Mul(ret, x);
x = Mul(x, x);
y >>= 1;
}
return ret;
}
void Init_ans(Matrix &x) {
x.size = n;
for (ll i = 1; i <= n; i++) {
for (ll j = head[i]; j; j = edge[j].pre) {
x.a[edge[j].to][i] = 1;
}
}
}
void deburg(Matrix x) {
for (ll i = 1; i <= x.size; i++) {
for (ll j = 1; j <= x.size; j++) {
cout << x.a[i][j] << " ";
} cout << endl;
} cout << endl;
}
ll read() {
ll ret = 0, w = 1;
char ch = getchar();
while (!isdigit(ch)) {
if (ch == '-') w = -1;
ch = getchar();
}
while (isdigit(ch)) {
ret = (ret << 1) + (ret << 3) + ch - '0';
ch = getchar();
}
return ret * w;
}
void write(ll x) {
if (x > 9) write(x / 10);
putchar(x % 10 + '0');
}
void print(ll x) {
if (x < 0) {
x = -x;
putchar('-');
}
write(x);
putchar('\n');
}
void add(ll u, ll v) {
edge[++tot].pre = head[u];
edge[tot].to = v;
head[u] = tot;
}
int main() {
// freopen("magic.in", "r", stdin);
// freopen("magic.out", "w", stdout);
n = read();
m = read();
q = read();
for (ll i = 1; i <= n; i++) f[i] = read();
for (ll i = 1, u, v; i <= m; i++) {
u = read(); v = read();
add(u, v);
add(v, u);
}
for (int i = 1; i <= q; i++) {
p[i] = read();
}
Init_ans(ans2[0]);
for (int i = 1; i <= 31; i++) {
ans2[i] = Mul(ans2[i - 1], ans2[i - 1]);
}
for (int t = 1; t <= q; t++) {
ans1.size = n;
init(ans1);
for (ll i = 1; i <= n; i++) {
for (ll j = 1; j <= n; j++) {
if (i == 1) ans1.a[i][j] = f[j];
else ans1.a[i][j] = 0;
}
}
for (int j = 31; j >= 0; j--) {
if (((p[t] >> j) & 1) == 1) {
ans1 = Mul(ans1, ans2[j]);
}
}
print(ans1.a[1][1]);
}
return 0;
}
5 和其它算法混用
例题:GT考试(kmp+dp+矩阵快速幂优化)
6 其他应用
广义的矩阵乘法。
如果我们把矩阵乘法的+*变为min+那有没有交换律呢(只有满足交换律才能用快速幂)。
其实是满足的,证明直接展开即可得证。
事实上只要满足分配律即有交换律(min(a,b)+c=min(a+c,b+c))。
应用:Cow Relays
非常的模板,直接套min+矩阵快速幂即可。
#include <bits/stdc++.h>
using namespace std;
const int inf = 0x3f3f3f3f;
struct Matrix{
int a[110][110], size;
}g;
int bl[1010];
int n, m, s, e;
int tot;
int u[110], v[110], l[110];
int read() {
int ret = 0, f = 1;
char ch = getchar();
while (!isdigit(ch)) {
if (ch == '-') {
f = -1;
}
ch = getchar();
}
while (isdigit(ch)) {
ret = (ret << 1) + (ret << 3) + ch - '0';
ch = getchar();
}
return ret * f;
}
void write(int x) {
if (x > 9) write(x / 10);
putchar(x % 10 + '0');
}
void print(int x) {
if (x < 0) {
x = -x;
putchar('-');
}
write(x);
}
Matrix init(Matrix a, int sz) {
a.size = sz;
for (int i = 1; i <= a.size; i++) {
for (int j = 1; j <= a.size; j++) {
a.a[i][j] = inf;
}
}
return a;
}
Matrix Mul(Matrix a, Matrix b) {
Matrix ret;
ret = init(ret, a.size);
for (int i = 1; i <= a.size; i++) {
for (int j = 1; j <= a.size; j++) {
for (int k = 1; k <= a.size; k++) {
ret.a[i][j] = min(ret.a[i][j], a.a[i][k] + b.a[k][j]);
}
}
}
return ret;
}
Matrix ksm(Matrix a, int b) {
Matrix ret;
ret = init(ret, a.size);
for (int i = 1; i <= ret.size; i++) {
ret.a[i][i] = 0;
}
while (b) {
if (b & 1) ret = Mul(ret, a);
a = Mul(a, a);
b >>= 1;
}
return ret;
}
void deburg(Matrix a) {
for (int i = 1; i <= a.size; i++) {
for (int j = 1; j <= a.size; j++) {
cout << a.a[i][j] << " ";
} cout << endl;
} cout << endl;
}
int main() {
n = read(), m = read(), s = read(), e = read();
for (int i = 1; i <= m; i++) {
l[i] = read();
u[i] = read();
v[i] = read();
if (!bl[u[i]]) bl[u[i]] = ++tot;
if (!bl[v[i]]) bl[v[i]] = ++tot;
}
g = init(g, tot);
for (int i = 1; i <= m; i++) {
g.a[bl[u[i]]][bl[v[i]]] = l[i];
g.a[bl[v[i]]][bl[u[i]]] = l[i];
}
g = ksm(g, n);
print(g.a[bl[s]][bl[e]] + 1);
return 0;
}
由此,我们还可以得出minmin,minmax等矩阵乘法。
7 一、扩展
yhc的论文,还在看
标签:int,ll,矩阵,times,某道题,Tutorial,alpha,向量 From: https://www.cnblogs.com/zcr-blog/p/16758828.html