首页 > 其他分享 >[Tutorial] 从某道题谈矩阵快速幂及其优化

[Tutorial] 从某道题谈矩阵快速幂及其优化

时间:2022-10-06 23:34:57浏览次数:49  
标签:int ll 矩阵 times 某道题 Tutorial alpha 向量

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

相关文章

  • python中的矩阵乘法
    1.np.multiply()函数 矩阵的对应位置相乘,如果其中一个矩阵的尺寸不够,会自动广播,但是尺寸不能广播就会报错2.np.dot()函数 矩阵的点积,又称数量积、标量积或内积,即一......
  • 最大子矩阵和 = 前缀和 + 最大子段和
    简单来说这道题就是求一个\(N\timesM\)的矩阵的最大子矩阵和。(因为求的是黑色石板与白色石板的数量差,所以代表白色石板的“0”可以看作-1,这样就将问题转化为了求最大......
  • 矩阵计算
    矩阵的二范数是他的最大特征值 https://www.zhihu.com/question/48945813/answer/113453186 矩阵的F范数等于矩阵的迹: ......
  • 关于Hessian矩阵的图像增强
    文章目录​​1.数字图像处理之尺度空间理论​​​​2.基于尺度理论的Hessian简化算法​​​​3.基于Hessian矩阵的图像增强​​本文是关于图像增强方面的知识。关于Ret......
  • 【机器学习中的矩阵求导】(八)标量函数f(x)的雅克比矩阵(迹函数)
    学习总结交换律:,需要满足、同维度行列式微分:文章目录​​学习总结​​​​一、标量函数的雅克比函数​​​​二、关于迹函数的性质​​​​2.1常用性质​​​​2.2迹函数的......
  • [CG从零开始] 5. 搞清 MVP 矩阵理论 + 实践
    在4中成功绘制了三角形以后,下面我们来加载一个fbx文件,然后构建MVP变换(model-view-projection)。简单介绍一下:从我们拿到模型(主要是网格信息)文件开始,模型网格(Mesh)里......
  • 位姿矩阵求逆(转)
    位姿矩阵(或者称为旋转平移矩阵)即若干旋转矩阵和平移矩阵的合成,可以用来描述物体的方位。位姿矩阵具有形式:   且其中3*3部分   是一个正交阵,表示合旋转。(......
  • 数据组中的矩阵两个对角线的和(优化版)
    #include<stdio.h>intmain(){ intMAP[4][4]={12,11,14,15, 21,22,23,24, 65,55,66,77, 87,98,51,32}; inti,j; intsu......
  • 矩阵乘法(快速幂)结合dp结合除法逆元的例题
    https://atcoder.jp/contests/abc271/tasks/abc271_g题目的思路为:构建dp矩阵,dp[i][j][k]表示开始前停在j,结束后停在k,且停下时恰好出现2^i次访问的概率则dp[i]=dp[i-1]*d......
  • ZenFS tutorial
    ZenFSgitclonehttps://github.com/facebook/rocksdb.gitcdrocksdbgitclonehttps://github.com/westerndigitalcorporation/zenfsplugin/zenfssudoDEBUG_LEVEL=......