Support for complex vector cones

[araujoms]

It would be nice to have the complex version of the vector cones defined in MOI, so that we have we nice a interface for the solvers that can handle them.

Specifically, the one I need is the infinity norm cone, but I think it also makes sense to add the 1-norm cone, as that is also directly supported by Hypatia, and the 2-norm cones, as they are so easy to do.

For solvers that cannot handle the cones directly maybe it's useful to add a bridge? One can rewrite t >= norm(x,Inf) in a horrendous way as

for i=1:lenght(x)
    [t, real(x[i]), imag(x[i])] in SecondOrderCone()
end

and t >= norm(x,1) in an even uglier way as

t >= sum(u)
for i=1:lenght(x)
    [u[i], real(x[i]), imag(x[i])] in SecondOrderCone()
end

where u is a vector of optimization variables of the same length as x.

For the 2-norm cones I think the only thing that makes sense is to add a bridge, as t >= norm(x) can be nicely written as [t;real.(x);imag.(x)] in SecondOrderCone(), and 2*t*u >= norm(x)^2 as [t;u;real.(x);imag.(x)] in RotatedSecondOrderCone()

[odow]

What is preventing Hypatia or a bridge from supporting VectorAffineFunction{ComplexF64}-in-NormOneCone at the moment?

Why do we need to define complex cones?

Perhaps this issue is an ask for a bridge?

[araujoms]

I don't know, I don't understand the architecture of MOI. I guessed that one needed to define the complex cone in MOI because it has both PositiveSemidefiniteConeTriangle and HermitianPositiveSemidefiniteConeTriangle. To be clear, I am already using the complex infinity norm cone in Hypatia; it's just that I'm doing it through a rather inconvenient low-level interface, and I'd like to be able to do simply

@variable(model, t)
@variable(model, x[1:d] in ComplexPlane())
@constraint(model, [t;x] in NormInfinityCone())

instead.

[odow]

@blegat should chime in here. He has more thoughts on the design of the Complex bridges.

[blegat]

In MOI, we assume that the components of cones are real because some bridges operations wouldn't be valid for complex numbers (e.g., SlackBridge creates real numbers). About the SOC cone, I assume you mean Hypatia.Cones.EpiNormEucl which is t >= |x|_2 where t is real and x is a vector of n complex numbers. I would create I complex version ComplexSecondOrderCone that has 2n + 1 entries, the first being t, then the n real parts and then the n imaginary parts. The, from JuMP, @constraint(model, v in SecondOrderCone()) would do the right thing depending on whether eltype(v) is an JuMP expression of real or complex coefficient type. For nonlinear expressions, it might be a bit tricky to determine whether it's real or complex though :thinking:

[araujoms]

I see, so the only way to get a nice interface is to define the complex cones in JuMP.

For the 2-norm cone it doesn't seem that a MOI version is necessary, as JuMP could directly translate it into a real 2-norm cone with 2n + 1 entries. Note that Hypatia doesn't support a complex 2-norm cone, presumably because there's no point.

For the other cones, though, a MOI layer is necessary (as far as I understand), as there is no direct translation into the real versions.

As for determining whether the expressions are real or complex, isn't it exactly what @odow did in jump-dev/JuMP.jl#3695?

[blegat]

As for determining whether the expressions are real or complex, isn't it exactly what @odow did in https://github.com/jump-dev/JuMP.jl/pull/3695?

Yes, you can see that it's considering that nonlinear expressions may be complex.

For the other cones, though, a MOI layer is necessary (as far as I understand), as there is no direct translation into the real versions.

What do you mean by "no direct translation" ? Can you give an example ?

[araujoms]

Well it's what I wrote in the first comment; if x is a complex vector, there's no way of rewriting t >= norm(x,Inf) in terms of the real infinity norm cone. At least as far as I can see. The only translation I could come with was that thing in terms of length(x) 2-norm cones.

[blegat]

So you could create I complex version ComplexNormInfinityCone that has 2n + 1 entries, the first being t, then the n real parts and then the n imaginary parts. The, from JuMP, @constraint(model, v in NormInfinityCone()) could try to do the right thing depending on whether eltype(v) is an JuMP expression of real or complex coefficient type (again, unclear for nonlinear expressions). For solvers not supporting ComplexNormInfinityCone, we could rewrite it by creating many real 2-norm cones using a bridge but Hypatia would support it directly. @chriscoey @lkapelevich have you already written such bridge ? Is there a big win in using your complex cone instead of the extended formulation ?

