Support for N-dimensional variables

[nicomignoni]

Hi all!
I noticed that cvxpy added support for n-dimensional variables relatively recently (v1.6.0), and I was wondering if something similar can be implemented in Convex.jl as well. There are a lot of applications where high-order tensors are pretty common (transports, energy systems, and warehousing are the ones I'm familiar with), so I think it would be a convenient addition.

I started going through the code, and I plan to begin tinkering with it myself, but I opened this issue for possible feedback on things that should require some attention. The first ones that pop to mind are:

• limiting positive semidefiniteness to matrices
• handling transposition, reshape, etc.

[odow]

I won't dissuade you from tinkering with the code to try and add this, and I will happily review a PR if you make one 😄

With that being said...

This will probably be much more complicated than you expect. Almost all of the code is written with the explicit assumption that atoms are a scalar, a vector, or a 2-dimensional matrix.

One issue is that Convex.jl was originally written in the early days of Julia. We've migrated the code to MOI and updated it, but the fundamental design decisions haven't changed, and are very inspired by CVX and Matlab (hence, everything is a matrix).

It's a lot more work, but I would be in favor of exploring a complete rewrite of Convex.jl that starts afresh and uses more modern Julia coding styles. It would be an experiment, and I don't know it it would be worth replacing Convex.jl with it, but we wouldn't know without trying.

The other approach would be to add more DCP-type features to JuMP, which supports n-dimenional variables.

[nicomignoni]

In my (still limited) experience with Julia, JuMP has all the features I need for optimization (convex and non-), specifically n-dimensional variables, easily enforceable sparsity, and index-based expressions.

However, what makes me go back to Convex.jl now and then is the ease of use of atomic functions. For instance, if I need to put a norm in the objective, adding a slack variable and using SecondOrderCone() is not a hassle. However, if part of the objective comes from an external function, being able to just write something like g(x) -> norm(x) would be nice. This is almost always the case in iterative methods for variational inequalities, something I've been working on lately.

I don't know what would be easier: adding a DCP side to JuMP (i.e., bringing some of Convex into JuMP) or writing a smaller JuMP-like package just for DCP (i.e., bringing some of JuMP into Convex).

[odow]

I don't know what would be easier: adding a DCP side to JuMP (i.e., bringing some of Convex into JuMP) or writing a smaller JuMP-like package just for DCP (i.e., bringing some of JuMP into Convex).

Me neither 😆

[langestefan]

The other approach would be to add more DCP-type features to JuMP, which supports n-dimenional variables.

I would be in favour of this approach.

I guess, what I would like to have is a Julia version of this: https://www.cvxpy.org/api_reference/cvxpy.reductions.back_end.html#dcp2cone

And I think the Julia version would export a set of MOI.bridges, which would then hook into the specific solver bridges? (I am just a user of these tools so I am a bit out of my depth here, if it makes no sense just say so)

But I completely understand this would be a ton of work and require quite some expertise that I currently don't have. Maybe the Stanford / cvxpy people would be able to help.

[ericphanson]

JuMP/MOI can do reformulations without an issue, the tricky bit is that these reformulations are only valid when the problem is DCP, see https://jump.dev/Convex.jl/stable/introduction/dcp/#An-invalid-formulation for a quick example of why. So the harder part is how to know when the reformulations can be used. The actual reformulations are well-understood and prominent in many JuMP tutorials (as well as in the Convex.jl source code).

Convex.jl, which was originally written by the CVX group, is an operator-overloading DSL that tracks DCPness throughout formulation and only supports DCP problems. So Convex always applies the reformulations and warns/errors if the DCP condition is violated.

JuMP does not track DCPness and moreover supports non-convex problems.

Both JuMP and Convex are frontends to MathOptInterface (MOI). Both of their variables lower to MOI variables and both populate a MOI model to be solved.

So the options as I understand it are:

• generalize (or rewrite) Convex to cover more cases, so you can use it as the MOI-front in more cases. That’s this issue to some extent
• add DCP tracking to the JuMP frontend. This needs to be non-intrusive for users who don’t need it, but accurate for when it is needed. Then add the reformulations as conditional MOI bridges (currently they are implemented on the Convex level rather than the MOI level)
• generalize the Convex.jl Optimizer interface. I don’t think this is best approach as it I believe requires rewriting the problem onto Convex variables in order to be lowered again back to MOI. This seems hard to do in the non-scalar case and indirect.

[langestefan]

Seems like SciML is working on this already: https://github.com/SciML/SymbolicAnalysis.jl

[odow]

Seems like SciML is working on this already: https://github.com/SciML/SymbolicAnalysis.jl

That's a long way from being a usable replacement to Convex.jl. There is definitely room to explore multiple approaches, so I wouldn't use the existence of that repo as a reason not to pursue improvements to Convex.jl, or a completely new approach again.

[langestefan]

Were you aware of it already? Would be nice to align some efforts where possible.