A+BBᵀ行列式计算中如何利用分块矩阵性质？

Question

在计算形如 $\det(A + BB^\top)$ 的行列式时，当 $A$ 为可逆矩阵，如何利用分块矩阵的行列式公式简化计算？常见问题在于：直接展开 $A + BB^\top$ 难以处理，尤其是当 $B$ 为 tall 矩阵（列数远小于行数）时。此时可构造分块矩阵 $$ M = \begin{bmatrix} I & B^\top \\ B & -A \end{bmatrix}, $$ 并利用分块矩阵行列式性质与Schur补技巧，推导出 $\det(A + BB^\top) = \det(A) \cdot \det(I + B^\top A^{-1} B)$。但实际应用中，若 $A$ 奇异或高维，该变换失效。如何在不显式求逆的前提下，有效结合矩阵分解（如LU或QR）与分块结构优化计算效率？

Answer 1

1. 问题背景与核心动机

在机器学习、高维统计推断和优化算法中，频繁出现形如 $\det(A + BB^\top)$ 的行列式计算。当 $A \in \mathbb{R}^{n\times n}$ 可逆且 $B \in \mathbb{R}^{n\times k}$（$k \ll n$）为 tall 矩阵时，直接展开 $A + BB^\top$ 将导致 $O(n^3)$ 计算复杂度，效率低下。

利用分块矩阵技巧可将原问题转化为低维子空间上的计算：

$$ M = \begin{bmatrix} I_k & B^\top \\ B & -A \end{bmatrix} \Rightarrow \det(M) = \det(-A)\cdot\det(I_k + B^\top A^{-1}B) $$

结合 Schur 补理论，可得重要恒等式：

$$ \det(A + BB^\top) = \det(A) \cdot \det(I_k + B^\top A^{-1} B) $$

该公式将 $n\times n$ 的行列式分解为一个 $k\times k$ 小矩阵的行列式，显著降低计算负担。然而，当 $A$ 奇异或维度极高时，显式求逆 $A^{-1}$ 不可行，需引入矩阵分解技术规避逆运算。

2. 分块矩阵与Schur补基础回顾

• Schur补定义：对于分块矩阵 $ \begin{bmatrix} P & Q \\ R & S \end{bmatrix} $，若 $P$ 可逆，则其Schur补为 $S - RP^{-1}Q$，且有：
• $$ \det\begin{bmatrix} P & Q \\ R & S \end{bmatrix} = \det(P)\cdot\det(S - RP^{-1}Q) $$
• 应用于构造矩阵 $M$，选择 $P=I_k$，则对应Schur补为 $-A - B I_k^{-1} B^\top = -(A + BB^\top)$
• 因此 $\det(M) = \det(I_k)\cdot\det(-(A + BB^\top)) = (-1)^n \det(A + BB^\top)$
• 另一方面，交换块顺序后使用另一方向的Schur补，可导出含 $A^{-1}$ 的表达式
• 关键洞察：通过不同分块策略，建立原始行列式与低秩更新之间的联系

3. 显式求逆的问题与挑战

场景	问题描述	影响
$A$ 奇异	无法计算 $A^{-1}$	公式失效
$n$ 极大（如 $n>10^4$）	存储和求逆成本过高	内存溢出或超时
条件数大	数值不稳定	结果不可靠
稀疏结构存在	直接求逆破坏稀疏性	丧失性能优势
需多次调用	重复求逆开销大	实时性差

4. 基于矩阵分解的替代方案设计

为避免显式求逆，采用如下分解策略：

1. LU分解：设 $A = PLU$，则 $A^{-1} = U^{-1}L^{-1}P^\top$，但无需显式构造，而是通过前代后代解线性系统 $Ax = b$ 实现 $A^{-1}B$ 的隐式计算
2. QR分解：对 $B$ 进行 $B = QR$，其中 $Q \in \mathbb{R}^{n\times k}, R \in \mathbb{R}^{k\times k}$，则：
3. $$ BB^\top = QRR^\top Q^\top \Rightarrow A + BB^\top = A + Q(RR^\top)Q^\top $$
4. 此时问题转为对 $A$ 添加一个低秩修正项，适合使用Woodbury恒等式处理逆矩阵，而行列式可通过下述引理计算：
5. $$ \det(A + UCV^\top) = \det(C^{-1} + V^\top A^{-1}U)\cdot\det(A)\cdot\det(C) $$
6. 特别地，当 $U=V=Q, C=RR^\top$，即得到：
7. $$ \det(A + BB^\top) = \det(RR^\top + Q^\top A^{-1} Q) \cdot \det(A) \cdot \det((RR^\top)^{-1}) \cdot \det(RR^\top) $$
8. 化简后仍归结为 $k\times k$ 矩阵行列式计算

5. 数值实现流程图与算法结构

def compute_logdet_A_plus_BBt(A, B):
    # 输入：A (n,n), B (n,k)，k << n
    # 输出：log(det(A + BB^T))

    if is_invertible(A):
        # 使用Cholesky或LU分解解多个右端项
        L = cholesky_decomp(A)  # 或 LU
        A_inv_B = solve_triangular_systems(L, B)  # 求 A^{-1}B
        C = I_k + B.T @ A_inv_B
        log_det_A = 2 * sum(log(diag(L)))  # 若Cholesky
        log_det_C = logdet(C)
        return log_det_A + log_det_C
    else:
        # 使用广义Schur补或正则化
        A_reg = A + eps * I
        return compute_logdet_A_plus_BBt(A_reg, B)

graph TD A[输入矩阵 A 和 B] --> B{A 是否可逆？} B -- 是 --> C[对A进行LU/Cholesky分解] B -- 否 --> D[添加小扰动正则化] D --> C C --> E[求解 A^{-1}B 隐式] E --> F[构造小矩阵 I + B^T A^{-1} B] F --> G[计算 det(A) 和 det(小矩阵)] G --> H[返回 det(A + BB^T)]

6. 高级优化技巧与扩展应用

• 低秩更新下的迭代更新：若 $B$ 动态变化，可用 Sherman-Morrison-Woodbury 公式快速更新逆和行列式
• 随机投影加速：对大型 $B$ 使用 Johnson-Lindenstrauss 投影预降维
• 分布式计算框架：将 $B = [B_1,\dots,B_m]$ 分块，利用 MapReduce 并行计算部分迹和行列式项
• GPU加速：利用 cuSOLVER 等库实现批量三角求解，提升 $A^{-1}B$ 计算效率
• 稀疏矩阵专用求解器：如 SuiteSparse 对稀疏 $A$ 使用左-右重排减少填充
• 自动微分兼容性：在 PyTorch/TensorFlow 中封装为可导操作，用于梯度回传
• 近似方法：采用随机迹估计（Hutchinson）估算 $\log\det(I + B^\top A^{-1}B)$
• 流形优化场景：在 Stiefel 流形上维护正交约束时，此类行列式常作为正则项出现