近似矩阵乘法

近似矩阵乘法

矩阵乘法是线性代数中的基本运算，也是众多科学计算和机器学习算法的核心操作。然而，随着数据规模的增长，精确矩阵乘法的计算复杂度变得难以承受。近似矩阵乘法应运而生，它通过牺牲一定的精度来换取计算效率的提升。

近似矩阵乘法的基本概念

在介绍近似矩阵乘法之前，我们先回顾一下精确矩阵乘法。给定两个矩阵 $A\in {\mathbb{R}}^{m\times n}$ 和 $B\in {\mathbb{R}}^{n\times p}$，它们的乘积 $C=A\times B\in {\mathbb{R}}^{m\times p}$ 定义为：

$$C_{ij}=\sum _{k=1}^n{A}_{ik}\cdot B_{kj}$$

对于维度较大的矩阵（如当 $m$、$n$、$p$ 很大时），这种计算的时间复杂度为 $O(mnp)$，在 $m=n=p$ 的情况下即为 $O(n^3)$。这种复杂度在大规模计算中往往难以接受。

近似矩阵乘法的核心思想是找到一种方法，使得计算结果 $\tilde{C}$ 与真实结果 $C$ 足够接近，同时大幅降低计算复杂度。形式上，我们希望：

$$\Vert C-\tilde{C}\Vert \le ϵ\Vert A\Vert \Vert B\Vert$$

其中 $ϵ$ 是一个小的误差容忍度，$\Vert \cdot \Vert$ 表示某种矩阵范数。

从数学上看，近似矩阵乘法可以视为一个优化问题：

$$\underset{\tilde{C}}{\min }\mathcal{C}(\tilde{C})\text{s.t.}\Vert C-\tilde{C}\Vert \le ϵ\Vert A\Vert \Vert B\Vert$$

其中 $\mathcal{C}(\tilde{C})$ 是计算 $\tilde{C}$ 的复杂度函数。

低秩近似方法

奇异值分解（SVD）

奇异值分解是近似矩阵乘法中最基础的方法之一。给定矩阵 $A\in {\mathbb{R}}^{m\times n}$，其SVD分解为：

$$A=U\Sigma V^T=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^T$$

其中 $U=[u_1,u_2,\dots ,u_m]\in {\mathbb{R}}^{m\times m}$ 和 $V=[v_1,v_2,\dots ,v_n]\in {\mathbb{R}}^{n\times n}$ 是正交矩阵，$\Sigma \in {\mathbb{R}}^{m\times n}$ 是对角矩阵，对角线上的元素 ${\sigma }_1\ge {\sigma }_2\ge ...\ge {\sigma }_r>0$ 称为奇异值，$r=\text{rank}(A)\le \min (m,n)$。

低秩近似的思想是只保留前 $k$ 个最大的奇异值及其对应的奇异向量，得到 $A$ 的秩为 $k$ 的最佳近似：

$$A_k=\sum _{i=1}^k{\sigma }_i{u}_i{v}_i^T=U_k{\Sigma }_k{V}_k^T$$

其中 $U_k\in {\mathbb{R}}^{m\times k}$，${\Sigma }_k\in {\mathbb{R}}^{k\times k}$，$V_k\in {\mathbb{R}}^{n\times k}$。

Eckart-Young定理表明，$A_k$ 是所有秩不超过 $k$ 的矩阵中，在Frobenius范数和谱范数下最接近 $A$ 的矩阵：

$$\Vert A-A_k{\Vert }_F=\underset{\text{rank}(X)\le k}{\min }\Vert A-X{\Vert }_F=\sqrt{\sum _{i=k+1}^r{\sigma }_i^2}$$

$$\Vert A-A_k{\Vert }_2=\underset{\text{rank}(X)\le k}{\min }\Vert A-X{\Vert }_2={\sigma }_{k+1}$$

类似地，我们可以对矩阵 $B$ 进行低秩近似：$B_l=P_l{\Lambda }_l{Q}_l^T$。然后，矩阵乘积可以近似为：

