草图绘制有效吗？

本文转载自我的个人网站。感谢lunarflu邀请我在此分享！

---

我很高兴地宣布，我的论文《线性最小二乘问题的快速且前向稳定的随机算法》已作为预印本发布在 arXiv 上。

随着这篇论文的发布，现在似乎是讨论一个我一直想写的主题的好时机：草图绘制。在过去的二十年里，草图绘制已成为矩阵计算中广泛使用的算法工具。尽管历史悠久，但关于草图绘制是否真的有效的问题似乎仍然存在。

草图绘制在实践中真的有效吗？我想知道为什么它没有真正在深度学习中使用。

— typedfemale (@typedfemale) 2022年9月22日

在这篇文章中，我想批判性地审视“草图绘制有效吗？”这个问题。回答这个问题需要回答两个基本问题：

1. 什么是草图绘制？
2. 草图绘制的“有效”意味着什么？

我认为关于草图绘制有效性的很大一部分分歧归结于对这些问题的不同回答。通过考虑这些问题的不同可能答案，我希望能够对草图绘制作为解决线性代数问题的算法原语的效用提供一个平衡的视角。

草图绘制

在矩阵计算中，草图绘制实际上是（线性）降维的同义词。假设我们正在解决一个涉及一个或多个高维向量 $b\in {\mathbb{R}}^{n}b\; \backslash in\; \backslash mathbb\{R\}^n$ 或许一个高矩阵 $A\in {\mathbb{R}}^{n\times k}A\backslash in\; \backslash mathbb\{R\}^\{n\backslash times\; k\}$ 的问题。草图矩阵是一个 $d\times nd\backslash times\; n$ 矩阵 $S\in {\mathbb{R}}^{d\times n}S\; \backslash in\; \backslash mathbb\{R\}^\{d\backslash times\; n\}$，其中 $d\ll nd\; \backslash ll\; n$。当与高维向量 $bb$ 或高矩阵 $AA$ 相乘时，草图矩阵 $SS$ 产生比原始向量 $bb$ 和矩阵 $AA$ 小得多的压缩或“草图化”版本 $SbSb$ 和 $SASA$。

设 $\mathsf{E}=\{x_1,\dots ,x_p\}\backslash mathsf\{E\}=\backslash \{x\_1,\backslash ldots,x\_p\backslash \}$ 为向量集合。为了使 $SS$ 成为 $\mathsf{E}\backslash mathsf\{E\}$ 的“好”草图矩阵，我们要求 $SS$ 能够将 $\mathsf{E}\backslash mathsf\{E\}$ 中每个向量的长度保持在失真参数 $\epsilon >0\backslash varepsilon>0$ 范围内：

$$\begin{array}{cccc}& (1-\epsilon )\mathrm{\parallel }x\mathrm{\parallel }\le \mathrm{\parallel }Sx\mathrm{\parallel }\le (1+\epsilon )\mathrm{\parallel }x\mathrm{\parallel }& & \text{(1)}\end{array}(1-\backslash varepsilon)\; \backslash |x\backslash |\backslash le\backslash |Sx\backslash |\backslash le(1+\backslash varepsilon)\backslash |x\backslash |\; \backslash tag\{1\}$$

对于 $\mathsf{E}\backslash mathsf\{E\}$ 中的每一个 $xx$。

对于线性代数问题，我们通常希望对矩阵 $AA$ 进行草图绘制。在这种情况下，我们希望草图“良好”的适当集合 $\mathsf{E}\backslash mathsf\{E\}$ 是矩阵 $AA$ 的列空间，定义为：

$$\mathrm{col}(A)\{Ax:x\in {\mathbb{R}}^k\}\mathrm{.}\backslash operatorname\{col\}(A)\; \backslash coloneqq\; \backslash \{\; Ax\; :\; x\; \backslash in\; \backslash mathbb\{R\}^k\; \backslash \}.$$

