矩阵运算

矩阵运算

OI-wiki Link

众所周知,数是可以进行加减乘除的,那矩阵为啥不可以呢?

假设现在我们有两个矩阵 \(A\)\(B\),矩阵大小分别为 \(n \times m\)\(x \times y\),矩阵元素对 \(mod\) 取模。

const int N = 110, mod = 1e9 + 7; // 矩阵大小和模数

struct matrix {
  int a[N][N], n, m; // 存储矩阵和矩阵大小
} ;

基本运算

矩阵加法

\(A + B = C\)

要求:\(n = x\) 并且 \(m = y\)

其实很简单,就是一一对应着加就行,即对于 \(1\leqslant i \leqslant n, 1\leqslant j \leqslant m\)\(C_{i, j} = A_{i, j} + B_{i, j}\)

所以 \(C\) 的大小为 \(n \times m\),矩阵减法同。

性质

  • 交换律:明显满足,即 \(A + B = B + A\)
  • 结合律:明显满足,即 \(A + B + C = A + (B + C)\)

单位矩阵

明显,就是一个大小为 \(n \times m\) 的全 \(0\) 矩阵。

矩阵加法 Code

时间复杂度:\(O(n \times m)\)

  matrix operator + (matrix y) { // 重载运算符 +
    matrix z;
    z.n = n, z.m = m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        z.a[i][j] = a[i][j] + y.a[i][j];
    return z;
  }

矩阵乘法

\(A \times B = C\)

要求:\(m = x\)

公式:对于 \(1 \leqslant i \leqslant n, 1 \leqslant j \leqslant y\):

\[C_{i, j} = \sum_{1 \leqslant k \leqslant m} A_{i, k} \times B_{k, j} \]

似乎看起来没啥用?因为它是在某些特定方面可以进行时间优化(例如用 \(O(8 \log n)\) 的时间复杂度求出斐波那契数列的第 \(n\) 项),所以一般就是把它和其他知识点捆绑起来(例如矩阵乘法优化 dp 等),单独拎出来出的题目很少(大多都是板子)。

性质

  • 交换律:不满足!
    • \(n = y\),交换以后还能乘法,但矩阵大小会变为 \(m \times x\),与原矩乘后答案不同。
    • 当然也有可能 \(n \neq y\),这样交换以后连矩乘的基本要求都无法满足。
  • 结合律:满足,设还有一个矩阵 \(C\),大小为 \(y \times z\)\(A \times B \times C = D\),则 \(D\) 大小为 \(n \times z\),而 \(A \times (B \times C) = A \times E(\texttt{大小为} x \times z) = D\)
  • 分配律:明显满足,设还有一个矩阵 \(C\),大小为 \(y \times z\),即 \(A \times (B + C) = A \times B + A \times C\)

单位矩阵

Link

提到矩乘,就不得不提到它的单位矩阵。

单位矩阵定义:若这个矩阵乘任意一个矩阵 \(C\),得到的结果都等于 \(C\),则这个矩阵为单位矩阵。

观察式子,容易推出:当一个大小为 \(n \times n\) 的矩阵除了主对角线(从左上到右下)以外的数为 \(1\),其他都为 \(0\) 时,它就是一个可以作为所有大小为 \(n \times m\) 的矩阵的单位矩阵。

时间复杂度:\(O(x \times y)\)

  void Clear (int x, int y, int f) { // 把矩阵大小设为 x * y,矩阵元素都为 f
    n = x, m = y;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        a[i][j] = f;
  }
  void Unit (int x, int y) { // 构造单位矩阵
    Clear(x, y, 0); // 先清空
    for (int i = 1; i <= n; i++)
      a[i][i] = 1; // 只有主对角线为 1
  }

矩阵乘法 Code