$$C\approx \tilde{C}=A_k{B}_l=(U_k{\Sigma }_k{V}_k^T)(P_l{\Lambda }_l{Q}_l^T)$$

这种计算的复杂度为 $O(mk+kl+lp)$，当 $k$ 和 $l$ 远小于原始维度时，可以显著降低计算量。

近似的误差分析可以通过矩阵范数的性质得到：

$$\Vert AB-A_k{B}_l{\Vert }_F\le \Vert A-A_k{\Vert }_F\Vert B{\Vert }_F+\Vert A_k{\Vert }_F\Vert B-B_l{\Vert }_F$$

$$\le \sqrt{\sum _{i=k+1}^{r_A}{\sigma }_i^2(A)}\Vert B{\Vert }_F+\sqrt{\sum _{i=1}^k{\sigma }_i^2(A)}\sqrt{\sum _{j=l+1}^{r_B}{\sigma }_j^2(B)}$$

其中 $r_A$ 和 $r_B$ 分别是 $A$ 和 $B$ 的秩。

随机化SVD

尽管SVD提供了最优的低秩近似，但计算完整的SVD本身就是一个复杂的过程。随机化SVD通过引入随机投影来加速这一过程。基本步骤如下：

1. 生成随机矩阵 $\Omega \in {\mathbb{R}}^{n\times (k+p)}$，其中 $p$ 是一个小的过采样参数
2. 计算 $Y=A\Omega$
3. 对 $Y$ 进行QR分解，得到 $Y=QR$
4. 计算 $B=Q^{T}A$
5. 对小矩阵 $B$ 进行SVD：$B=\tilde{U}\Sigma V^T$
6. 计算 $U=Q\tilde{U}$

随机化SVD的理论基础是Johnson-Lindenstrauss引理和随机投影的性质。可以证明，当 $\Omega$ 的元素是独立同分布的高斯随机变量时，存在常数 $C$ 和 $c$，使得以概率至少 $1-\delta$：

$$\Vert A-Q{Q}^{T}A{\Vert }_F\le (1+C\sqrt{(k+p)/p})\Vert A-A_k{\Vert }_F+\frac{e\sqrt{k+p}}{\sqrt{p}}\Vert A{\Vert }_F\delta$$

当过采样参数 $p$ 适当选择时，这个误差界接近于最优低秩近似的误差界。

进一步地，可以证明随机化SVD得到的近似 ${\hat{A}}_k=U\Sigma V^T$ 满足：

$$\mathbb{E}[\Vert A-{\hat{A}}_k{\Vert }_F^2]\le (1+\frac{k}{p-1})\Vert A-A_k{\Vert }_F^2$$

这表明随机化SVD在期望意义上接近于最优低秩近似。

随机投影方法

Johnson-Lindenstrauss变换

Johnson-Lindenstrauss引理告诉我们，包含 $n$ 个点的欧几里得空间中的点集 $S$，可以嵌入到 $d=O(\log n/{ϵ}^2)$ 维空间，同时保持点之间距离的相对误差在 $(1\pm ϵ)$ 之内。

基于这一原理，我们可以设计随机投影矩阵 $R\in {\mathbb{R}}^{n\times d}$ ($d\ll n$)，并将矩阵乘法近似为：

$$C\approx \tilde{C}=(AR)(BR{)}^T$$

其中 $AR\in {\mathbb{R}}^{m\times d}$，$BR\in {\mathbb{R}}^{p\times d}$。这样计算的复杂度降为 $O(mnd+npd+mpd)$，当 $d$ 远小于 $n$ 时，效率提升显著。

常用的随机投影矩阵包括：

