Support for vector dual

[yemine0x01]

For constraints such as second-order cones or the exponential cone, the dual variable is a vector. Currently, however, the return type for the dual value is always a double. Would it be useful to add support for returning vector-valued duals in these cases?

[metab0t]

The dual variable should have the same length as the dimension of cone ($N$ for SOC and 3 for exponential cone), but only solvers with the genuine conic formulation supports it (like Mosek and Clarabel). For example, SOC is treated as a single quadratic constraint in Gurobi.

If we decide to adopt vector dual for conic constraint, it is important to implement the transformation from the scalar dual in some solvers.

[yemine0x01]

I agree. I know how to do it on KNITRO and Gurobi, but unfortunately, I don't know how to do it for other solvers. Can we, for example, raise that it is not supported on them?

[metab0t]

Yes, it would be great if you can provide implementation on KNITRO and Gurobi.

I think IPOPT can be done similarly like Gurobi because it also treats conic constraint like a single constraint.

[metab0t]

For $\mu$ as the scalar dual for SOC constraint $t^2 - x_1^2 - x_2^2 \geq 0$, the dual vector can be written as:

$\lambda_\mu = 2 \mu t$
$\lambda_{x_1} = -2 \mu x_1$
$\lambda_{x_2} = -2 \mu x_2$

I am not sure whether my understanding is correct.

A small example shows it holds.

Python part:

import pyoptinterface as poi
from pyoptinterface import gurobi


def test_soc(model: gurobi.Model):
    x = model.add_variable(lb=0.0, name="x")
    y = model.add_variable(lb=3.0, name="y")
    z = model.add_variable(lb=4.0, name="z")

    obj = x + y + z
    model.set_objective(obj, poi.ObjectiveSense.Minimize)

    soc_con = model.add_second_order_cone_constraint([x, y, z])
    model.optimize()
    status = model.get_model_attribute(poi.ModelAttribute.TerminationStatus)
    assert status == poi.TerminationStatusCode.OPTIMAL
    x_val = model.get_value(x)
    y_val = model.get_value(y)
    z_val = model.get_value(z)
    obj_val = model.get_value(obj)

    mu = model.get_constraint_attribute(soc_con, poi.ConstraintAttribute.Dual)
    print(f"Dual value of soc constraint: {mu}")

    vector_dual = [2 * mu * x_val, -2 * mu * y_val, -2 * mu * z_val]
    print(f"Vector dual: {vector_dual}")


if __name__ == "__main__":
    model = gurobi.Model()
    model.set_raw_parameter("QCPDual", 1)
    test_soc(model)

Julia part:

using JuMP
using Clarabel

function test_soc()
    model = Model(Clarabel.Optimizer)

    @variable(model, x >= 0)
    @variable(model, y >= 3)
    @variable(model, z >= 4)

    @objective(model, Min, x + y + z)

    @constraint(model, soc_con, [x, y, z] in SecondOrderCone())

    optimize!(model)
    @assert termination_status(model) == OPTIMAL

    soc_con_dual = dual(soc_con)
    println("Dual value of soc constraint: $soc_con_dual")
end

test_soc()

Python result:

Dual value of soc constraint: 0.10000000000659735
Vector dual: [1.0000000000659734, -0.6000000000395841, -0.8000000000527788]

Julia result:

Dual value of soc constraint: [1.0000000002709364, -0.6000013885393909, -0.7999989586922533]

[yemine0x01]

I’m not sure that’s correct. My understanding is that you need to derive the dual variables corresponding to each primal variable $x_i$; let’s denote them by $y_i$. Then the dual vector would be $(\mu, y_1, y_2)$

[yemine0x01]

@metab0t you’re probably right. How would you approach implementing this? Would you extend the C interface to support returning a vector, or do you see a better alternative?

[metab0t]

Write a C++ function where the return value is std::vector<double> and nanobind will convert it to a list automatically.

I think a good solution will be implementing a get_vector_dual-like function in C++ side and dispatch on the Python side according to the constraint type. C++ is not good at handling different types of return values.

Besides, vector dual should only be implemented for solvers with native conic support (KNITRO and Mosek) now as a temporary fix.

Other solvers need the transformation I mentioned above, so additional bookkeeping must be added to record the variables in cones. This problem should be left for future work to implement a uniform API to support more powerful affine conic constraint ($Ax+b \in \mathscr{K}$) and their duals. Mosek and COPT have supported this form. (But COPT currently does not return duals other than linear constraints AFAIK)

[yemine0x01]

Totally agree. I will try to add this as soon as possible then.

[yemine0x01]

By the way, would it make more sense to define SOC and Exp as constraint types?

[metab0t]

Yes, go ahead.

The constraint type enum needs to be cleaned up. As you can see, it has contained too much solver-specific stuff. Luckily, we do not need to maintain compatibility because it is mostly an internal implementation detail.

I think this would be a good start:

enum class ConstraintType
{
    Linear,
    Quadratic,
    SOS,
    SecondOrderCone,
    ExponentialCone,
    SolverDefined_NL,
};