值得注意的是，存在许多草图矩阵能够为 $\mathsf{E}=\mathrm{col}(A)\backslash mathsf\{E\}=\backslash operatorname\{col\}(A)$ 实现失真 $\epsilon \backslash varepsilon$，且输出维度大约为 $d\approx k\mathrm{/}{\epsilon }^{2}d\; \backslash approx\; k\; /\; \backslash varepsilon^2$。特别是，草图维度 $dd$ 与 $AA$ 的列数 $kk$ 成正比。这非常棒！我们可以设计一个单一的草图矩阵 $SS$，它能保留 $AA$ 的列空间中所有无限多个向量 $AxAx$ 的长度。

草图矩阵

草图矩阵有许多类型，每种都有不同的优缺点。许多草图矩阵基于随机化构造，即 $SS$ 的条目被选择为随机数。广义上，草图矩阵可以分为两种类型：

• 数据相关草图。草图矩阵 $SS$ 构造用于特定输入向量集 $\mathsf{E}\backslash mathsf\{E\}$。
• 无知草图。草图矩阵 $SS$ 旨在适用于给定大小（即 $\mathsf{E}\backslash mathsf\{E\}$ 具有 $pp$ 个元素）或维度（即 $\mathsf{E}\backslash mathsf\{E\}$ 是一个 $kk$ 维线性子空间）的任意输入向量集 $\mathsf{E}\backslash mathsf\{E\}$。

本篇文章将只讨论无知草图。我们将研究三种类型的草图矩阵：高斯嵌入、欠采样随机三角变换和稀疏符号嵌入。

这些草图矩阵的构建细节及其优缺点可能有些技术性。这三种构造都独立于本文的其余部分，可以在初次阅读时跳过。主要观点是，好的草图矩阵是存在的，并且应用速度很快：将 $b\in {\mathbb{R}}^{n}b\backslash in\backslash mathbb\{R\}^n$ 缩减为 $Sb\in {\mathbb{R}}^{d}Sb\backslash in\backslash mathbb\{R\}^\{d\}$ 大致需要 $\mathcal{O}(n\log n)\backslash mathcal\{O\}(n\backslash log\; n)$ 次操作，而不是我们预期将 $d\times nd\backslash times\; n$ 矩阵乘以长度为 $nn$ 的向量所需的 $\mathcal{O}(dn)\backslash mathcal\{O\}(dn)$ 次操作。这里，$\mathcal{O}(\cdot )\backslash mathcal\{O\}(\backslash cdot)$ 是大 O 符号。

高斯嵌入

最简单的草图矩阵 $S\in {\mathbb{R}}^{d\times n}S\backslash in\backslash mathbb\{R\}^\{d\backslash times\; n\}$ 是通过（独立地）将 $SS$ 的每个条目设置为零均值和 $1\mathrm{/}d1/d$ 方差的高斯随机数来获得的。这样的草图矩阵被称为高斯嵌入。这里的嵌入是草图矩阵的同义词。

优点。高斯嵌入易于编写代码，只需一个标准的矩阵乘法即可应用于向量 $SbSb$ 或矩阵 $SASA$。高斯嵌入具有清晰的理论分析，其数学特性已得到充分理解。

缺点。对于高斯嵌入，计算 $SbSb$ 需要 $\mathcal{O}(dn)\backslash mathcal\{O\}(dn)$ 次操作，明显慢于我们将在下面讨论的其他草图矩阵。此外，生成和存储高斯嵌入可能在计算上非常昂贵。

欠采样随机三角变换

欠采样随机三角变换（SRTT）草图矩阵的形式更为复杂。草图矩阵定义为三个矩阵的比例积：

$$S=\sqrt{\frac{n}{d}}\cdot R\cdot F\cdot D\mathrm{.}S\; =\; \backslash sqrt\{\backslash frac\{n\}\{d\}\}\; \backslash cdot\; R\; \backslash cdot\; F\; \backslash cdot\; D.$$

这些矩阵的定义如下：

