跳到主要内容

特殊多项式

特征的这部分只研究方阵，即矩阵 AA 对应的线性变换将 nn 个向量映射到 nn 个向量。

由于在实际问题中，经常要考虑连续进行重复的变换，如果只用「矩阵 AA 对应的线性变换将单位阵 II 变换为 AA」的描述，就会很抽象。此时最好的办法是找「不动点」，即变换当中不动的部分。

然而事实上，矩阵 AA 对应的线性变换很可能没有不动点，于是退而求其次，寻找共线或者类似于简单变形的部分。

特征值与特征向量​

在矩阵 AA 对应的线性变换作用下，一些向量的方向不改变，只是伸缩了。

设 VV 是 FF 上的线性空间，TT 是 VV 上的线性变换。若存在 FF 中的 λ\lambda 与 VV 中的 非零向量 ξ\xi，使得：

Tξ=λξT\xi=\lambda\xi

则称 λ\lambda 为 TT 的一个 特征值，而 ξ\xi 为 TT 的 属于特征值 λ\lambda 的一个特征向量。

特征向量在同一直线上，在线性变换作用下保持方向不改变（压缩到零也认为是方向不改变）。特征向量不唯一，与特征向量共线的向量都是特征向量，但是规定零向量不是特征向量，拥有方向的向量自然是非零向量。特征向量的特征值就是它伸缩的倍数。

在实际应用中，一般对于拥有相同特征值的特征向量，会选取一组基作为它们全体的代表。

设 α1,α2,⋯ ,αn\alpha_1,\alpha_2,\cdots,\alpha_n 是 VV 的一组基，TT 在这组基下的矩阵为 AA，即：

T(α1,α2,⋯ ,αn)=(α1,α2,⋯ ,αn)AT(\alpha_1,\alpha_2,\cdots,\alpha_n)=(\alpha_1,\alpha_2,\cdots,\alpha_n)A

设 λ0\lambda_0 是 TT 的一个特征值，ξ\xi 为 TT 的属于特征值 λ0\lambda_0 的一个特征向量，且有非零向量 XX 满足：

ξ=(α1,α2,⋯ ,αn)X\xi=(\alpha_1,\alpha_2,\cdots,\alpha_n)X

于是有：

Tξ=λ0ξT\xi=\lambda_0\xi T(α1,α2,⋯ ,αn)X=λ0(α1,α2,⋯ ,αn)XT(\alpha_1,\alpha_2,\cdots,\alpha_n)X=\lambda_0(\alpha_1,\alpha_2,\cdots,\alpha_n)X (α1,α2,⋯ ,αn)AX=λ0(α1,α2,⋯ ,αn)X(\alpha_1,\alpha_2,\cdots,\alpha_n)AX=\lambda_0(\alpha_1,\alpha_2,\cdots,\alpha_n)X AX=λ0XAX=\lambda_0X (A−λ0I)X=0(A-\lambda_0I)X=0

所以相应的行列式也为 00。

特征多项式​

考虑一个 n×nn\times n 的矩阵 AA，其中 n≥0∧n∈Zn\geq 0\land n\in\mathbb{Z}。设 λ\lambda 为一个参量，矩阵 λI−A\lambda I-A 称为 AA 的 特征矩阵。

特征矩阵的行列式称为 AA 的 特征多项式，展开为一个 nn 次多项式，根为 AA 的特征值，记为 pA(λ)p_A(\lambda)：

pA(λ)=det⁡(λIn−A)=∣λ−a11−a12⋯−a1n−a21λ−a22⋯−a2n⋮⋮⋮−an1−an2⋯λ−ann∣p_A(\lambda)=\det(\lambda I_n-A)=\begin{vmatrix} \lambda-a_{11} & -a_{12} & \cdots & -a_{1n} \\ -a_{21} & \lambda-a_{22} & \cdots & -a_{2n} \\ \vdots & \vdots & & \vdots \\ -a_{n1} & -a_{n2} & \cdots & \lambda-a_{nn} \\ \end{vmatrix}

其中 InI_n 为一个 n×nn\times n 的单位矩阵。一些地方会定义为 pA(λ)=det⁡(A−λIn)p_A(\lambda)=\det(A-\lambda I_n) 与我们的定义仅相差了一个符号 (−1)n(-1)^n，但采用这种定义得到的 pA(λ)p_A(\lambda) 一定为首一多项式，而另外的定义则仅当 nn 为偶数时才是首一多项式。需要注意的是 0×00\times 0 的矩阵行列式为 11 是良定义的。