时间复杂度:\(O(n \times m \times y)\)

  matrix operator * (matrix y) { // 重载运算符 *
    matrix z;
    z.n = n, z.m = y.m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= y.m; j++) {
        z.a[i][j] = 0;
        for (int k = 1; k <= m; k++)
          z.a[i][j] = (z.a[i][j] + 1ll * a[i][k] * y.a[k][j] % mod) % mod; // 套公式
      }
    return z;
  }

矩阵除法

\(\frac{A}{B} = C\)

要求:同矩阵乘法。

由于浮点数不可以取模,数里面除法取模等于乘上除数的逆,所以矩阵除法同矩阵除法,不过是把 \(B_{i, j}\) 变成 \(B_{i, j}\) 在模 \(mod\) 意义下的逆。

矩乘与它的特殊运算和应用

P3390 矩乘快速幂

给定一个矩阵 \(A\),大小为 \(n \times m\),一个指数 \(k\) 和一个模数 \(mod\),求 \(A^k\ \%\ mod\)

要求:矩阵的大小必须为 \(n \times n\)

既然矩乘有结合律,那么就可以考虑使用类似于整数快速幂的方法以 \(O(\log k)\) 的方式求出 \(A^k\),由于一次矩阵乘法需要 \(O(n^3)\),所以总时间复杂度为 \(O(n ^ 3 \log k)\)

矩乘快速幂 Code

时间复杂度:\(O(n ^ 3 \log k)\)

matrix MATqmod (matrix x, ll y) { // 矩阵乘法快速幂
  matrix sum;
  sum.Unit(x.n, x.m); // 初始化为单位矩阵
  while (y) sum = (y & 1 ? sum * x : sum), x = x * x, y >>= 1;
  return sum;
}

P1349 广义斐波那契数列

P1962 斐波那契数列

广义的斐波那契数列是指形如 \(a_n=p\times a_{n-1}+q\times a_{n-2}\) 的数列。

今给定数列的两系数 \(p\)\(q\),以及数列的最前两项 \(a_1\) 和 $ a_2$,另给出两个整数 \(n\)\(m\),试求数列的第 \(n\)\(a_n \bmod m\)

\(F_i\) 表示当前情况下斐波那契数列第 \(i\) 项的值。

观察一下,把相邻两项构成一个矩阵,即要从 \(f_i = \begin{bmatrix}F_i,F_{i - 1}\end{bmatrix}\) 通过操作变为 \(f_{i + 1} = \begin{bmatrix}F_{i + 1},F_i\end{bmatrix} = \begin{bmatrix}p \times F_i + q \times F_{i - 1},F_i\end{bmatrix}\),考虑用矩乘完成。

要能与 \(\begin{bmatrix}F_i,F_{i - 1}\end{bmatrix}\) 做矩乘,就需要一个大小为 \(2 \times k\) 的矩阵,又因为矩乘之后的矩阵仍然大小为 \(1 \times 2\),所以 \(k = 2\),即要求一个矩阵 \(A = \begin{bmatrix}a,b\\c,d\end{bmatrix}\),使得 \(A \times f_i = f_{i + 1}\)

套入公式,就是要求 \(\begin{cases}a\times F_i + c \times F_{i - 1} = p \times F_i + q \times F_{i - 1}\\b\times F_i + d \times F_{i - 1}=F_i\end{cases}\),可以推出 \(\begin{cases}a=p\\b=1\\c=q\\d=0\end{cases}\),现在要求 \(F_n\),初始矩阵为\(\begin{bmatrix}a_2,a_1\end{bmatrix}\),每次乘矩阵 \(A = \begin{bmatrix}p,1\\q,0\end{bmatrix}\) 后都会使矩阵两项下标加 \(1\),特判 \(n = 1\) 的情况(答案为 \(a_1\)),那么答案就是 \(\begin{bmatrix}a_2,a_1\end{bmatrix} \times A^{n - 2}\) 取第一个元素即可,需要矩阵乘法快速幂。

普通的其实就相当于 \(a_1=a_2=p=q=1\)