• $D\in {\mathbb{R}}^{n\times n}D\backslash in\backslash mathbb\{R\}^\{n\backslash times\; n\}$ 是一个对角矩阵，其每个条目都是一个随机的 $\pm 1\backslash pm\; 1$（以等概率独立选择）。
• $F\in {\mathbb{R}}^{n\times n}F\backslash in\backslash mathbb\{R\}^\{n\backslash times\; n\}$ 是一个快速三角变换，例如快速离散余弦变换。也可以使用普通的快速傅里叶变换，但这会产生复数值草图。
• $R\in {\mathbb{R}}^{d\times n}R\backslash in\backslash mathbb\{R\}^\{d\backslash times\; n\}$ 是一个选择矩阵。为了生成 $RR$，设 $i_1,\dots ,i_{d}i\_1,\backslash ldots,i\_d$ 是从 $1,\dots ,n\backslash \backslash \{1,\backslash ldots,n\backslash \backslash \}$ 中随机选取的子集，不重复。 $RR$ 定义为一个矩阵，对于每个向量 $bb$，都有 $Rb=(b_{i_1},\dots ,b_{i_d})Rb\; =\; (b\_\{i\_1\},\backslash ldots,b\_\{i\_d\})$。

要在计算机上存储 $SS$，只需存储 $DD$ 的 $nn$ 个对角线元素以及定义 $RR$ 的 $dd$ 个选定坐标 $i_1,\dots ,i_{d}i\_1,\backslash ldots,i\_d$。将 $SS$ 乘以向量 $bb$ 的操作应按顺序应用矩阵 $RR$、$FF$ 和 $DD$，如下面的 MATLAB 代码所示：

% Generate randomness for S
signs = 2*randi(2,m,1)-3; % diagonal entries of D
idx = randsample(m,d); % indices i_1,...,i_d defining R

% Multiply S against b
c = signs .* b % multiply by D
c = dct(c) % multiply by F
c = c(idx) % multiply by R
c = sqrt(n/d) * c % scale

优点。$SS$ 可以非常快速地应用于向量 $bb$（根据参数选择，为 $\mathcal{O}(n)\backslash mathcal\{O\}(n)$ 或 $\mathcal{O}(n\log k)\backslash mathcal\{O\}(n\backslash log\; k)$ 次操作）。借助良好的稀疏矩阵库，稀疏符号嵌入通常是速度最快的草图矩阵。

缺点。将 $SS$ 应用于向量需要良好的快速三角变换实现。即使使用高质量的三角变换，SRTT 仍可能比稀疏符号嵌入（下文定义）慢得多。例如，请参见本文图 2。SRTT 难以并行化。然而，块 SRTT 更具并行性。理论上，草图维度应选择为 $d\approx (k\log k)\mathrm{/}{\epsilon }^{2}d\; \backslash approx\; (k\backslash log\; k)/\backslash varepsilon^2$，比高斯草图更大。

稀疏符号嵌入

稀疏符号嵌入的形式为：

此处，每列 $s_{i}s\_i$ 都是一个独立生成的随机向量，具有精确的 $\zeta \backslash zeta$ 个非零项，且在均匀随机位置上具有随机 $\pm 1\backslash pm\; 1$ 值。结果是一个 $d\times nd\backslash times\; n$ 矩阵，其中只有 $\zeta \cdot n\backslash zeta\; \backslash cdot\; n$ 个非零项。参数 $\zeta \backslash zeta$ 在实践中通常设置为一个小的常数，例如 $88$。此建议来自这篇论文，我在我自己的工作中成功使用了此参数设置。

优点。通过使用专用稀疏矩阵库，$SS$ 可以非常快速地应用于向量 $bb$（根据参数选择，为 $\mathcal{O}(n)\backslash mathcal\{O\}(n)$ 或 $\mathcal{O}(n\log k)\backslash mathcal\{O\}(n\backslash log\; k)$ 次操作）。借助良好的稀疏矩阵库，稀疏符号嵌入通常是速度最快的草图矩阵，且优势显著。

