From 0d03ddbb2206e23d5bda0482f632e00a845aaf7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 30 Jul 2024 10:30:37 +0200 Subject: [PATCH 01/13] Add support for rank-1 constraints --- src/MOI_wrapper.jl | 126 +++++++++++++++++++++++++++++---------------- src/sdpcone.jl | 48 ++--------------- 2 files changed, 86 insertions(+), 88 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index c895e8e..c0ace1c 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -15,6 +15,10 @@ mutable struct Optimizer <: MOI.AbstractOptimizer # * `-d` means a diagonal block with diagonal of length `d` # * `d` means a symmetric `d x d` block blockdims::Vector{Int} + # MOI variable index -> rank 1 matrix it corresponds to + rank_one::Vector{Union{Nothing,MOI.LowRankMatrix{Cdouble}}} + # To avoid it being free'd + cached_ind::Vector{Vector{Cint}} varmap::Vector{Tuple{Int,Int,Int}} # Variable Index vi -> blk, i, j # If `blockdims[i] < 0`, `blk[i]` is the offset in `lpdvars`. # That is the **sum of length** of diagonal block before @@ -32,6 +36,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer sdpdinds::Vector{Vector{Vector{Cint}}} sdpdcoefs::Vector{Vector{Vector{Cdouble}}} y::Vector{Cdouble} + solve_time::Float64 silent::Bool options::Dict{Symbol,Any} @@ -43,6 +48,8 @@ mutable struct Optimizer <: MOI.AbstractOptimizer 1, Cdouble[], Int[], + Union{Nothing,MOI.LowRankMatrix{Cdouble}}[], + Vector{Cint}[], Tuple{Int,Int,Int}[], Int[], C_NULL, @@ -53,6 +60,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer Vector{Int}[], Vector{Cdouble}[], Cdouble[], + NaN, false, Dict{Symbol,Any}(), ) @@ -80,6 +88,8 @@ function MOI.empty!(optimizer::Optimizer) optimizer.objective_sign = 1 empty!(optimizer.b) empty!(optimizer.blockdims) + empty!(optimizer.rank_one) + empty!(optimizer.cached_ind) empty!(optimizer.varmap) empty!(optimizer.blk) optimizer.nlpdrows = 0 @@ -89,6 +99,7 @@ function MOI.empty!(optimizer::Optimizer) empty!(optimizer.sdpdinds) empty!(optimizer.sdpdcoefs) empty!(optimizer.y) + optimizer.solve_time = NaN return end @@ -218,8 +229,16 @@ end MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false -const SupportedSets = - Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle} +const _SetWithDotProd = MOI.SetWithDotProducts{ + MOI.PositiveSemidefiniteConeTriangle, + MOI.TriangleVectorization{Cdouble,MOI.LowRankMatrix{Cdouble}}, +} + +const SupportedSets = Union{ + MOI.Nonnegatives, + MOI.PositiveSemidefiniteConeTriangle, + _SetWithDotProd, +} function MOI.supports_add_constrained_variables( ::Optimizer, @@ -241,6 +260,7 @@ function new_block(optimizer::Optimizer, set::MOI.Nonnegatives) blk = length(optimizer.blockdims) for i in 1:MOI.dimension(set) push!(optimizer.varmap, (blk, i, i)) + push!(optimizer.rank_one, nothing) end return end @@ -254,14 +274,27 @@ function new_block( for j in 1:set.side_dimension for i in 1:j push!(optimizer.varmap, (blk, i, j)) + push!(optimizer.rank_one, nothing) end end return end +function new_block(model::Optimizer, set::_SetWithDotProd) + println("______________________Low-Rank") + blk = length(model.blockdims) + 1 + for i in eachindex(set.vectors) + push!(model.varmap, (blk, 0, 0)) + push!(model.rank_one, set.vectors[i].matrix) + end + new_block(model, set.set) +end + function _add_constrained_variables(optimizer::Optimizer, set::SupportedSets) offset = length(optimizer.varmap) new_block(optimizer, set) + @assert length(optimizer.varmap) == offset + MOI.dimension(set) + @assert length(optimizer.rank_one) == offset + MOI.dimension(set) ci = MOI.ConstraintIndex{MOI.VectorOfVariables,typeof(set)}(offset + 1) return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci end @@ -306,57 +339,54 @@ function constrain_variables_on_creation( return end -function _setcoefficient!( - m::Optimizer, - coef, - constr::Integer, - blk::Integer, - i::Integer, - j::Integer, -) - if m.blockdims[blk] < 0 - @assert i == j - push!(m.lpdvars, constr + 1) - push!(m.lpdrows, m.blk[blk] + i - 1) # -1 because indexing starts at 0 in DSDP - push!(m.lpcoefs, coef) - else - sdp = m.blk[blk] - push!(m.sdpdinds[end][sdp], i + (j - 1) * m.blockdims[blk] - 1) - if i != j - coef /= 2 +function _setcoefficient!(dest::Optimizer, coef, constr::Integer, vi::MOI.VariableIndex) + blk, i, j = varmap(dest, vi) + rank_one = dest.rank_one[vi.value] + if isnothing(rank_one) + if dest.blockdims[blk] < 0 + @assert i == j + push!(dest.lpdvars, constr + 1) + push!(dest.lpdrows, dest.blk[blk] + i - 1) # -1 because indexing starts at 0 in DSDP + push!(dest.lpcoefs, coef) + else + sdp = dest.blk[blk] + push!(dest.sdpdinds[end][sdp], i + (j - 1) * dest.blockdims[blk] - 1) + if i != j + coef /= 2 + end + push!(dest.sdpdcoefs[end][sdp], coef) end - push!(m.sdpdcoefs[end][sdp], coef) + else + d = Cint(dest.blockdims[blk]) + push!(dest.cached_ind, collect(Cint(0):(d - 1))) + # We use `Add` and not `Set` because I think (if I interpret the name correctly) that would allow mixing with sparse matrices for the same block and constraint + DSDP.SDPCone.SetARankOneMat( + dest.sdpcone, + dest.blk[blk] - 1, + constr, + d, + coef * rank_one.diagonal[], + 0, + last(dest.cached_ind), + collect(eachcol(rank_one.factor)[]), + d, + ) end return end # Loads objective coefficient α * vi -function load_objective_term!( - optimizer::Optimizer, - index_map, - α, - vi::MOI.VariableIndex, -) - blk, i, j = varmap(optimizer, vi) +function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex) coef = optimizer.objective_sign * α - _setcoefficient!(optimizer, coef, 0, blk, i, j) + _setcoefficient!(optimizer, coef, 0, vi) return end function _set_A_matrices(m::Optimizer, i) for (blk, blkdim) in zip(m.blk, m.blockdims) - if blkdim > 0 - SDPCone.SetASparseVecMat( - m.sdpcone, - blk - 1, - i, - blkdim, - 1.0, - 0, - m.sdpdinds[end][blk], - m.sdpdcoefs[end][blk], - length(m.sdpdcoefs[end][blk]), - ) + if blkdim > 0 && !isempty(m.sdpdcoefs[end][blk]) + @show (blk - 1, i, blkdim, 1.0, 0, m.sdpdinds[end][blk], m.sdpdcoefs[end][blk], length(m.sdpdcoefs[end][blk])) + SDPCone.SetASparseVecMat(m.sdpcone, blk - 1, i, blkdim, 1.0, 0, m.sdpdinds[end][blk], m.sdpdcoefs[end][blk], length(m.sdpdcoefs[end][blk])) end end return @@ -386,6 +416,12 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) index_map, MOI.PositiveSemidefiniteConeTriangle, ) + constrain_variables_on_creation( + dest, + src, + index_map, + _SetWithDotProd, + ) vis_src = MOI.get(src, MOI.ListOfVariableIndices()) if length(vis_src) < length(index_map.var_map) _error( @@ -464,8 +500,7 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) _new_A_matrix(dest) for t in func.terms if !iszero(t.coefficient) - blk, i, j = varmap(dest, index_map[t.variable]) - _setcoefficient!(dest, t.coefficient, k, blk, i, j) + _setcoefficient!(dest, t.coefficient, k, index_map[t.variable]) end end _set_A_matrices(dest, k) @@ -508,7 +543,6 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) if !iszero(term.coefficient) load_objective_term!( dest, - index_map, term.coefficient, index_map[term.variable], ) @@ -533,7 +567,9 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) end function MOI.optimize!(m::Optimizer) + start_time = time() Solve(m.dsdp) + m.solve_time = time() - start_time # Calling `ComputeX` not right after `Solve` seems to sometime cause segfaults or weird Heisenbug's # let's call it directly what `DSDP/examples/readsdpa.c` does ComputeX(m.dsdp) @@ -543,6 +579,8 @@ function MOI.optimize!(m::Optimizer) return end +MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) = optimizer.solve_time + function MOI.get(m::Optimizer, ::MOI.RawStatusString) if m.dsdp == C_NULL return "`optimize!` not called" diff --git a/src/sdpcone.jl b/src/sdpcone.jl index 8500452..e60d059 100644 --- a/src/sdpcone.jl +++ b/src/sdpcone.jl @@ -95,28 +95,8 @@ function SetADenseVecMat( ) sdpcone arg2 arg3 arg4 arg5 arg6 arg7 end -function SetARankOneMat( - sdpcone::SDPConeT, - arg2::Integer, - arg3::Integer, - arg4::Integer, - arg5::Cdouble, - arg6::Integer, - arg7::Vector{Cint}, - arg8::Vector{Cdouble}, - arg9::Integer, -) - @dsdp_ccall SDPConeSetARankOneMat ( - SDPConeT, - Cint, - Cint, - Cint, - Cdouble, - Cint, - Ptr{Cint}, - Ptr{Cdouble}, - Cint, - ) sdpcone arg2 arg3 arg4 arg5 arg6 arg7 arg8 arg9 +function SetARankOneMat(sdpcone::SDPConeT, blockj::Integer, vari::Integer, n::Integer, alpha::Cdouble, ishift::Integer, ind::Vector{Cint}, val::Vector{Cdouble}, nnz::Integer) + @dsdp_ccall SDPConeSetARankOneMat (SDPConeT, Cint, Cint, Cint, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint) sdpcone blockj vari n alpha ishift ind val nnz end function SetConstantMat( @@ -220,28 +200,8 @@ function AddIdentity( @dsdp_ccall SDPConeAddIdentity (SDPConeT, Cint, Cint, Cint, Cdouble) sdpcone arg2 arg3 arg4 arg5 end -function AddARankOneMat( - sdpcone::SDPConeT, - arg2::Integer, - arg3::Integer, - arg4::Integer, - arg5::Cdouble, - arg6::Integer, - arg7::Vector{Cint}, - arg8::Vector{Cdouble}, - arg9::Integer, -) - @dsdp_ccall SDPConeAddARankOneMat ( - SDPConeT, - Cint, - Cint, - Cint, - Cdouble, - Cint, - Ptr{Cint}, - Ptr{Cdouble}, - Cint, - ) sdpcone arg2 arg3 arg4 arg5 arg6 arg7 arg8 arg9 +function AddARankOneMat(sdpcone::SDPConeT, blockj::Integer, vari::Integer, n::Integer, alpha::Cdouble, ishift::Integer, ind::Vector{Cint}, val::Vector{Cdouble}, nnz::Integer) + @dsdp_ccall SDPConeAddARankOneMat (SDPConeT, Cint, Cint, Cint, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint) sdpcone blockj vari n alpha ishift ind val nnz end function AddSparseVecMat( From bce1aa3e5ef28033b33c9a76577fc1d9e14adfd6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 15 Apr 2025 10:57:28 +0200 Subject: [PATCH 02/13] Checkout LRO --- .github/workflows/ci.yml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8f4c9c2..3f4b12e 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,6 +25,13 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 + - name: LRO + shell: julia --project=@. {0} + run: | + using Pkg + Pkg.add([ + PackageSpec(url="https://github.com/blegat/LowRankOpt.jl/"), + ]) - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From 825cf9f8f856d06dc198abdb59b0168c98d8f60a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 15 Apr 2025 10:58:59 +0200 Subject: [PATCH 03/13] Fix format --- src/DSDP.jl | 2 +- src/MOI_wrapper.jl | 53 +++++++++++++++++++++++++++++++--------------- src/sdpcone.jl | 48 +++++++++++++++++++++++++++++++++++++---- test/maxcut.jl | 2 +- 4 files changed, 82 insertions(+), 23 deletions(-) diff --git a/src/DSDP.jl b/src/DSDP.jl index 41399b9..94728d7 100644 --- a/src/DSDP.jl +++ b/src/DSDP.jl @@ -10,7 +10,7 @@ import DSDP_jll using LinearAlgebra macro dsdp_ccall(f, args...) - quote + return quote # QuoteNode prevents the interpretion of the symbol # and leave it as a symbol info = diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index c0ace1c..a2e6a26 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -234,11 +234,8 @@ const _SetWithDotProd = MOI.SetWithDotProducts{ MOI.TriangleVectorization{Cdouble,MOI.LowRankMatrix{Cdouble}}, } -const SupportedSets = Union{ - MOI.Nonnegatives, - MOI.PositiveSemidefiniteConeTriangle, - _SetWithDotProd, -} +const SupportedSets = + Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle,_SetWithDotProd} function MOI.supports_add_constrained_variables( ::Optimizer, @@ -287,7 +284,7 @@ function new_block(model::Optimizer, set::_SetWithDotProd) push!(model.varmap, (blk, 0, 0)) push!(model.rank_one, set.vectors[i].matrix) end - new_block(model, set.set) + return new_block(model, set.set) end function _add_constrained_variables(optimizer::Optimizer, set::SupportedSets) @@ -339,7 +336,12 @@ function constrain_variables_on_creation( return end -function _setcoefficient!(dest::Optimizer, coef, constr::Integer, vi::MOI.VariableIndex) +function _setcoefficient!( + dest::Optimizer, + coef, + constr::Integer, + vi::MOI.VariableIndex, +) blk, i, j = varmap(dest, vi) rank_one = dest.rank_one[vi.value] if isnothing(rank_one) @@ -350,7 +352,10 @@ function _setcoefficient!(dest::Optimizer, coef, constr::Integer, vi::MOI.Variab push!(dest.lpcoefs, coef) else sdp = dest.blk[blk] - push!(dest.sdpdinds[end][sdp], i + (j - 1) * dest.blockdims[blk] - 1) + push!( + dest.sdpdinds[end][sdp], + i + (j - 1) * dest.blockdims[blk] - 1, + ) if i != j coef /= 2 end @@ -358,7 +363,7 @@ function _setcoefficient!(dest::Optimizer, coef, constr::Integer, vi::MOI.Variab end else d = Cint(dest.blockdims[blk]) - push!(dest.cached_ind, collect(Cint(0):(d - 1))) + push!(dest.cached_ind, collect(Cint(0):(d-1))) # We use `Add` and not `Set` because I think (if I interpret the name correctly) that would allow mixing with sparse matrices for the same block and constraint DSDP.SDPCone.SetARankOneMat( dest.sdpcone, @@ -385,8 +390,27 @@ end function _set_A_matrices(m::Optimizer, i) for (blk, blkdim) in zip(m.blk, m.blockdims) if blkdim > 0 && !isempty(m.sdpdcoefs[end][blk]) - @show (blk - 1, i, blkdim, 1.0, 0, m.sdpdinds[end][blk], m.sdpdcoefs[end][blk], length(m.sdpdcoefs[end][blk])) - SDPCone.SetASparseVecMat(m.sdpcone, blk - 1, i, blkdim, 1.0, 0, m.sdpdinds[end][blk], m.sdpdcoefs[end][blk], length(m.sdpdcoefs[end][blk])) + @show ( + blk - 1, + i, + blkdim, + 1.0, + 0, + m.sdpdinds[end][blk], + m.sdpdcoefs[end][blk], + length(m.sdpdcoefs[end][blk]), + ) + SDPCone.SetASparseVecMat( + m.sdpcone, + blk - 1, + i, + blkdim, + 1.0, + 0, + m.sdpdinds[end][blk], + m.sdpdcoefs[end][blk], + length(m.sdpdcoefs[end][blk]), + ) end end return @@ -416,12 +440,7 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) index_map, MOI.PositiveSemidefiniteConeTriangle, ) - constrain_variables_on_creation( - dest, - src, - index_map, - _SetWithDotProd, - ) + constrain_variables_on_creation(dest, src, index_map, _SetWithDotProd) vis_src = MOI.get(src, MOI.ListOfVariableIndices()) if length(vis_src) < length(index_map.var_map) _error( diff --git a/src/sdpcone.jl b/src/sdpcone.jl index e60d059..da18462 100644 --- a/src/sdpcone.jl +++ b/src/sdpcone.jl @@ -95,8 +95,28 @@ function SetADenseVecMat( ) sdpcone arg2 arg3 arg4 arg5 arg6 arg7 end -function SetARankOneMat(sdpcone::SDPConeT, blockj::Integer, vari::Integer, n::Integer, alpha::Cdouble, ishift::Integer, ind::Vector{Cint}, val::Vector{Cdouble}, nnz::Integer) - @dsdp_ccall SDPConeSetARankOneMat (SDPConeT, Cint, Cint, Cint, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint) sdpcone blockj vari n alpha ishift ind val nnz +function SetARankOneMat( + sdpcone::SDPConeT, + blockj::Integer, + vari::Integer, + n::Integer, + alpha::Cdouble, + ishift::Integer, + ind::Vector{Cint}, + val::Vector{Cdouble}, + nnz::Integer, +) + @dsdp_ccall SDPConeSetARankOneMat ( + SDPConeT, + Cint, + Cint, + Cint, + Cdouble, + Cint, + Ptr{Cint}, + Ptr{Cdouble}, + Cint, + ) sdpcone blockj vari n alpha ishift ind val nnz end function SetConstantMat( @@ -200,8 +220,28 @@ function AddIdentity( @dsdp_ccall SDPConeAddIdentity (SDPConeT, Cint, Cint, Cint, Cdouble) sdpcone arg2 arg3 arg4 arg5 end -function AddARankOneMat(sdpcone::SDPConeT, blockj::Integer, vari::Integer, n::Integer, alpha::Cdouble, ishift::Integer, ind::Vector{Cint}, val::Vector{Cdouble}, nnz::Integer) - @dsdp_ccall SDPConeAddARankOneMat (SDPConeT, Cint, Cint, Cint, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}, Cint) sdpcone blockj vari n alpha ishift ind val nnz +function AddARankOneMat( + sdpcone::SDPConeT, + blockj::Integer, + vari::Integer, + n::Integer, + alpha::Cdouble, + ishift::Integer, + ind::Vector{Cint}, + val::Vector{Cdouble}, + nnz::Integer, +) + @dsdp_ccall SDPConeAddARankOneMat ( + SDPConeT, + Cint, + Cint, + Cint, + Cdouble, + Cint, + Ptr{Cint}, + Ptr{Cdouble}, + Cint, + ) sdpcone blockj vari n alpha ishift ind val nnz end function AddSparseVecMat( diff --git a/test/maxcut.jl b/test/maxcut.jl index cd496a2..aa56785 100644 --- a/test/maxcut.jl +++ b/test/maxcut.jl @@ -70,7 +70,7 @@ function maxcut(nnodes, edges) yy = zeros(nnodes) indd = zeros(Cint, nnodes + nedges) val = zeros(nnodes + nedges) - indd[nedges.+(1:nnodes)] = iptr + indd[nedges .+ (1:nnodes)] = iptr tval = 0.0 for (i, (u, v, w)) in enumerate(edges) indd[i] = di(u, v) From 14ceb6d0237f44ee93222f7b4a9ad1620842baba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 16 Apr 2025 11:38:40 +0200 Subject: [PATCH 04/13] Fix format --- test/maxcut.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/maxcut.jl b/test/maxcut.jl index aa56785..cd496a2 100644 --- a/test/maxcut.jl +++ b/test/maxcut.jl @@ -70,7 +70,7 @@ function maxcut(nnodes, edges) yy = zeros(nnodes) indd = zeros(Cint, nnodes + nedges) val = zeros(nnodes + nedges) - indd[nedges .+ (1:nnodes)] = iptr + indd[nedges.+(1:nnodes)] = iptr tval = 0.0 for (i, (u, v, w)) in enumerate(edges) indd[i] = di(u, v) From 7435eb0f42eb6e469458aabbc27cb57a44616f94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 16 Apr 2025 11:51:19 +0200 Subject: [PATCH 05/13] Fixes --- Project.toml | 1 + src/MOI_wrapper.jl | 37 ++++++++++++++++--------------------- 2 files changed, 17 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 8bc5e08..3bdbd1b 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "0.2.1" [deps] DSDP_jll = "1065e140-e56c-5613-be8b-7480bf7138df" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [compat] diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index a2e6a26..6d6cca8 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -4,6 +4,17 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. import MathOptInterface as MOI +import LowRankOpt as LRO + +const _RankOneMatrix{F<:AbstractVector{Cdouble},D<:AbstractArray{Cdouble,0}} = + LRO.Factorization{Cdouble,F,D} + +const _SetDotProd{F<:AbstractMatrix{Cdouble},D<:AbstractVector{Cdouble}} = + LRO.SetDotProducts{ + LRO.WITH_SET, + MOI.PositiveSemidefiniteConeTriangle, + LRO.TriangleVectorization{Cdouble,_RankOneMatrix{F,D}}, + } mutable struct Optimizer <: MOI.AbstractOptimizer dsdp::DSDPT @@ -16,7 +27,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer # * `d` means a symmetric `d x d` block blockdims::Vector{Int} # MOI variable index -> rank 1 matrix it corresponds to - rank_one::Vector{Union{Nothing,MOI.LowRankMatrix{Cdouble}}} + rank_one::Vector{Union{Nothing,_RankOneMatrix}} # To avoid it being free'd cached_ind::Vector{Vector{Cint}} varmap::Vector{Tuple{Int,Int,Int}} # Variable Index vi -> blk, i, j @@ -48,7 +59,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer 1, Cdouble[], Int[], - Union{Nothing,MOI.LowRankMatrix{Cdouble}}[], + Union{Nothing,_RankOneMatrix}[], Vector{Cint}[], Tuple{Int,Int,Int}[], Int[], @@ -229,13 +240,8 @@ end MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false -const _SetWithDotProd = MOI.SetWithDotProducts{ - MOI.PositiveSemidefiniteConeTriangle, - MOI.TriangleVectorization{Cdouble,MOI.LowRankMatrix{Cdouble}}, -} - const SupportedSets = - Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle,_SetWithDotProd} + Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle,_SetDotProd} function MOI.supports_add_constrained_variables( ::Optimizer, @@ -277,8 +283,7 @@ function new_block( return end -function new_block(model::Optimizer, set::_SetWithDotProd) - println("______________________Low-Rank") +function new_block(model::Optimizer, set::_SetDotProd) blk = length(model.blockdims) + 1 for i in eachindex(set.vectors) push!(model.varmap, (blk, 0, 0)) @@ -390,16 +395,6 @@ end function _set_A_matrices(m::Optimizer, i) for (blk, blkdim) in zip(m.blk, m.blockdims) if blkdim > 0 && !isempty(m.sdpdcoefs[end][blk]) - @show ( - blk - 1, - i, - blkdim, - 1.0, - 0, - m.sdpdinds[end][blk], - m.sdpdcoefs[end][blk], - length(m.sdpdcoefs[end][blk]), - ) SDPCone.SetASparseVecMat( m.sdpcone, blk - 1, @@ -440,7 +435,7 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) index_map, MOI.PositiveSemidefiniteConeTriangle, ) - constrain_variables_on_creation(dest, src, index_map, _SetWithDotProd) + constrain_variables_on_creation(dest, src, index_map, _SetDotProd) vis_src = MOI.get(src, MOI.ListOfVariableIndices()) if length(vis_src) < length(index_map.var_map) _error( From 82635739de682d7a7ca10cdaa4c4b836b5c39c45 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 16 Apr 2025 11:51:31 +0200 Subject: [PATCH 06/13] Add test/Project.toml --- test/Project.toml | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 test/Project.toml diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 0000000..27b18d3 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,4 @@ +[deps] +DSDP = "2714ae6b-e930-5b4e-9c21-d0bacf577842" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 9c830904f898129250396dc80fb61ed4b4790159 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 16 Apr 2025 21:05:42 +0200 Subject: [PATCH 07/13] Add test --- test/MOI_wrapper.jl | 19 +++++++++++++++++++ test/Project.toml | 1 + 2 files changed, 20 insertions(+) diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 6f90d76..209c824 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -7,6 +7,7 @@ module TestDSDP using Test import MathOptInterface as MOI +import LowRankOpt as LRO import DSDP function runtests() @@ -108,6 +109,24 @@ function test_runtests() return end +function test_LRO_runtests() + T = Float64 + model = MOI.instantiate( + DSDP.Optimizer, + with_bridge_type = T, + with_cache_type = T, + ) + LRO.Bridges.add_all_bridges(model, T) + MOI.set(model, MOI.Silent(), true) + config = MOI.Test.Config( + rtol = 1e-2, + atol = 1e-2, + ) + MOI.Test.runtests(model, config, test_module = LRO.Test) + return +end + + end # module TestDSDP.runtests() diff --git a/test/Project.toml b/test/Project.toml index 27b18d3..729396d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] DSDP = "2714ae6b-e930-5b4e-9c21-d0bacf577842" +LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 03d5ba822afc295e56d41d4be8967335cf88abcc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 28 Apr 2026 17:14:31 +0200 Subject: [PATCH 08/13] Fix --- Project.toml | 4 +++- src/MOI_wrapper.jl | 12 +++++++++--- test/Project.toml | 1 + 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 3bdbd1b..235e5b1 100644 --- a/Project.toml +++ b/Project.toml @@ -11,11 +11,13 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [compat] DSDP_jll = "0.0.1" +LowRankOpt = "0.2" MathOptInterface = "1" julia = "1.10" [extras] +Dualization = "191a621a-6537-11e9-281d-650236a99e60" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test"] +test = ["Dualization", "Test"] diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 6d6cca8..575f7c1 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -9,7 +9,7 @@ import LowRankOpt as LRO const _RankOneMatrix{F<:AbstractVector{Cdouble},D<:AbstractArray{Cdouble,0}} = LRO.Factorization{Cdouble,F,D} -const _SetDotProd{F<:AbstractMatrix{Cdouble},D<:AbstractVector{Cdouble}} = +const _SetDotProd{F<:AbstractVector{Cdouble},D<:AbstractArray{Cdouble,0}} = LRO.SetDotProducts{ LRO.WITH_SET, MOI.PositiveSemidefiniteConeTriangle, @@ -435,9 +435,15 @@ function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) index_map, MOI.PositiveSemidefiniteConeTriangle, ) - constrain_variables_on_creation(dest, src, index_map, _SetDotProd) + # `_SetDotProd` is a parametric UnionAll type, so we need to find the + # concrete type from the model's actual constraint types. + for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent()) + if F == MOI.VectorOfVariables && S <: _SetDotProd + constrain_variables_on_creation(dest, src, index_map, S) + end + end vis_src = MOI.get(src, MOI.ListOfVariableIndices()) - if length(vis_src) < length(index_map.var_map) + if length(vis_src) > length(index_map.var_map) _error( "Free variables are not supported by DSDP", "to bridge free variables into `x - y` where `x` and `y` are nonnegative.", diff --git a/test/Project.toml b/test/Project.toml index 729396d..a495482 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] DSDP = "2714ae6b-e930-5b4e-9c21-d0bacf577842" +Dualization = "191a621a-6537-11e9-281d-650236a99e60" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 11094d7250cbd79a75db02bf07e11ed197ddef69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 29 Apr 2026 10:12:10 +0200 Subject: [PATCH 09/13] Use our own test --- Project.toml | 3 +-- src/MOI_wrapper.jl | 9 +++++---- test/MOI_wrapper.jl | 37 +++++++++++++++++++++++++++++++++---- test/Project.toml | 1 - 4 files changed, 39 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 235e5b1..6116ec1 100644 --- a/Project.toml +++ b/Project.toml @@ -16,8 +16,7 @@ MathOptInterface = "1" julia = "1.10" [extras] -Dualization = "191a621a-6537-11e9-281d-650236a99e60" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Dualization", "Test"] +test = ["Test"] diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 575f7c1..af4fc7c 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -369,16 +369,17 @@ function _setcoefficient!( else d = Cint(dest.blockdims[blk]) push!(dest.cached_ind, collect(Cint(0):(d-1))) - # We use `Add` and not `Set` because I think (if I interpret the name correctly) that would allow mixing with sparse matrices for the same block and constraint - DSDP.SDPCone.SetARankOneMat( + # Use `Add` (not `Set`) so that multiple rank-one contributions + # to the same block and constraint accumulate instead of overwriting. + DSDP.SDPCone.AddARankOneMat( dest.sdpcone, dest.blk[blk] - 1, constr, d, - coef * rank_one.diagonal[], + coef * rank_one.scaling[], 0, last(dest.cached_ind), - collect(eachcol(rank_one.factor)[]), + rank_one.factor, d, ) end diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 209c824..1667526 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -110,7 +110,26 @@ function test_runtests() end function test_LRO_runtests() + # The LRO.Test tests (polynomial and moment) formulate problems with a + # free variable γ: + # max γ s.t. [3-γ, -1-γ] ∈ SetDotProducts(PSD(2), [v₁v₁ᵀ, v₂v₂ᵀ]) + # Since DSDP is a dual-only interior-point method, free variables + # (bridged as z⁺ - z⁻ with z⁺, z⁻ ≥ 0) create a degenerate dual with + # no strict interior, causing poor accuracy. See Anjos & Burer (2007), + # "On handling free variables in interior-point methods for conic linear + # optimization", SIAM J. Optim. + # + # Instead, we hardcode the equivalent problem with γ eliminated: + # min ⟨v₁v₁ᵀ, X⟩ s.t. ⟨v₁v₁ᵀ - v₂v₂ᵀ, X⟩ = 4, X ≽ 0 + # expressed via SetDotProducts{WITH_SET} as a variable constraint + # (no free variables needed). T = Float64 + v1 = LRO.positive_semidefinite_factorization(T[1, -1]) + v2 = LRO.positive_semidefinite_factorization(T[1, 1]) + set = LRO.SetDotProducts{LRO.WITH_SET}( + MOI.PositiveSemidefiniteConeTriangle(2), + LRO.TriangleVectorization.([v1, v2]), + ) model = MOI.instantiate( DSDP.Optimizer, with_bridge_type = T, @@ -118,11 +137,21 @@ function test_LRO_runtests() ) LRO.Bridges.add_all_bridges(model, T) MOI.set(model, MOI.Silent(), true) - config = MOI.Test.Config( - rtol = 1e-2, - atol = 1e-2, + vars, cv = MOI.add_constrained_variables(model, set) + y1, y2 = vars[1], vars[2] + # y1 = ⟨v₁v₁ᵀ, X⟩, y2 = ⟨v₂v₂ᵀ, X⟩ + # Constraint: y1 - y2 = 4 (from eliminating γ in y1 + γ = 3, y2 + γ = -1) + MOI.add_constraint(model, T(1) * y1 - T(1) * y2, MOI.EqualTo(T(4))) + MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.set( + model, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), + T(1) * y1, ) - MOI.Test.runtests(model, config, test_module = LRO.Test) + MOI.optimize!(model) + @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL + # Optimal: y1 = 4 (= 3 - γ* = 3 - (-1)), y2 = 0 (= -1 - γ* = -1 - (-1)) + @test MOI.get(model, MOI.ObjectiveValue()) ≈ T(4) atol = 1e-2 return end diff --git a/test/Project.toml b/test/Project.toml index a495482..729396d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,5 @@ [deps] DSDP = "2714ae6b-e930-5b4e-9c21-d0bacf577842" -Dualization = "191a621a-6537-11e9-281d-650236a99e60" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 305415e76f810d3a98f85086d9f1e875772b7020 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 29 Apr 2026 14:03:15 +0200 Subject: [PATCH 10/13] Fix format --- test/MOI_wrapper.jl | 3 +-- test/maxcut.jl | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 1667526..6a6d709 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -131,7 +131,7 @@ function test_LRO_runtests() LRO.TriangleVectorization.([v1, v2]), ) model = MOI.instantiate( - DSDP.Optimizer, + DSDP.Optimizer; with_bridge_type = T, with_cache_type = T, ) @@ -155,7 +155,6 @@ function test_LRO_runtests() return end - end # module TestDSDP.runtests() diff --git a/test/maxcut.jl b/test/maxcut.jl index cd496a2..aa56785 100644 --- a/test/maxcut.jl +++ b/test/maxcut.jl @@ -70,7 +70,7 @@ function maxcut(nnodes, edges) yy = zeros(nnodes) indd = zeros(Cint, nnodes + nedges) val = zeros(nnodes + nedges) - indd[nedges.+(1:nnodes)] = iptr + indd[nedges .+ (1:nnodes)] = iptr tval = 0.0 for (i, (u, v, w)) in enumerate(edges) indd[i] = di(u, v) From 3511c0708716d41ee02e541b579f5add6dfb0107 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 29 Apr 2026 14:07:29 +0200 Subject: [PATCH 11/13] Fix format --- test/maxcut.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/maxcut.jl b/test/maxcut.jl index aa56785..cd496a2 100644 --- a/test/maxcut.jl +++ b/test/maxcut.jl @@ -70,7 +70,7 @@ function maxcut(nnodes, edges) yy = zeros(nnodes) indd = zeros(Cint, nnodes + nedges) val = zeros(nnodes + nedges) - indd[nedges .+ (1:nnodes)] = iptr + indd[nedges.+(1:nnodes)] = iptr tval = 0.0 for (i, (u, v, w)) in enumerate(edges) indd[i] = di(u, v) From 505785da5f57e4b391a6d5f77e1e17034a681267 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 29 Apr 2026 14:19:50 +0200 Subject: [PATCH 12/13] Fix --- .github/workflows/ci.yml | 7 ------- 1 file changed, 7 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3f4b12e..8f4c9c2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -25,13 +25,6 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v2 - - name: LRO - shell: julia --project=@. {0} - run: | - using Pkg - Pkg.add([ - PackageSpec(url="https://github.com/blegat/LowRankOpt.jl/"), - ]) - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From 3669edcb8741594c3b32eda56b06334498cdc2fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 29 Apr 2026 14:21:34 +0200 Subject: [PATCH 13/13] Remove dup --- test/MOI_wrapper.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 6a6d709..6660380 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -86,7 +86,6 @@ function test_runtests() r"test_solve_TerminationStatus_DUAL_INFEASIBLE$", r"test_DualObjectiveValue_Max_VariableIndex_LessThan$", r"test_DualObjectiveValue_Min_VariableIndex_GreaterThan$", - # ArgumentError: DSDP does not support problems with no constraint. r"test_conic_SecondOrderCone_negative_initial_bound$", r"test_conic_SecondOrderCone_negative_post_bound$", r"test_conic_SecondOrderCone_nonnegative_initial_bound$",