Code

  matrix c, d;
  int n, mod;
  c.Clear(2, 2, 0), d.Clear(1, 2, 0);
  cin >> c.a[1][1] >> c.a[2][1] >> d.a[1][2] >> d.a[1][1] >> n >> mod;
  if (n == 1) {
    cout << d.a[1][2];
    return 0;
  }
  c.a[1][2] = 1, d = d * qmod(c, n - 2);
  cout << d.a[1][1];

完整代码

点击查看最新代码
#include <bits/stdc++.h>
#define _1 (__int128)1

using namespace std;
using ll = long long;

void FileIO (string s) {
  freopen((s + ".in").c_str(), "r", stdin);
  freopen((s + ".out").c_str(), "w", stdout);
}

const int mod = 1e9 + 7, MATsize = 110; // 矩阵大小和模数

struct matrix {
  int a[MATsize][MATsize], n, m; // 存储矩阵和矩阵大小
  void Clear (int x, int y, int f) { // 把矩阵大小设为 x * y,矩阵元素都为 f
    n = x, m = y;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        a[i][j] = f;
  }
  void Unit (int x, int y) { // 构造单位矩阵
    Clear(x, y, 0); // 先清空
    for (int i = 1; i <= n; i++)
      a[i][i] = 1; // 只有主对角线为 1
  }
  matrix operator + (matrix y) { // 重载运算符 +
    matrix z;
    z.n = n, z.m = m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        z.a[i][j] = a[i][j] + y.a[i][j];
    return z;
  }
  matrix operator * (matrix y) { // 重载运算符 *
    matrix z;
    z.n = n, z.m = y.m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= y.m; j++) {
        z.a[i][j] = 0;
        for (int k = 1; k <= m; k++)
          z.a[i][j] = (z.a[i][j] + 1ll * a[i][k] * y.a[k][j] % mod) % mod; // 套公式
      }
    return z;
  }
} ;

matrix MATqmod (matrix x, ll y) { // 矩阵乘法快速幂
  matrix sum;
  sum.Unit(x.n, x.m); // 初始化为单位矩阵
  while (y) sum = (y & 1 ? sum * x : sum), x = x * x, y >>= 1;
  return sum;
}

signed main () {
  ios::sync_with_stdio(0), cin.tie(0);
  // FileIO("");
  return 0;
}

点击查看无注释代码
#include <bits/stdc++.h>
#define _1 (__int128)1

using namespace std;
using ll = long long;

void FileIO (string s) {
  freopen((s + ".in").c_str(), "r", stdin);
  freopen((s + ".out").c_str(), "w", stdout);
}

const int mod = 1e9 + 7, MATsize = 110;

struct matrix {
  int a[MATsize][MATsize], n, m;
  void Clear (int x, int y, int f) {
    n = x, m = y;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        a[i][j] = f;
  }
  void Unit (int x, int y) {
    Clear(x, y, 0);
    for (int i = 1; i <= n; i++)
      a[i][i] = 1;
  }
  matrix operator + (matrix y) {
    matrix z;
    z.n = n, z.m = m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= m; j++)
        z.a[i][j] = a[i][j] + y.a[i][j];
    return z;
  }
  matrix operator * (matrix y) {
    matrix z;
    z.n = n, z.m = y.m;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j <= y.m; j++) {
        z.a[i][j] = 0;
        for (int k = 1; k <= m; k++)
          z.a[i][j] = (z.a[i][j] + 1ll * a[i][k] * y.a[k][j] % mod) % mod;
      }
    return z;
  }
} ;

matrix MATqmod (matrix x, ll y) {
  matrix sum;
  sum.Unit(x.n, x.m);
  while (y) sum = (y & 1 ? sum * x : sum), x = x * x, y >>= 1;
  return sum;
}

signed main () {
  ios::sync_with_stdio(0), cin.tie(0);
  // FileIO("");
  return 0;
}

posted @ 2023-08-05 11:51  wnsyou  阅读(58)  评论(0编辑  收藏  举报