Add Support for Sparse Matrices in quadform()

src/reformulations/quadform.jl

@@ -121,10 +121,19 @@ 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))
+    # Calculates 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 _square_root(A::SparseArrays.SparseMatrixCSC)
+    chol = LinearAlgebra.cholesky(A; shift = 1e-10)
+    L = SparseArrays.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...)

test/test_atoms.jl

@@ -12,6 +12,7 @@ using Test
 import LinearAlgebra
 
 import MathOptInterface as MOI
+import SparseArrays
 
 function runtests()
     for name in names(@__MODULE__; all = true)
@@ -2459,6 +2460,15 @@ function test_quadform()
         quadform(Variable(2), [1 0; -2 1]),
     )
     @test quadform(constant([1, 2]), constant([1 2; 2 3])) == 21
+    target = """
+    variables: u, t, x1, x2
+    minobjective: 1.0 * u
+    [t, 2.000000000025 * x1, 3.0000000000166667*x2] in SecondOrderCone(3)
+    [u, 0.5, t] in RotatedSecondOrderCone(3)
+    """
+    _test_reformulation(target) do context
+        return quadform(Variable(2), A)
+    end
     return
 end