缺点。为了快速实现，稀疏符号嵌入需要一个好的稀疏矩阵库。它们需要生成和存储大约 $\zeta n\backslash zeta\; n$ 个随机数，这比 SRTT（大约 $nn$ 个数）要高，但远低于高斯嵌入（正好 $dndn$ 个数）。理论上，草图维度应选择为 $d\approx (k\log k)\mathrm{/}{\epsilon }^{2}d\; \backslash approx\; (k\backslash log\; k)/\backslash varepsilon^2$，稀疏度应设置为 $\zeta \approx (\log k)\mathrm{/}\epsilon \backslash zeta\; \backslash approx\; (\backslash log\; k)/\backslash varepsilon$；理论上批准的草图维度（至少根据现有理论）比高斯草图大。在实践中，我们通常可以使用 $d\approx k\mathrm{/}{\epsilon }^{2}d\; \backslash approx\; k/\backslash varepsilon^2$ 和 $\zeta =8\backslash zeta=8$。

摘要

使用 SRTT 或稀疏映射，将长度为 $nn$ 的向量草图化到 $dd$ 维仅需要 $\mathcal{O}(n)\backslash mathcal\{O\}(n)$ 到 $\mathcal{O}(n\log n)\backslash mathcal\{O\}(n\backslash log\; n)$ 次操作。因此，将草图应用于整个 $n\times kn\backslash times\; k$ 矩阵 $AA$ 大约需要 $\mathcal{O}(nk)\backslash mathcal\{O\}(nk)$ 次操作。因此，草图有望加速涉及 $AA$ 的线性代数计算，这些计算通常需要 $\mathcal{O}(n{k}^2)\backslash mathcal\{O\}(nk^2)$ 次操作。

如何使用草图？

使用草图最简单的方法是首先将草图应用于所有数据进行降维，然后使用降维后的数据应用标准算法解决问题。这种使用草图的方法称为 草图-求解。

例如，我们将草图-求解应用于最小二乘问题

$$\begin{array}{cccc}& \underset{x\in {\mathbb{R}}^k}{\mathrm{minimize}}\mathrm{\parallel }Ax-b\mathrm{\parallel }\mathrm{.}& & \text{(2)}\end{array}\backslash operatorname*\{minimize\}\_\{x\backslash in\backslash mathbb\{R\}^k\}\; \backslash |Ax\; -\; b\backslash |.\; \backslash tag\{2\}$$

我们假设这个问题是一个高度超定的问题，其中 $AA$ 的行数 $nn$ 远多于列数 $mm$。

为了用草图-求解解决这个问题，我们为集合 生成一个好的草图矩阵 $SS$。将 $SS$ 应用于我们的数据 $AA$ 和 $bb$ 后，我们得到了一个降维的最小二乘问题

$$\begin{array}{cccc}& \underset{\widehat{x}\in {\mathbb{R}}^k}{\mathrm{minimize}}\mathrm{\parallel }(SA)\widehat{x}-Sb\mathrm{\parallel }\mathrm{.}& & \text{(3)}\end{array}\backslash operatorname*\{minimize\}\_\{\backslash hat\{x\}\backslash in\backslash mathbb\{R\}^k\}\; \backslash |(SA)\backslash hat\{x\}\; -\; Sb\backslash |.\; \backslash tag\{3\}$$

解 $\widehat{x}\backslash hat\{x\}$ 是最小二乘问题的草图-求解解决方案，我们可以将其用作原始最小二乘问题的近似解决方案。

最小二乘法只是草图-求解范式的一个例子。我们还可以使用草图来加速其他算法。例如，我们可以将草图-求解应用于聚类。为了对数据点 $x_1,\dots ,x_{p}x\_1,\backslash ldots,x\_p$ 进行聚类，首先应用草图得到 $S{x}_1,\dots ,S{x}_{p}Sx\_1,\backslash ldots,Sx\_p$，然后对草图化后的数据点应用开箱即用的聚类算法，如k-means。

草图是否有效？

通常，当草图批评者说“草图无效”时，他们的意思是“草图-求解无效”。

为了在一个更具体的背景下解决这个问题，让我们回到最小二乘问题（2）。设 $x_{\star }x\_\backslash star$ 表示最优最小二乘解，设 $\widehat{x}\backslash hat\{x\}$ 为草图-求解解决方案（3）。然后，利用失真条件（1），可以证明

