稀疏矩阵运算与Cholesky分解性能优化问询（LCPsolve.jl）

问题描述 #

我正在求解含变分不等式的微分方程组，当前使用LCPsolve.jl中的LCP求解器。经代码性能分析，发现瓶颈位于指定位置的第216-219行，具体为：

H = Jk'*Jk + μ*I

其中Jk是大小不超过150×150的稀疏矩阵，类型信息：

typeof(Jk) = SparseMatrixCSC{Float64, Int64}
typeof(μ) = Float64

该步骤存在内存分配问题。此外，性能分析还高亮了如下操作：

-(cholesky(H)\Jphi)

类型信息：

typeof(H) = SparseMatrixCSC{Float64, Int64}
typeof(Jphi) = Vector{Float64}

请问是否有方法优化这些步骤的性能及内存分配？

优化方案 #

针对H = Jk'*Jk + μ*I的内存与性能优化 #

• 预分配与原地计算：避免临时矩阵分配，先预分配对称矩阵存储空间，用mul!原地计算Jk'*Jk，再直接给对角线加μ，跳过构造μ*I的临时矩阵：

  # 预分配稀疏矩阵存储空间
  H = spzeros(size(Jk,2), size(Jk,2))
  mul!(H, Jk', Jk)  # 原地执行Jk'*Jk计算
  H.diag .+= μ  # 直接修改对角线值
  

• 利用对称矩阵特性：Jk'*Jk天然对称，用Symmetric包装后，后续分解操作能自动利用对称特性减少计算量与存储：

  H = Symmetric(Jk'*Jk)
  H.data.diag .+= μ
  

针对-(cholesky(H)\Jphi)的优化 #

• 原地求解减少向量分配：预分配结果向量，用ldiv!原地执行求解，避免生成新向量：

  res = similar(Jphi)
  # 利用对称特性构造Cholesky分解
  chol_H = cholesky(Symmetric(H), check=false)
  ldiv!(res, chol_H, Jphi)  # 求解结果直接存入res
  res .*= -1  # 原地取负
  

• 跳过H的显式构造，用迭代法求解：由于H = Jk'*Jk + μ*I，可直接定义矩阵-向量乘操作，用共轭梯度法（CG）求解，完全省去构造H的内存开销：

  using IterativeSolvers
  # 定义H与向量的乘法操作
  function mul_H!(Hv, v, Jk, μ)
      temp = Jk * v
      mul!(Hv, Jk', temp)
      Hv .+= μ .* v
  end
  # 用CG求解Hx = Jphi，结果直接取负
  x = cg((Hv, v) -> mul_H!(Hv, v, Jk, μ), Jphi)
  x .*= -1
  

• 复用分解存储空间：若Jk和μ在多轮迭代中仅小幅变化，可预先分配Cholesky分解的存储空间，避免重复初始化的开销。

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