相应于 (λ0I−A)X=0(\lambda_0 I-A)X=0 的非零解向量 XX，称为 AA 的属于 λ0\lambda_0 的特征向量。

线性变换 TT 有特征值 λ0\lambda_0 等价于矩阵 AA 有特征值 λ0\lambda_0。

线性变换 TT 有特征向量 ξ\xi 等价于矩阵 AA 有特征向量 XX，其中有：

ξ=(α1,⋯ ,αn)X\xi=(\alpha_1,\cdots,\alpha_n)X

根据代数基本定理，特征多项式可以分解为：

f(λ)=∣λI−A∣=(λ−λ1)d1⋯(λ−λm)dmf(\lambda)=|\lambda I-A|={(\lambda-\lambda_1)}^{d_1}\cdots{(\lambda-\lambda_m)}^{d_m}

称 did_i 为特征值 λi\lambda_i 的 代数重数。全体代数重数的和为空间维数 nn。

求解矩阵的全部特征值及特征向量​

分为以下步骤：

  • 计算行列式 ∣λI−A∣|\lambda I-A|。
  • 求出多项式 f(λ)=∣λI−A∣f(\lambda)=|\lambda I-A| 在域 FF 中的全部根，即 AA 的特征值。
  • 对 AA 的每个特征值 λ\lambda，解齐次线性方程组 (λI−A)X=0(\lambda I-A)X=0，求出它的一组基础解系 X1,⋯ ,XtX_1,\cdots,X_t，则 AA 的属于 λ\lambda 的全部特征向量为：
k1X1+k2X2+⋯+ktXtk_1X_1+k_2X_2+\cdots+k_tX_t

该表达式中的 kik_i 不全为零。

  • 线性变换 TT 的属于 λ\lambda 的特征向量为：
ξi=(α1,⋯ ,αn)Xi\xi_i=(\alpha_1,\cdots,\alpha_n)X_i

因此，属于 λ\lambda 的全部特征向量为：

k1ξ1+k2ξ2+⋯+ktξtk_1\xi_1+k_2\xi_2+\cdots+k_t\xi_t

该表达式中的 kik_i 不全为零。

特征值与特征向量是否存在，依赖于 VV 所在的域。

相似变换​

引入​

若 n×nn\times n 的矩阵 AA 为上三角矩阵如

A=[a1,1a1,2⋯a1,na2,2⋯a2,n⋱⋮an,n]A= \begin{bmatrix} a_{1,1}&a_{1,2}&\cdots &a_{1,n}\\ &a_{2,2}&\cdots &a_{2,n}\\ &&\ddots &\vdots \\ &&&a_{n,n} \end{bmatrix}

那么

pA(x)=det⁡(xIn−A)=[x−a1,1−a1,2⋯−a1,nx−a2,2⋯−a2,n⋱⋮x−an,n]=∏i=1n(x−ai,i)\begin{aligned} p_A(x)&=\det(xI_n-A)\\ &= \begin{bmatrix} x-a_{1,1}&-a_{1,2}&\cdots &-a_{1,n}\\ &x-a_{2,2}&\cdots &-a_{2,n}\\ &&\ddots &\vdots \\ &&&x-a_{n,n} \end{bmatrix} \\ &=\prod_{i=1}^n(x-a_{i,i}) \end{aligned}

可轻松求得，下三角矩阵也是类似的。但如果 AA 不属于这两种矩阵，则需要使用相似变换，使得矩阵变为容易求得特征多项式的形式。

定义​

对于 n×nn\times n 的矩阵 AA 和 BB，当存在 n×nn\times n 的可逆矩阵 PP 满足

B=P−1APB=P^{-1}AP

则矩阵 AA 和 BB 相似，记变换 A↦P−1APA\mapsto P^{-1}AP 为相似变换。且 AA 和 P−1APP^{-1}AP 有相同的特征多项式。

考虑