$$\mathrm{\parallel }A\widehat{x}-b\mathrm{\parallel }\le \frac{1+\epsilon }{1-\epsilon }\mathrm{\parallel }Ax-b\mathrm{\parallel }\mathrm{.}\backslash |A\backslash hat\{x\}\; -\; b\backslash |\; \backslash le\; \backslash frac\{1+\backslash varepsilon\}\{1-\backslash varepsilon\}\; \backslash |Ax\; -\; b\backslash |.$$

如果我们使用失真为 $\epsilon =1\mathrm{/}3\backslash varepsilon\; =\; 1/3$ 的草图矩阵，那么这个边界告诉我们

$$\begin{array}{cccc}& \mathrm{\parallel }A\widehat{x}-b\mathrm{\parallel }\le 2\mathrm{\parallel }A{x}_{\star }-b\mathrm{\parallel }\mathrm{.}& & \text{(4)}\end{array}\backslash |A\backslash hat\{x\}\; -\; b\backslash |\; \backslash le\; 2\backslash |Ax\_\backslash star\; -\; b\backslash |.\; \backslash tag\{4\}$$

这是一个好结果还是坏结果？最终，这取决于具体情况。在某些应用中，一个假设的最小二乘解 $\widehat{x}\backslash hat\{x\}$ 的质量可以通过残差范数 $\mathrm{\parallel }A\widehat{x}-b\mathrm{\parallel }\backslash |A\backslash hat\{x\}\; -\; b\backslash |$ 来评估。对于此类应用，边界（4）确保 $\mathrm{\parallel }A\widehat{x}-b\mathrm{\parallel }\backslash |A\backslash hat\{x\}\; -\; b\backslash |$ 最多是 $\mathrm{\parallel }A{x}_{\star }-b\mathrm{\parallel }\backslash |Ax\_\backslash star-b\backslash |$ 的两倍。通常，这意味着 $\widehat{x}\backslash hat\{x\}$ 是最小二乘问题的一个相当不错的近似解。

对于其他问题，适当的精度度量是所谓的前向误差 $\mathrm{\parallel }\widehat{x}-x_{\star }\mathrm{\parallel }\backslash |\backslash hat\{x\}\; -\; x\_\backslash star\backslash |$，衡量 $\widehat{x}\backslash hat\{x\}$ 与 $x_{\star }x\_\backslash star$ 的接近程度。对于这些情况，即使残差具有可比性（4），$\mathrm{\parallel }\widehat{x}-x_{\star }\mathrm{\parallel }\backslash |\backslash hat\{x\}\; -\; x\_\backslash star\backslash |$ 仍然可能很大。

让我们看一个例子，使用我论文中的 MATLAB 代码

[A, b, x, r] = random_ls_problem(1e4, 1e2, 1e8, 1e-4); % Random LS problem
S = sparsesign(4e2, 1e4, 8); % Sparse sign embedding
sketch_and_solve = (S*A) \ (S*b); % Sketch-and-solve
direct = A \ b; % MATLAB mldivide

在这里，我们生成一个大小为 10,000x100 的随机最小二乘问题（条件数为 ${10}^{8}10^8$，残差范数为 ${10}^{-4}10^\{-4\}$）。然后，我们生成一个维度为 $d=400d\; =\; 400$ 的稀疏符号嵌入（对应于大约 $\epsilon \approx 1\mathrm{/}2\backslash varepsilon\; \backslash approx\; 1/2$ 的失真）。然后，我们计算草图-求解解决方案，并作为参考，通过 MATLAB 的 \ 计算一个“直接”解决方案。

我们使用残差和前向误差比较草图-求解解决方案与直接解决方案的质量

fprintf('Residuals: sketch-and-solve %.2e, direct %.2e, optimal %.2e\n',...
           norm(b-A*sketch_and_solve), norm(b-A*direct), norm(r))
fprintf('Forward errors: sketch-and-solve %.2e, direct %.2e\n',...
           norm(x-sketch_and_solve), norm(x-direct))