[araujoms]

Yep, that's it, although I think the most convenient definition would be [t, real(x[1]), imag(x[1]), real(x[2]), imag(x[2]), ...] instead of [t; real.(x); imag.(x)].

[blegat]

Why would it be more convenient ? This isn't consistent with HermitianPositiveSemidefiniteConeTriangle

[araujoms]

Because that's how Hypatia does it, so if you define it the same way no reshuffling will be necessary.

I'm not a fan of the indexing in HermitianPositiveSemidefiniteConeTriangle, it drove me crazy when I was trying to support it in SeDuMi. I believe @chriscoey was also unhappy about it.

[blegat]

For the hermitian PSD cone, we can't change before MOI v2 since that would be breaking and I don't remember getting any comment on the choice of ordering (see https://github.com/jump-dev/MathOptInterface.jl/pull/1962). About the norm infinity cone, I would find it convenient to have first all the real part and then all the complex part so that we can build it with [t; real.(x); imag.(x)] and recover each part with by taking sub-vectors but no strong opinion. Having it consistent with the Hermitian cone would be nice as well. However, we could mix like the Hypatia cone so that the bridge is trivial.

[araujoms]

I was referring to this comment and this bug. In any case I wasn't suggesting that the indexing should be changed, I was just saying it is not an example we should follow (although to be fair Hypatia is the only external package that is using the indexing, so it wouldn't really be breaking).

But as for the infinity norm cone I don't think it really makes a difference, the reshuffling is simple enough that it won't cause any bugs.

[chriscoey]

@araujoms is right, it was a surprisingly big pain actually to implement the various transformations between Hypatia's format and MOI's format. It is my fault though for not noticing that the proposed MOI format did not match the format Hypatia was already using at that point. I would love for the format to change to match Hypatia's in future (where complex off-diag entries are represented with consecutive real then imaginary parts).

[odow]

I don't really understand the problem of supporting VectorAffineFunction{ComplexF64}-in-NormOneCone. We do it for ScalarAffineFunction{ComplexF64}-in-EqualTo{ComplexF64}.

If some bridges are invalid, then that is a bug and they should be fixed to not apply.

I don't want to start adding a bunch of ComplexXXX sets.

[odow]

See https://github.com/jump-dev/MathOptInterface.jl/pull/2451

[blegat]

I don't really understand the problem of supporting VectorAffineFunction{ComplexF64}-in-NormOneCone. We do it for ScalarAffineFunction{ComplexF64}-in-EqualTo{ComplexF64}.

What if you add VectorAffineFunction{ComplexF64}-in-NormOneCone and it gets bridge by the SlackBridge which adds VectorOfVariables-in-NormOneCone and then equal these real variables to the entries of the VectorAffineFunction. This wouldn't be correct. Indeed we have the same issue with ScalarAffineFunction{ComplexF64}-in-EqualTo{ComplexF64} so maybe we could fix it for both of them. The difference though is that we see in MOI.EqualTo{ComplexF64} that this is a cone over complex numbers while NormOneCone is always a cone over real numbers. What does it mean to create constrained variables in NormOneCone ? It will be a vector of real variables. I feel having this ambiguity for these real cones might be dangerous but I might be convinced otherwise

[araujoms]

:shrug: If it is possible to support complex variables in the regular cones it would be more convenient for the users as well, I'd rather not have to switch between x in PSDCone() and x in HermitianPSDCone() depending on whether x is real or complex. Incidentally, this would also make it possible to change the indexing of the complex version without it being even in principle breaking.

Some syntax would be needed to constrain variables upon creation, perhaps @variable(model, x[1:d,1:d] in PSDCone(Complex)) and @variable(model, v[1:d] in NormOneCone(Complex))?

[odow]

and it gets bridge by the SlackBridge which

My point is that it should NOT get bridged by slack bridge, because slack cannot apply to complex functions. That is a bug.

[odow]

See https://github.com/jump-dev/MathOptInterface.jl/pull/2468 for the bug fix removing Complex support from some bridges

[blegat]

I don't really understand the is_maybe_real function added in https://github.com/jump-dev/MathOptInterface.jl/pull/2468. Don't we instead want may_be_complex ? The bridges need to have a guarantee that the function is always real right so why is_maybe_real is not really the right naming. And actually, you return true for functions for which it's always real so it's just a naming issue I guess. There is also the ambiguity of the type vs the actual value is a complex type with zero imaginary considered real ? For instance