1. 高斯随机投影：$R_{ij}\sim N(0,1/d)$
2. 稀疏随机投影：$R_{ij}=\{\begin{array}{ll}+\sqrt{s/d}& \text{概率为}1/2s\\ 0& \text{概率为}1-1/s\\ -\sqrt{s/d}& \text{概率为}1/2s\end{array}$

其中 $s$ 是稀疏度参数。

对于高斯随机投影，可以证明对于任意固定向量 $x,y\in {\mathbb{R}}^n$，当 $d=O(\log (1/\delta )/{ϵ}^2)$ 时，以概率至少 $1-\delta$：

$$(1-ϵ)\Vert x{\Vert }_2\Vert y{\Vert }_2\le d\cdot (Rx{)}^T(Ry)\le (1+ϵ)\Vert x{\Vert }_2\Vert y{\Vert }_2$$

这意味着随机投影近似保留了向量内积，从而保留了矩阵乘法的结构。

更精确地，对于矩阵乘法 $C=AB$，可以证明：

$$\mathbb{E}[\Vert (AR)(BR{)}^T-AB{\Vert }_F^2]=\frac{\Vert A{\Vert }_F^2\Vert B{\Vert }_F^2}{d}$$

这表明当 $d$ 增大时，近似的误差按 $1/d$ 的速率减小。

伪随机变换

完全随机的矩阵仍需生成和存储大量随机数，为了进一步提高效率，可以使用结构化随机矩阵，如：

$$R=\frac{1}{\sqrt{d}}DHS$$

其中 $D$ 是随机对角矩阵，对角元素为独立的 $\pm 1$ 随机变量；$H$ 是Hadamard矩阵，定义为 $H_n=\left[\begin{array}{cc}{H}_{n/2}& H_{n/2}\\ H_{n/2}& -H_{n/2}\end{array}\right]$，$H_1=[1]$；$S$ 是随机抽样矩阵，从单位矩阵 $I_n$ 中随机选择 $d$ 列组成。使用这种结构化随机矩阵，矩阵-向量乘法 $Rx$ 的复杂度可以从 $O(nd)$ 降低到 $O(n\log d)$，进一步提高了效率。

Ailon和Chazelle证明了这种伪随机变换仍然满足Johnson-Lindenstrauss引理，即对于任意固定向量 $x,y\in {\mathbb{R}}^n$，当 $d=O(\log (1/\delta )/{ϵ}^2)$ 时，以概率至少 $1-\delta$：

$$(1-ϵ)\Vert x{\Vert }_2\Vert y{\Vert }_2\le d\cdot (Rx{)}^T(Ry)\le (1+ϵ)\Vert x{\Vert }_2\Vert y{\Vert }_2$$

蒙特卡洛采样方法

基本蒙特卡洛方法

蒙特卡洛方法的核心思想是通过随机采样来近似矩阵乘法。考虑矩阵乘积 $C=AB$，我们可以将其表示为：

$$C_{ij}=\sum _{k=1}^n{A}_{ik}{B}_{kj}=n\cdot {\mathbb{E}}_K[A_{iK}{B}_{Kj}]$$

其中 $K$ 是在 { 1 , 2 , . . . , n } \{1, 2, ..., n\} {1,2,...,n} 上均匀分布的随机变量。通过对 $K$ 进行多次采样，我们可以得到 $C_{ij}$ 的无偏估计：

$${\tilde{C}}_{ij}=\frac{n}{s}\sum _{t=1}^s{A}_{i{K}_t}{B}_{K_{t}j}$$

其中 $s$ 是采样次数，$K_1,K_2,...,K_s$ 是独立同分布的随机下标。当 $s\ll n$ 时，计算复杂度可以从 $O(mnp)$ 降低到 $O(mps)$。

该估计的方差为：

$$\text{Var}({\tilde{C}}_{ij})=\frac{n^2}{s}\cdot \text{Var}(A_{iK}{B}_{Kj})=\frac{n^2}{s}\left(\frac{1}{n}\sum _{k=1}^n(A_{ik}{B}_{kj}{)}^2-{\left(\frac{1}{n}\sum _{k=1}^n{A}_{ik}{B}_{kj}\right)}^2\right)$$

根据Hoeffding不等式，对于任意固定的 $i,j$，当 $s\ge \frac{2n^2\log (2/\delta )}{{ϵ}^2\Vert A_{i:}{\Vert }_2^2\Vert B_{:j}{\Vert }_2^2}$ 时，以概率至少 $1-\delta$：

$$|{\tilde{C}}_{ij}-C_{ij}|\le ϵ\Vert A_{i:}{\Vert }_2\Vert B_{:j}{\Vert }_2$$

其中 $A_{i:}$ 是 $A$ 的第 $i$ 行，$B_{:j}$ 是 $B$ 的第 $j$ 列。

重要性采样

基本蒙特卡洛方法的问题在于方差可能较大，导致估计不稳定。重要性采样通过引入非均匀采样概率来降低方差。具体来说，我们选择概率分布：

$$p(k)=\frac{\Vert A_{:,k}{\Vert }_2\Vert B_{k,:}{\Vert }_2}{{\sum }_{l=1}^n\Vert A_{:,l}{\Vert }_2\Vert B_{l,:}{\Vert }_2}$$

其中 $A_{:,k}$ 表示 $A$ 的第 $k$ 列，$B_{k,:}$ 表示 $B$ 的第 $k$ 行。相应的估计为：

$${\tilde{C}}_{ij}=\sum _{t=1}^s\frac{1}{s\cdot p(K_t)}{A}_{i{K}_t}{B}_{K_{t}j}$$

可以证明这个估计是无偏的，即 $\mathbb{E}[{\tilde{C}}_{ij}]=C_{ij}$。更重要的是，它的方差比均匀采样小：

$$\text{Var}({\tilde{C}}_{ij})=\frac{1}{s}\sum _{k=1}^n\frac{(A_{ik}{B}_{kj}{)}^2}{p(k)}-\frac{1}{s}{C}_{ij}^2$$

通过适当选择概率分布 $p(k)$，可以显著降低估计的方差。理论上最优的采样分布是：

$$p^{\ast }(k)=\frac{|A_{ik}{B}_{kj}|}{{\sum }_{l=1}^n|A_{il}{B}_{lj}|}$$

然而，这个分布取决于特定的 $i,j$，难以在实践中直接应用。因此，我们通常使用上述基于列范数和行范数的分布，它在整体上提供了良好的方差减小效果。

Drineas等人证明了，当使用重要性采样时，对于 $\delta \in (0,1)$ 和 $ϵ>0$，如果采样数

$$s\ge \frac{68\Vert A{\Vert }_F^2\Vert B{\Vert }_F^2}{{ϵ}^2\Vert AB{\Vert }_F^2}\ln \left(\frac{2}{\delta }\right)$$

则以概率至少 $1-\delta$：

$$\Vert AB-\tilde{C}{\Vert }_F^2\le {ϵ}^2\Vert AB{\Vert }_F^2$$

张量分解方法

CP分解

对于高维张量的计算，可以通过张量分解方法实现近似。CP（CANDECOMP/PARAFAC）分解将张量表示为多个秩为1的张量之和：

$$\mathcal{T}\approx \sum _{r=1}^R{a}_r\otimes b_r\otimes c_r\otimes \cdots$$

其中 $\otimes$ 表示张量积，$R$ 是分解的秩。将矩阵视为二阶张量，CP分解等价于矩阵的低秩分解。

形式上，对于三阶张量 $\mathcal{T}\in {\mathbb{R}}^{I\times J\times K}$，其CP分解可以写为：

$${\mathcal{T}}_{ijk}\approx \sum _{r=1}^R{a}_{ir}{b}_{jr}{c}_{kr}$$

其中 $a_r\in {\mathbb{R}}^I$, $b_r\in {\mathbb{R}}^J$, $c_r\in {\mathbb{R}}^K$。

CP分解可以通过交替最小二乘（ALS）方法求解。给定一个张量 $\mathcal{T}$，CP-ALS算法迭代地更新因子矩阵 $A\in {\mathbb{R}}^{I\times R}$，$B\in {\mathbb{R}}^{J\times R}$，$C\in {\mathbb{R}}^{K\times R}$，使得

$$\underset{A,B,C}{\min }\Vert \mathcal{T}-\sum _{r=1}^R{a}_r\otimes b_r\otimes c_r{\Vert }_F^2$$

在每一步中，固定两个因子矩阵，更新第三个矩阵。例如，固定 $B$ 和 $C$，更新 $A$ 的问题变为：

$$\underset{A}{\min }\Vert {\mathcal{T}}_{(1)}-A(C\odot B{)}^T{\Vert }_F^2$$

其中 ${\mathcal{T}}_{(1)}\in {\mathbb{R}}^{I\times JK}$ 是 $\mathcal{T}$ 的模式-1展开，$\odot$ 表示Khatri-Rao积（列到列的Kronecker积）。该问题的解为：

$$A={\mathcal{T}}_{(1)}[(C\odot B{)}^T{]}^{†}={\mathcal{T}}_{(1)}(C\odot B)(C^{T}C\ast B^{T}B{)}^{-1}$$

其中 $\ast$ 表示Hadamard积（元素对应的乘积）。

Tucker分解

Tucker分解是另一种常用的张量分解方法：

$$\mathcal{T}\approx \mathcal{G}{\times }_{1}A{\times }_{2}B{\times }_{3}C{\times }_4\cdots$$

其中 $\mathcal{G}$ 是核心张量，$A$, $B$, $C$, … 是因子矩阵，${\times }_i$ 表示沿第 $i$ 维度的张量-矩阵乘积。对于三阶张量，Tucker分解可以写为：

$${\mathcal{T}}_{ijk}\approx \sum _{p=1}^P\sum _{q=1}^Q\sum _{r=1}^R{\mathcal{G}}_{pqr}{A}_{ip}{B}_{jq}{C}_{kr}$$

其中 $\mathcal{G}\in {\mathbb{R}}^{P\times Q\times R}$ 是核心张量，$A\in {\mathbb{R}}^{I\times P}$, $B\in {\mathbb{R}}^{J\times Q}$, $C\in {\mathbb{R}}^{K\times R}$ 是因子矩阵。

Tucker分解可以通过高阶奇异值分解（HOSVD）或交替最小二乘（ALS）方法求解。HOSVD方法的步骤如下：

1. 计算 $\mathcal{T}$ 的模式-1、模式-2、模式-3展开的SVD：
   • ${\mathcal{T}}_{(1)}\approx U_A{\Sigma }_A{V}_A^T$
   • ${\mathcal{T}}_{(2)}\approx U_B{\Sigma }_B{V}_B^T$
   • ${\mathcal{T}}_{(3)}\approx U_C{\Sigma }_C{V}_C^T$
2. 取 $A=U_A(:,1:P)$, $B=U_B(:,1:Q)$, $C=U_C(:,1:R)$
3. 计算核心张量 $\mathcal{G}=\mathcal{T}{\times }_1{A}^T{\times }_2{B}^T{\times }_3{C}^T$

在矩阵乘法的上下文中，Tucker分解可以用于近似高维数据的运算，将原本的大规模矩阵乘法转化为多个小规模矩阵乘法和张量收缩操作。

快速矩阵乘法算法

Strassen算法

Strassen算法是精确矩阵乘法的优化算法，其计算复杂度为 $O(n^{{\log }_{2}7})\approx O(n^{2.81})$，优于传统的 $O(n^3)$ 复杂度。该算法的基本思想是将 $n\times n$ 矩阵分解为四个 $n/2\times n/2$ 的子矩阵，然后通过7次子矩阵乘法（而非直接的8次）和多次加减法得到结果。

具体来说，假设我们有两个 $n\times n$ 矩阵 $A$ 和 $B$，将它们分块为：

$$A=\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right],B=\left[\begin{array}{cc}{B}_{11}& B_{12}\\ B_{21}& B_{22}\end{array}\right]$$

