Decomp — 矩阵分解¶

概述¶

两种 QR 分解实现：修正 Gram-Schmidt (MGS) 和 Householder QR (HQR)。TinyMath 依赖链中的第 3 层。

头文件： #include "tiny_decomp.h"

如何选择 QR 方法？¶

维度	MGS	HQR
速度	快 (\(O(m \cdot k \cdot n)\))	慢 (\(O(m \cdot n^2)\))
正交性	好（可重正交化）	优秀（后向稳定）
矩阵形状	\(m \geq k\)（高瘦）	\(m \geq n\)（通用）
原地操作	是（\(A\) 被 \(Q\) 覆盖）	否（\(Q\)、\(R\) 单独输出）
适用场景	ERA / ITD / 最小二乘	SSI / 通过双对角化的 SVD

经验法则： ERA 中速度重要且 \(m \gg n\) 时，使用 reorthogonalize = 1 的 MGS。SSI 中正交性至关重要时，使用 HQR。

MGS QR¶

修正 Gram-Schmidt 计算 \(A_{m \times k} = Q_{m \times k} \cdot R_{k \times k}\)。

工作空间大小¶

size_t tiny_mgs_f32_workspace_size(int m, int k);

返回所需工作空间的字节数。当前返回 0（不需要外部工作空间——实现使用局部数组）。

核心分解¶

tiny_error_t tiny_mgs_f32(float *A, int m, int k,
                           float *R,
                           int reorthogonalize,
                           void *ws, size_t ws_size);

参数：

参数	类型	说明
`A`	`float *`	输入/输出。入口：\(m \times k\) 矩阵。出口：被 \(Q\) 覆盖（行优先，列步长 = \(k\)）
`m`	`int`	行数。必须满足 \(m \geq k\)
`k`	`int`	列数（同时也是 \(R\) 的秩）
`R`	`float *`	输出。 \(k \times k\) 上三角矩阵（行优先，步长 = \(k\)）。下三角被清零
`reorthogonalize`	`int`	`0` = 单次，`1` = 双次（推荐）。将正交性从 \(\sim 10^{-4}\) 提升到 \(\sim 10^{-7}\)
`ws`	`void *`	工作空间缓冲区（当前未使用，可为 `NULL`）
`ws_size`	`size_t`	工作空间字节数（当前未使用）

返回： 成功返回 TINY_OK，无效输入返回错误码。

务必使用重正交化

不开启重正交化时，float32 的 MGS 在列条件数超过 \(10^3\) 时会丧失正交性。ERA 中 \(H_0\) 的条件数可能超过 \(10^4\)——始终设置 reorthogonalize = 1。

正交性检查¶

float tiny_mgs_check_orthogonality_f32(const float *Q, int m, int k);

计算 \(\|Q^T Q - I\|_F\)——接近 0 的值表示良好的正交性。分解后用于验证数值质量。

Householder QR¶

Householder QR 使用反射矩阵计算 \(A_{m \times n} = Q_{m \times m} \cdot R_{m \times n}\)。

工作空间大小¶

size_t tiny_hqr_f32_workspace_size(int m, int n);

返回所需工作空间的字节数。包括中间反射向量和累积 \(Q\) 矩阵的内存。

核心分解¶

tiny_error_t tiny_hqr_f32(const float *A, int m, int n,
                           float *Q, float *R,
                           void *ws, size_t ws_size);

参数：

参数	类型	说明
`A`	`const float *`	仅输入。 \(m \times n\) 矩阵（不修改）
`m`	`int`	行数
`n`	`int`	列数
`Q`	`float *`	输出。 \(m \times m\) 正交矩阵（行优先）
`R`	`float *`	输出。 \(m \times n\) 上三角矩阵（行优先）。对角线以下的值被清零
`ws`	`void *`	工作空间缓冲区。至少为 `tiny_hqr_f32_workspace_size(m, n)` 字节
`ws_size`	`size_t`	工作空间字节数

返回： 成功返回 TINY_OK，失败返回错误码。

工程说明¶

每种方法的工作原理¶

MGS 逐列处理： 1. 复制 \(A\) 的第 \(j\) 列 2. 将其投影到之前正交化的每一列 \(q_0 \dots q_{j-1}\) 3. 减去投影——与所有前列正交 4. 归一化——成为 \(q_j\)

HQR 使用一次反射即可消去整个次对角线块——每次反射更高效，但总计需要 \(O(m \cdot n^2)\)。

性能指引¶

问题规模	建议
\(n \leq 8\)	两者均可——速度差异可忽略
\(n \leq 50\)、\(m \gg n\)	使用 `reorthogonalize = 1` 的 MGS（ERA/ITD）
\(n \leq 50\)、通用形状	HQR（SSI）
精度关键	HQR（后向稳定，\(Q^T Q \approx I\) 达到机器精度）

内存布局¶

两种算法都假设 行优先、连续 存储（无 padding，步长 = 列数）。非连续矩阵必须在调用前复制到连续缓冲区。