Avoid eigendecomposition in check for positive semidefiniteness

[nrontsis]

Uses pivoted Cholesky instead of eigendecomposition for checking positive/negative semidefiniteness. Advantages of the approach include that Cholesky is faster than eigen-decomposition, scales much better with parallelism and can exploit sparsity.

The positive semidefiniteness checks can be become computationally relevant when using approximate projections, as, in this case, we are trying to avoid big eigendecompositions as much as possible.

References [1] [2] [3]

[codecov-io]

Codecov Report

Merging #94 into master will increase coverage by 0.34%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master      #94      +/-   ##
==========================================
+ Coverage   89.89%   90.23%   +0.34%     
==========================================
  Files          22       22              
  Lines        1712     1803      +91     
==========================================
+ Hits         1539     1627      +88     
- Misses        173      176       +3

Impacted Files	Coverage Δ
src/algebra.jl	98.83% <100%> (+0.11%)	⬆️
src/convexset.jl	94.25% <100%> (-0.52%)	⬇️
src/infeasibility.jl	100% <100%> (ø)	⬆️
src/solver.jl	98.66% <0%> (+0.01%)	⬆️
src/chordal_decomposition.jl	93.63% <0%> (+1.83%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 340ecde...4b8396d. Read the comment docs.

[nrontsis]

After reading LAPACK Working Note #161 (Section 7) I think that using pivoting in Cholesky does not help with determining indefiniteness. Pivoting is useful when we know (a priori) that the matrix is positive semidefinite, and we want e.g. to estimate its rank.

The non-pivoted Cholesky decomposition fails "iff the matrix is positive definite". I think this covers our needs because we always check definiteness with tolerances much higher than machine precision. However, to better emphasise the fact that we cannot (always) check if a matrix is positive semidefinite, I renamed the functions to:

is_pos_sem_def -> is_pos_def
is_neg_sem_def -> is_neg_def

[migarstka]

Thanks Nikitas. I squashed your commits into one and added some changes.

• I got an error for the sparse version of is_pos_def, saying Unsupported keyword split=. For now the method is only used in the infeasibility checks anyway and the input variable will be dense. I removed the sparse version for now.

• I added an in-place version of is_neg_def

• Your proposed methods pass some tests that I wrote:

using BenchmarkTools, Random, SparseArrays, LinearAlgebra, COSMO, Test


rng = Random.MersenneTwister(213213);

# utility function to generate matrix with specific eigenvalues
 function generate_matrix(dim::Int64, rng, eigs)
      X = rand(rng, dim, dim)
      Q, R = qr(X)
      X = Q*Diagonal(eigs)*Q'
      X = 0.5*(X+X')
      return X
  end

# current implementation
function is_pos_sem_def_old!(X, tol)
    # set option 'N' to only compute eigenvalues, s is ordered from min to max
   s, U = LAPACK.syevr!('N', 'A', 'U', X, 0.0, 0.0, 0, 0, -1.0);
   return s[1] >= -tol
end

dim = 100
tol = 1e-4

######
# Unit tests
######

# negative definite matrix
eigs = -100 .+ rand(rng, dim) .* (-1 - (-100));
Xn = generate_matrix(dim, rng, eigs);

@test COSMO.is_pos_def(Xn, tol) == false

# positive definite matrix
@. eigs *= -1;
Xp = generate_matrix(dim, rng, eigs);
@test COSMO.is_pos_def(Xp, tol) == true

# positive semi-definite matrix (edge case negative)
eigs[1] = -1e-5;
Xe1 = generate_matrix(dim, rng, eigs);
@test COSMO.is_pos_def(Xe1, tol) == true

• Furthermore, your version is a lot faster and the in-place version allocates less memory than the current one:

dim = 400
Xtest = generate_matrix(dim, rng, rand(rng, dim));
Xtest2 = copy(Xtest);
Xtest3 = copy(Xtest);

# with copy
println("With copy of matrix:")
@btime COSMO.is_pos_def($Xtest, 1e-4)

# without allocation
println("In place version:")
@btime COSMO.is_pos_def!($Xtest2, 1e-4)
@btime is_pos_sem_def_old!($Xtest3, 1e-4)

Benchmarking yields the following results:

  617.922 μs (6 allocations: 1.22 MiB) ( for COSMO.is_pos_def ) 🥈
  500.471 ns (4 allocations: 112 bytes) ( for COSMO.is_pos_def! ) 🥇
  6.431 ms (10 allocations: 147.75 KiB) ( for current implementation ) 🥉

• Currently, we are using non-modifying versions of the functions is_pos_def and is_neg_def and they are only used in the infeasibility detection routines:
  • is_primal_infeasible -> support_function -> in_dual -> is_pos_def
  • is_dual_infeasible -> in_pol_recc -> is_neg_def
    However, since the functions operate on temporary variables unrelated to the ADMM loop iterates, we can use the faster in-place versions is_pos_def! and is_neg_def. I will do that in a separate commit.

[nrontsis]

Great, thanks for the detailed feedback!