传统方法需要8次子矩阵乘法：

$$C_{11}=A_{11}{B}_{11}+A_{12}{B}_{21}$$
$$C_{12}=A_{11}{B}_{12}+A_{12}{B}_{22}$$
$$C_{21}=A_{21}{B}_{11}+A_{22}{B}_{21}$$
$$C_{22}=A_{21}{B}_{12}+A_{22}{B}_{22}$$

Strassen算法引入以下7个中间矩阵：

$$M_1=(A_{11}+A_{22})(B_{11}+B_{22})$$
$$M_2=(A_{21}+A_{22})B_{11}$$
$$M_3=A_{11}(B_{12}-B_{22})$$
$$M_4=A_{22}(B_{21}-B_{11})$$
$$M_5=(A_{11}+A_{12})B_{22}$$
$$M_6=(A_{21}-A_{11})(B_{11}+B_{12})$$
$$M_7=(A_{12}-A_{22})(B_{21}+B_{22})$$

然后，结果矩阵 $C=AB$ 的子块为：

$$C_{11}=M_1+M_4-M_5+M_7$$
$$C_{12}=M_3+M_5$$
$$C_{21}=M_2+M_4$$
$$C_{22}=M_1-M_2+M_3+M_6$$

可以验证这两种计算方法给出相同的结果，但Strassen算法只需要7次子矩阵乘法，而不是8次。通过递归应用这一技术，可以将 $n\times n$ 矩阵乘法的复杂度从 $O(n^3)$ 降低到 $O(n^{{\log }_{2}7})\approx O(n^{2.81})$。

