简介面向数值线性代数研究人员与工程开发者的正交矩阵正弦余弦分解CSD工具包可对 m×m 分区正交矩阵执行分解计算。分区结构意味着矩阵被划分为多个子块并支持单独处理适合大规模矩阵运算、并行计算以及控制系统设计、图像处理、信号分析等应用场景。压缩包共 5 个文件包含 C 源码实现的 CSD 核心算法、头文件、MATLAB 构建与调用脚本、接口文件及开源许可说明整体仅 6KB结构紧凑便于抽取集成或二次开发。已有 329 人学习浏览适合已具备基本矩阵理论基础、正在做算法落地或课程设计的读者。借助该工具使用者可通过 MATLAB 或 C 语言环境直接完成分解运算免去手写 LAPACK 底层调用的复杂度有效缩短开发周期同时提升数值计算稳定性与可并行性。1. 正交矩阵的正弦余弦分解CSD到底在算什么先说结论CSD即 Cosine-Sine Decomposition中文常译为正弦余弦分解是专门处理“分块正交矩阵”的一套工具。你大概熟悉对称矩阵的对角化也熟悉一般矩阵的 SVD但当你面对的是一个带预设分块的正交矩阵比如控制论中的传递矩阵、信号处理中的子空间变换直接套 SVD 会把分块边界打散得到的奇异向量不再对应你关心的行块和列块。CSD 的核心思路是用两个块对角正交矩阵去夹原始矩阵把块之间的信息压缩成一组余弦值和一组正弦值。这组角度就是两个互补子空间之间的主夹角。下面从数学形式开始给你一条能直接落地的链路定义、代码、参数、验证技巧。2. 从正交矩阵到 C 和 SCSD 的数学形式与存在条件CSD 是广义 SVD 的一种特殊形态但它比 SVD 多了一层“分块约束”。要理解这一点建议先看最干净的平衡分块。2.1 平衡分块2p×2p 的正交矩阵怎么拆设 Q 是 2p×2p 的正交矩阵按行和列各切一道得到四个 p×p 块Q [ Q11 Q12 ] [ Q21 Q22 ]CSD 断言存在 p×p 的正交矩阵 U1、U2、V1、V2使得[ U1 0 ]^T [ Q11 Q12 ] [ V1 0 ] [ C S ] [ 0 U2 ] [ Q21 Q22 ] [ 0 V2 ] [ -S C ]其中 C 和 S 都是 p×p 对角矩阵C diag(cos θ1, ..., cos θp) S diag(sin θ1, ..., sin θp)θi 的取值范围是 [0, π/2]。因为正交矩阵保证 C² S² I所以每一对 (cos θi, sin θi) 对应一个平面旋转角。这个形式非常直观四个块被化简成一个“旋转块”左右两个块对角正交矩阵只负责调整每个子空间内部的坐标真正跨块的信息全部落在角度上。2.2 C 和 S 从哪来对称矩阵特征分解是关键很多人会问为什么偏偏是余弦和正弦答案来自正交矩阵的 Gram 矩阵约束。由 QᵀQ I 可以直接推出Q11ᵀ Q11 Q21ᵀ Q21 I Q12ᵀ Q12 Q22ᵀ Q22 I这两个式子说明Q11 的 Gram 矩阵和 Q21 的 Gram 矩阵互为补余。取 Q11 的 SVDQ11 U1 C V1ᵀ代回去会发现 V1 恰好也能对角化 Q21ᵀ Q21V1ᵀ Q21ᵀ Q21 V1 I - C² S²也就是说C 和 S 本质上来自两个互补对称半正定矩阵的同步对角化。这里的 C² 和 S² 就是特征值分解结果而 U1、V1 是特征向量矩阵U2 再由 Q21 V1 补全。这种“互补分解”和初等矩阵的关系也很直接当 p1 时CSD 退化为一个初等平面旋转矩阵也就是 Givens 旋转变换。所以你可以把 CSD 理解为“把一组 Givens 旋转叠加到分块正交矩阵上”。这也就解释了为什么对称矩阵在 CSD 里地位特殊你不需要对 Q 整体做特征分解只需要对 Q11 的对称 Gram 矩阵做特征分解就能拿到角度。这也是后续代码实现的数学基础。2.3 非平衡分块p≠q 时的可复现处理实际使用时左上块不一定是方阵。设 Q11 是 p×q那么最终块矩阵需要用单位阵和零矩阵补足才能保证整个变换仍然正交。常见做法是当 p q 时在 C 左边补 I在 S 右边补 0使块矩阵的行列尺寸匹配当 p q 时在 C 下边补 I在 S 下边补 0并对符号做对应调整。具体排列顺序和符号约定在不同文献里略有差异但用户不需要手推。你只要向 SciPy 的cossin函数传对 p、q它会自动处理非平衡分块。理解平衡分块已经足够掌握 CSD 的核心行为非平衡情况只是用 I/0 做维度对齐。3. 用 Python 把 CSD 跑起来最小可复现代码理论立住之后最关心的问题一定是给我一段能跑的代码。这里给出基于 NumPy 和 SciPy 的最小实现直接用scipy.linalg.cossin计算 CSD。这个函数是数值稳定的标准实现比手写两次 SVD 再拼装要可靠得多。3.1 安装依赖并构造一个正交矩阵先保证环境里有 NumPy 和 SciPy然后生成一个随机正交矩阵用于演示。生成正交矩阵的常见做法是对高斯随机矩阵做 QR 分解import numpy as np from scipy.linalg import cossin n 6 p q 3 rng np.random.default_rng(42) A rng.standard_normal((n, n)) Q, _ np.linalg.qr(A) print(np.round(Q Q.T, 6))这段代码先用standard_normal生成 6×6 随机矩阵再通过qr分解得到正交矩阵 Q。最后打印 Q Qᵀ结果应该是单位阵用来确认构造没有出错。3.2 调用 cossin 并解读返回值对平衡分块直接传 p 和 q 即可u, cs, vh cossin(Q, p, q) print(np.round(cs, 4)) print(重建误差, np.linalg.norm(Q - u cs vh))cossin返回三个对象u是左块对角正交矩阵vh是右块对角正交矩阵的共轭转置cs是中间那个带 C/S 的对角块矩阵。cs的左上块接近对角阵右上和左下接近正弦对角块右下接近负余弦对角块。重建误差在 1e-14 量级说明分解正确。3.3 用 separateTrue 直接拿角度很多时候你只想拿主角度不想自己从cs里拆。cossin提供了separate参数u1, u2, theta, v1h, v2h cossin(Q, p, q, separateTrue) print(theta (弧度), theta) print(theta (角度), np.degrees(theta))separateTrue时返回五个对象u1、u2是两个对角块中的左正交矩阵theta是主角度数组v1h、v2h是对应的右正交矩阵共轭转置。theta的单位是弧度范围在 [0, π/2]按降序排列。拿到theta后可以用cos(theta)和sin(theta)重建 C 和 S。3.4 验证 cos² sin² 1CSD 最容易被忽略的一点是角度必须满足互补关系。以下代码验证数值结果是否满足这个恒等式cos_vals np.cos(theta) sin_vals np.sin(theta) print(np.max(np.abs(cos_vals**2 sin_vals**2 - 1)))输出应该非常接近 0。如果这个值偏大说明你传入的矩阵不是正交矩阵或者分块尺寸与矩阵结构不匹配。这是排查 CSD 结果是否可信的第一道门。4. CSD 的参数怎么设四个绕不开的细节CSD 的 API 看起来简单真正用起来容易在四个地方翻车分块尺寸、角度顺序、病态角度和语义选型。4.1 p 和 q 的维度匹配你必须知道自己在切哪里的子空间cossin(Q, p, q)中p 是左上块的行数q 是左上块的列数。它必须满足 0 ≤ p ≤ n0 ≤ q ≤ n。这里的 p 和 q 决定了主角度描述的是“前 p 行子空间与前 q 列子空间之间的夹角”。如果你误以为 p 和 q 必须满足 pqn就会把本来无关的维度塞进分解里。实际应用中p 和 q 可以独立选择但一定要先想清楚你关心的两个子空间分别对应哪些行、哪些列。比如在控制系统里行块对应输入通道分解列块对应状态通道分解p、q 分别取输入维数和状态维数。4.2 角度顺序的陷阱SciPy 返回降序theta按降序排列这与np.linalg.svd的奇异值排序一致。但如果你自己用np.linalg.eigh对 Q11ᵀ Q11 做特征分解默认返回的是升序特征值得到的特征向量顺序刚好相反。此时如果不做np.flip就直接拼 C 和 S分块结构会乱掉。建议统一走cossin的返回值避免手动排序。若必须自己写教学版记得先确认顺序。4.3 角度趋近 0 或 π/2特征值模糊问题当某个 θi 接近 0 时C 的对角元素接近 1S 接近 0当 θi 接近 π/2 时情况反过来。数值上Q11ᵀ Q11 的特征值会非常接近 1 或 0特征分解可能产生微小负特征值导致sqrt报 NaN。一个常见补救是eigvals np.clip(eigvals, 0.0, 1.0) c_vals np.sqrt(eigvals) s_vals np.sqrt(1.0 - c_vals**2)这个 clip 操作能避免负特征值带来的 NaN但也会掩盖真实的病态。更稳妥的做法是检查奇异值序列是否存在多个数量级差异。如果某些角度相差超过 12 个数量级建议减小分块规模或者改用cossin内部的均衡化算法不要手工处理。4.4 别拿 CSD 当 SVD 用有一个常见误用把 CSD 输出的角度当作普通奇异值直接用阈值截断来降维。CSD 的目标是保留分块结构而 SVD 的目标是全局最优低秩近似。两者的目标函数不同维度CSDSVD适用矩阵方阵正交矩阵任意矩阵输出结构块对角正交矩阵 C/S 对角块三个满矩阵核心信息子空间主夹角矩阵的奇异值降维能力不直接支持低秩近似直接支持截断如果你需要对正交矩阵做低秩近似应该先提取矩阵的线性部分再做普通 SVD。CSD 更适合分析两个子空间之间的几何关系而不是压缩矩阵本身。5. 进阶技巧用 CSD 构造带指定主夹角的正交矩阵理解了 CSD 之后你可以反向利用它给定一组角度构造一个正交矩阵使其主夹角精确等于这些角度。这个技巧非常适合做单元测试也可以用于生成仿真数据。5.1 构造公式核心公式就是第 2 章的分解逆过程。给定角度 theta随机生成四个正交矩阵 U1、U2、V1、V2再拼装中间旋转块def make_orth_with_angles(theta, seed0): rng np.random.default_rng(seed) p len(theta) n 2 * p def rand_orth(): A rng.standard_normal((p, p)) return np.linalg.qr(A)[0] U1, U2 rand_orth(), rand_orth() V1, V2 rand_orth(), rand_orth() C np.diag(np.cos(theta)) S np.diag(np.sin(theta)) mid np.block([[C, S], [-S, C]]) U np.block([[U1, np.zeros((p, p))], [np.zeros((p, p)), U2]]) V np.block([[V1, np.zeros((p, p))], [np.zeros((p, p)), V2]]) return U mid V.T这段代码先把 θ 转成 C 和 S再组装中间分块矩阵最后用左右块对角正交矩阵做相似变换。由于 U、V 和 mid 都正交乘积 Q 一定正交。5.2 用 cossin 验证构造结果构造完成后马上用cossin反向提取角度对比输入值theta_in np.array([0.0, 0.3, 1.2]) Q make_orth_with_angles(theta_in, seed7) _, _, theta_out, _, _ cossin(Q, len(theta_in), len(theta_in)) print(np.degrees(theta_out))由于 CSD 的角度排序是降序而输入的 theta 不一定有序所以先对输入做降序排序再比较。这一步能同时验证两件事构造公式是否正确以及cossin的结果是否稳定。如果残差接近零说明你的正交矩阵生成流程可以作为测试夹具复用。5.3 一个小技巧用主夹角做正交矩阵的可视化拿到 theta 之后可以直接用平面旋转的视角观察一个正交矩阵。把 theta 画成条形图每个条的长度代表一个子空间夹角角度越接近 π/2表示对应分量子空间越正交角度接近 0则表示两个子空间有较大重合。这种可视化比看 6×6 的矩阵数值直观得多也更容易向非数值背景的人解释 CSD 的物理含义。对于高维正交矩阵CSD 把几何信息压缩成 p 个标量本身就是一种降维表达。你在调试控制系统或分析高维信号时不妨把 CSD 角度序列作为第一诊断指标。本文还有配套的精品资源点击获取