Quad Form Fails with Sparse Matrix

[RoyiAvital]

It seems that Convex.quadform() fails with sparse matrices even with assume_psd = true:

using ECOS;
using SparseArrays;

function SolveQuadForm( mP :: AbstractMatrix{T}, vQ :: Vector{T} ) where {T <: AbstractFloat}

    dataDim = size(mP, 1);

    vX = Convex.Variable(dataDim);
    sConvProb = minimize( T(0.5) * Convex.quadform(vX, mP; assume_psd = true) + Convex.dot(vQ, vX) );
    solve!(sConvProb, ECOS.Optimizer; silent = true);
    
    return vec(vX.value);
    
end

dataDim = 10;

mP = sprand(dataDim, dataDim, 0.5);
mP = (mP' * mP) + 0.5 * I;
mP = 0.5 * (mP + mP');
vQ = randn(dataDim);

vX = SolveQuadForm(mP, vQ);

I get: ERROR: MethodError: no method matching eigen!(::Hermitian{Float64, SparseMatrixCSC{Float64, Int64}}; sortby::Nothing).
Which is from:

Convex.jl/src/reformulations/quadform.jl

Line 124 in 0843b60

P = sqrt(LinearAlgebra.Hermitian(factor * A))

Speculating it tries to do Eigen Decomposition for the root of the matrix which fails for sparse matrix.
Any idea about a temporary bypass?

[RoyiAvital]

Current bypass:

function SolveQuadForm( mP :: AbstractMatrix{T}, vQ :: Vector{T} ) where {T <: AbstractFloat}

    dataDim = size(mP, 1);

    sC = cholesky(mP; check = false, perm = 1:dataDim);
    mU = sparse(sparse(sC.L)');

    vX = Convex.Variable(dataDim);
    # sConvProb = minimize( T(0.5) * Convex.quadform(vX, mP; assume_psd = true) + Convex.dot(vQ, vX) );
    sConvProb = minimize( T(0.5) * Convex.square(Convex.norm2(mU * vX)) + Convex.dot(vQ, vX) );
    solve!(sConvProb, ECOS.Optimizer; silent = true);
    
    return vec(vX.value);
    
end

Works only if mP is Symmetric Positive Definite (SPD).

[odow]

The MethodError isn't coming from Convex, so your only option is to convert the matrix into a dense matrix.

[RoyiAvital]

What do you mean?
Convex.jl tries to compute the matrix square root which is not feasible for sparse matrices unless they are SPD.
Or do I miss something?

[odow]

Use Matrix(mP):

julia> using Convex, ECOS, LinearAlgebra, SparseArrays

julia> function SolveQuadForm(mP::AbstractMatrix{T}, vQ::Vector{T}) where {T<:AbstractFloat}
           dataDim = size(mP, 1)
           vX = Convex.Variable(dataDim)
           sConvProb = minimize(T(0.5) * quadform(vX, mP; assume_psd = true) + dot(vQ, vX))
           solve!(sConvProb, ECOS.Optimizer; silent = true)
           return vec(vX.value)
       end
SolveQuadForm (generic function with 1 method)

julia> dataDim = 10
10

julia> mP = sprand(dataDim, dataDim, 0.5)
10×10 SparseMatrixCSC{Float64, Int64} with 46 stored entries:
  ⋅         ⋅          ⋅         ⋅           ⋅        …  0.246577   ⋅         0.59069   0.0372561
  ⋅         ⋅         0.356732  0.321219    0.773096      ⋅         ⋅          ⋅        0.208567
 0.168875  0.295742    ⋅         ⋅           ⋅            ⋅         ⋅         0.879747  0.276595
 0.432293   ⋅         0.940476  0.00724681   ⋅            ⋅        0.360426   0.52251   0.285745
 0.840612  0.0522051  0.793347  0.345789     ⋅            ⋅        0.375546    ⋅        0.168664
  ⋅        0.285009   0.512001   ⋅           ⋅        …   ⋅        0.519182   0.236298  0.68705
  ⋅         ⋅          ⋅        0.954531     ⋅            ⋅        0.0525299   ⋅        0.824198
  ⋅         ⋅         0.681335   ⋅          0.727378      ⋅         ⋅         0.831776  0.845443
 0.853151   ⋅          ⋅         ⋅           ⋅            ⋅         ⋅         0.685721   ⋅ 
  ⋅         ⋅         0.189244  0.697526    0.696164      ⋅        0.495016    ⋅        0.159639

julia> mP = (mP' * mP) + 0.5 * I
10×10 SparseMatrixCSC{Float64, Int64} with 80 stored entries:
 2.14989    0.0938275  1.07346   0.293807     ⋅        …   ⋅         0.471498  0.959468    0.312017
 0.0938275  0.671418   0.187342  0.0180519    ⋅            ⋅         0.167577  0.327525    0.286421
 1.07346    0.187342   2.90333   0.527738    0.903122      ⋅         0.996411  1.17911     1.43496
 0.293807   0.0180519  0.527738  2.12048     0.733925      ⋅         0.5279    0.00378653  1.02546
  ⋅          ⋅         0.903122  0.733925    2.1114        ⋅         0.344612  0.605015    0.887333
 0.403849    ⋅         1.40071   0.827486    0.597638  …   ⋅         0.380063  1.08086     1.57101
  ⋅          ⋅          ⋅         ⋅           ⋅           0.5608      ⋅        0.145651    0.0091865
 0.471498   0.167577   0.996411  0.5279      0.344612      ⋅         1.28829   0.311008    0.645354
 0.959468   0.327525   1.17911   0.00378653  0.605015     0.145651   0.311008  3.11379     1.28021
 0.312017   0.286421   1.43496   1.02546     0.887333     0.0091865  0.645354  1.28021     2.62309

julia> mP = 0.5 * (mP + mP')
10×10 SparseMatrixCSC{Float64, Int64} with 80 stored entries:
 2.14989    0.0938275  1.07346   0.293807     ⋅        …   ⋅         0.471498  0.959468    0.312017
 0.0938275  0.671418   0.187342  0.0180519    ⋅            ⋅         0.167577  0.327525    0.286421
 1.07346    0.187342   2.90333   0.527738    0.903122      ⋅         0.996411  1.17911     1.43496
 0.293807   0.0180519  0.527738  2.12048     0.733925      ⋅         0.5279    0.00378653  1.02546
  ⋅          ⋅         0.903122  0.733925    2.1114        ⋅         0.344612  0.605015    0.887333
 0.403849    ⋅         1.40071   0.827486    0.597638  …   ⋅         0.380063  1.08086     1.57101
  ⋅          ⋅          ⋅         ⋅           ⋅           0.5608      ⋅        0.145651    0.0091865
 0.471498   0.167577   0.996411  0.5279      0.344612      ⋅         1.28829   0.311008    0.645354
 0.959468   0.327525   1.17911   0.00378653  0.605015     0.145651   0.311008  3.11379     1.28021
 0.312017   0.286421   1.43496   1.02546     0.887333     0.0091865  0.645354  1.28021     2.62309

julia> vQ = randn(dataDim)
10-element Vector{Float64}:
 -0.44713347675937803
 -0.12165519515276513
  1.1691996905022073
  1.013833594370483
  1.0678462387137522
 -0.8630538619227617
  1.2072577636715998
 -0.5999336039316524
  0.33712681339992073
 -0.054263156946659535

julia> # vX = SolveQuadForm(mP, vQ)
       vX = SolveQuadForm(Matrix(mP), vQ)
10-element Vector{Float64}:
  0.6904015459888632
  0.20526385720942314
 -1.5229487620225548
 -1.1499484068228385
 -0.06364484978257556
  1.1622043558074542
 -2.068957917672532
  1.4168613491392508
 -0.3446816727614091
  0.35118489110979656

[RoyiAvital]

I know it works with a matrix.
I am reporting about the sparse case.

The issue the code of Convex.jl uses Julia's sqrtm() which does not have a path for sparse matrix.
I think there is no such case in general. So maybe there should be a special path for sparse matrices in quadform()?
Wouldn't Pivoted Cholesky work?

[odow]

Not sure what the best approach is. PRs accepted.

[RoyiAvital]

Would you be OK with making the current method available with _quadform() and a PR for the Pivoted Cholesky?

[odow]

Would you be OK with making the current method available with _quadform()

No, don't break user code.

Don't you just need a new method for sparse arrays?

function quadform(x::AbstractExpr, A::SparseArrays.SparseMatrixCSC; assume_psd = false)
    ....
end

[RoyiAvital]

I thought about something that should benefit (Speed) work for both the dense and sparse cases.
I can do a path only for sparse as well. What would you prefer?

[odow]

I don't have a preference so long as that any code that currently calls quadform continues to work.