需要注意的是，Strassen算法引入的额外加减法操作会增加数值误差，因此在某些情况下可能需要权衡精度和效率。

Coppersmith-Winograd算法

Coppersmith-Winograd算法及其改进版本是理论上最快的矩阵乘法算法，其计算复杂度可达 $O(n^{2.373})$。该算法是对Strassen方法的进一步泛化，通过更复杂的代数结构来减少递归中的乘法次数。Coppersmith-Winograd算法的核心是一个基于三线性形式的结构。给定三个向量 $x,y,z$，我们考虑所有可能的三重积 $x_i{y}_j{z}_k$。传统方法需要计算所有 $n^3$ 个三重积，而Coppersmith-Winograd算法通过巧妙的代数变换，将其转化为更少的乘法操作。

具体来说，算法基于以下恒等式：

$$\left(\sum _{i=0}^{n-1}{x}_i\right)\left(\sum _{j=0}^{n-1}{y}_j\right)\left(\sum _{k=0}^{n-1}{z}_k\right)=\sum _{i,j,k=0}^{n-1}{x}_i{y}_j{z}_k$$

通过分组和重新组合，可以用少于 $n^3$ 的乘法来计算所有 $x_i{y}_j{z}_k$ 的值。这一思想可以推广到矩阵乘法，从而实现复杂度的降低。Coppersmith-Winograd算法的理论复杂度为 $O(n^{\omega })$，其中 $\omega \approx 2.373$ 是当前已知的最小指数。然而，由于算法常数项较大且实现复杂，它在实践中较少使用。

