如何在Numpy/Scipy中高效稳定计算Gram矩阵的对数行列式？

直接计算J.dot(J.T)再求行列式的对数存在两个明显问题：一是额外的矩阵乘法会增加计算量，二是当矩阵条件数较大时，数值稳定性差，容易出现溢出或精度损失。以下是两种更高效稳定的解决方案：

方法一：利用奇异值分解（SVD） #

从数学性质出发，Gram矩阵J*J^T的行列式等于J的所有非零奇异值的平方乘积，因此：
$$\log(\det(JJ^T)) = 2 \times \sum \log(s_i)$$
其中$s_i$是J的奇异值。

Numpy中可以通过经济型SVD直接计算奇异值，无需生成冗余的正交矩阵，大幅提升效率：

import numpy as np
m, n = 100, 150
J = np.random.randn(m, n)
# 仅计算奇异值，跳过U和V矩阵的生成
s = np.linalg.svd(J, full_matrices=False, compute_uv=False)
log_det = 2 * np.sum(np.log(s))

SVD的数值稳定性极强，能轻松应对高条件数矩阵，完全避免了构造Gram矩阵的中间步骤带来的精度风险。

方法二：利用QR分解 #

对J做QR分解（$J=QR$，Q为正交矩阵，R为上三角矩阵），由于正交矩阵的行列式绝对值为1，因此：
$$\det(JJ^T) = \det(RR^T) = \left(\det(R_{m\times m})\right)^2$$
其中$R_{m\times m}$是R矩阵前m列构成的上三角方阵，其行列式等于对角线元素的乘积。由此推导：
$$\log(\det(JJ^T)) = 2 \times \sum \log(|r_{ii}|)$$
$r_{ii}$是R的对角线元素。

通过Scipy可以仅计算R矩阵，进一步压缩计算开销：

import numpy as np
from scipy.linalg import qr

m, n = 100, 150
J = np.random.randn(m, n)
# 仅生成R矩阵，跳过Q矩阵的构造
R = qr(J, mode='r')
log_det = 2 * np.sum(np.log(np.abs(np.diag(R))))

QR分解的计算速度通常比SVD更快，数值稳定性也远优于直接构造Gram矩阵，适合对性能要求较高的场景。

底层LAPACK例程说明 #

• Numpy的np.linalg.svd底层调用LAPACK的dgesvd（双精度）或cgesvd（单精度）例程；
• scipy.linalg.qr的mode='r'模式调用LAPACK的dgeqrf例程，仅计算R矩阵，无需构造Q矩阵，进一步降低计算成本。

内容的提问来源于stack exchange，提问作者Euler_Salter