det⁡(xIn−P−1AP)=det⁡(xP−1InP−P−1AP)=det⁡(P−1xInP−P−1AP)=det⁡(P−1)⋅det⁡(P)⋅det⁡(xIn−A)=det⁡(xIn−A)=pA(x)\begin{aligned} \det(xI_n-P^{-1}AP)&=\det(xP^{-1}I_nP-P^{-1}AP)\\ &=\det(P^{-1}xI_nP-P^{-1}AP)\\ &=\det(P^{-1})\cdot \det(P)\cdot \det(xI_n-A)\\ &=\det(xI_n-A)\\ &=p_A(x) \end{aligned}

得证，对于 A↦PAP−1A\mapsto PAP^{-1} 也是一样的。另外 pA(0)=(−1)n⋅det⁡(A)p_A(0)=(-1)^n\cdot \det(A)，因为 pA(0)=det⁡(−1⋅InA)=det⁡(−1⋅In)⋅det⁡(A)p_A(0)=\det(-1\cdot I_nA)=\det(-1\cdot I_n)\cdot \det(A) 故 det⁡(A)=det⁡(P−1AP)\det(A)=\det(P^{-1}AP)。

定理：相似矩阵有相同的特征多项式及特征值，反之不然。

定理表明，线性变换的矩阵的特征多项式与基的选取无关，而直接由线性变换决定，故可称之为线性变换的特征多项式。

矩阵 AA 的特征多项式 f(λ)=∣λI−A∣f(\lambda)=|\lambda I-A| 是一个首一的多项式。根据韦达定理，它的 n−1n-1 次系数为：

−(λ1+⋯+λn)=−(a11+⋯+ann)=−trA-(\lambda_1+\cdots+\lambda_n)=-(a_{11}+\cdots+a_{nn})=-tr A

其中 trAtr A 称为 AA 的迹，为 AA 的主对角线元素之和。

根据韦达定理，特征多项式的常数项为：

(−1)n∣A∣=(−1)n(λ1⋯λn){(-1)}^n|A|={(-1)}^n(\lambda_1\cdots\lambda_n)

定理：相似的矩阵有相同的迹。

换位公式​

定理：无论矩阵 AA 和矩阵 BB 是否方阵，只要乘法能进行，则矩阵 ABAB 的迹等于矩阵 BABA 的迹。

一种证法是直接展开，即证毕。另一种证法用到换位公式。

定理：设 AA 为 mm 行 nn 列矩阵，设 BB 为 nn 行 mm 列矩阵，则有：

λn∣λIm−AB∣=λm∣λIn−BA∣\lambda^n|\lambda I_m-AB|=\lambda^m|\lambda I_n-BA|

该公式表明 ABAB 与 BABA 有相同的非零特征值。

舒尔（Schur）引理​

任意的 nn 阶矩阵 AA 都相似于一个上三角阵，即存在满秩阵 PP，使得 P−1APP^{-1}AP 为上三角阵，它的主对角线上元素为 AA 的全部特征值。

推论：设 AA 的 nn 个特征值为 λ1,⋯ ,λn\lambda_1,\cdots,\lambda_n，ϕ(x)\phi(x) 为任一多项式，则矩阵多项式 ϕ(A)\phi(A) 的 nn 个特征值为：

ϕ(λ1),⋯ ,ϕ(λn)\phi(\lambda_1),\cdots,\phi(\lambda_n)

特别地，kAkA 的特征值为 kλ1,⋯ ,kλnk\lambda_1,\cdots,k\lambda_n，AmA^m 的特征值为 λ1m,⋯ ,λnm{\lambda_1}^m,\cdots,{\lambda_n}^m。

使用高斯消元进行相似变换​

对 n×nn\times n 的矩阵 BB 可以进行高斯消元，其基本操作为初等行变换。

在对矩阵使用上述操作（左乘初等矩阵）后再右乘其逆矩阵即相似变换，左乘为行变换，易发现右乘即列变换。

若能将矩阵通过相似变换变为上三角或下三角的形式，那么可以轻松求出其特征多项式。但若对主对角线上的元素应用变换 A↦Tij(k)ATij(−k)A\mapsto T_{ij}(k)AT_{ij}(-k) 后会导致原本通过 A↦Tij(k)AA\mapsto T_{ij}(k)A 将第 ii 行第 jj 列的元素消为零后右乘 Tij(−k)T_{ij}(-k) 即将 AA 的第 ii 列的 −k-k 倍加到第 jj 列这一操作使得之前消为零的元素现在可能不为零，可能不能将其变为上三角或下三角形式。

