Add conic to ridge sentitivity (#231)

* Quadratic functions with conic constraints

* Fix format

* Remove extra section

* Move compat entries

* Add docstring

Author: blegat

Project.toml

@@ -17,22 +17,9 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 [compat]
 BlockDiagonals = "0.1"
 ChainRulesCore = "1"
-HiGHS = "1"
-Ipopt = "1"
 IterativeSolvers = "0.9"
 JuMP = "1"
 LazyArrays = "0.21, 0.22, 1"
-MathOptInterface = "1"
-MathOptSetDistances = "0.2"
-SCS = "1"
+MathOptInterface = "1.14.1"
+MathOptSetDistances = "0.2.7"
 julia = "1.6"
-
-[extras]
-DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b"
-Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
-SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["DelimitedFiles", "HiGHS", "Ipopt", "SCS", "Test"]

src/ConicProgram/ConicProgram.jl

@@ -50,6 +50,12 @@ const Form{T} = MOI.Utilities.GenericModel{
         DiffOpt.ProductOfSets{T},
     },
 }
+function MOI.supports(
+    ::Form{T},
+    ::MOI.ObjectiveFunction{F},
+) where {T,F<:MOI.AbstractFunction}
+    return F === MOI.ScalarAffineFunction{T}
+end
 
 """
     Diffopt.ConicProgram.Model <: DiffOpt.AbstractModel
@@ -199,6 +205,8 @@ function _gradient_cache(model::Model)
     n = A.n
     N = m + n + 1
     # NOTE: w = 1.0 systematically since we asserted the primal-dual pair is optimal
+    # `inv(M)((x, y, 1), (0, s, 0)) = (x, y, 1) - (0, s, 0)`,
+    # see Minty parametrization in https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
     (u, v, w) = (model.x, model.y - model.s, 1.0)
 
     # find gradient of projections on dual of the cones
@@ -209,16 +217,16 @@ function _gradient_cache(model::Model)
     #     -A   0    b;
     #     -c' -b'   0;
     # ]
-    # M = (Q- I) * B + I
+    # M = (Q - I) * B + I
     # with B =
     # [
-    #  I    .   .    # Πx = x because x is a solution and hence satistfies the constraints
-    #  .  Dπv   .
-    #  .    .   1    # w >= 0, but in the solution x = 1
+    #  I    .   .    # Πx = x because it projects on R^n
+    #  .  Dπv   .    # Derivative of the projection onto the dual cones
+    #  .    .   1    # Projection onto R_+ but w is 1 so the derivative is 1
     # ]
     # see: https://stanford.edu/~boyd/papers/pdf/cone_prog_refine.pdf
     # for the definition of Π and why we get I and 1 for x and w respectectively
-    # K is defined in (5), Π in sect 2, and projection sin sect 3
+    # K is defined in (5), Π in sect 2, and projections in sect 3
 
     M = [
         SparseArrays.spzeros(n, n) (A'*Dπv) c

src/DiffOpt.jl

@@ -20,6 +20,7 @@ include("diff_opt.jl")
 include("moi_wrapper.jl")
 include("jump_moi_overloads.jl")
 
+include("copy_dual.jl")
 include("bridges.jl")
 
 include("QuadraticProgram/QuadraticProgram.jl")

src/bridges.jl

@@ -3,6 +3,23 @@
 # Use of this source code is governed by an MIT-style license that can be found
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
+function MOI.get(
+    model::MOI.ModelLike,
+    ::ObjectiveFunctionAttribute{ReverseObjectiveFunction},
+    bridge::MOI.Bridges.Objective.SlackBridge,
+)
+    return MOI.get(model, ReverseConstraintFunction(), bridge.constraint)
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    ::ObjectiveFunctionAttribute{ForwardObjectiveFunction},
+    bridge::MOI.Bridges.Objective.SlackBridge,
+    value,
+)
+    return MOI.set(model, ForwardConstraintFunction(), bridge.constraint, value)
+end
+
 function MOI.set(
     model::MOI.ModelLike,
     attr::ForwardConstraintFunction,
@@ -44,6 +61,23 @@ function MOI.set(
     return MOI.set(model, attr, bridge.constraint, mapped_func)
 end
 
+function _variable_to_index_map(bridge)
+    return Dict{MOI.VariableIndex,MOI.VariableIndex}(
+        v => MOI.VariableIndex(i) for
+        (i, v) in enumerate(bridge.index_to_variable_map)
+    )
+end
+
+function _U(func::MOI.VectorAffineFunction, n, variable_to_index_map)
+    # x'Qx/2 + a'x + β is bridged into
+    # (1, -a'x - β, U*x) in RSOC
+    # where Q = U' * U
+    # The linear part of `func` is `[0; -a'; U]`
+    Ux = MOI.Utilities.eachscalar(func)[3:end]
+    func_array = sparse_array_representation(Ux, n, variable_to_index_map)
+    return func_array.terms
+end
+
 function MOI.get(
     model::MOI.ModelLike,
     attr::DiffOpt.ReverseConstraintFunction,
@@ -53,55 +87,167 @@ function MOI.get(
     return MOI.Bridges.adjoint_map_function(typeof(bridge), func)
 end
 
-function MOI.set(
+function MOI.get(
     model::MOI.ModelLike,
-    attr::DiffOpt.ForwardConstraintFunction,
+    attr::DiffOpt.ReverseConstraintFunction,
     bridge::MOI.Bridges.Constraint.QuadtoSOCBridge{T},
-    diff::MOI.ScalarQuadraticFunction,
 ) where {T}
-    func = MOI.get(model, MOI.ConstraintFunction(), bridge.soc)
-    index_map = index_map_to_oneto(func)
-    index_map_to_oneto!(index_map, diff)
-    n = length(index_map.var_map)
-    func_array = sparse_array_representation(func, n, index_map)
-    # x'Qx/2 + a'x + β is bridged into
-    # (1, -a'x - β, U*x) in RSOC
-    # where Q = U' * U
-    # `func_array.terms` is `[0; -a'; U]`
-    U = func_array.terms[3:end, :]
-    diff_array = sparse_array_representation(diff, n, index_map)
-    dU = Matrix(diff_array.quadratic_terms)
-    dU = dU_from_dQ!(Matrix(diff_array.quadratic_terms), U)
-    x = Vector{MOI.VariableIndex}(undef, n)
-    for v in keys(index_map.var_map)
-        x[index_map[v].value] = v
+    variable_to_index_map = _variable_to_index_map(bridge)
+    n = length(bridge.index_to_variable_map)
+    U = _U(
+        MOI.get(model, MOI.ConstraintFunction(), bridge.soc),
+        n,
+        variable_to_index_map,
+    )
+    Δfunc =
+        convert(MOI.VectorAffineFunction{T}, MOI.get(model, attr, bridge.soc))
+    filter!(Δfunc.terms) do t
+        return haskey(variable_to_index_map, t.scalar_term.variable)
     end
-    diff_aff = MOI.VectorAffineFunction{T}(
-        MOI.VectorAffineTerm{T}[],
-        zeros(T, size(dU, 1) + 2),
+    aff = sparse_array_representation(
+        -MOI.Utilities.eachscalar(Δfunc)[2],
+        n,
+        variable_to_index_map,
+    )
+    ΔU = _U(Δfunc, n, variable_to_index_map)
+    ΔQ = ΔQ_from_ΔU!(ΔU, U)
+    func = MatrixScalarQuadraticFunction(
+        convert(VectorScalarAffineFunction{T,Vector{T}}, aff),
+        ΔQ,
     )
-    for t in diff.affine_terms
+    index_map = MOI.Utilities.IndexMap()
+    for (i, vi) in enumerate(bridge.index_to_variable_map)
+        index_map[MOI.VariableIndex(i)] = vi
+    end
+    Δ = standard_form(IndexMappedFunction(func, index_map))
+    if !bridge.less_than
+        Δ = -Δ
+    end
+    return Δ
+end
+
+function _quad_to_soc_diff(
+    ::MOI.ModelLike,
+    bridge::MOI.Bridges.Constraint.QuadtoSOCBridge{T},
+    diff::MOI.ScalarAffineFunction{T},
+) where {T}
+    n = length(bridge.index_to_variable_map)
+    diff_soc =
+        MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], zeros(T, n + 2))
+    for t in diff.terms
         push!(
-            diff_aff.terms,
+            diff_soc.terms,
             MOI.VectorAffineTerm(
                 2,
                 MOI.ScalarAffineTerm(-t.coefficient, t.variable),
             ),
         )
     end
-    diff_aff.constants[2] = -diff.constant
+    diff_soc.constants[2] = -diff.constant
+    return diff_soc
+end
+
+function _quad_to_soc_diff(
+    model::MOI.ModelLike,
+    bridge::MOI.Bridges.Constraint.QuadtoSOCBridge{T},
+    diff::MOI.ScalarQuadraticFunction{T},
+) where {T}
+    variable_to_index_map = _variable_to_index_map(bridge)
+    n = length(bridge.index_to_variable_map)
+    U = _U(
+        MOI.get(model, MOI.ConstraintFunction(), bridge.soc),
+        n,
+        variable_to_index_map,
+    )
+    # We assume `diff.quadratic_terms` only contains quadratic terms with
+    # variables already in a quadratic term of the initial constraint
+    # Otherwise, the `U` matrix will not be square so it's a TODO
+    # We remove the affine terms here as they might have other variables not
+    # in `variable_to_index_map`
+    diff_quad = MOI.ScalarQuadraticFunction(
+        diff.quadratic_terms,
+        MOI.ScalarAffineTerm{T}[],
+        zero(T),
+    )
+    diff_aff = MOI.ScalarAffineFunction(diff.affine_terms, diff.constant)
+    diff_array =
+        sparse_array_representation(diff_quad, n, variable_to_index_map)
+    dU = dU_from_dQ!(Matrix(diff_array.quadratic_terms), U)
+    diff_soc = _quad_to_soc_diff(model, bridge, diff_aff)
     for i in axes(dU, 1)
         for j in axes(dU, 2)
             if !iszero(dU[i, j])
-                scalar = MOI.ScalarAffineTerm(dU[i, j], x[j])
-                push!(diff_aff.terms, MOI.VectorAffineTerm(2 + i, scalar))
+                scalar = MOI.ScalarAffineTerm(
+                    dU[i, j],
+                    bridge.index_to_variable_map[j],
+                )
+                push!(diff_soc.terms, MOI.VectorAffineTerm(2 + i, scalar))
             end
         end
     end
-    MOI.set(model, attr, bridge.soc, diff_aff)
+    return diff_soc
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::DiffOpt.ForwardConstraintFunction,
+    bridge::MOI.Bridges.Constraint.QuadtoSOCBridge{T},
+    diff::MOI.AbstractScalarFunction,
+) where {T}
+    if !bridge.less_than
+        diff = -diff
+    end
+    diff_soc = _quad_to_soc_diff(model, bridge, diff)
+    MOI.set(model, attr, bridge.soc, diff_soc)
     return
 end
 
+"""
+    ΔQ_from_ΔU!(ΔU, U)
+
+Return the symmetric solution `ΔQ` of the matrix equation
+`triu(ΔU) = 2triu(U * ΔQ)`
+where `ΔU` and `U` are the two argument of the function.
+
+This function overwrites the first argument `ΔU` to store the solution.
+The matrix `U` is not however modified.
+
+The matrix `U` is assumed to be upper triangular.
+
+We can exploit the structure of `U` here:
+
+* If the factorization was obtained from SVD, `U` would be orthogonal
+* If the factorization was obtained from Cholesky, `U` would be upper triangular.
+
+The MOI bridge uses Cholesky in order to exploit sparsity so we are in the
+second case.
+
+We can find each column of `ΔQ` by solving a triangular linear system.
+"""
+function ΔQ_from_ΔU!(ΔU, U)
+    n = LinearAlgebra.checksquare(ΔU)
+    LinearAlgebra.rdiv!(ΔU, 2)
+    for j in n:-1:1
+        # FIXME MA does not support mutating subarray
+        #MA.operate!(MA.sub_mul, view(ΔU, 1:j, j), view(U, 1:j, (j + 1):n), view(ΔU, j, (j+1):n))
+        for row in 1:j
+            for col in (j+1):n
+                ΔU[row, j] -= U[row, col] * ΔU[j, col]
+            end
+        end
+        _U = LinearAlgebra.UpperTriangular(view(U, 1:j, 1:j))
+        LinearAlgebra.ldiv!(_U, view(ΔU, 1:j, j))
+    end
+    for j in 1:n
+        for i in 1:(j-1)
+            if i != j
+                ΔU[j, i] = ΔU[i, j]
+            end
+        end
+    end
+    return ΔU
+end
+
 """
     dU_from_dQ!(dQ, U)