这是输出

Residuals: sketch-and-solve 1.13e-04, direct 1.00e-04, optimal 1.00e-04
Forward errors: sketch-and-solve 1.06e+03, direct 8.08e-07

草图-求解解决方案的残差范数为 $1.13\times {10}^{-4}$，接近直接方法的残差范数 $1.00\times {10}^{-4}1.00\backslash times\; 10^\{-4\}$。然而，草图-求解的前向误差为 $1\times {10}^{3}1\backslash times\; 10^3$，比直接方法的前向误差 $8\times {10}^{-7}8\backslash times\; 10^\{-7\}$ 大九个数量级。

草图-求解是否有效？最终，这取决于您的应用程序需要什么样的精度。如果只需要足够小的残差，那么草图-求解就完全足够了。如果需要小的前向误差，草图-求解可能会相当糟糕。

改进草图-求解的一种方法是增加草图维度 $dd$ 并降低失真 $\epsilon \backslash varepsilon$。不幸的是，改善草图的失真代价高昂。由于 $d\approx k\mathrm{/}{\epsilon }^{2}d\; \backslash approx\; k\; /\backslash varepsilon^2$ 的关系，将失真降低十倍需要将草图维度 $dd$ 增加一百倍！因此，草图-求解只适用于需要低失真 $\epsilon \backslash varepsilon$ 的情况。

迭代草图：结合草图与迭代

草图-求解是快速获得最小二乘问题低精度解的方法。但这并不是唯一一种将草图用于最小二乘的方法。我们还可以通过将草图与迭代方法相结合来获得高精度解。

最小二乘问题有许多迭代方法。迭代方法生成一系列近似解 $x_1,x_2,\dots x\_1,x\_2,\backslash ldots$，我们希望它们能快速收敛到真正的最小二乘解 $x_{\star }x\_\backslash star$。

为了使用草图迭代地解决最小二乘问题，我们可以利用以下观察结果

如果 $SS$ 是 $\mathsf{E}=\mathrm{col}(A)\backslash mathsf\{E\}\; =\; \backslash operatorname\{col\}(A)$ 的草图矩阵，那么 $(SA{)}^{\mathrm{\top }}SA\approx A^{\mathrm{\top }}A(SA)^\backslash top\; SA\; \backslash approx\; A^\backslash top\; A$。

因此，如果我们计算 QR 分解

$$SA=QR,SA\; =\; QR,$$

那么

$$A^{\mathrm{\top }}A\approx (SA{)}^{\mathrm{\top }}(SA)=R^{\mathrm{\top }}{Q}^{\mathrm{\top }}QR=R^{\mathrm{\top }}R\mathrm{.}A^\backslash top\; A\; \backslash approx\; (SA)^\backslash top\; (SA)\; =\; R^\backslash top\; Q^\backslash top\; Q\; R\; =\; R^\backslash top\; R.$$

请注意，我们使用了 $Q^{\mathrm{\top }}Q=IQ^\backslash top\; Q\; =\; I$ 的事实，因为 $QQ$ 具有正交列。结论是 $R^{\mathrm{\top }}R\approx A^{\mathrm{\top }}AR^\backslash top\; R\; \backslash approx\; A^\backslash top\; A$。

让我们使用近似值 $R^{\mathrm{\top }}R\approx A^{\mathrm{\top }}AR^\backslash top\; R\; \backslash approx\; A^\backslash top\; A$ 来迭代地解决最小二乘问题。从正规方程 [脚注 1] 开始

$$\begin{array}{cccc}& (A^{\mathrm{\top }}A)x=A^{\mathrm{\top }}b\mathrm{.}& & \text{(5)}\end{array}(A^\backslash top\; A)x\; =\; A^\backslash top\; b.\; \backslash tag\{5\}$$

我们可以通过替换 (5) 中的 $A^{\mathrm{\top }}AA^\backslash top\; A$ 为 $R^{\mathrm{\top }}RR^\backslash top\; R$ 并求解，从而得到最小二乘问题的近似解。所得解为

