Solving a linear system with a symmetric positive definite matrix

[loiseaujc]

Motivation

I'm currently doing some more modernization work on the QuadProg solver for solving strictly convex quadratic programs. One of the critical steps is solving a linear system P x = q where P is a symmetric positive definite matrix. The matrix is being factorized using its Cholesky decomposition (originally using dpotrf from lapack, now replaced with cholesky from stdlib). The linear solve currently makes an explicit call to dpotrs from lapack but, for the sake of code clarity, I'd like to replace it with an equivalent chol_solve which does not currently exist in stdlib.

As of today, stdlib has a solve_lu subroutine whose interface reads

  interface solve_lu        
     #:for nd,ndsuf,nde in ALL_RHS
     #:for rk,rt,ri in RC_KINDS_TYPES
     pure module subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$(a,b,x,pivot,overwrite_a,err)     
         !> Input matrix a[n,n]
         ${rt}$, intent(inout), target :: a(:,:)
         !> Right hand side vector or array, b[n] or b[n,nrhs]
         ${rt}$, intent(in) :: b${nd}$
         !> Result array/matrix x[n] or x[n,nrhs]     
         ${rt}$, intent(inout), contiguous, target :: x${nd}$
         !> [optional] Storage array for the diagonal pivot indices
         integer(ilp), optional, intent(inout), target :: pivot(:)
         !> [optional] Can A data be overwritten and destroyed?
         logical(lk), optional, intent(in) :: overwrite_a
         !> [optional] state return flag. On error if not requested, the code will stop
         type(linalg_state_type), optional, intent(out) :: err
     end subroutine stdlib_linalg_${ri}$_solve_lu_${ndsuf}$
     #:endfor
     #:endfor    
  end interface solve_lu     

We could write an equivalent routine for chol_solve (or solve_chol if we want to be consistent with the LU call). Almost all of the necessary routines are already present in stdlib_lapack_solve_chol.fypp and stdlib_lapack_solve_chol_comp.fypp. Solving a linear system with a symmetric positive definite matrix could then be done using two successives calls, e.g. (omitting optional arguments)

!> Factorize the matrix (in-place).
call cholesky(P)
!> Solve the linear system.
call chol_solve(P, q, x)

Prior Art

• scipy.linalg has the function cho_solve whose interface reads cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True) (see here).
• In Julia, you would do F = cholesky(P) followed by x = F \ b where F is the equivalent of a Fortran derived-type specialized for the Cholesky factorization (see Derived type for matrix factorizations #1040 for a discussion with the QR example).

Additional Information

ping: @perazz, @jvdp1, @jalvesz

[manvendra5345]

Hello! I’m interested in contributing to this project as part of GSoC.
Could you please assign this issue to me? I’d love to work on it and get started.
Thanks!

[loiseaujc]

Be my guest. Pick on of the two issues you're interested, fork the repo, create your branch, implement the stuff and send a PR. We'll iterate from there to polish your code.

[manvendra5345]

Sure sir I will start working on it.

[aamrindersingh]

Hi @loiseaujc
I have implemented solve_chol and cholesky_solve for this. PR is ready for review.
Thanks

[Nitin962dev]

linalg: add solve_chol for SPD matrices

Add a new solve_chol routine to solve P x = b when P is symmetric
positive definite and already factorized using cholesky.

This avoids using LU decomposition (which is slower and needs
pivoting) and instead performs two simple triangular solves.
It makes the code cleaner and better suited for SPD systems.

linalg: add solve_chol for SPD systems

Provide a Cholesky-based solver for P x = b.
Cleaner and more efficient than LU for SPD matrices.