julia> isreal(1 + im - im)
true

So may_be_complex makes more sense, it's true for functions for which, even when evaluated with real values (MOI variables are always real), their value might have an imaginary part.

[odow]

So I didn't want to break existing code (like Polynomials) for which we don't know. So we had to default to false positives.

So we'd need to define is_maybe_complex(::Any) = false or is_maybe_real(::Any) = true. We really want to disallow some bridges for which we definitely know they don't apply.

In MOI 2.0, I'd like to switch to making bridges opt-in instead of opt-out.

[blegat]

So I didn't want to break existing code (like Polynomials) for which we don't know. So we had to default to false positives.

MOI 2.0 is still far in the future so I'd rather have a cleaner solution in the meantime. If we add the function and let PolyJuMP implement it before we use the function, I don't think we go for the function that default to assuming functions may be false (which is closer to the opt-in behavior you want).

In MOI 2.0, I'd like to switch to making bridges opt-in instead of opt-out.

A bridge is always equivalent to a mathematical relation between sets so if the function satisfies a few mathematical properties (e.g., being real), there is no reason we cannot opt-in by default.

[odow]

(e.g., being real)

But this was the entire point. We have no way of knowing if an arbitrary MOI.AbstractFunction satisfies the required properties, or the required promote_op etc. At the moment, if you use an unsupported function you can hit a MethodError because of a bridge that isn't relevant. That isn't good.

[blegat]

At the moment, if you use an unsupported function you can hit a MethodError because of a bridge that isn't relevant. That isn't good.

I think it's sufficiently rare to define new MOI functions that we can have a list of functions that are required to be implemented for any AbstractFunction, including the one added in https://github.com/jump-dev/MathOptInterface.jl/pull/2468 and the other ones added in the bridges. If one of these function is not implemented, having a bridge throw an MethodError is expected behavior.

[odow]

Kinda agree. I see both arguments for MOI 2.0.

But for now, disallowing bridges unless functions implement is_maybe_complex is a breaking change, so it needs to be is_maybe_real.

[odow]

I'm also okay to make this function private for now.

[blegat]

If the default of is_maybe_complex is false then it's non-breaking (we say it's mandatory to implement it for new functions). We can add a note in the docstring that in MOI 2.0 it will default to true. So it's basically the same behavior, just another name. We can also consider that this is a bug-fix that arbitrary functions were used in these bridges that assume realness in which case we can use true as default without it being breaking.

[odow]

Is there anything left to do here?

[araujoms]

Could you please clarify what you have decided? So in the end you're not defining new complex sets, but allowing complex vectors to be added to the existing ones? So in order to get support for the complex norm one cone in Hypatia I need to add support for VectorAffineFunction{ComplexF64}-in-NormOneCone there?

[blegat]

The variables are real but the cones may have complex entries. So NormOneCone changes meaning if the function is complex or real. So Dualization.jl is now incorrect with complex vector in front, I think Dualization.jl should say that it does not support VectorAffineFunction{ComplexF64}-in-NormInfinityCone so that it's bridged into VectorAffineFunction{Float64}-in-ComplexNormInfinityCone where ComplexNormInfinityCone has real and imaginary parts are separate entries so that the dual can be real variables in that cones. So we'll still have to add complex version of the cones to be able to use them with add_constrained_variables. We can allow VectorAffineFunction{ComplexF64}-in-FooCone but if we add ComplexFooCone, we should bridge VectorAffineFunction{ComplexF64}-in-FooCone to VectorAffineFunction{Float64}-in-ComplexFooCone. Long story short, since we don't have ComplexNormOneCone, Hypatia should add support for VectorAffineFunction{ComplexF64}-in-NormOneCone

[araujoms]

Thanks for the information. I can't do this now, but I'll get around to it eventually.

[odow]

So Dualization.jl is now incorrect with complex vector in front

Yip. But this is just a straight up bug that has existed since we started allowing, for example, ScalarAffineFunction{<:Complex}-in-EqualTo{<:Complex}.

[odow]

Is there anything left to do here?

Sooo...

[odow]

Closing because I don't know if there is anything left to do here. Current suggestion is to add VectorAffineFunction{Complex}-in-Cone constraints. There might be bugs with various bridges etc because the Complex codepath is not well tested, but they can be opened and addressed as separate issues.