$$x_0=R^{-1}(R^{-\mathrm{\top }}(A^{\mathrm{\top }}b))\mathrm{.}x\_0\; =\; R^\{-1\}\; (R^\{-\backslash top\}(A^\backslash top\; b)).$$

这个解 $x_{0}x\_0$ 通常不是最小二乘问题 (2) 的一个好解，因此我们需要迭代。为此，我们将尝试求解误差 $x-x_{0}x\; -\; x\_0$。为了推导出误差方程，从法方程 (5) 的两边减去 $A^{\mathrm{\top }}A{x}_{0}A^\backslash top\; A\; x\_0$，得到

$$(A^{\mathrm{\top }}A)(x-x_0)=A^{\mathrm{\top }}(b-A{x}_0)\mathrm{.}(A^\backslash top\; A)(x-x\_0)\; =\; A^\backslash top\; (b-Ax\_0).$$

现在，为了求解误差，再次用 $R^{\mathrm{\top }}RR^\backslash top\; R$ 替换 $A^{\mathrm{\top }}AA^\backslash top\; A$，并求解 $xx$，得到新的近似解 $x_{1}x\_1$

$$x\approx x_1{x}_0+R^{-\mathrm{\top }}(R^{-1}(A^{\mathrm{\top }}(b-A{x}_0)))\mathrm{.}x\backslash approx\; x\_1\; \backslash coloneqq\; x\_0\; +\; R^\{-\backslash top\}(R^\{-1\}(A^\backslash top(b-Ax\_0))).$$

我们现在可以再进一步：推导出误差 $x-x_{1}x-x\_1$ 的方程，近似 $A^{\mathrm{\top }}A\approx R^{\mathrm{\top }}RA^\backslash top\; A\; \backslash approx\; R^\backslash top\; R$，并得到一个新的近似解 $x_{2}x\_2$。继续这个过程，我们得到一个迭代式

$$\begin{array}{cccc}& x_{i+1}=x_i+R^{-\mathrm{\top }}(R^{-1}(A^{\mathrm{\top }}(b-A{x}_i)))\mathrm{.}& & \text{(6)}\end{array}x\_\{i+1\}\; =\; x\_i\; +\; R^\{-\backslash top\}(R^\{-1\}(A^\backslash top(b-Ax\_i))).\backslash tag\{6\}$$

这种迭代方法被称为“*迭代速写*”方法。[脚注2]

让我们将迭代速写应用于我们上面考虑的例子。我们将速写即解法和直接方法的正向误差显示为水平虚线，分别为紫色和红色。迭代速写从速写即解法的正向误差大致开始，误差以指数速率下降，直到在十四次迭代过程中达到直接方法的误差。至少对于这个问题，迭代速写为最小二乘问题提供了高精度解！

总而言之，我们现在已经看到了两种截然不同的速写使用方式

• 速写即解法。 速写数据 $AA$ 和 $bb$，并求解速写最小二乘问题 (3)。得到的解 $\widehat{x}\backslash hat\{x\}$ 成本低廉，但精度可能不高。
• 迭代速写。 速写矩阵 $AA$，并得到 $A^{\mathrm{\top }}AA^\backslash top\; A$ 的近似 $R^{\mathrm{\top }}R=(SA{)}^{\mathrm{\top }}(SA)R^\backslash top\; R\; =\; (SA)^\backslash top\; (SA)$。使用近似 $R^{\mathrm{\top }}RR^\backslash top\; R$ 通过迭代 (6) 生成一系列越来越好的最小二乘解 $x_{i}x\_i$。如果迭代次数 $qq$ 足够多，迭代速写解 $x_{q}x\_q$ 的精度可以相当高。

通过将速写与更复杂的迭代方法（例如共轭梯度法和LSQR）相结合，我们可以得到一个收敛速度更快的最小二乘算法，称为*速写预处理*。下面是上面添加了速写预处理的相同图表；我们看到速写预处理的收敛速度甚至比迭代速写更快！

“速写有效吗？”即使对于像最小二乘这样简单的问题，答案也很复杂