量化方法

低精度计算

量化是一种通过降低数值精度来加速计算的技术。例如，将32位浮点数转换为8位定点数，可以大幅减少存储空间和计算开销。在矩阵乘法中，我们可以先将矩阵 $A$ 和 $B$ 量化为低精度表示，计算乘积后再转回高精度：

$$\tilde{C}=\text{dequantize}(\text{quantize}(A)\times \text{quantize}(B))$$

形式上，量化可以表示为：

$$\text{quantize}(x)=\text{round}\left(\frac{x-\min (x)}{\max (x)-\min (x)}\cdot (2^b-1)\right)$$

其中 $b$ 是目标位宽。相应的反量化操作为：

$$\text{dequantize}(y)=\min (x)+\frac{y}{2^b-1}\cdot (\max (x)-\min (x))$$

量化误差的理论分析可以基于量化噪声模型。假设量化误差是独立同分布的均匀随机变量，那么量化误差的方差为：

$${\sigma }_q^2=\frac{(\max (x)-\min (x){)}^2}{(2^b-1{)}^2}\cdot \frac{1}{12}$$

对于量化后的矩阵乘法，误差可以建模为：

$$E_{AB}=(A+\Delta A)(B+\Delta B)-AB=\Delta A\cdot B+A\cdot \Delta B+\Delta A\cdot \Delta B$$

