HermitianPositiveSemidefiniteConeTriangle with convenient indexing

[araujoms]

As we've discussed before, I find the indexing of HermitianPositiveSemidefiniteConeTriangle rather inconvenient. I'm about to write the JuMP interface for the complex PSD cone of SCS, and I'd rather not use it again (as I already did for SeDuMi and COSMO). Changing the existing set would break backward compatibility, though, so I suggest adding a new set, say ComplexPositiveSemidefiniteConeTriangle, and deprecating HermitianPositiveSemidefiniteConeTriangle. What do you think? I can do it myself if you're interested, I just don't want to waste my time writing a PR that won't be welcomed.

[odow]

I don't really think we need yet-another-cone. I understand that the indexing is a pain, but it's a couple of lines of code in a couple of places. To me, that doesn't justify adding another cone. @blegat needs to decide this though.

I'm about to write the JuMP interface for the complex PSD cone of SCS

I don't understand this part.

[araujoms]

I don't understand this part.

I contributed a complex PSD cone to the solver SCS. I'm also going to make it accessible through MOI, by adding the appropriate bridge to the SCS.jl. Clear now?

[odow]

Are you talking about this? cvxgrp/scs#301

Since that PR is not merged yet, you could always choose to use MOI's vectorization 😄

At minimum, we should hold off on any decision here until there SCS PR is actually merged and a new version is tagged...

[araujoms]

Yes, that's the PR. There's no way I'm going to use the MOI vectorization in SCS. Not only it's a pain in the ass to program, it's also rather inefficient due to keeping the real and imaginary parts far apart. LAPACK requires a complex matrix with the real and imaginary parts interleaved, so the vectorization needs to match. And the efficiency does matter, as the vectorization needs to be done hundreds, if not thousands of times.

So no, nothing depends on the PR being merged. Also, I have time to program now that I'm on vacations, after I get back to work it will be harder.

[odow]

My smile was slightly sarcastic 😄

The underlying point is that the user's input format does not need to match the solver's internal computation format. It's just a convention with a re-order.

[araujoms]

I know that. Note that SCS does not use Hypatia's format, because it is written in C, and thus uses a vectorization that's more appropriate for a row-major language, whereas Hypatia's is more appropriate for a col-major language. I'm proposing here to use Hypatia's format in MOI.

What I'm talking about is how convenient it is to write the bridge. I think it's worth rewriting the cone in order to spare myself the pain of writing yet another bridge with the MOI vectorization. It's also a gift for the future programmers.

[odow]

But if we add a cone that uses Hypatia's format, how does that help SCS?

[blegat]

