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

General Usage Performance
, ,
1

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.

2

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())

The answer tyo your questions are available in detail in this video by @ChrisRackauckas : https://www.youtube.com/watch?v=JWI34_w-yYw

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
3

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

4

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

1 Like
5

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.

6

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
7

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

3 Likes
8

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

4 Likes