# Solving linear system BC = X where B is a symmetric positive semi-definite matrix

I am trying to solve for X in BC=X where B is a symmetric positive semi-definite matrix and X is simply a vector with 1 in all of it’s entries. As far as I understand, this can be solved using `C=B\X` , but I was wondering that how can I take advantage of the positive semi-definite nature of B in order to make the solving faster or more numerically stable. I was using the below code

``````    B = zeros(p, p)
for i in 1:p
for j in 1:i
B[i,j] = dot(R_iter_storage[i],R_iter_storage[j])
B[j,i] = B[i,j]
end
end
B = Symmetric(B)
X = ones(p)
C = B\X
``````

I know that as I have declared `B` as a symmetric matrix, Julia would automatically use the optimized algorithm for symmetric matrices instead of the generalized solving algorithm. However, I have an additional constraint on B which is that it is a positive semi-definite matrix. Any ideas on what functions to use to solve the system such that the later property is also exploited ? Note - Here `R_iter_storage` is a vector whose elements are tensors of rank 4.

Please look into linearsolve. This package lets you solve system of linear equations very efficiently. Since your matrix is symm-positive definite, the method is recomended: LinearSolve.KrylovJL_CG()

Basically the command is :
using LinearAlgebra, LinearSolve, SparseArrays (if you have sparse vectors/matrices)
prob = LinearProblem(B, X);
sol = solve(prob,LinearSolve.KrylovJL_CG())

Also here is a matlab help center link with a flowchart on which method you have to use and when (just to give you an idea): Iterative Methods for Linear Systems - MATLAB & Simulink - MathWorks India

2 Likes

Perhaps you could compute a Cholesky decomposition and use that in the linear solve?

Cholesky or CG can both be used. Which is better depends on size (and if you have a good preconditioner).

1 Like

The matrices I would be working with are in the order of `6x6`. Using the `LinearSolve.KrylovJL_CG()` seems to be causing numerical instability at the later stages of the DIIS solver for Coupled Cluster Equations I am trying to implement. The matrices I am trying to work with have small values of entries like

``````6×6 Symmetric{Float64, Matrix{Float64}}:
1.72458e-16  1.06023e-16  1.01555e-16  9.45545e-17  8.19558e-17  7.61703e-17
1.06023e-16  8.87195e-17  8.45981e-17  7.75149e-17  6.32949e-17  5.68955e-17
1.01555e-16  8.45981e-17  8.3458e-17   7.63146e-17  6.20451e-17  5.55372e-17
9.45545e-17  7.75149e-17  7.63146e-17  7.54038e-17  5.98022e-17  5.33082e-17
8.19558e-17  6.32949e-17  6.20451e-17  5.98022e-17  5.52861e-17  4.88773e-17
7.61703e-17  5.68955e-17  5.55372e-17  5.33082e-17  4.88773e-17  4.73348e-17
``````

As of now I have not provided the exact code where the `LinearSolve` is used, but would be happy to provide if necessary.

I would start with just `cholesky(Hermitian(B)) \ X` (you can use `Symmetric` instead of `Hermitian`, it makes no difference for real matrices). Using the Cholesky factorization will bestow the numerical and efficiency benefits you’re after.

Although if your matrix is singular (not strictly positive definite, or sufficiently close to singular that numerical issues arise) then `cholesky` will fail. In that situation, you will need to use `pinv(B) * X` or one of these iterative methods others have proposed.

3 Likes

Have you considerd using StaticArrays.jl, perhaps in conjunction with the custom BLAS methods mentioned here?

3 Likes

No need anymore when StaticArrays.jl 1.9.5 is released! The code was recently added to `main`

4 Likes