Iterative — Krylov 方法与随机 SVD¶

概述¶

Krylov 子空间方法（Arnoldi、Lanczos）、随机 SVD 和线性算子接口。TinyMath 依赖链中的第 5 层——最高层。

头文件： #include "tiny_iterative.h"

算子接口¶

算子抽象使迭代方法可以与 任何线性变换 配合使用，无需存储完整的矩阵。

typedef struct {
    int rows, cols; void *ctx;
    tiny_error_t (*matvec)(const float *x, float *y, void *ctx);
    tiny_error_t (*rmatvec)(const float *x, float *y, void *ctx);
} tiny_operator_f32_t;

预置算子¶

初始化函数	矩阵结构	存储	用途
`tiny_dense_operator_init`	通用连续矩阵	\(O(m \cdot n)\)	通用
`tiny_hankel_operator_init`	信号 \(s(t)\) 的 Hankel 矩阵	\(O(N)\)	ERA/SSI 相关（隐式）
`tiny_block_hankel_operator_init`	\(R_{lag}[ch][ch]\) 的块 Hankel 矩阵	\(O(n_{ch}^2 \cdot N_{lags})\)	SSI Toeplitz-free
`tiny_covariance_operator_init`	特征值的协方差矩阵	\(O(n_{ch} \cdot i_{max})\)	SSI 协方差

块 Hankel 内存节省

无需构建 \(T_{dim}^2\) 的 Toeplitz 矩阵（\(O(n^2)\) 内存），算子即时计算 matvec（\(O(N \cdot n)\) 时间）。\(T_{dim} = 100\) 时节省 40 KB RAM。

Arnoldi 迭代¶

size_t       tiny_arnoldi_f32_workspace_size(int n, int k);
tiny_error_t tiny_arnoldi_f32(const tiny_operator_f32_t *op,
                               const float *v0, int n, int k,
                               float *V, float *H,
                               void *ws, size_t ws_size);

Arnoldi = 基于 MGS 的 Krylov 子空间构建。 给定算子 \(A\) 和起始向量 \(v_0\)，计算 \(V_k = [v_1, ..., v_k]\) 张成 \(\mathcal{K}_k(A, v_0)\)，以及 Hessenberg 矩阵 \(H_k\)。

参数：

参数	类型	说明
`op`	`const tiny_operator_f32_t *`	线性算子（matvec）
`v0`	`const float *`	起始向量，长度 \(n\)
`n`	`int`	问题维度
`k`	`int`	Krylov 子空间维度（迭代次数）
`V`	`float *`	输出。 \(n \times k\) 基，列优先（`V[i + j*n]` = 第 \(j\) 列）
`H`	`float *`	输出。 \(k \times k\) 上 Hessenberg（行优先）
`ws`	`void *`	工作空间缓冲区
`ws_size`	`size_t`	工作空间字节数

返回： 成功返回 TINY_OK。

\(H_k\) 的特征值近似 \(A\) 的极端特征值（Ritz 值）。

用途： 前 \(k\) 个特征值的子空间迭代、稀疏特征值问题。

Lanczos 迭代¶

size_t       tiny_lanczos_f32_workspace_size(int n, int k);
tiny_error_t tiny_lanczos_f32(const tiny_operator_f32_t *op,
                               const float *v0, int n, int k,
                               float *alpha, float *beta, float *V,
                               void *ws, size_t ws_size);

Lanczos = 对称算子的 Arnoldi（\(A = A^T\)）。生成三对角矩阵而非完整的 Hessenberg。

参数：

参数	类型	说明
`op`	`const tiny_operator_f32_t *`	线性算子（必须对称）
`v0`	`const float *`	起始向量，长度 \(n\)
`n`	`int`	问题维度
`k`	`int`	Krylov 子空间维度
`alpha`	`float *`	输出。对角线元素（\(k\) 个值）
`beta`	`float *`	输出。次对角线元素（\(k\) 个值）
`V`	`float *`	输出。 \(n \times k\) 基，列优先
`ws`	`void *`	工作空间缓冲区
`ws_size`	`size_t`	工作空间字节数

性能： 比 Arnoldi 快 2 倍（每步只需一次 MGS 正交化）。

