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

**URL:** <https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280>\
**Category:** Performance\
**Tags:** linearalgebra, matrices, linearsolve\
**Created:** [June 6, 2024, 12:52pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280 "2024-06-06T12:52:10Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [June 6, 2024, 12:52pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/1 "2024-06-06T12:52:10Z")

</div>

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

```julia
    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.

---

<div class="post-metadata">

**Author:** ![balaji\_sriram](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balaji_sriram/32/208000_2.png) [@balaji\_sriram](https://discourse.julialang.org/u/balaji_sriram)\
**Post date:** [June 6, 2024, 1:15pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/2 "2024-06-06T13:15:23Z")

</div>

Please look into [linearsolve](https://docs.sciml.ai/LinearSolve/stable/). 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](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](https://in.mathworks.com/help/matlab/math/iterative-methods-for-linear-systems.html)

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [June 6, 2024, 1:21pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/3 "2024-06-06T13:21:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 6, 2024, 1:22pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/4 "2024-06-06T13:22:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [June 6, 2024, 2:23pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/5 "2024-06-06T14:23:13Z")

</div>

> [@ChrisRackauckas](#):
>
> Which is better depends on size

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](https://discourse.julialang.org/t/how-to-employ-nlsolve-to-a-function-which-has-parameters-which-are-difficult-to-calculate/115126) I am trying to implement. The matrices I am trying to work with have small values of entries like

```julia
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.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [June 6, 2024, 2:40pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/6 "2024-06-06T14:40:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [June 6, 2024, 3:06pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/7 "2024-06-06T15:06:41Z")

</div>

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

> [@Package to avoid allocations in some functions when using StaticArrays.jl](https://discourse.julialang.org/t/package-to-avoid-allocations-in-some-functions-when-using-staticarrays-jl/114539):
>
> Hi! As I mentioned before, I am trying to build a full satellite attitude control subsystem (ACS) using Julia, which maybe one day will hopefully fly in a CubeSat. One very important requirement is to avoid runtime allocations. The entire system is built using StaticArrays from StaticArrays.jl and I manage to create a mockup of an entire ACS without allocations except for one single case: the Kalman filters. We need to use at least two Kalman filters in our project: one for the attitude and o…

---

<div class="post-metadata">

**Author:** ![Ronis\_BR](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ronis_br/32/50999_2.png) [@Ronis\_BR](https://discourse.julialang.org/u/Ronis_BR)\
**Post date:** [June 6, 2024, 4:19pm UTC](https://discourse.julialang.org/t/solving-linear-system-bc-x-where-b-is-a-symmetric-positive-semi-definite-matrix/115280/8 "2024-06-06T16:19:18Z")

</div>

> [@gdalle](#):
>
> Have you considerd using [StaticArrays.jl](https://juliahub.com/ui/Packages/General/StaticArrays), perhaps in conjunction with the custom BLAS methods mentioned here?

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

> <https://github.com/JuliaArrays/StaticArrays.jl/pull/1259>
>
> This PR implements direct calls to BLAS library for computing SVDs, reducing the… allocations if the input types are \`SMatrix\` ou \`MMatrix\` with real numbers.
> 
> The approach is the convert the input matrix to \`MMatrix\`, call BLAS function, and convert it back to the input type. Since \`MMatrix\` does not leave the function scope, the compiler is clever enough to avoid allocations.
> 
> With this PR, we can call \`svdvals\`, \`svd\`, and \`pinv\` with 0 allocations.
> 
> The current implementation of \`svd\` and \`svdvals\` converts the static matrix to \`Matrix\` and call \`LinearAlgebra.svd\`. Here, we are calling the same BLAS function as in LinearAlgebra stdlib. Hence, the performance gain here seems "free" except for maintaining new code inside StaticArrays.
> 
> This PR does not add any noticeable delay when importing StaticArrays:
> 
> Before: 0.126202 seconds (253.52 k allocations: 13.714 MiB)
> After: 0.128506 seconds (253.58 k allocations: 13.751 MiB)
> 
> Here are some benchmarks comparing the new and old implementation:
> 
> \`svd\`
> 
> | \*\*Dimension\*\* | \*\*SMatrix (old)\*\* | \*\*SMatrix (new)\*\* | \*\*MMatrix (old)\*\* | \*\*MMatrix (new)\*\* |
> |--------------:|-------------------------:|-------------------------:|--------------------------:|-------------------------:|
> | 2 | 0.911 us (8 allocations) | 0.441 us (0 allocations) | 0.907 us (11 allocations) | 0.450 us (3 allocations) |
> | 3 | 1.328 us (8 allocations) | 0.875 us (0 allocations) | 1.504 us (11 allocations) | 1.042 us (3 allocations) |
> | 4 | 2.102 us (8 allocations) | 1.562 us (0 allocations) | 2.250 us (11 allocations) | 1.721 us (3 allocations) |
> | 5 | 2.893 us (8 allocations) | 2.343 us (0 allocations) | 2.847 us (11 allocations) | 2.292 us (3 allocations) |
> | 6 | 3.651 us (8 allocations) | 3.026 us (0 allocations) | 4.054 us (11 allocations) | 3.432 us (3 allocations) |
> 
> We also have a performance boost when computing \`pinv\` because it basically calls \`svd\`. For example, computing \`pinv\` in a 3 x 3 matrix takes 0.995 us in this commit versus 1.462 us in \`master\`.
> 
> Closes #1255