Add Support for Sparse Matrices in quadform()

src/reformulations/quadform.jl

@@ -121,10 +121,24 @@ function quadform(x::AbstractExpr, A::Value; assume_psd = false)
     else
         error("Quadratic forms supported only for semidefinite matrices")
     end
-    P = sqrt(LinearAlgebra.Hermitian(factor * A))
+    # Caluclates matrix `P` such that `A = P' * P`
+    P = _square_root(factor * A)
     return factor * square(norm2(P * x))
 end
 
+_square_root(A::Value) = sqrt(LinearAlgebra.Hermitian(A))
+
+function _squareroot(A::SparseArrays.SparseMatrixCSC; shift = 1e-10)
+    chol = cholesky(A; shift, check = false)
+    if !issuccess(chol)
+        error(
+            "The Shifted Cholesky decomposition failed. The matrix may not be an SPSD.",
+        )
+    end
+    L = sparse(chol.L)
+    return L'[:, invperm(chol.p)]
+end
+
 # These `evaluate` calls are safe since they are not a `fix!`ed variable
 # (must be a constant):
 function quadform(x::Constant, A::AbstractExpr; kwargs...)