正交性丧失

在 float32 中，Lanczos 约 \(\sim 10\) 次迭代后因舍入丧失正交性。未内置重正交化；若正交性至关重要，请使用 Arnoldi。

随机 SVD¶

使用随机投影近似 \(A \approx U \cdot \Sigma \cdot V^T\)（Halko, 2011）。

算法¶

1. Ω = 随机 n × (k+p) 高斯矩阵
2. Y = A · Ω              — 草图 m × (k+p)
3. Y = Q · R              — 草图的 QR
4. B = Q^T · A            — 投影到 (k+p) × n
5. B = Ũ · Σ · V^T       — 小型 B 的完整 SVD
6. U = Q · Ũ,  Σ,  V     — 重构

稠密矩阵形式¶

size_t       tiny_rsvd_dense_f32_workspace_size(int m, int n, int k, int p, int q);
tiny_error_t tiny_rsvd_dense_f32(const float *A, int m, int n, int k, int p, int q,
                                   float *U, float *S, float *Vt,
                                   void *ws, size_t ws_size);

算子形式¶

size_t       tiny_rsvd_f32_workspace_size(int m, int n, int k, int p, int q);
tiny_error_t tiny_rsvd_f32(const tiny_operator_f32_t *op, int m, int n, int k,
                            float *U, float *S, float *Vt,
                            void *ws, size_t ws_size);

参数：

参数	类型	说明
`A` / `op`	—	输入矩阵或算子
`m`, `n`	`int`	矩阵维度
`k`	`int`	目标秩（要计算的奇异值/向量数）
`p`	`int`	过采样（默认 5）。\(k\) 之外的额外列
`q`	`int`	幂迭代次数（默认 1）。每次提升慢衰减奇异值的精度
`U`	`float *`	输出。 \(m \times k\) 左奇异向量（列优先）
`S`	`float *`	输出。 \(k\) 个奇异值（降序）
`Vt`	`float *`	输出。 \(k \times n\) 右奇异向量（行优先）
`ws`	`void *`	工作空间缓冲区
`ws_size`	`size_t`	工作空间字节数

工作空间： \(\sim (m \cdot r + n \cdot r + r^2)\) floats，\(r = k + p\)。\(m=75, n=50, k=10\) 时：\(\sim 8\) KB。

RSVD 性能

\(k=10, p=5, q=1\) 的 RSVD 对 \(50 \times 50\) 矩阵提供与完整 Jacobi 相当的 SVD 精度，速度快 10 倍。作为配对 ERA 实现中的主要 SVD 使用。

SVD 工具函数¶

void  tiny_svd_cumulative_energy_f32(const float *S, int len, float *energy_out);
int   tiny_svd_rank_by_energy_f32(const float *S, int len, float threshold);
int   tiny_svd_rank_by_threshold_f32(const float *S, int len, float tol);
float tiny_svd_condition_estimate_f32(const float *S, int len);

函数	用途	示例
`cumulative_energy`	\(\text{energy}[i] = \frac{\sum_{j=0}^i S_j^2}{\sum S^2}\)	跟踪多少 SV 解释 99% 能量
`rank_by_energy`	累积能量 \(\geq\) 阈值的第一个索引	`rank_by_energy(S, len, 0.99)` → 模型阶数
`rank_by_threshold`	计数 \(S_k > \text{tol} \cdot S_0\)	`rank_by_threshold(S, len, 1e-3)` → 数值秩
`condition_estimate`	\(S_0 / S_{n-1}\)	条件数 \(\kappa\) — 若 \(> 10^6\)，病态

工程说明¶

随机种子： RSVD 使用 Box-Muller 高斯 RNG，静态状态。引导周期之间是确定性的。
幂迭代： 每 \(q\) 次放大接近快速衰减的奇异值。对于大多数呈指数衰减的 SHM 数据，\(q=1\) 足够。
算子 vs 稠密： 如果矩阵适合 RAM（\(< 100 \times 100\)），使用 tiny_rsvd_dense_f32。算子形式仅对太大而无法显式存储的矩阵有帮助（例如作为算子的 Toeplitz）。