后文将说明对次对角线上的元素应用变换后得到的矩阵依然可以轻松得到其特征多项式。

上 Hessenberg 矩阵​

对于 n>2n\gt 2 的形如

H=[α1h12……h1nβ2α2h23…⋮⋱⋱⋱⋮⋱⋱h(n−1)nβnαn]H= \begin{bmatrix} \alpha_{1}&h_{12}&\dots&\dots&h_{1n}\\ \beta_{2}&\alpha_{2}&h_{23}&\dots &\vdots \\ &\ddots &\ddots & \ddots &\vdots \\ & &\ddots &\ddots & h_{(n-1)n}\\ &&& \beta_{n}& \alpha_{n} \end{bmatrix}

的矩阵我们称为上 Hessenberg 矩阵，其中 β\beta 为次对角线。

我们使用相似变换将次对角线以下的元素消为零后即能得到上 Hessenberg 矩阵，而求出一个 n×nn\times n 上 Hessenberg 矩阵的特征多项式则可在 O(n3)O(n^3) 时间完成。

我们记 HiH_i 为只保留 HH 的前 ii 行和前 ii 列的矩阵，记 pi(x)=det⁡(xIi−Hi)p_i(x)=\det(xI_i-H_i) 那么

H0=[],p0(x)=1H_0= \begin{bmatrix} \end{bmatrix},\quad p_0(x)=1 H1=[α1],p1(x)=det⁡(xI1−H1)=x−α1H_1= \begin{bmatrix} \alpha_1 \end{bmatrix},\quad p_1(x)=\det(x I_1-H_1)=x -\alpha_1 H2=[α1h12β2α2],p2(x)=det⁡(xI2−H2)=(x−α2)p1(x)−β2h12p0(x)H_2= \begin{bmatrix} \alpha_1&h_{12}\\ \beta_2&\alpha_2 \end{bmatrix},\quad p_2(x)=\det(xI_2-H_2)=(x-\alpha_2)p_1(x)-\beta_2h_{12}p_0(x)

在计算行列式时我们一般选择按零最多的行或列余子式展开，余子式即删除了当前选择的元素所在行和列之后的矩阵，在这里我们选择按最后一行进行展开，有

p3(x)=det⁡(xI3−H3)=∣x−α1−h12−h13−β2x−α2−h23−β3x−α3∣=(x−α3)⋅(−1)3+3p2(x)−β3⋅(−1)3+2∣x−α1−h13−β2−h23∣=(x−α3)p2(x)−β3(h23p1(x)+β2h13p0(x))\begin{aligned} p_3(x)&= \det(xI_3-H_3)\\ &=\begin{vmatrix} x-\alpha_1&-h_{12}&-h_{13}\\ -\beta_2&x-\alpha_2&-h_{23}\\ &-\beta_3&x-\alpha_3 \end{vmatrix}\\ &=(x-\alpha_3)\cdot (-1)^{3+3}p_2(x)-\beta_3\cdot (-1)^{3+2} \begin{vmatrix} x-\alpha_1&-h_{13}\\ -\beta_2&-h_{23} \end{vmatrix}\\ &=(x-\alpha_3)p_2(x)-\beta_3(h_{23}p_1(x)+\beta_2h_{13}p_0(x)) \end{aligned}

观察并归纳，对 2≤i≤n2\leq i\leq n 有

pi(x)=(x−αi)pi−1(x)−∑m=1i−1hi−m,i(∏j=i−m+1iβj)pi−m−1(x)p_i(x)=(x-\alpha_i)p_{i-1}(x)- \sum_{m=1}^{i-1}h_{i-m,i} \left( \prod_{j=i-m+1}^{i}\beta_j \right) p_{i-m-1}(x)

至此完成了整个算法，该算法一般被称为 Hessenberg 算法。

Cayley–Hamilton 定理​

对于任意的 nn 阶矩阵 AA，特征多项式为 f(λ)=∣λI−A∣f(\lambda)=|\lambda I-A|，则必有 f(A)=0f(A)=0。

对于线性变换 TT 有平行的结果：如果 f(λ)f(\lambda) 为 TT 的特征多项式，则 f(T)f(T) 为零变换。