其中 $\Delta A$ 和 $\Delta B$ 是量化误差矩阵。

混合精度计算

混合精度计算是量化的一种变体，它在计算过程中使用不同的精度级别。例如，可以用低精度进行主要计算，但用高精度累加结果，以减少舍入误差：

$${\tilde{C}}_{ij}=\sum _{k=1}^n\text{float32}(\text{float16}(A_{ik})\cdot \text{float16}(B_{kj}))$$

混合精度计算的理论基础是数值分析中的舍入误差累积理论。给定舍入误差界 ${ϵ}_{16}$ 和 ${ϵ}_{32}$（分别对应16位和32位浮点数），混合精度矩阵乘法的相对误差界为：

$$\frac{\Vert \tilde{C}-C{\Vert }_F}{\Vert C{\Vert }_F}\le (n+2){ϵ}_{16}+O(n{ϵ}_{16}^2)$$

相比之下，纯16位计算的误差界为：

$$\frac{\Vert \tilde{C}-C{\Vert }_F}{\Vert C{\Vert }_F}\le 2n{ϵ}_{16}+O(n^2{ϵ}_{16}^2)$$

可以看出，混合精度方法有效地减小了高阶误差项，特别是当 $n$ 较大时。

近似矩阵乘法的误差分析

确定性误差界

对于基于低秩近似的方法，可以通过奇异值给出精确的误差界。假设 $A_k$ 和 $B_l$ 分别是 $A$ 和 $B$ 的秩为 $k$ 和 $l$ 的最佳近似，则：

$$\Vert A-A_k{\Vert }_F=\sqrt{\sum _{i=k+1}^{\min (m,n)}{\sigma }_i^2(A)}$$

$$\Vert B-B_l{\Vert }_F=\sqrt{\sum _{i=l+1}^{\min (n,p)}{\sigma }_i^2(B)}$$

矩阵乘积的误差可以通过三角不等式得到：

$$\Vert AB-A_k{B}_l{\Vert }_F\le \Vert AB-A_{k}B{\Vert }_F+\Vert A_{k}B-A_k{B}_l{\Vert }_F\le \Vert A-A_k{\Vert }_F\cdot \Vert B{\Vert }_F+\Vert A_k{\Vert }_F\cdot \Vert B-B_l{\Vert }_F$$

利用奇异值分解的性质，$\Vert A_k{\Vert }_F=\sqrt{{\sum }_{i=1}^k{\sigma }_i^2(A)}$，可以进一步得到：

$$\Vert AB-A_k{B}_l{\Vert }_F\le \sqrt{\sum _{i=k+1}^{\min (m,n)}{\sigma }_i^2(A)}\cdot \Vert B{\Vert }_F+\sqrt{\sum _{i=1}^k{\sigma }_i^2(A)}\cdot \sqrt{\sum _{i=l+1}^{\min (n,p)}{\sigma }_i^2(B)}$$

这个界限是确定性的，不依赖于任何概率假设。

概率误差界

对于基于随机方法的近似，通常给出概率误差界。例如，对于随机投影方法，假设 $R\in {\mathbb{R}}^{n\times d}$ 是高斯随机矩阵，元素服从 $N(0,1/d)$ 分布，则对于任意 $ϵ,\delta >0$，当 $d\ge \frac{c\log (1/\delta )}{{ϵ}^2}$ 时，以概率至少 $1-\delta$：