直接使用速写（即速写即解法）可以快速、低精度地解决最小二乘问题。但通过将速写与迭代方法（迭代速写和速写预处理）相结合，速写可以实现更高的最小二乘问题精度。

本节我们重点讨论了最小二乘问题，但这些结论可能更普遍适用。如果“速写不起作用”在您的应用程序中，那么如果将其与迭代方法结合使用，它可能就会起作用。

速写能有多精确？

我们上一节关于速写加迭代方法的讨论以积极的基调结束，但还有一个悬而未决的问题有待回答。我们曾说过迭代速写（和速写预处理）以指数速率收敛。但我们的计算机只能存储有限精度的数字；在实践中，迭代方法的精度必须在某个点达到饱和。

如果一个（迭代）最小二乘求解器在运行足够次数的迭代 $qq$ 后，最终精度 $\mathrm{\parallel }{x}_q-x_{\star }\mathrm{\parallel }\backslash |x\_q\; -\; x\_\backslash star\backslash |$ 与标准*直接方法*（例如 MATLAB 的 \ 命令 或 Python 的 *scipy.linalg.lstsq*）的精度相当，则称其为*前向稳定*。前向稳定性不是关于*速度*或*收敛速率*，而是关于*可达到的最大精度*。

速写预处理的稳定性在Meier、Nakatsukasa、Townsend 和 Webb 的一篇近期论文中进行了研究。他们证明，当初始迭代 $x_0=0x\_0\; =\; 0$ 时，速写预处理不是前向稳定的。可达到的最大精度比标准求解器差几个数量级！也许速写毕竟不起作用？

幸运的是，有好消息

• 迭代速写方法被证明是前向稳定的。这个结果在我的新发布论文中有所展示；如果您感兴趣，可以查看一下！
• 如果我们将速写即解法用作速写预处理的*初始迭代 $x_0=\widehat{x}x\_0\; =\; \backslash hat\{x\}$*，那么速写预处理在实践中似乎是前向稳定的。目前尚无支持这一发现的理论分析。对于感兴趣的人来说，迭代速写和速写预处理都不是后向稳定的，后向稳定性是比前向稳定性更强的稳定性保证。幸运的是，前向稳定性对于许多（但不是所有）应用来说，是一种完全充分的稳定性保证。

这些结论相当微妙。要了解发生了什么，查看图表可能会有所帮助：对于另一个具有相同大小的随机生成的最小二乘问题，其条件数${10}^{10}10^\{10\}$，残差${10}^{-6}10^\{-6\}$。

不同方法的性能可以总结如下：速写即解法可能具有非常差的前向误差。零初始化 $x_0=0x\_0\; =\; 0$ 的速写预处理更好，但仍远不如直接方法。迭代速写和 $x_0=\widehat{x}x\_0\; =\; \backslash hat\{x\}$ 的速写预处理表现更好，最终达到与直接方法相当的精度。

更简单地说，如果实现得当，速写终究是有效的！

结论

速写是一种计算工具，就像快速傅里叶变换或随机 SVD 一样。速写可以有效地解决一些问题。但是，像任何计算工具一样，速写也不是万能药。速写允许您对矩阵和向量进行降维，但它会对其进行相当大的失真。这种失真是否可以接受取决于您的问题（需要多少精度？）以及您使用速写的方式（速写即解法或迭代方法）。

---

脚注

1. 正如我在上一篇文章中描述的那样，通常不建议使用法方程来解决最小二乘问题。在这里，我们只是将法方程作为推导最小二乘问题算法的概念工具。
2. “*迭代速写*”这个名字是历史原因。该过程的原始版本涉及在每次迭代中获取新的速写 $S_{i}A=Q_i{R}_{i}S\_iA\; =\; Q\_iR\_i$。后来，人们意识到只需要一次速写 $SASA$ 就足够了，尽管收敛速度会变慢。通常，只需进行一次速写和 QR 分解就值得收敛速度的降低。如果失真足够小，该方法将以指数速率收敛，在几次迭代后产生高精度的最小二乘解。