由本定理可知，对于任意的矩阵 AA，必有可以使其零化的多项式。

最小多项式​

设 VV 是一个 nn 维向量空间，由于线性变换对应的矩阵有 n2n^2 个元素，一切线性变换构成 n2n^2 维线性空间。

对于一个特定的线性变换 TT，从作用 00 次到作用 nn 次，总共 n2+1n^2+1 个线性变换，它们对应的矩阵一定线性相关。于是存在非零多项式 ff，使得 f(T)f(T) 为零变换，称变换 TT 满足多项式 ff。在 TT 满足的所有多项式 ff 中，存在次数最低的。

可以将矩阵 AA 零化的最小次数的首一多项式称为 AA 的最小多项式，记为 mA(λ)m_A(\lambda)。

根据多项式的辗转相除法，最小多项式是唯一的，且可整除任一 AA 的零化多项式。特别地，最小多项式整除特征多项式。

定理：在不计重数的情况下，矩阵 AA 的特征多项式 f(λ)f(\lambda) 与最小多项式 mA(λ)m_A(\lambda) 有相同的根。

定理：矩阵 AA 的属于不同特征值的特征向量线性无关。

应用​

在信息学中我们一般考虑 (Z/mZ)n×n(\mathbb{Z}/m\mathbb{Z})^{n\times n} 上的矩阵，通常 mm 为素数，进行上述相似变换是简单的，当 mm 为合数时，我们可以考虑类似辗转相除的方法来进行。

实现代码
实现代码
#include <cassert>
#include <iostream>
#include <random>
#include <vector>

using Matrix = std::vector<std::vector<int>>;
using i64 = int64_t;

Matrix to_upper_Hessenberg(const Matrix &M, int mod) {
  Matrix H(M);
  int n = H.size();
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < n; ++j) {
      if ((H[i][j] %= mod) < 0) H[i][j] += mod;
    }
  }
  for (int i = 0; i < n - 1; ++i) {
    int pivot = i + 1;
    for (; pivot < n; ++pivot) {
      if (H[pivot][i] != 0) break;
    }
    if (pivot == n) continue;
    if (pivot != i + 1) {
      for (int j = i; j < n; ++j) std::swap(H[i + 1][j], H[pivot][j]);
      for (int j = 0; j < n; ++j) std::swap(H[j][i + 1], H[j][pivot]);
    }
    for (int j = i + 2; j < n; ++j) {
      for (;;) {
        if (H[j][i] == 0) break;
        if (H[i + 1][i] == 0) {
          for (int k = i; k < n; ++k) std::swap(H[i + 1][k], H[j][k]);
          for (int k = 0; k < n; ++k) std::swap(H[k][i + 1], H[k][j]);
          break;
        }
        if (H[j][i] >= H[i + 1][i]) {
          int q = H[j][i] / H[i + 1][i], mq = mod - q;
          for (int k = i; k < n; ++k)
            H[j][k] = (H[j][k] + i64(mq) * H[i + 1][k]) % mod;
          for (int k = 0; k < n; ++k)
            H[k][i + 1] = (H[k][i + 1] + i64(q) * H[k][j]) % mod;
        } else {
          int q = H[i + 1][i] / H[j][i], mq = mod - q;
          for (int k = i; k < n; ++k)
            H[i + 1][k] = (H[i + 1][k] + i64(mq) * H[j][k]) % mod;
          for (int k = 0; k < n; ++k)
            H[k][j] = (H[k][j] + i64(q) * H[k][i + 1]) % mod;
        }
      }
    }
  }
  return H;
}

std::vector<int> get_charpoly(const Matrix &M, int mod) {
  Matrix H(to_upper_Hessenberg(M, mod));
  int n = H.size();
  std::vector<std::vector<int>> p(n + 1);
  p[0] = {1 % mod};
  for (int i = 1; i <= n; ++i) {
    const std::vector<int> &pi_1 = p[i - 1];
    std::vector<int> &pi = p[i];
    pi.resize(i + 1, 0);
    int v = mod - H[i - 1][i - 1];
    if (v == mod) v -= mod;
    for (int j = 0; j < i; ++j) {
      pi[j] = (pi[j] + i64(v) * pi_1[j]) % mod;
      if ((pi[j + 1] += pi_1[j]) >= mod) pi[j + 1] -= mod;
    }
    int t = 1;
    for (int j = 1; j < i; ++j) {
      t = i64(t) * H[i - j][i - j - 1] % mod;
      int prod = i64(t) * H[i - j - 1][i - 1] % mod;
      if (prod == 0) continue;
      prod = mod - prod;
      for (int k = 0; k <= i - j - 1; ++k)
        pi[k] = (pi[k] + i64(prod) * p[i - j - 1][k]) % mod;
    }
  }
  return p[n];
}