$$\Vert AB-(AR)(BR{)}^T{\Vert }_F\le ϵ\Vert A{\Vert }_F\Vert B{\Vert }_F$$

其中 $c$ 是一个适当的常数。

这个界限可以通过矩阵集中不等式来证明。首先，注意到 $(AR)(BR{)}^T$ 可以写为：

$$(AR)(BR{)}^T=\sum _{j=1}^d(AR{)}_{:,j}(BR{)}_{:,j}^T$$

这是 $d$ 个独立随机矩阵的和。定义 $X_j=(AR{)}_{:,j}(BR{)}_{:,j}^T$，则 $\mathbb{E}[X_j]=\frac{1}{d}AB$，且 $(AR)(BR{)}^T={\sum }_{j=1}^d{X}_j$。

应用矩阵Bernstein不等式，当 $d$ 足够大时，可以得到上述概率界限。

类似地，对于蒙特卡洛方法，当采样数 $s$ 足够大时，可以证明：

$$P\left(\Vert AB-\tilde{C}{\Vert }_F\le ϵ\Vert A{\Vert }_F\Vert B{\Vert }_F\right)\ge 1-\delta$$

其中 $s\ge c\cdot \log (mp/\delta )/{ϵ}^2$，$c$ 是一个适当的常数。

近似矩阵乘法的实例

为了更具体地说明近似矩阵乘法的工作原理，我们通过一个简单的例子来演示。假设我们有两个矩阵：

$$A=\left[\begin{array}{cccc}2& 1& 3& 0\\ 1& 2& 0& 1\\ 0& 1& 1& 2\end{array}\right],B=\left[\begin{array}{ccc}1& 0& 2\\ 2& 1& 0\\ 0& 3& 1\\ 1& 0& 1\end{array}\right]$$

精确的矩阵乘积为：

$$C=AB=\left[\begin{array}{ccc}5& 10& 7\\ 6& 2& 3\\ 5& 1& 3\end{array}\right]$$

现在，我们应用随机采样方法，假设采样次数 $s=2$，随机选择的下标为 $K_1=1$ 和 $K_2=3$，则近似矩阵乘积为：

$${\tilde{C}}_{ij}=\frac{4}{2}\left(A_{i1}{B}_{1j}+A_{i3}{B}_{3j}\right)$$

计算得到：

$$\tilde{C}=\left[\begin{array}{ccc}2\cdot (2\cdot 1+3\cdot 0)& 2\cdot (2\cdot 0+3\cdot 3)& 2\cdot (2\cdot 2+3\cdot 1)\\ 2\cdot (1\cdot 1+0\cdot 0)& 2\cdot (1\cdot 0+0\cdot 3)& 2\cdot (1\cdot 2+0\cdot 1)\\ 2\cdot (0\cdot 1+1\cdot 0)& 2\cdot (0\cdot 0+1\cdot 3)& 2\cdot (0\cdot 2+1\cdot 1)\end{array}\right]=\left[\begin{array}{ccc}4& 18& 10\\ 2& 0& 4\\ 0& 6& 2\end{array}\right]$$

计算误差为：

$$\Vert C-\tilde{C}{\Vert }_F=\sqrt{\sum _{i,j}(C_{ij}-{\tilde{C}}_{ij}{)}^2}=\sqrt{(5-4{)}^2+(10-18{)}^2+\cdots +(3-2{)}^2}\approx 9.85$$

相对误差为：

$$\frac{\Vert C-\tilde{C}{\Vert }_F}{\Vert C{\Vert }_F}=\frac{9.85}{\sqrt{5^2+10^2+\cdots +3^2}}\approx 0.68$$

这个误差较大，主要是因为采样次数太少。在实际应用中，我们通常会使用更多的采样点或结合其他技术来提高近似精度。