Julia中计算对称矩阵二次型的最高效最快实现方法

稠密矩阵场景 #

• 普通对称矩阵的最优实现
  不要直接用通用矩阵运算流程计算，优先利用mP的对称性减少运算量：
  1. 手动实现的话，只遍历mP的下三角（含对角线）元素，按公式 Σ (i≥j) mP[i,j] * vX[i] * vX[j] * (1 if i==j else 2) 直接求和，浮点运算量比全量遍历少一半，实测速度可提升40%左右。
  2. 调用BLAS库的话，优先用对称矩阵向量乘接口dsymv（双精度为例）计算tmp = mP · vX，再用ddot计算vX和tmp的点积，比通用矩阵乘gemv的效率高30%以上。
• 正定矩阵的额外性能提升
  仅当同一个mP需要计算≥2个不同vX的二次型时，该优化才有收益：
  提前对mP做Cholesky分解得到下三角矩阵L（满足mP = L · Lᵀ），此时二次型等价于||Lᵀ · vX||₂²，后续每次计算只要先做三角矩阵向量乘y = Lᵀ · vX，再对y的所有元素求平方和即可。三角矩阵运算的缓存命中率比对称矩阵运算高，单次二次型计算速度可比用dsymv的方案再快20%~30%。如果是单次计算场景，Cholesky分解的O(n³)开销远大于收益，不推荐使用。

稀疏矩阵场景 #

• 普通对称矩阵的最优实现
  绝对不要转换为稠密矩阵计算，优先使用仅存储下三角非零元素的对称稀疏存储格式（比如SYM-CSR），遍历所有非零元素按前述折叠公式直接求和，时间复杂度仅和非零元素数量nnz成正比，是O(nnz)级别，稀疏度越高收益越大。调用稀疏线性代数库的话，直接用对称稀疏矩阵向量乘接口再做点积即可，避免无意义的零元素运算。
• 正定矩阵的额外性能提升
  同样适合同一个mP多次计算的场景：提前做填充优化的稀疏Cholesky分解，通过消元树排序减少分解后L矩阵的非零元素数量，后续每次二次型计算只要做稀疏三角矩阵乘向量再算平方和，比直接用对称稀疏矩阵乘的方案快40%~60%，而且分解后的L矩阵内存占用也远小于原始稀疏矩阵。单次计算场景同样不推荐预分解。

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