If I understand correctly, you would like the chain of bridges:
MOI.HermitianPSDTriangle -> ComplexPSDTriangle -> MOI.Scaled{ComplexPSDTriangle} -> SCS.ScaledRowMajorComplexPSD.
I agree it's a clean way to do it and you can start implementing it in SCS.jl like:
https://github.com/jump-dev/SCS.jl/blob/d5b0b09758be35fa3193dcf478d44508ca75bbbd/src/MOI_wrapper/scaled_psd_cone_bridge.jl
by adding two cones : ComplexPSDTriangle and ScaledRowMajorComplexPSD.
Once it is added in SCS and working, it's easier to have the discussion on adding it to MOI.
The issue with adding things to MOI is that MOI has reached v1 so you better think twice before merging a PR (HermitianPSDTriangle is a good example, we first started as an extension https://github.com/jump-dev/ComplexOptInterface.jl and then merged it to MOI but even then it was not enough to have your useful feedback about the ordering of real and complex part ;) ).
So I wouldn't be against adding it, if HermitianPSD was a bad decision, it's better to fix it now that we don't have too many solvers using this cone than later but I advice you to start with everything in a SCS.jl PR.

I know that. Note that SCS does not use Hypatia's format, because it is written in C, and thus uses a vectorization that's more appropriate for a row-major language, whereas Hypatia's is more appropriate for a col-major language. I'm proposing here to use Hypatia's format in MOI.

Why is row-major more appropriate in C ? You have an option to choose between row-major and col-major in lapack if I'm not mistaken. Julia and SCS just use different options. Note that col-major indexing is much easier as the mapping from (i, j) to the vectorized index does not depend on the dimension of the matrix which is the reason it was chosen.

[araujoms]

If I have to write the bridge from MOI.HermitianPSDTriangle anyway I'd rather just write MOI.Scaled{MOI.HermitianPSDTriangle} -> SCS.ScaledRowMajorComplexPSD directly and be done with it. The goal is less work, not more. I understand your reticence about changing MOI, but I don't see what could possible change, there is no other indexing under consideration.

Why is row-major more appropriate in C ? You have an option to choose between row-major and col-major in lapack if I'm not mistaken.

Because C is row-major. I don't know about this LAPACK option, and it doesn't sound like a good idea anyway, as LAPACK is written in FORTRAN, which is col-major, and telling it to work in row-major format would be rather inefficient.

Julia and SCS just use different options. Note that col-major indexing is much easier as the mapping from (i, j) to the vectorized index does not depend on the dimension of the matrix which is the reason it was chosen.

You can also have this nice property with a row-major ordering by taking the rows of the lower diagonal, instead of the upper one. It's not what SCS does it, though, and I can't fix that without breaking backwards compatibility. So I kept taking the rows of the upper diagonal.

[blegat]

If I have to write the bridge from MOI.HermitianPSDTriangle anyway I'd rather just write MOI.Scaled{MOI.HermitianPSDTriangle} -> SCS.ScaledRowMajorComplexPSD directly and be done with it. The goal is less work, not more.

I don't understand, you want the ComplexPSDCone to be col-major so it will be different to SCS so you'll still need to write a bridge anyway. Your suggestion is precisely to add that ComplexPSD cone to MOI and a bridge from the Hermitian one. I'm just suggesting to do exactly what you suggested but first in a SCS PR

I understand your reticence about changing MOI, but I don't see what could possible change, there is no other indexing under consideration

I'm just saying in general that you can always define your cone and bridge outside of MOI and then do a PR to MOI once we see that several solvers would benefit from it. Here it seems to be the case for Hypatia and SCS so I don't have any objection but it's still easier to start in an SCS PR.

You can also have this nice property with a row-major ordering by taking the rows of the lower diagonal, instead of the upper one. It's not what SCS does it, though, and I can't fix that without breaking backwards compatibility. So I kept taking the rows of the upper diagonal

I don't understand. MOI is both col-major upper-triangular and row-major lower- triangular. There two are the same by symmetry of the matrix and these are the ones with that property that the index mapping does not depend on the dimension. SCS is the other one so it does not have that property.

Because C is row-major

What does that mean ? C has no concept of matrix.

[araujoms]

I don't understand, you want the ComplexPSDCone to be col-major so it will be different to SCS so you'll still need to write a bridge anyway. Your suggestion is precisely to add that ComplexPSD cone to MOI and a bridge from the Hermitian one. I'm just suggesting to do exactly what you suggested but first in a SCS PR

Yes, but the bridge I want to write is MOI.Scaled{MOI.ComplexPSDTriangle} -> SCS.ScaledRowMajorComplexPSD, which would make a bridge MOI.HermitianPSDTriangle -> SCS.ComplexPSDTriangle unnecessary.

I don't understand. MOI is both col-major upper-triangular and row-major lower- triangular. There two are the same by symmetry of the matrix and these are the ones with that property that the index mapping does not depend on the dimension. SCS is the other one so it does not have that property.

The underlying matrix is not necessarily Hermitian (and in fact never is when working with LAPACK), so col-major upper-triangular and row-major lower-triangular are not the same thing.

What does that mean ? C has no concept of matrix.

If you declare int matrix[2][2] = {{1, 2}, {3, 4}}; and read the matrix elements in the order they are in memory, you'll get 1,2,3,4. For this reason even when working with one-dimensional arrays and doing index-juggling manually C programs are usually row-major.

[blegat]

The underlying matrix is not necessarily Hermitian (and in fact never is when working with LAPACK), so col-major upper-triangular and row-major lower-triangular are not the same thing.

At the MOI level, we only store the triangle part for the ...Triangle sets so it can only be symmetry/hermitian. Col-major upper and row-major lower are therefore the same. The solver can interpret it either way.

Yes, but the bridge I want to write is MOI.Scaled{MOI.ComplexPSDTriangle} -> SCS.ScaledRowMajorComplexPSD, which would make a bridge MOI.HermitianPSDTriangle -> SCS.ComplexPSDTriangle unnecessary.

JuMP still creates MOI.HermitianPSDTriangle so we'll have to bridge it to MOI.ComplexPSDTriangle at some point. If we decide to deprecate MOI.HermitianPSDTriangle then JuMP will use MOI.ComplexPSDTriangle and we'll need a bridge from MOI.ComplexPSDTriangle to MOI.HermitianPSDTriangle.

But you can indeed start with the bridge fromMOI.Scaled{MOI.ComplexPSDTriangle} to SCS.ScaledRowMajorComplexPSD and should should still be able to use it from JuMP by doing @constraint(model, x in ComplexPSDTriangle(n)).

[araujoms]

At the MOI level, we only store the triangle part for the ...Triangle sets so it can only be symmetry/hermitian. Col-major upper and row-major lower are therefore the same. The solver can interpret it either way.

Well, yes, except for complex conjugation. Btw @chriscoey why are you taking the complex conjugate of the upper triangular in Hypatia?

JuMP still creates MOI.HermitianPSDTriangle so we'll have to bridge it to MOI.ComplexPSDTriangle at some point. If we decide to deprecate MOI.HermitianPSDTriangle then JuMP will use MOI.ComplexPSDTriangle and we'll need a bridge from MOI.ComplexPSDTriangle to MOI.HermitianPSDTriangle.

But you can indeed start with the bridge fromMOI.Scaled{MOI.ComplexPSDTriangle} to SCS.ScaledRowMajorComplexPSD and should should still be able to use it from JuMP by doing @constraint(model, x in ComplexPSDTriangle(n)).

Fair enough, had forgotten about JuMP; I'll write the chain of bridges then.