bool verify(const Matrix &M, const std::vector<int> &charpoly, int mod) {
  if (mod == 1) return true;
  int n = M.size();
  std::vector<int> randvec(n), sum(n, 0);
  std::mt19937 gen(std::random_device{}());
  std::uniform_int_distribution<int> dis(1, mod - 1);
  for (int i = 0; i < n; ++i) randvec[i] = dis(gen);
  for (int i = 0; i <= n; ++i) {
    int v = charpoly[i];
    for (int j = 0; j < n; ++j) sum[j] = (sum[j] + i64(v) * randvec[j]) % mod;
    std::vector<int> prod(n, 0);
    for (int j = 0; j < n; ++j) {
      for (int k = 0; k < n; ++k) {
        prod[j] = (prod[j] + i64(M[j][k]) * randvec[k]) % mod;
      }
    }
    randvec.swap(prod);
  }
  for (int i = 0; i < n; ++i)
    if (sum[i] != 0) return false;
  return true;
}

int main() {
  std::ios::sync_with_stdio(false);
  std::cin.tie(nullptr);
  int n, mod;
  std::cin >> n >> mod;
  Matrix M(n, std::vector<int>(n));
  for (int i = 0; i < n; ++i)
    for (int j = 0; j < n; ++j) std::cin >> M[i][j];
  std::vector<int> charpoly(get_charpoly(M, mod));
  for (int i = 0; i <= n; ++i) std::cout << charpoly[i] << ' ';
  assert(verify(M, charpoly, mod));
  return 0;
}

上述 Hessenberg 算法不具有数值的稳定性，所以 Rn×n\mathbb{R}^{n\times n} 上的矩阵在使用前需要其他算法进行调整或改用其他具有数值稳定性的算法。

我们可以将特征多项式与常系数齐次线性递推联系起来，也可结合 Cayley–Hamilton 定理、多项式取模加速一些域上求矩阵幂次的算法。

Cayley–Hamilton 定理指出

pA(A)=An+c1An−1+⋯+cn−1A+cnI=O\begin{aligned} p_A(A)&=A^n+c_1A^{n-1}+\cdots +c_{n-1}A+c_nI\\ &=O \end{aligned}

其中 OO 为 n×nn\times n 的零矩阵，A∈Cn×nA\in\mathbb{C}^{n\times n} 且 pA(x)=xn+∑i=1ncixn−i∈C[x]p_A(x)=x^n+\sum_{i=1}^nc_ix^{n-i}\in\mathbb{C}[x] 为 AA 的特征多项式。

若我们要求 AKA^K 其中 KK 较大，那么可以求出 f(x)=xK mod pA(x)f(x)=x^K\bmod{p_A(x)} 后利用 f(A)=AKf(A)=A^K。

而 deg⁡(f(x))<n\deg(f(x))\lt n 显然。我们令 f(x)=∑i=0n−1fixif(x)=\sum_{i=0}^{n-1}f_ix^i 且 n=kmn=km 那么

fkm−1xkm−1+⋯+f1x+f0=(⋯(fkm−1xk−1+⋯+fk(m−1))xk+fk(m−1)−1xk−1+⋯+fk(m−2))xk+⋯+fk−1xk−1+⋯+f1x+f0\begin{aligned} f_{km-1}x^{km-1}+\cdots +f_1x+f_0&=(\cdots (f_{km-1}x^{k-1}+\cdots +f_{k(m-1)})x^k\\ &+f_{k(m-1)-1}x^{k-1}+\cdots +f_{k(m-2)})x^k\\ &+\cdots\\ &+f_{k-1}x^{k-1}+\cdots +f_1x+f_0 \end{aligned}

令 k=nk=\sqrt{n} 可以发现计算 f(A)f(A) 大约需要 O(n)O(\sqrt{n}) 次矩阵与矩阵的乘法。