From f61f012ffbf785a686d2eaf65c2b9fc988ba64e4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 28 May 2025 20:01:36 +0200 Subject: [PATCH 01/48] Add BurerMonteiro formulation --- src/burer_monteiro.jl | 37 ++++++++ src/model.jl | 201 ++++++++++++++++++++++++++++++++++++++++++ src/schur.jl | 173 ++++++++++++++++++++++++++++++++++++ 3 files changed, 411 insertions(+) create mode 100644 src/burer_monteiro.jl create mode 100644 src/model.jl create mode 100644 src/schur.jl diff --git a/src/burer_monteiro.jl b/src/burer_monteiro.jl new file mode 100644 index 0000000..34efe40 --- /dev/null +++ b/src/burer_monteiro.jl @@ -0,0 +1,37 @@ +import NLPModels + +struct BurerMonteiro{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} + model::Model{T} + meta::NLPModels.NLPModelMeta{T,Vector{T}} + counters::NLPModels.Counters + function BurerMonteiro(model::Model{T}) where {T} + n = num_scalars(model) + sum(side_dimension(model, i) for i in matrix_indices(model); init = 0) + ncon = num_constraints(model) + return new( + ad, + NLPModels.NLPModelMeta( + n, #nvar + ncon = ncon, + nnzj = 0, + nnzh = 0, + x0 = rand(n), + y0 = rand(ncon), + lvar = fill(-Inf, n), + uvar = fill(Inf, n), + lcon = cons_constant(model), + ucon = cons_constant(model), + minimize = true, + ), + NLPModels.Counters(), + ) + end +end + +function NLPModels.obj(model::BurerMonteiro, x::AbstractVector) + return obj(model.model, x) +end + +function NLPModels.grad!(model::BurerMonteiro, x::AbstractVector, g::AbstractVector) + grad!(model.model, x, g) + return g +end diff --git a/src/model.jl b/src/model.jl new file mode 100644 index 0000000..b9f34de --- /dev/null +++ b/src/model.jl @@ -0,0 +1,201 @@ +# Adapted from Loraine.jl + +import SparseArrays +import LinearAlgebra +import MutableArithmetics as MA +import MathOptInterface as MOI + +""" + MyModel + +Model representing the problem: +```math +\\begin{aligned} +\\max {} & b^\\top y - b_\\text{const} +\\\\ +& \\sum_{j=1}^n y_j A_{i,j} \\preceq C_i +\\qquad +\\forall i \\in \\{1,\\ldots,\\text{nlmi}\\} +\\\\ +& C_\\text{lin}^\\top y \\le d_\\text{lin} +\\end{aligned} +``` +The fields of the `struct` as related to the arrays of the above formulation as follows: + +* The ``i``th PSD constraint is of size `msize[i] × msisze[i]` +* The matrix ``C_i`` is given by `C[i]`. +* The matrix ``A_{i,j}`` is given by `-A[i,j]`. +""" +mutable struct MyModel{T,A<:AbstractMatrix{T}} + C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} + A::Matrix{A} + b::Vector{T} + b_const::T + d_lin::SparseArrays.SparseVector{T, Int64} + C_lin::SparseArrays.SparseMatrixCSC{T, Int64} + msizes::Vector{Int64} + + function MyModel( + C::Vector{SparseArrays.SparseMatrixCSC{T,Int}}, + A::Matrix{AT}, + b::Vector{T}, + b_const::T, + d_lin::SparseArrays.SparseVector{T, Int64}, + C_lin::SparseArrays.SparseMatrixCSC{T, Int64}, + msizes::Vector{Int64}, + ) where {T,AT<:AbstractMatrix{T}} + + model = new{T,AT}() + model.C = C + model.A = A + model.b = b + model.b_const = b_const + model.d_lin = d_lin + model.C_lin = C_lin + model.msizes = msizes + return model + end +end + +struct ScalarIndex + value::Int64 +end + +num_scalars(model::MyModel) = length(model.d_lin) + +function scalar_indices(model::MyModel) + return MOI.Utilities.LazyMap{ScalarIndex}(ScalarIndex, Base.OneTo(num_scalars(model))) +end + +struct MatrixIndex + value::Int64 +end + +num_matrices(model::MyModel) = length(model.C) + +function matrix_indices(model::MyModel) + return MOI.Utilities.LazyMap{MatrixIndex}(MatrixIndex, Base.OneTo(num_matrices(model))) +end + +side_dimension(model::MyModel, i::MatrixIndex) = model.msizes[i.value] + +struct ConstraintIndex + value::Int64 +end +num_constraints(model::MyModel) = length(model.b) +function constraint_indices(model::MyModel) + return MOI.Utilities.LazyMap{ConstraintIndex}(ConstraintIndex, Base.OneTo(num_constraints(model))) +end + +# Should be only used with `norm` +jac(model::MyModel, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] +function norm_jac(model::MyModel{T}, i::MatrixIndex) where {T} + if isempty(model.A) + return zero(T) + end + return norm(model.A[i.value, :]) +end + +function obj(model::MyModel, X, i::MatrixIndex) + return -dot(model.C[i.value], X) +end + +function obj(model::MyModel, X, ::Type{MatrixIndex}) + result = zero(eltype(eltype(X))) + for mat_idx in matrix_indices(model) + result += obj(model, X[mat_idx.value], mat_idx) + end + return result +end + +function obj(model::MyModel, X_lin, ::Type{ScalarIndex}) + return -dot(model.d_lin, X_lin) +end + +function obj(model::MyModel, X_lin, X) + return model.b_const + obj(model, X, MatrixIndex) - dot(model.d_lin, X_lin) +end + +dual_obj(model::MyModel, y) = -dot(model.b, y) + model.b_const + +function jtprod(model::MyModel, ::Type{ScalarIndex}, y) + return -model.C_lin' * y +end + +function dual_cons(model::MyModel, ::Type{ScalarIndex}, y, S) + return model.d_lin - S + jtprod(model, ScalarIndex, y) +end + +function buffer_for_jtprod(model::MyModel) + if iszero(num_matrices(model)) + return + end + return map(Base.Fix1(buffer_for_jtprod, model), matrix_indices(model)) +end + +function buffer_for_jtprod(model::MyModel, mat_idx::MatrixIndex) + if iszero(num_constraints(model)) + return + end + # FIXME: at some point, switch to dense + return sum( + abs.(model.A[mat_idx.value, j]) + for j in 1:num_constraints(model) + ) +end + +function _add_mul!(A::SparseMatrixCSC, B::SparseMatrixCSC, α) + for col in axes(A, 2) + range_A = SparseArrays.nzrange(A, col) + it_A = iterate(range_A) + for k in SparseArrays.nzrange(B, col) + row_B = SparseArrays.rowvals(B)[k] + while SparseArrays.rowvals(A)[it_A[1]] < row_B + it_A = iterate(range_A, it_A[2]) + end + @assert row_B == SparseArrays.rowvals(A)[it_A[1]] + SparseArrays.nonzeros(A)[it_A[1]] += SparseArrays.nonzeros(B)[k] * α + end + end +end + +_zero!(A::SparseMatrixCSC) = fill!(SparseArrays.nonzeros(A), 0.0) + +function jtprod!(buffer, model::MyModel, mat_idx::MatrixIndex, y) + if iszero(num_constraints(model)) + return MA.Zero() + end + _zero!(buffer) + for j in eachindex(y) + _add_mul!(buffer, model.A[mat_idx.value, j], y[j]) + end + return buffer +end + +function dual_cons!(buffer, model::MyModel, mat_idx::MatrixIndex, y, S) + i = mat_idx.value + return jtprod!(buffer[i], model, mat_idx, y) + model.C[i] - S[i] +end + +objgrad(model::MyModel, ::Type{ScalarIndex}) = model.d_lin +objgrad(model::MyModel, i::MatrixIndex) = model.C[i.value] + +cons_constant(model::MyModel) = model.b + +function cons(model::MyModel, x, X) + return model.b - jprod(model, x, X) +end + +function jprod(model::MyModel, i::MatrixIndex, W) + return eltype(W)[ + -dot(model.A[i.value, j], W) for j in 1:num_constraints(model) + ] +end + +function jprod(model::MyModel, w, W) + h = model.C_lin * w + for i in matrix_indices(model) + h += jprod(model, i, W[i.value]) + end + return h +end diff --git a/src/schur.jl b/src/schur.jl new file mode 100644 index 0000000..6f401b1 --- /dev/null +++ b/src/schur.jl @@ -0,0 +1,173 @@ +# Adapted from Loraine.jl + +# Computes `⟨A * W, W * B⟩` for symmetric sparse matrices `A` and `B` +function _dot(A::SparseMatrixCSC, B::SparseMatrixCSC, W::AbstractMatrix) + @assert LinearAlgebra.checksquare(W) == LinearAlgebra.checksquare(A) == LinearAlgebra.checksquare(B) + # After these asserts, we know that `A`, `B` and `W` are square and + # have the same sizes so we can safely use `@inbounds` + result = zero(eltype(A)) + @inbounds for i in axes(A, 2) + nzA = nzrange(A, i) + if !isempty(nzA) + for j in axes(B, 2) + nzB = nzrange(B, j) + if !isempty(nzB) + AW = zero(result) + for k in nzA + AW += nonzeros(A)[k] * W[rowvals(A)[k], j] + end + WB = zero(result) + for k in nzB + WB += W[i, rowvals(B)[k]] * nonzeros(B)[k] + end + result += AW * WB + end + end + end + end + return result +end + +function buffer_for_schur_complement(model::MyModel, κ) + n = num_constraints(model) + σ = zeros(Int64, n, num_matrices(model)) + last_dense = zeros(Int64, num_matrices(model)) + + for mat_idx in matrix_indices(model) + i = mat_idx.value + nzA = [nnz(model.A[i, j]) for j in 1:n] + σ[:,i] = sortperm(nzA, rev = true) + sorted = nzA[σ[:,i]] + + last_dense[i] = something(findlast(Base.Fix1(isless, κ), sorted), 0) + end + + return σ, last_dense +end + +function makeBBBB_rank1(n,nlmi,B,G) + tmp = zeros(Float64, n, n) + BBBB = zeros(Float64, n, n) + for ilmi = 1:nlmi + BB = transpose(B[ilmi] * G[ilmi]) + mul!(tmp,BB',BB) + if ilmi == 1 + BBBB = tmp .^ 2 + else + BBBB += tmp .^ 2 + end + end + return BBBB +end + +######################### + +function schur_complement(buffer, model::MyModel, W, ::Type{MatrixIndex}) + n = num_constraints(model) + BBBB = zeros(eltype(eltype(W)), n, n) + for mat_idx in matrix_indices(model) + BBBB += schur_complement(buffer, model, mat_idx, W[mat_idx.value]) + end + return BBBB +end + +##### +function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T} + σ, last_dense = buffer + ilmi = mat_idx.value + n = num_constraints(model) + BBBB = zeros(T, n, n) + dim = side_dimension(model, mat_idx) + @assert dim == size(W, 1) == size(W, 2) + tmp1 = Matrix{T}(undef, size(W, 2), dim) + tmp = zeros(T, size(W, 2), dim) + + for ii = 1:n + i = σ[ii,ilmi] + Ai = model.A[ilmi, i] + if nnz(Ai) > 0 + if ii <= last_dense[ilmi] + mul!(tmp1, W, Ai) + mul!(tmp, tmp1, W) + tmp2 = jprod(model, mat_idx, tmp) + indi = σ[ii:end,ilmi] + BBBB[indi,i] .= -tmp2[indi] + BBBB[i,indi] .= -tmp2[indi] + else + if !iszero(nnz(Ai)) + if nnz(Ai) > 1 + @inbounds for jj = ii:n + j = σ[jj,ilmi] + Aj = model.A[ilmi, j] + if !iszero(nnz(Aj)) + ttt = _dot(Ai, Aj, W) + if i >= j + BBBB[i,j] = ttt + else + BBBB[j,i] = ttt + end + end + end + else + # A is symmetric + iiiiAi = jjjiAi = only(rowvals(Ai)) + vvvi = only(nonzeros(Ai)) + @inbounds for jj = ii:n + j = σ[jj,ilmi] + Ajjj = model.A[ilmi, j] + # As we sort the matrices in decreasing `nnz` order, + # the rest of matrices is either zero or have only + # one entry + if !iszero(nnz(Ajjj)) + iiijAj = jjjjAj = only(rowvals(Ajjj)) + vvvj = only(nonzeros(Ajjj)) + ttt = vvvi * W[iiiiAi,iiijAj] * W[jjjiAi,jjjjAj] * vvvj + if i >= j + BBBB[i,j] = ttt + else + BBBB[j,i] = ttt + end + end + end + end + end + end + end + end + return BBBB +end + +# [HKS24, (5b)] +# Returns the matrix equal to the sum, for each equation, of +# ⟨A_i, WA_jW⟩ +function schur_complement(buffer, model::MyModel, w, W::AbstractVector) + H = MA.Zero() + if num_matrices(model) > 0 + H = MA.add!!(H, schur_complement(buffer, model, W, MatrixIndex)) + end + if num_scalars(model) > 0 + H = MA.add!!(H, schur_complement(model, w, ScalarIndex)) + end + if H isa MA.Zero + n = num_constraints(model) + H = zeros(eltype(w), n, n) + end + return Hermitian(H, :L) +end + +function schur_complement(model::MyModel, w, ::Type{ScalarIndex}) + return model.C_lin * spdiagm(w) * model.C_lin' +end + +# [HKS24, (5b)] +# Returns the matrix equal to the sum, for each equation, of +# ⟨A_i, WA(y)W⟩ +function eval_schur_complement!(buffer, result, model::MyModel, w, W, y) + result .= 0.0 + for mat_idx in matrix_indices(model) + i = mat_idx.value + result .-= jprod(model, mat_idx, W[i] * jtprod!(buffer[i], model, mat_idx, y) * W[i]) + end + result .+= model.C_lin * (w .* (model.C_lin' * y)) + return result +end From ec17b21d79b4b18d70b3c363ba0d71a11b5fefce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 29 May 2025 23:34:18 +0200 Subject: [PATCH 02/48] Rename --- src/{burer_monteiro.jl => BurerMonteiro.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/{burer_monteiro.jl => BurerMonteiro.jl} (100%) diff --git a/src/burer_monteiro.jl b/src/BurerMonteiro.jl similarity index 100% rename from src/burer_monteiro.jl rename to src/BurerMonteiro.jl From be132e932dba6850507e42e2a172a1dd45f7e689 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 30 May 2025 11:01:38 +0200 Subject: [PATCH 03/48] Fixes --- Project.toml | 4 + examples/maxcut.jl | 9 +- .../Variable/bridges/DotProductsBridge.jl | 2 +- src/BurerMonteiro.jl | 62 +++- src/LowRankOpt.jl | 4 + src/MOI_wrapper.jl | 328 ++++++++++++++++++ src/Test/Test.jl | 10 + src/model.jl | 62 ++-- src/schur.jl | 12 +- test/BurerMonteiro.jl | 6 + 10 files changed, 450 insertions(+), 49 deletions(-) create mode 100644 src/MOI_wrapper.jl create mode 100644 test/BurerMonteiro.jl diff --git a/Project.toml b/Project.toml index ee62264..3c8511e 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,8 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] @@ -17,5 +19,7 @@ KrylovKit = "0.9.5" LinearAlgebra = "1.10" MathOptInterface = "1.40" MutableArithmetics = "1.6.4" +NLPModelsJuMP = "0.13.2" +SparseArrays = "1.11.0" Test = "1.10" julia = "1.10" diff --git a/examples/maxcut.jl b/examples/maxcut.jl index 4f97207..dc61d1e 100644 --- a/examples/maxcut.jl +++ b/examples/maxcut.jl @@ -7,9 +7,14 @@ function e_i(i, n) return ei end -function maxcut(weights, solver) +function maxcut_objective(weights) N = LinearAlgebra.checksquare(weights) L = Diagonal(weights * ones(N)) - weights + return L / 4 +end + +function maxcut(weights, solver) + N = LinearAlgebra.checksquare(weights) model = Model(solver) LRO.Bridges.add_all_bridges(backend(model).optimizer, Float64) cone = MOI.PositiveSemidefiniteConeTriangle(N) @@ -26,7 +31,7 @@ function maxcut(weights, solver) dot_prod_set[length(factors) .+ (1:MOI.dimension(cone))], SymmetricMatrixShape(N), ) - @objective(model, Max, dot(L, X) / 4) + @objective(model, Max, dot(maxcut_objective(weights), X)) @constraint(model, dot_prod .== 1) return model end diff --git a/src/Bridges/Variable/bridges/DotProductsBridge.jl b/src/Bridges/Variable/bridges/DotProductsBridge.jl index c22dfad..966c747 100644 --- a/src/Bridges/Variable/bridges/DotProductsBridge.jl +++ b/src/Bridges/Variable/bridges/DotProductsBridge.jl @@ -63,7 +63,7 @@ function MOI.Bridges.map_function( bridge.set.set, ) else - return scalars[i-length(bridge.set.vectors)] + return scalars[i.value-length(bridge.set.vectors)] end end diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 34efe40..b5fd4bd 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -1,19 +1,42 @@ +module BurerMonteiro + import NLPModels +import LowRankOpt as LRO + +struct Dimensions + num_scalars::Int64 + side_dimensions::Vector{Int64} + ranks::Vector{Int64} + offsets::Vector{Int64} +end -struct BurerMonteiro{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} +function Dimensions(model::LRO.Model, ranks) + side_dimensions = [LRO.side_dimension(model, i) for i in LRO.matrix_indices(model)] + return Dimensions( + LRO.num_scalars(model), + side_dimensions, + ranks, + [0; cumsum(side_dimensions .* ranks)], + ) +end + +Base.length(d::Dimensions) = d.side_dimensions[end] + +struct Model{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} model::Model{T} + d::Dimension meta::NLPModels.NLPModelMeta{T,Vector{T}} counters::NLPModels.Counters - function BurerMonteiro(model::Model{T}) where {T} - n = num_scalars(model) + sum(side_dimension(model, i) for i in matrix_indices(model); init = 0) + function Model(model::Model{T}, ranks) where {T} + d = Dimensions(model, ranks) + n = length(d) ncon = num_constraints(model) return new( ad, + d, NLPModels.NLPModelMeta( n, #nvar ncon = ncon, - nnzj = 0, - nnzh = 0, x0 = rand(n), y0 = rand(ncon), lvar = fill(-Inf, n), @@ -27,11 +50,32 @@ struct BurerMonteiro{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} end end -function NLPModels.obj(model::BurerMonteiro, x::AbstractVector) - return obj(model.model, x) +struct Solution{T,VT<:AbstractVector{T}} + x::VT + d::Dimensions +end + +function Base.getindex(s::Solution, ::Type{ScalarIndex}) + return view(s.x, Base.OneTo(s.d.num_scalars)) +end + +function Base.getindex(s::Solution, mi::MatrixIndex) + i = mi.value + U = reshape( + view(s.x, s.d.offsets[i]:(s.d.offsets[i+1] - 1)), + s.d.side_dimensions[i], + s.d.ranks[i], + ) + return LRO.positive_semidefinite_factorization(U) end -function NLPModels.grad!(model::BurerMonteiro, x::AbstractVector, g::AbstractVector) - grad!(model.model, x, g) +function NLPModels.obj(model::Model, x::AbstractVector) + return obj(model.model, Solution(x, model.d)) +end + +function NLPModels.grad!(model::Model, x::AbstractVector, g::AbstractVector) + grad!(model.model, Solution(x, model.d), Solution(g, model.d)) return g end + +end diff --git a/src/LowRankOpt.jl b/src/LowRankOpt.jl index b2e1e45..2547616 100644 --- a/src/LowRankOpt.jl +++ b/src/LowRankOpt.jl @@ -15,4 +15,8 @@ include("distance_to_set.jl") include("Test/Test.jl") include("Bridges/Bridges.jl") +include("model.jl") +include("schur.jl") +include("MOI_wrapper.jl") + end # module LowRankOpt diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl new file mode 100644 index 0000000..4deeaa0 --- /dev/null +++ b/src/MOI_wrapper.jl @@ -0,0 +1,328 @@ +import NLPModelsJuMP + +const VAF{T} = MOI.VectorAffineFunction{T} +const PSD = MOI.PositiveSemidefiniteConeTriangle +const NNG = MOI.Nonnegatives + +MOI.Utilities.@product_of_sets(NNGCones, NNG) + +MOI.Utilities.@product_of_sets(PSDCones, PSD) + +MOI.Utilities.@struct_of_constraints_by_set_types(PSDOrNot, PSD, NNG) + +const OptimizerCache{T} = MOI.Utilities.GenericModel{ + T, + MOI.Utilities.ObjectiveContainer{T}, + MOI.Utilities.VariablesContainer{T}, + PSDOrNot{T}{ + MOI.Utilities.MatrixOfConstraints{ + T, + MOI.Utilities.MutableSparseMatrixCSC{ + T, + Int64, + MOI.Utilities.OneBasedIndexing, + }, + Vector{T}, + PSDCones{T}, + }, + MOI.Utilities.MatrixOfConstraints{ + T, + MOI.Utilities.MutableSparseMatrixCSC{ + T, + Int64, + MOI.Utilities.OneBasedIndexing, + }, + Vector{T}, + NNGCones{T}, + }, + }, +} + +mutable struct Optimizer{T} <: MOI.AbstractOptimizer + model::Union{Nothing,Model{T}} + lmi_id::Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64} + lin_cones::Union{Nothing,NNGCones{T}} + max_sense::Bool + objective_constant::T + silent::Bool + options::Dict{String,Any} + + function Optimizer{T}() where {T} + return new{T}( + nothing, + Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64}(), + nothing, + false, + 0.0, + false, + Dict{String,Any}(), + ) + end +end + +Optimizer() = Optimizer{Float64}() + +function MOI.default_cache(::Optimizer, ::Type{T}) where {T} + return MOI.Utilities.UniversalFallback(OptimizerCache{T}()) +end + +MOI.is_empty(optimizer::Optimizer) = isnothing(optimizer.model) + +function MOI.empty!(optimizer::Optimizer) + optimizer.model = nothing + optimizer.lin_cones = nothing + return +end + +MOI.get(::Optimizer, ::MOI.SolverName) = "Loraine" + +# MOI.RawOptimizerAttribute + +function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute) + return haskey(Solvers.DEFAULT_OPTIONS, param.name) +end + +function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value) + if !MOI.supports(optimizer, param) + throw(MOI.UnsupportedAttribute(param)) + end + optimizer.options[param.name] = value + if !isnothing(optimizer.solver) + setproperty!(optimizer.solver, Symbol(param.name), value) + end + return +end + +function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute) + if !MOI.supports(optimizer, param) + throw(MOI.UnsupportedAttribute(param)) + end + return optimizer.options[param.name] +end + +# MOI.Silent + +MOI.supports(::Optimizer, ::MOI.Silent) = true + +function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool) + optimizer.silent = value + return +end + +MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent + +function MOI.set(optimizer::Optimizer, ::MOI.ObjectiveSense, value::Bool) + optimizer.max_sense = value + return +end + +# MOI.supports + +function MOI.supports( + ::Optimizer, + ::Union{MOI.ObjectiveSense,MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}}, +) where {T} + return true +end + +const SUPPORTED_CONES = Union{NNG,PSD} + +function MOI.supports_constraint(::Optimizer{T}, ::Type{VAF{T}}, ::Type{<:SUPPORTED_CONES}) where {T} + return true +end + +function MOI.optimize!(model::Optimizer) + options = Dict{Symbol, Any}( + Symbol(key) => model.options[key] for key in keys(model.options) if key != "solver" + ) + if model.silent + options[:verbose] = 0 + else + options[:verbose] = 1 + end + model.stats = SolverCore.GenericExecutionStats(model.nlp) + SolverCore.solve!(model.solver, model.model, model.stats; options...) + return +end + +function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} + MOI.empty!(dest) + psd_AC = MOI.Utilities.constraints(src.constraints, VAF{T}, PSD) + Cd_lin = MOI.Utilities.constraints(src.constraints, VAF{T}, NNG) + SM = SparseArrays.SparseMatrixCSC{T,Int64} + psd_A = convert(SM, psd_AC.coefficients) + C_lin = convert(SM, Cd_lin.coefficients) + C_lin = -convert(SM, C_lin') + n = MOI.get(src, MOI.NumberOfVariables()) + nlmi = MOI.get(src, MOI.NumberOfConstraints{VAF{T},PSD}()) + A = Matrix{Tuple{Vector{Int64},Vector{Int64},Vector{T},Int64,Int64}}(undef, nlmi, n + 1) + back = Vector{Tuple{Int64,Int64,Int64}}(undef, size(psd_A, 1)) + empty!(dest.lmi_id) + row = 0 + msizes = Int64[] + for (lmi_id, ci) in enumerate(MOI.get(src, MOI.ListOfConstraintIndices{VAF{T},PSD}())) + dest.lmi_id[ci] = lmi_id + set = MOI.get(src, MOI.ConstraintSet(), ci) + d = set.side_dimension + push!(msizes, d) + for k = 1:(n+1) + A[lmi_id, k] = (Int64[], Int64[], T[], d, d) + end + for j = 1:d + for i = 1:j + row += 1 + back[row] = (lmi_id, i, j) + end + end + end + function __add(lmi_id, k, i, j, v) + I, J, V, _, _ = A[lmi_id, k] + push!(I, i) + push!(J, j) + push!(V, v) + return + end + function _add(lmi_id, k, i, j, coef) + __add(lmi_id, k, i, j, coef) + if i != j + __add(lmi_id, k, j, i, coef) + end + return + end + for row in eachindex(back) + lmi_id, i, j = back[row] + _add(lmi_id, 1, i, j, -psd_AC.constants[row]) + end + for var = 1:n + for k in SparseArrays.nzrange(psd_A, var) + lmi_id, i, j = back[SparseArrays.rowvals(psd_A)[k]] + col = 1 + var + _add(lmi_id, col, i, j, SparseArrays.nonzeros(psd_A)[k]) + end + end + dest.max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE + obj = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}()) + # objective_constant = MOI.constant(obj) # TODO # MK: done(?) + b_const = obj.constant + b_const = dest.max_sense ? -b_const : b_const + b0 = zeros(T, n) + for term in obj.terms + b0[term.variable.value] += term.coefficient + end + b = dest.max_sense ? b0 : -b0 + # b = max_sense ? -b0 : b0 + + AA = SparseArrays.SparseMatrixCSC{T,Int}[SparseArrays.sparse(IJV...) for IJV in A] + dest.model = Model( + -AA[:,1], + AA[:,2:end], + b, + b_const, + convert(SparseArrays.SparseVector{T,Int64}, SparseArrays.sparsevec(Cd_lin.constants)), + C_lin, + msizes, + ) + # FIXME this does not work if an option is changed between `MOI.copy_to` and `MOI.optimize!` + options = copy(dest.options) + if dest.silent + options["verb"] = 0 + end + dest.lin_cones = Cd_lin.sets + dest.solver = dest.options["solver"](dest.nlp) + return MOI.Utilities.identity_index_map(src) +end + +function MOI.copy_to(dest::Optimizer{T}, src::MOI.ModelLike) where {T} + cache = OptimizerCache{T}() + index_map = MOI.copy_to(cache, src) + MOI.copy_to(dest, cache) + return index_map +end + +function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) + return optimizer.stats.elapsed_time +end + +function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString) + return SolverCore.STATUSES[optimizer.stats.status] +end + +struct RawStatus <: MOI.AbstractModelAttribute + name::Symbol +end + +MOI.is_set_by_optimize(::RawStatus) = true + +function MOI.get(optimizer::Optimizer, attr::RawStatus) + return getfield(optimizer.stats, attr.name) +end + +function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) + if isnothing(optimizer.stats) + return MOI.OPTIMIZE_NOT_CALLED + end + return NLPModelsJuMP.TERMINATION_STATUS[optimizer.stats.status] +end + +function MOI.get(model::Optimizer, ::MOI.ResultCount) + if MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED + return 0 + else + return 1 + end +end + +function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) + if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) + return MOI.NO_SOLUTION + elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED + return MOI.FEASIBLE_POINT + else + # TODO + return MOI.UNKNOWN_RESULT_STATUS + end +end + +function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} + MOI.check_result_index_bounds(optimizer, attr) + val = optimizer.stats.objective + return optimizer.max_sense ? -val : val +end + +function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) + MOI.check_result_index_bounds(optimizer, attr) + return optimizer.stats.solution[vi.value] +end + +function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} + MOI.check_result_index_bounds(optimizer, attr) + val = Solvers.obj(optimizer.solver.model, optimizer.solver.X_lin, optimizer.solver.X)::T + return optimizer.max_sense ? -val : val +end + +function MOI.get(::Optimizer, ::MOI.DualStatus) + # TODO + return MOI.NO_SOLUTION +end + +function MOI.get( + optimizer::Optimizer{T}, + attr::MOI.ConstraintDual, + ci::MOI.ConstraintIndex{VAF{T},PSD}, +) where {T} + MOI.check_result_index_bounds(optimizer, attr) + lmi_id = optimizer.lmi_id[ci] + X = optimizer.solver.X[lmi_id] + n = optimizer.solver.model.msizes[lmi_id] + return [X[i, j] for j = 1:n for i = 1:j]::Vector{T} +end + +function MOI.get( + optimizer::Optimizer{T}, + attr::MOI.ConstraintDual, + ci::MOI.ConstraintIndex{VAF{T},NNG}, +) where {T} + MOI.check_result_index_bounds(optimizer, attr) + rows = MOI.Utilities.rows(optimizer.lin_cones, ci) + return optimizer.solver.X_lin[rows]::Vector{T} +end diff --git a/src/Test/Test.jl b/src/Test/Test.jl index afaacca..52d5440 100644 --- a/src/Test/Test.jl +++ b/src/Test/Test.jl @@ -43,6 +43,11 @@ function test_conic_PositiveSemidefinite_RankOne_polynomial( LRO.positive_semidefinite_factorization(T[1, 1]), ]), ) + @show MOI.supports_constraint( + model, + MOI.VectorAffineFunction{T}, + typeof(set), + ) MOI.Test.@requires MOI.supports_constraint( model, MOI.VectorAffineFunction{T}, @@ -70,6 +75,7 @@ function test_conic_PositiveSemidefinite_RankOne_polynomial( if MOI.Test._supports(config, MOI.ConstraintDual) @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT end + @show MOI.get(model, MOI.ObjectiveValue()) @test ≈(MOI.get(model, MOI.ObjectiveValue()), T(-1), config) if MOI.Test._supports(config, MOI.DualObjectiveValue) @test ≈(MOI.get(model, MOI.DualObjectiveValue()), T(-1), config) @@ -128,6 +134,10 @@ function test_conic_PositiveSemidefinite_RankOne_moment( LRO.positive_semidefinite_factorization(T[1, 1]), ]), ) + @show MOI.supports_add_constrained_variables( + model, + typeof(set), + ) MOI.Test.@requires MOI.supports_add_constrained_variables( model, typeof(set), diff --git a/src/model.jl b/src/model.jl index b9f34de..0351635 100644 --- a/src/model.jl +++ b/src/model.jl @@ -6,7 +6,7 @@ import MutableArithmetics as MA import MathOptInterface as MOI """ - MyModel + Model Model representing the problem: ```math @@ -26,7 +26,7 @@ The fields of the `struct` as related to the arrays of the above formulation as * The matrix ``C_i`` is given by `C[i]`. * The matrix ``A_{i,j}`` is given by `-A[i,j]`. """ -mutable struct MyModel{T,A<:AbstractMatrix{T}} +mutable struct Model{T,A<:AbstractMatrix{T}} C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} A::Matrix{A} b::Vector{T} @@ -35,7 +35,7 @@ mutable struct MyModel{T,A<:AbstractMatrix{T}} C_lin::SparseArrays.SparseMatrixCSC{T, Int64} msizes::Vector{Int64} - function MyModel( + function Model( C::Vector{SparseArrays.SparseMatrixCSC{T,Int}}, A::Matrix{AT}, b::Vector{T}, @@ -61,9 +61,9 @@ struct ScalarIndex value::Int64 end -num_scalars(model::MyModel) = length(model.d_lin) +num_scalars(model::Model) = length(model.d_lin) -function scalar_indices(model::MyModel) +function scalar_indices(model::Model) return MOI.Utilities.LazyMap{ScalarIndex}(ScalarIndex, Base.OneTo(num_scalars(model))) end @@ -71,36 +71,36 @@ struct MatrixIndex value::Int64 end -num_matrices(model::MyModel) = length(model.C) +num_matrices(model::Model) = length(model.C) -function matrix_indices(model::MyModel) +function matrix_indices(model::Model) return MOI.Utilities.LazyMap{MatrixIndex}(MatrixIndex, Base.OneTo(num_matrices(model))) end -side_dimension(model::MyModel, i::MatrixIndex) = model.msizes[i.value] +side_dimension(model::Model, i::MatrixIndex) = model.msizes[i.value] struct ConstraintIndex value::Int64 end -num_constraints(model::MyModel) = length(model.b) -function constraint_indices(model::MyModel) +num_constraints(model::Model) = length(model.b) +function constraint_indices(model::Model) return MOI.Utilities.LazyMap{ConstraintIndex}(ConstraintIndex, Base.OneTo(num_constraints(model))) end # Should be only used with `norm` -jac(model::MyModel, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] -function norm_jac(model::MyModel{T}, i::MatrixIndex) where {T} +jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] +function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) return zero(T) end return norm(model.A[i.value, :]) end -function obj(model::MyModel, X, i::MatrixIndex) +function obj(model::Model, X, i::MatrixIndex) return -dot(model.C[i.value], X) end -function obj(model::MyModel, X, ::Type{MatrixIndex}) +function obj(model::Model, X, ::Type{MatrixIndex}) result = zero(eltype(eltype(X))) for mat_idx in matrix_indices(model) result += obj(model, X[mat_idx.value], mat_idx) @@ -108,32 +108,32 @@ function obj(model::MyModel, X, ::Type{MatrixIndex}) return result end -function obj(model::MyModel, X_lin, ::Type{ScalarIndex}) +function obj(model::Model, X_lin, ::Type{ScalarIndex}) return -dot(model.d_lin, X_lin) end -function obj(model::MyModel, X_lin, X) +function obj(model::Model, X_lin, X) return model.b_const + obj(model, X, MatrixIndex) - dot(model.d_lin, X_lin) end -dual_obj(model::MyModel, y) = -dot(model.b, y) + model.b_const +dual_obj(model::Model, y) = -dot(model.b, y) + model.b_const -function jtprod(model::MyModel, ::Type{ScalarIndex}, y) +function jtprod(model::Model, ::Type{ScalarIndex}, y) return -model.C_lin' * y end -function dual_cons(model::MyModel, ::Type{ScalarIndex}, y, S) +function dual_cons(model::Model, ::Type{ScalarIndex}, y, S) return model.d_lin - S + jtprod(model, ScalarIndex, y) end -function buffer_for_jtprod(model::MyModel) +function buffer_for_jtprod(model::Model) if iszero(num_matrices(model)) return end return map(Base.Fix1(buffer_for_jtprod, model), matrix_indices(model)) end -function buffer_for_jtprod(model::MyModel, mat_idx::MatrixIndex) +function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) if iszero(num_constraints(model)) return end @@ -144,7 +144,7 @@ function buffer_for_jtprod(model::MyModel, mat_idx::MatrixIndex) ) end -function _add_mul!(A::SparseMatrixCSC, B::SparseMatrixCSC, α) +function _add_mul!(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, α) for col in axes(A, 2) range_A = SparseArrays.nzrange(A, col) it_A = iterate(range_A) @@ -159,9 +159,9 @@ function _add_mul!(A::SparseMatrixCSC, B::SparseMatrixCSC, α) end end -_zero!(A::SparseMatrixCSC) = fill!(SparseArrays.nonzeros(A), 0.0) +_zero!(A::SparseArrays.SparseMatrixCSC) = fill!(SparseArrays.nonzeros(A), 0.0) -function jtprod!(buffer, model::MyModel, mat_idx::MatrixIndex, y) +function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) if iszero(num_constraints(model)) return MA.Zero() end @@ -172,27 +172,27 @@ function jtprod!(buffer, model::MyModel, mat_idx::MatrixIndex, y) return buffer end -function dual_cons!(buffer, model::MyModel, mat_idx::MatrixIndex, y, S) +function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y, S) i = mat_idx.value return jtprod!(buffer[i], model, mat_idx, y) + model.C[i] - S[i] end -objgrad(model::MyModel, ::Type{ScalarIndex}) = model.d_lin -objgrad(model::MyModel, i::MatrixIndex) = model.C[i.value] +objgrad(model::Model, ::Type{ScalarIndex}) = model.d_lin +objgrad(model::Model, i::MatrixIndex) = model.C[i.value] -cons_constant(model::MyModel) = model.b +cons_constant(model::Model) = model.b -function cons(model::MyModel, x, X) +function cons(model::Model, x, X) return model.b - jprod(model, x, X) end -function jprod(model::MyModel, i::MatrixIndex, W) +function jprod(model::Model, i::MatrixIndex, W) return eltype(W)[ -dot(model.A[i.value, j], W) for j in 1:num_constraints(model) ] end -function jprod(model::MyModel, w, W) +function jprod(model::Model, w, W) h = model.C_lin * w for i in matrix_indices(model) h += jprod(model, i, W[i.value]) diff --git a/src/schur.jl b/src/schur.jl index 6f401b1..e9ba5c8 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -1,7 +1,7 @@ # Adapted from Loraine.jl # Computes `⟨A * W, W * B⟩` for symmetric sparse matrices `A` and `B` -function _dot(A::SparseMatrixCSC, B::SparseMatrixCSC, W::AbstractMatrix) +function _dot(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, W::AbstractMatrix) @assert LinearAlgebra.checksquare(W) == LinearAlgebra.checksquare(A) == LinearAlgebra.checksquare(B) # After these asserts, we know that `A`, `B` and `W` are square and # have the same sizes so we can safely use `@inbounds` @@ -28,7 +28,7 @@ function _dot(A::SparseMatrixCSC, B::SparseMatrixCSC, W::AbstractMatrix) return result end -function buffer_for_schur_complement(model::MyModel, κ) +function buffer_for_schur_complement(model::Model, κ) n = num_constraints(model) σ = zeros(Int64, n, num_matrices(model)) last_dense = zeros(Int64, num_matrices(model)) @@ -62,7 +62,7 @@ end ######################### -function schur_complement(buffer, model::MyModel, W, ::Type{MatrixIndex}) +function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) n = num_constraints(model) BBBB = zeros(eltype(eltype(W)), n, n) for mat_idx in matrix_indices(model) @@ -140,7 +140,7 @@ end # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA_jW⟩ -function schur_complement(buffer, model::MyModel, w, W::AbstractVector) +function schur_complement(buffer, model::Model, w, W::AbstractVector) H = MA.Zero() if num_matrices(model) > 0 H = MA.add!!(H, schur_complement(buffer, model, W, MatrixIndex)) @@ -155,14 +155,14 @@ function schur_complement(buffer, model::MyModel, w, W::AbstractVector) return Hermitian(H, :L) end -function schur_complement(model::MyModel, w, ::Type{ScalarIndex}) +function schur_complement(model::Model, w, ::Type{ScalarIndex}) return model.C_lin * spdiagm(w) * model.C_lin' end # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA(y)W⟩ -function eval_schur_complement!(buffer, result, model::MyModel, w, W, y) +function eval_schur_complement!(buffer, result, model::Model, w, W, y) result .= 0.0 for mat_idx in matrix_indices(model) i = mat_idx.value diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl new file mode 100644 index 0000000..335a580 --- /dev/null +++ b/test/BurerMonteiro.jl @@ -0,0 +1,6 @@ +using Test +using LowRankOpt +include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) +weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; +model = maxcut(weights, LowRankOpt.Optimizer) +optimize!(model) From 69353e19ec1af0e9df668f5af0f80d46733a54b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Jun 2025 15:58:19 +0200 Subject: [PATCH 04/48] Fixes --- Project.toml | 4 ++ src/BurerMonteiro.jl | 111 +++++++++++++++++++++++++++++++----------- src/LowRankOpt.jl | 1 + src/MOI_wrapper.jl | 21 ++++++-- src/Test/Test.jl | 10 ---- src/model.jl | 87 ++++++++++++++++++++++----------- test/BurerMonteiro.jl | 5 ++ test/Project.toml | 1 + 8 files changed, 169 insertions(+), 71 deletions(-) diff --git a/Project.toml b/Project.toml index 3c8511e..9bb671f 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,9 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" @@ -19,7 +21,9 @@ KrylovKit = "0.9.5" LinearAlgebra = "1.10" MathOptInterface = "1.40" MutableArithmetics = "1.6.4" +NLPModels = "0.21.5" NLPModelsJuMP = "0.13.2" +SolverCore = "0.3.8" SparseArrays = "1.11.0" Test = "1.10" julia = "1.10" diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index b5fd4bd..40f229b 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -1,5 +1,7 @@ module BurerMonteiro +import LinearAlgebra +import SolverCore import NLPModels import LowRankOpt as LRO @@ -12,28 +14,25 @@ end function Dimensions(model::LRO.Model, ranks) side_dimensions = [LRO.side_dimension(model, i) for i in LRO.matrix_indices(model)] - return Dimensions( - LRO.num_scalars(model), - side_dimensions, - ranks, - [0; cumsum(side_dimensions .* ranks)], - ) + num_scalars = LRO.num_scalars(model) + offsets = num_scalars .+ [0; cumsum(side_dimensions .* ranks)] + return Dimensions(num_scalars, side_dimensions, ranks, offsets) end -Base.length(d::Dimensions) = d.side_dimensions[end] +Base.length(d::Dimensions) = d.offsets[end] -struct Model{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} - model::Model{T} - d::Dimension +struct Model{T,AT} <: NLPModels.AbstractNLPModel{T,Vector{T}} + model::LRO.Model{T,AT} + dim::Dimensions meta::NLPModels.NLPModelMeta{T,Vector{T}} counters::NLPModels.Counters - function Model(model::Model{T}, ranks) where {T} - d = Dimensions(model, ranks) - n = length(d) - ncon = num_constraints(model) - return new( - ad, - d, + function Model(model::LRO.Model{T,AT}, ranks) where {T,AT} + dim = Dimensions(model, ranks) + n = length(dim) + ncon = LRO.num_constraints(model) + return new{T,AT}( + model, + dim, NLPModels.NLPModelMeta( n, #nvar ncon = ncon, @@ -41,8 +40,8 @@ struct Model{T} <: NLPModels.AbstractNLPModel{T,Vector{T}} y0 = rand(ncon), lvar = fill(-Inf, n), uvar = fill(Inf, n), - lcon = cons_constant(model), - ucon = cons_constant(model), + lcon = LRO.cons_constant(model), + ucon = LRO.cons_constant(model), minimize = true, ), NLPModels.Counters(), @@ -52,30 +51,86 @@ end struct Solution{T,VT<:AbstractVector{T}} x::VT - d::Dimensions + dim::Dimensions end -function Base.getindex(s::Solution, ::Type{ScalarIndex}) - return view(s.x, Base.OneTo(s.d.num_scalars)) +Base.eltype(::Type{<:Solution{T}}) where {T} = T +Base.eltype(x::Solution) = eltype(typeof(x)) + +function Base.getindex(s::Solution, ::Type{LRO.ScalarIndex}) + return view(s.x, Base.OneTo(s.dim.num_scalars)) end -function Base.getindex(s::Solution, mi::MatrixIndex) +function Base.getindex(s::Solution, mi::LRO.MatrixIndex) i = mi.value U = reshape( - view(s.x, s.d.offsets[i]:(s.d.offsets[i+1] - 1)), - s.d.side_dimensions[i], - s.d.ranks[i], + view(s.x, (1 + s.dim.offsets[i]):s.dim.offsets[i+1]), + s.dim.side_dimensions[i], + s.dim.ranks[i], ) return LRO.positive_semidefinite_factorization(U) end function NLPModels.obj(model::Model, x::AbstractVector) - return obj(model.model, Solution(x, model.d)) + return NLPModels.obj(model.model, Solution(x, model.dim)) end function NLPModels.grad!(model::Model, x::AbstractVector, g::AbstractVector) - grad!(model.model, Solution(x, model.d), Solution(g, model.d)) + X = Solution(x, model.dim) + G = Solution(g, model.dim) + copyto!(G[LRO.ScalarIndex], NLPModels.grad(model.model, LRO.ScalarIndex)) + for i in LRO.matrix_indices(model.model) + C = NLPModels.grad(model.model, i) + LinearAlgebra.mul!(G[i].factor, C, X[i].factor) + G[i].factor ./= 2 + end return g end +function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) + NLPModels.cons!(model.model, Solution(x, model.dim), cx) +end + +function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, Jtv::AbstractVector) + X = Solution(x, model.dim) + JtV = Solution(Jtv, model.dim) + LinearAlgebra.mul!( + JtV[LRO.ScalarIndex], + NLPModels.jac(model.model, LRO.ScalarIndex)', + y, + ) + for i in LRO.matrix_indices(model.model) + U = JtV[i].factor + for j in eachindex(y) + A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) + LinearAlgebra.mul!(U, A, X[i].factor) + U ./= 2 + end + end + return Jtv +end + +struct Solver{T,ST} <: SolverCore.AbstractOptimizationSolver + model::Model{T} + solver::ST + stats::SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any} +end + +function Solver(src::LRO.Model; sub_solver, ranks, kws...) + model = Model(src, ranks) + solver = sub_solver(model; kws...) + stats = SolverCore.GenericExecutionStats(model) + return Solver(model, solver, stats) +end + +function SolverCore.solve!( + solver::Solver, + model::NLPModels.AbstractNLPModel, # Same as `solver.model.model` + stats::SolverCore.GenericExecutionStats; + kws..., +) + SolverCore.solve!(solver.solver, solver.model, solver.stats; kws...) +end + + end diff --git a/src/LowRankOpt.jl b/src/LowRankOpt.jl index 2547616..3f1e257 100644 --- a/src/LowRankOpt.jl +++ b/src/LowRankOpt.jl @@ -18,5 +18,6 @@ include("Bridges/Bridges.jl") include("model.jl") include("schur.jl") include("MOI_wrapper.jl") +include("BurerMonteiro.jl") end # module LowRankOpt diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 4deeaa0..63a860a 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -1,3 +1,4 @@ +import SolverCore import NLPModelsJuMP const VAF{T} = MOI.VectorAffineFunction{T} @@ -39,7 +40,12 @@ const OptimizerCache{T} = MOI.Utilities.GenericModel{ } mutable struct Optimizer{T} <: MOI.AbstractOptimizer + solver::Union{Nothing,SolverCore.AbstractOptimizationSolver} model::Union{Nothing,Model{T}} + stats::Union{ + Nothing, + SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any}, + } lmi_id::Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64} lin_cones::Union{Nothing,NNGCones{T}} max_sense::Bool @@ -49,6 +55,8 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer function Optimizer{T}() where {T} return new{T}( + nothing, + nothing, nothing, Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64}(), nothing, @@ -79,7 +87,7 @@ MOI.get(::Optimizer, ::MOI.SolverName) = "Loraine" # MOI.RawOptimizerAttribute function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute) - return haskey(Solvers.DEFAULT_OPTIONS, param.name) + return true end function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value) @@ -131,16 +139,18 @@ function MOI.supports_constraint(::Optimizer{T}, ::Type{VAF{T}}, ::Type{<:SUPPOR return true end +SOLVER_OPTIONS = ["solver", "sub_solver", "ranks"] + function MOI.optimize!(model::Optimizer) options = Dict{Symbol, Any}( - Symbol(key) => model.options[key] for key in keys(model.options) if key != "solver" + Symbol(key) => model.options[key] for key in keys(model.options) if !(key in SOLVER_OPTIONS) ) if model.silent options[:verbose] = 0 else options[:verbose] = 1 end - model.stats = SolverCore.GenericExecutionStats(model.nlp) + model.stats = SolverCore.GenericExecutionStats(model.model) SolverCore.solve!(model.solver, model.model, model.stats; options...) return end @@ -228,7 +238,10 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} options["verb"] = 0 end dest.lin_cones = Cd_lin.sets - dest.solver = dest.options["solver"](dest.nlp) + options = Dict{Symbol, Any}( + Symbol(key) => dest.options[key] for key in keys(dest.options) if key in SOLVER_OPTIONS && key != "solver" + ) + dest.solver = dest.options["solver"](dest.model; options...) return MOI.Utilities.identity_index_map(src) end diff --git a/src/Test/Test.jl b/src/Test/Test.jl index 52d5440..afaacca 100644 --- a/src/Test/Test.jl +++ b/src/Test/Test.jl @@ -43,11 +43,6 @@ function test_conic_PositiveSemidefinite_RankOne_polynomial( LRO.positive_semidefinite_factorization(T[1, 1]), ]), ) - @show MOI.supports_constraint( - model, - MOI.VectorAffineFunction{T}, - typeof(set), - ) MOI.Test.@requires MOI.supports_constraint( model, MOI.VectorAffineFunction{T}, @@ -75,7 +70,6 @@ function test_conic_PositiveSemidefinite_RankOne_polynomial( if MOI.Test._supports(config, MOI.ConstraintDual) @test MOI.get(model, MOI.DualStatus()) == MOI.FEASIBLE_POINT end - @show MOI.get(model, MOI.ObjectiveValue()) @test ≈(MOI.get(model, MOI.ObjectiveValue()), T(-1), config) if MOI.Test._supports(config, MOI.DualObjectiveValue) @test ≈(MOI.get(model, MOI.DualObjectiveValue()), T(-1), config) @@ -134,10 +128,6 @@ function test_conic_PositiveSemidefinite_RankOne_moment( LRO.positive_semidefinite_factorization(T[1, 1]), ]), ) - @show MOI.supports_add_constrained_variables( - model, - typeof(set), - ) MOI.Test.@requires MOI.supports_add_constrained_variables( model, typeof(set), diff --git a/src/model.jl b/src/model.jl index 0351635..80de35d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -4,6 +4,7 @@ import SparseArrays import LinearAlgebra import MutableArithmetics as MA import MathOptInterface as MOI +import NLPModels """ Model @@ -26,13 +27,14 @@ The fields of the `struct` as related to the arrays of the above formulation as * The matrix ``C_i`` is given by `C[i]`. * The matrix ``A_{i,j}`` is given by `-A[i,j]`. """ -mutable struct Model{T,A<:AbstractMatrix{T}} +mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vector{T}} + meta::NLPModels.NLPModelMeta{T,Vector{T}} C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} A::Matrix{A} b::Vector{T} b_const::T - d_lin::SparseArrays.SparseVector{T, Int64} - C_lin::SparseArrays.SparseMatrixCSC{T, Int64} + d_lin::SparseArrays.SparseVector{T,Int64} + C_lin::SparseArrays.SparseMatrixCSC{T,Int64} msizes::Vector{Int64} function Model( @@ -40,11 +42,10 @@ mutable struct Model{T,A<:AbstractMatrix{T}} A::Matrix{AT}, b::Vector{T}, b_const::T, - d_lin::SparseArrays.SparseVector{T, Int64}, - C_lin::SparseArrays.SparseMatrixCSC{T, Int64}, + d_lin::SparseArrays.SparseVector{T,Int64}, + C_lin::SparseArrays.SparseMatrixCSC{T,Int64}, msizes::Vector{Int64}, ) where {T,AT<:AbstractMatrix{T}} - model = new{T,AT}() model.C = C model.A = A @@ -53,10 +54,25 @@ mutable struct Model{T,A<:AbstractMatrix{T}} model.d_lin = d_lin model.C_lin = C_lin model.msizes = msizes + model.meta = NLPModels.NLPModelMeta( + num_scalars(model) + sum( + Base.Fix1(side_dimension, model), + matrix_indices(model); + init = 0 + ), + ) return model end end +function NLPModels.unconstrained(model::Model) + return iszero(num_constraints(model)) +end + +# TODO the scalar actually have lower bounds and the SDP variables too +# but these are not box constraints +NLPModels.has_bounds(::Model) = false + struct ScalarIndex value::Int64 end @@ -88,7 +104,9 @@ function constraint_indices(model::Model) end # Should be only used with `norm` -jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] +NLPModels.jac(model::Model, ::Type{ScalarIndex}) = model.C_lin +NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = model.A[j.value, i.value] +NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) return zero(T) @@ -96,24 +114,32 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} return norm(model.A[i.value, :]) end -function obj(model::Model, X, i::MatrixIndex) - return -dot(model.C[i.value], X) +function NLPModels.obj(model::Model, X, i::MatrixIndex) + return -LinearAlgebra.dot(model.C[i.value], X) end -function obj(model::Model, X, ::Type{MatrixIndex}) - result = zero(eltype(eltype(X))) - for mat_idx in matrix_indices(model) - result += obj(model, X[mat_idx.value], mat_idx) +function NLPModels.obj(model::Model, x, ::Type{MatrixIndex}) + result = zero(eltype(x)) + for i in matrix_indices(model) + result += NLPModels.obj(model, x[i], i) end return result end -function obj(model::Model, X_lin, ::Type{ScalarIndex}) - return -dot(model.d_lin, X_lin) +function NLPModels.obj(model::Model, x, ::Type{ScalarIndex}) + return -LinearAlgebra.dot(model.d_lin, x[ScalarIndex]) end -function obj(model::Model, X_lin, X) - return model.b_const + obj(model, X, MatrixIndex) - dot(model.d_lin, X_lin) +function NLPModels.obj(model::Model, x) + return model.b_const + NLPModels.obj(model, x, MatrixIndex) + NLPModels.obj(model, x, ScalarIndex) +end + +function NLPModels.grad!(model::Model, _, g) + copyto!(g[ScalarIndex], model.d_lin) + for i in matrix_indices(model) + copyto!(g[i], model.C[i.value]) + end + return g end dual_obj(model::Model, y) = -dot(model.b, y) + model.b_const @@ -177,25 +203,28 @@ function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y, S) return jtprod!(buffer[i], model, mat_idx, y) + model.C[i] - S[i] end -objgrad(model::Model, ::Type{ScalarIndex}) = model.d_lin -objgrad(model::Model, i::MatrixIndex) = model.C[i.value] +NLPModels.grad(model::Model, ::Type{ScalarIndex}) = model.d_lin +NLPModels.grad(model::Model, i::MatrixIndex) = model.C[i.value] cons_constant(model::Model) = model.b -function cons(model::Model, x, X) - return model.b - jprod(model, x, X) +function NLPModels.cons!(model::Model, x, cx) + NLPModels.jprod!(model, x, x, cx) + cx .*= -1 + cx .+= model.b + return cx end -function jprod(model::Model, i::MatrixIndex, W) - return eltype(W)[ - -dot(model.A[i.value, j], W) for j in 1:num_constraints(model) - ] +function add_jprod!(model::Model, i::MatrixIndex, V, Jv) + for j in 1:num_constraints(model) + Jv[j] -= LinearAlgebra.dot(model.A[i.value, j], V) + end end -function jprod(model::Model, w, W) - h = model.C_lin * w +function NLPModels.jprod!(model::Model, _, v, Jv) + LinearAlgebra.mul!(Jv, model.C_lin, v[ScalarIndex]) for i in matrix_indices(model) - h += jprod(model, i, W[i.value]) + add_jprod!(model, i, v[i], Jv) end - return h + return Jv end diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 335a580..6cea41a 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -3,4 +3,9 @@ using LowRankOpt include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; model = maxcut(weights, LowRankOpt.Optimizer) + +import Percival +set_attribute(model, "solver", LRO.BurerMonteiro.Solver) +set_attribute(model, "sub_solver", Percival.PercivalSolver) +set_attribute(model, "ranks", [1]) optimize!(model) diff --git a/test/Project.toml b/test/Project.toml index f6e56aa..311b17b 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,4 +2,5 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From 6d80835c95f9f5c4f181f1d776f5d8b688005f8c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Jun 2025 17:43:47 +0200 Subject: [PATCH 05/48] Define jprod --- src/BurerMonteiro.jl | 30 ++++- src/LowRankOpt.jl | 1 + src/factorization.jl | 256 +++++++++++++++++++++++++++++++++++++++++++ src/sets.jl | 194 -------------------------------- 4 files changed, 284 insertions(+), 197 deletions(-) create mode 100644 src/factorization.jl diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 40f229b..b53d732 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -1,6 +1,7 @@ module BurerMonteiro import LinearAlgebra +import FillArrays import SolverCore import NLPModels import LowRankOpt as LRO @@ -54,13 +55,22 @@ struct Solution{T,VT<:AbstractVector{T}} dim::Dimensions end -Base.eltype(::Type{<:Solution{T}}) where {T} = T -Base.eltype(x::Solution) = eltype(typeof(x)) +struct _OuterProduct{T,UT<:AbstractVector{T},VT<:AbstractVector{T}} + x::Solution{T,VT} + v::Solution{T,UT} +end + +Base.eltype(::Type{<:Union{Solution{T},_OuterProduct{T}}}) where {T} = T +Base.eltype(x::Union{Solution,_OuterProduct}) = eltype(typeof(x)) function Base.getindex(s::Solution, ::Type{LRO.ScalarIndex}) return view(s.x, Base.OneTo(s.dim.num_scalars)) end +function Base.getindex(s::_OuterProduct, ::Type{LRO.ScalarIndex}) + return getindex(s.v, LRO.ScalarIndex) +end + function Base.getindex(s::Solution, mi::LRO.MatrixIndex) i = mi.value U = reshape( @@ -71,6 +81,12 @@ function Base.getindex(s::Solution, mi::LRO.MatrixIndex) return LRO.positive_semidefinite_factorization(U) end +function Base.getindex(s::_OuterProduct{T}, i::LRO.MatrixIndex) where {T} + U = s.x[i] + V = s.v[i] + return LRO.AsymmetricFactorization(U, V, FillArrays.Fill(T(2), size(U, 2))) +end + function NLPModels.obj(model::Model, x::AbstractVector) return NLPModels.obj(model.model, Solution(x, model.dim)) end @@ -91,6 +107,14 @@ function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) NLPModels.cons!(model.model, Solution(x, model.dim), cx) end +function NLPModels.jprod!(model::Model, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) + X = Solution(x, model.dim) + V = Solution(v, model.dim) + # The second argument is ignored as it is linear so it does + # not matter that we give `x` + NLPModels.jprod!(model.model, X, _OuterProduct(X, V), Jv) +end + function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, Jtv::AbstractVector) X = Solution(x, model.dim) JtV = Solution(Jtv, model.dim) @@ -104,7 +128,7 @@ function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, J for j in eachindex(y) A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) LinearAlgebra.mul!(U, A, X[i].factor) - U ./= 2 + U .*= 2 end end return Jtv diff --git a/src/LowRankOpt.jl b/src/LowRankOpt.jl index 3f1e257..ea44268 100644 --- a/src/LowRankOpt.jl +++ b/src/LowRankOpt.jl @@ -10,6 +10,7 @@ import KrylovKit import MathOptInterface as MOI include("sets.jl") +include("factorization.jl") include("attributes.jl") include("distance_to_set.jl") include("Test/Test.jl") diff --git a/src/factorization.jl b/src/factorization.jl new file mode 100644 index 0000000..74cdb1c --- /dev/null +++ b/src/factorization.jl @@ -0,0 +1,256 @@ +# Copyright (c) 2024: Benoît Legat and contributors +# +# 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. + +abstract type AbstractFactorization{T,F} <: AbstractMatrix{T} end + +function Base.size(m::AbstractFactorization) + n = size(left_factor(m), 1) + return (n, n) +end + +function Base.getindex(m::AbstractFactorization, i::Int, j::Int) + left = left_factor(m) + right = left_factor(m) + return sum( + left[i, k] * m.scaling[k] * right[j, k]' for + k in eachindex(m.scaling) + ) +end + +""" + struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, + } <: AbstractMatrix{T} + factor::F + scaling::D + end + +Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. +If `factor` is a vector and `diagonal` is a scalar, this corresponds to +the matrix `diagonal * factor * factor'`. +If `factor` is a matrix and `diagonal` is a vector, this corresponds to +the matrix `factor * Diagonal(scaling) * factor'`. +""" +struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{AbstractArray{T,0},AbstractVector{T}}, +} <: AbstractFactorization{T,F} + factor::F + scaling::D + function Factorization{T,F,S}( + factor::AbstractMatrix{T}, + scaling::AbstractVector{T}, + ) where {T,F<:AbstractMatrix{T},S<:AbstractVector{T}} + if length(scaling) != size(factor, 2) + error( + "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", + ) + end + return new{T,F,S}(factor, scaling) + end + function Factorization{T,F,S}( + factor::AbstractVector{T}, + scaling::AbstractArray{T,0}, + ) where {T,F<:AbstractVector{T},S<:AbstractArray{T,0}} + return new{T,F,S}(factor, scaling) + end +end + +function Factorization( + factor::AbstractMatrix{T}, + scaling::AbstractVector{T}, +) where {T} + return Factorization{T,typeof(factor),typeof(scaling)}(factor, scaling) +end + +function Factorization( + factor::AbstractVector{T}, + scaling::AbstractArray{T,0}, +) where {T} + return Factorization{T,typeof(factor),typeof(scaling)}(factor, scaling) +end + +function Factorization(factor::AbstractVector{T}, scaling::T) where {T} + return Factorization(factor, fill(scaling, tuple())) +end + +left_factor(m::Factorization) = m.factor +right_factor(m::Factorization) = m.factor + +function Base.promote_rule( + ::Type{Factorization{T,M,S1}}, + ::Type{Factorization{T,V,S2}}, +) where {T,M<:AbstractMatrix{T},V<:AbstractVector{T},S1,S2} + return Factorization{T,M,S1} +end + +function Base.convert( + ::Type{Factorization{T,F,V}}, + f::Factorization{S,<:AbstractVector{S},<:AbstractArray{S,0}}, +) where {T,S,F<:AbstractMatrix{T},V<:AbstractVector{T}} + return Factorization{T,F,V}( + reshape(f.factor, length(f.factor), 1), + reshape(f.scaling, 1), + ) +end + +function Base.convert( + ::Type{<:Factorization{T,F,V}}, + f::Factorization{S,<:AbstractMatrix{S},<:AbstractVector{S}}, +) where {T,S,F<:AbstractVector{T},V<:AbstractArray{T,0}} + return Factorization{T,F,V}( + reshape(f.factor, size(f.factor, 1)), + reshape(f.scaling, tuple()), + ) +end + +function MOI.Bridges.Constraint.conversion_cost( + ::Type{<:Factorization{T,<:AbstractMatrix{T},<:AbstractVector{T}}}, + ::Type{<:Factorization{T,<:AbstractVector{T},<:AbstractArray{T,0}}}, +) where {T} + return 1.0 +end + +function MOI.Bridges.Constraint.conversion_cost( + ::Type{<:AbstractMatrix}, + ::Type{<:AbstractMatrix}, +) + return Inf +end + +# Solvers are recommented to use this constant instead of hardcoding this +# `FillArrays` type so that the solver does not have to explicitly `import` +# `FillArrays` nor explicitly add it to its dependency so that it remains +# a detail that's internal to LowRankOpt that we can easily change later +# The rest of the code of LowRankOpt should also use these two constants +# and not `FillArrays` directly. + +import FillArrays + +const One{T} = FillArrays.Ones{T,0,Tuple{}} +const Ones{T} = FillArrays.Ones{T,1,Tuple{Base.OneTo{Int}}} + +function positive_semidefinite_factorization( + factor::AbstractVector{T}, +) where {T} + return Factorization(factor, One{T}(tuple())) +end + +function positive_semidefinite_factorization( + factor::AbstractMatrix{T}, +) where {T} + return Factorization(factor, Ones{T}(Base.OneTo(size(factor, 2)))) +end + +struct AsymmetricFactorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{AbstractArray{T,0},AbstractVector{T}}, +} <: AbstractFactorization{T,F} + left::F + right::F + scaling::D + function AsymmetricFactorization{T,F,S}( + left::AbstractMatrix{T}, + right::AbstractMatrix{T}, + scaling::AbstractVector{T}, + ) where {T,F<:AbstractMatrix{T},S<:AbstractVector{T}} + if size(left) != size(right) + error( + "Size `$(size(left))` of left factor does not match size `$(size(right))` of right factor", + ) + end + if length(scaling) != size(left, 2) + error( + "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", + ) + end + return new{T,F,S}(left, right, scaling) + end + function AsymmetricFactorization{T,F,S}( + left::AbstractVector{T}, + right::AbstractVector{T}, + scaling::AbstractArray{T,0}, + ) where {T,F<:AbstractVector{T},S<:AbstractArray{T,0}} + return new{T,F,S}(left, right, scaling) + end +end + +function AsymmetricFactorization( + left::AbstractMatrix{T}, + right::AbstractMatrix{T}, + scaling::AbstractVector{T}, +) where {T} + return AsymmetricFactorization{T,typeof(left),typeof(scaling)}(left, right, scaling) +end + +function AsymmetricFactorization( + left::AbstractVector{T}, + right::AbstractVector{T}, + scaling::AbstractArray{T,0}, +) where {T} + return AsymmetricFactorization{T,typeof(left),typeof(scaling)}(left, right, scaling) +end + +left_factor(m::AsymmetricFactorization) = m.left +right_factor(m::AsymmetricFactorization) = m.right + +""" + symmetrize_factorization(L, R) + +Factorization corresponding to the symmetrization `(L * R' + R * L') / 2` of `L * R'`. + +## Example + +```jldoctest +julia> LowRankOpt.symmetrize_factorization([1, 0], [0, 1]) +2×2 LowRankOpt.Factorization{Float64, Matrix{Float64}, Vector{Float64}}: + 0.0 0.5 + 0.5 0.0 +``` +""" +function symmetrize_factorization(L, R; use_krylov = true) + sym = LinearAlgebra.Symmetric((L * R' + R * L') / 2) + r = 2size(L, 2) + if use_krylov + eigvals, factors = KrylovKit.eigsolve(sym, r) + factor = reduce(hcat, factors) + else + eigvals, factor = LinearAlgebra.eigen(sym) + end + σ = sortperm(abs.(eigvals), rev = true) + keep = σ[1:r] + return Factorization(factor[:, keep], eigvals[keep]) +end + +struct TriangleVectorization{T,M<:AbstractMatrix{T}} <: AbstractVector{T} + matrix::M +end + +function Base.convert( + ::Type{TriangleVectorization{T,M}}, + t::TriangleVectorization, +) where {T,M} + return TriangleVectorization{T,M}(t.matrix) +end + +function MOI.Bridges.Constraint.conversion_cost( + ::Type{TriangleVectorization{T,M1}}, + ::Type{TriangleVectorization{T,M2}}, +) where {T,M1,M2} + return MOI.Bridges.Constraint.conversion_cost(M1, M2) +end + +function Base.size(v::TriangleVectorization) + n = size(v.matrix, 1) + return (MOI.Utilities.trimap(n, n),) +end + +function Base.getindex(v::TriangleVectorization, k::Int) + return getindex(v.matrix, MOI.Utilities.inverse_trimap(k)...) +end diff --git a/src/sets.jl b/src/sets.jl index af41b58..1fcc1d7 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -139,197 +139,3 @@ function MOI.Utilities.dot_coefficients( c[(n+1):end] = MOI.Utilities.dot_coefficients(x[(n+1):end], set.set) return c end - -abstract type AbstractFactorization{T,F} <: AbstractMatrix{T} end - -function Base.size(m::AbstractFactorization) - n = size(m.factor, 1) - return (n, n) -end - -""" - struct Factorization{ - T, - F<:Union{AbstractVector{T},AbstractMatrix{T}}, - D<:Union{T,AbstractVector{T}}, - } <: AbstractMatrix{T} - factor::F - scaling::D - end - -Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. -If `factor` is a vector and `diagonal` is a scalar, this corresponds to -the matrix `diagonal * factor * factor'`. -If `factor` is a matrix and `diagonal` is a vector, this corresponds to -the matrix `factor * Diagonal(scaling) * factor'`. -""" -struct Factorization{ - T, - F<:Union{AbstractVector{T},AbstractMatrix{T}}, - D<:Union{AbstractArray{T,0},AbstractVector{T}}, -} <: AbstractFactorization{T,F} - factor::F - scaling::D - function Factorization{T,F,S}( - factor::AbstractMatrix{T}, - scaling::AbstractVector{T}, - ) where {T,F<:AbstractMatrix{T},S<:AbstractVector{T}} - if length(scaling) != size(factor, 2) - error( - "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", - ) - end - return new{T,F,S}(factor, scaling) - end - function Factorization{T,F,S}( - factor::AbstractVector{T}, - scaling::AbstractArray{T,0}, - ) where {T,F<:AbstractVector{T},S<:AbstractArray{T,0}} - return new{T,F,S}(factor, scaling) - end -end - -function Factorization( - factor::AbstractMatrix{T}, - scaling::AbstractVector{T}, -) where {T} - return Factorization{T,typeof(factor),typeof(scaling)}(factor, scaling) -end - -function Factorization( - factor::AbstractVector{T}, - scaling::AbstractArray{T,0}, -) where {T} - return Factorization{T,typeof(factor),typeof(scaling)}(factor, scaling) -end - -function Factorization(factor::AbstractVector{T}, scaling::T) where {T} - return Factorization(factor, fill(scaling, tuple())) -end - -function Base.getindex(m::Factorization, i::Int, j::Int) - return sum( - m.factor[i, k] * m.scaling[k] * m.factor[j, k]' for - k in eachindex(m.scaling) - ) -end - -function Base.promote_rule( - ::Type{Factorization{T,M,S1}}, - ::Type{Factorization{T,V,S2}}, -) where {T,M<:AbstractMatrix{T},V<:AbstractVector{T},S1,S2} - return Factorization{T,M,S1} -end - -function Base.convert( - ::Type{Factorization{T,F,V}}, - f::Factorization{S,<:AbstractVector{S},<:AbstractArray{S,0}}, -) where {T,S,F<:AbstractMatrix{T},V<:AbstractVector{T}} - return Factorization{T,F,V}( - reshape(f.factor, length(f.factor), 1), - reshape(f.scaling, 1), - ) -end - -function Base.convert( - ::Type{<:Factorization{T,F,V}}, - f::Factorization{S,<:AbstractMatrix{S},<:AbstractVector{S}}, -) where {T,S,F<:AbstractVector{T},V<:AbstractArray{T,0}} - return Factorization{T,F,V}( - reshape(f.factor, size(f.factor, 1)), - reshape(f.scaling, tuple()), - ) -end - -function MOI.Bridges.Constraint.conversion_cost( - ::Type{<:Factorization{T,<:AbstractMatrix{T},<:AbstractVector{T}}}, - ::Type{<:Factorization{T,<:AbstractVector{T},<:AbstractArray{T,0}}}, -) where {T} - return 1.0 -end - -function MOI.Bridges.Constraint.conversion_cost( - ::Type{<:AbstractMatrix}, - ::Type{<:AbstractMatrix}, -) - return Inf -end - -# Solvers are recommented to use this constant instead of hardcoding this -# `FillArrays` type so that the solver does not have to explicitly `import` -# `FillArrays` nor explicitly add it to its dependency so that it remains -# a detail that's internal to LowRankOpt that we can easily change later -# The rest of the code of LowRankOpt should also use these two constants -# and not `FillArrays` directly. - -import FillArrays - -const One{T} = FillArrays.Ones{T,0,Tuple{}} -const Ones{T} = FillArrays.Ones{T,1,Tuple{Base.OneTo{Int}}} - -function positive_semidefinite_factorization( - factor::AbstractVector{T}, -) where {T} - return Factorization(factor, One{T}(tuple())) -end - -function positive_semidefinite_factorization( - factor::AbstractMatrix{T}, -) where {T} - return Factorization(factor, Ones{T}(Base.OneTo(size(factor, 2)))) -end - -""" - symmetrize_factorization(L, R) - -Factorization corresponding to the symmetrization `(L * R' + R * L') / 2` of `L * R'`. - -## Example - -```jldoctest -julia> LowRankOpt.symmetrize_factorization([1, 0], [0, 1]) -2×2 LowRankOpt.Factorization{Float64, Matrix{Float64}, Vector{Float64}}: - 0.0 0.5 - 0.5 0.0 -``` -""" -function symmetrize_factorization(L, R; use_krylov = true) - sym = LinearAlgebra.Symmetric((L * R' + R * L') / 2) - r = 2size(L, 2) - if use_krylov - eigvals, factors = KrylovKit.eigsolve(sym, r) - factor = reduce(hcat, factors) - else - eigvals, factor = LinearAlgebra.eigen(sym) - end - σ = sortperm(abs.(eigvals), rev = true) - keep = σ[1:r] - return Factorization(factor[:, keep], eigvals[keep]) -end - -struct TriangleVectorization{T,M<:AbstractMatrix{T}} <: AbstractVector{T} - matrix::M -end - -function Base.convert( - ::Type{TriangleVectorization{T,M}}, - t::TriangleVectorization, -) where {T,M} - return TriangleVectorization{T,M}(t.matrix) -end - -function MOI.Bridges.Constraint.conversion_cost( - ::Type{TriangleVectorization{T,M1}}, - ::Type{TriangleVectorization{T,M2}}, -) where {T,M1,M2} - return MOI.Bridges.Constraint.conversion_cost(M1, M2) -end - -function Base.size(v::TriangleVectorization) - n = size(v.matrix, 1) - return (MOI.Utilities.trimap(n, n),) -end - -function Base.getindex(v::TriangleVectorization, k::Int) - return getindex(v.matrix, MOI.Utilities.inverse_trimap(k)...) -end From 0957df908d5369a4eb29f16d09d9ec74d78c7ea3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Jun 2025 21:16:28 +0200 Subject: [PATCH 06/48] Add tests --- src/BurerMonteiro.jl | 24 ++++++++++++++++++++---- src/MOI_wrapper.jl | 24 +++++++++--------------- src/model.jl | 2 +- test/BurerMonteiro.jl | 24 +++++++++++++++++++++++- test/Project.toml | 3 +++ 5 files changed, 56 insertions(+), 21 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index b53d732..04f80e4 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -98,7 +98,7 @@ function NLPModels.grad!(model::Model, x::AbstractVector, g::AbstractVector) for i in LRO.matrix_indices(model.model) C = NLPModels.grad(model.model, i) LinearAlgebra.mul!(G[i].factor, C, X[i].factor) - G[i].factor ./= 2 + G[i].factor .*= -2 end return g end @@ -128,12 +128,29 @@ function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, J for j in eachindex(y) A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) LinearAlgebra.mul!(U, A, X[i].factor) - U .*= 2 + U .*= (2y[j]) end end return Jtv end +function NLPModels.hprod!(model::Model{T}, ::AbstractVector, y, v::AbstractVector, Hv::AbstractVector; obj_weight = one(T)) where {T} + V = Solution(v, model.dim) + HV = Solution(Hv, model.dim) + fill!(Hv, zero(eltype(Hv))) + for i in LRO.matrix_indices(model.model) + Vi = V[i].factor + C = NLPModels.grad(model.model, i) + Hvi = HV[i].factor + Hvi .+= C * Vi + Hvi .*= 2obj_weight + for j in 1:model.meta.ncon + A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) + Hvi .+= A * Vi .* (2y[j]) + end + end +end + struct Solver{T,ST} <: SolverCore.AbstractOptimizationSolver model::Model{T} solver::ST @@ -149,8 +166,7 @@ end function SolverCore.solve!( solver::Solver, - model::NLPModels.AbstractNLPModel, # Same as `solver.model.model` - stats::SolverCore.GenericExecutionStats; + model::NLPModels.AbstractNLPModel; # Same as `solver.model.model` kws..., ) SolverCore.solve!(solver.solver, solver.model, solver.stats; kws...) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 63a860a..83d9029 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -42,10 +42,6 @@ const OptimizerCache{T} = MOI.Utilities.GenericModel{ mutable struct Optimizer{T} <: MOI.AbstractOptimizer solver::Union{Nothing,SolverCore.AbstractOptimizationSolver} model::Union{Nothing,Model{T}} - stats::Union{ - Nothing, - SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any}, - } lmi_id::Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64} lin_cones::Union{Nothing,NNGCones{T}} max_sense::Bool @@ -55,7 +51,6 @@ mutable struct Optimizer{T} <: MOI.AbstractOptimizer function Optimizer{T}() where {T} return new{T}( - nothing, nothing, nothing, Dict{MOI.ConstraintIndex{VAF{T},PSD},Int64}(), @@ -82,7 +77,7 @@ function MOI.empty!(optimizer::Optimizer) return end -MOI.get(::Optimizer, ::MOI.SolverName) = "Loraine" +MOI.get(::Optimizer, ::MOI.SolverName) = "LowRankOpt" # MOI.RawOptimizerAttribute @@ -150,8 +145,7 @@ function MOI.optimize!(model::Optimizer) else options[:verbose] = 1 end - model.stats = SolverCore.GenericExecutionStats(model.model) - SolverCore.solve!(model.solver, model.model, model.stats; options...) + SolverCore.solve!(model.solver, model.model; options...) return end @@ -253,11 +247,11 @@ function MOI.copy_to(dest::Optimizer{T}, src::MOI.ModelLike) where {T} end function MOI.get(optimizer::Optimizer, ::MOI.SolveTimeSec) - return optimizer.stats.elapsed_time + return optimizer.solver.stats.elapsed_time end function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString) - return SolverCore.STATUSES[optimizer.stats.status] + return SolverCore.STATUSES[optimizer.solver.stats.status] end struct RawStatus <: MOI.AbstractModelAttribute @@ -267,14 +261,14 @@ end MOI.is_set_by_optimize(::RawStatus) = true function MOI.get(optimizer::Optimizer, attr::RawStatus) - return getfield(optimizer.stats, attr.name) + return getfield(optimizer.solver.stats, attr.name) end function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) - if isnothing(optimizer.stats) + if isnothing(optimizer.solver.stats) return MOI.OPTIMIZE_NOT_CALLED end - return NLPModelsJuMP.TERMINATION_STATUS[optimizer.stats.status] + return NLPModelsJuMP.TERMINATION_STATUS[optimizer.solver.stats.status] end function MOI.get(model::Optimizer, ::MOI.ResultCount) @@ -298,13 +292,13 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) - val = optimizer.stats.objective + val = optimizer.solver.stats.objective return optimizer.max_sense ? -val : val end function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) MOI.check_result_index_bounds(optimizer, attr) - return optimizer.stats.solution[vi.value] + return optimizer.solver.stats.solution[vi.value] end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} diff --git a/src/model.jl b/src/model.jl index 80de35d..b1d75e9 100644 --- a/src/model.jl +++ b/src/model.jl @@ -203,7 +203,7 @@ function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y, S) return jtprod!(buffer[i], model, mat_idx, y) + model.C[i] - S[i] end -NLPModels.grad(model::Model, ::Type{ScalarIndex}) = model.d_lin +NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin NLPModels.grad(model::Model, i::MatrixIndex) = model.C[i.value] cons_constant(model::Model) = model.b diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 6cea41a..fa135b9 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -1,11 +1,33 @@ using Test using LowRankOpt +import Percival + include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; model = maxcut(weights, LowRankOpt.Optimizer) -import Percival set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) set_attribute(model, "ranks", [1]) optimize!(model) +solution_summary(model) + +import NLPModels, FiniteDiff +function jac_check(model, x; tol = 1e-6) + f(x) = NLPModels.cons(model, x) + J = FiniteDiff.finite_difference_jacobian(f, x) + v = rand(model.meta.nvar) + @test NLPModels.jprod(model, x, v) ≈ J * v rtol = tol atol = tol + v = rand(model.meta.ncon) + @test NLPModels.jtprod(model, x, v) ≈ J' * v rtol = tol atol = tol +end + + +@testset "Diff check" begin + b = unsafe_backend(model) + using NLPModelsTest + bm = b.solver.model + x = rand(bm.meta.nvar) + @test isempty(NLPModelsTest.gradient_check(bm; x)) + jac_check(bm, x) +end diff --git a/test/Project.toml b/test/Project.toml index 311b17b..3ea19cf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,9 @@ [deps] +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From a0eea2b5fdc93fc5226084ea50e5532d8b4b1445 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 08:00:39 +0200 Subject: [PATCH 07/48] Fixes --- src/model.jl | 28 ++++++++++++++++++++++++++-- src/schur.jl | 2 +- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/src/model.jl b/src/model.jl index b1d75e9..86b26e6 100644 --- a/src/model.jl +++ b/src/model.jl @@ -6,6 +6,26 @@ import MutableArithmetics as MA import MathOptInterface as MOI import NLPModels +struct Solution{T} <: AbstractVector{T} + scalars::Vector{T} + matrices::Vector{LinearAlgebra.Symmetric{T,Matrix{T}}} +end + +function Base.zero(::Type{Solution{T}}, num_scalars::Integer, side_dimensions) where {T} + return Solution{T}( + zeros(T, num_scalars), + [LinearAlgebra.Symmetric(zeros(T, d, d)) for d in side_dimensions] + ) +end + +struct Meta{T} <: NLPModels.AbstractNLPModelMeta{T,Solution{T}} + nvar::Int + x0::Solution{T} + ncon::Int + y0::Vector{T} + minimize::Bool +end + """ Model @@ -28,7 +48,7 @@ The fields of the `struct` as related to the arrays of the above formulation as * The matrix ``A_{i,j}`` is given by `-A[i,j]`. """ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vector{T}} - meta::NLPModels.NLPModelMeta{T,Vector{T}} + meta::Meta{T} C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} A::Matrix{A} b::Vector{T} @@ -54,12 +74,16 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vec model.d_lin = d_lin model.C_lin = C_lin model.msizes = msizes - model.meta = NLPModels.NLPModelMeta( + model.meta = Meta{T}( num_scalars(model) + sum( Base.Fix1(side_dimension, model), matrix_indices(model); init = 0 ), + zero(Solution{T}, num_scalars(model), msizes), + length(b), + zero(b), + true, ) return model end diff --git a/src/schur.jl b/src/schur.jl index e9ba5c8..ed92f8b 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -35,7 +35,7 @@ function buffer_for_schur_complement(model::Model, κ) for mat_idx in matrix_indices(model) i = mat_idx.value - nzA = [nnz(model.A[i, j]) for j in 1:n] + nzA = [SparseArrays.nnz(model.A[i, j]) for j in 1:n] σ[:,i] = sortperm(nzA, rev = true) sorted = nzA[σ[:,i]] From ce58ccf3f8de1474f05d65137145975ed39cd2b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 15:33:17 +0200 Subject: [PATCH 08/48] Fixes for Loraine --- src/MOI_wrapper.jl | 59 +++++++++++++++++--------- src/model.jl | 103 +++++++++++++++++++++++++++++---------------- src/schur.jl | 63 +++++++++++++-------------- 3 files changed, 138 insertions(+), 87 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 83d9029..871942e 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -72,6 +72,7 @@ end MOI.is_empty(optimizer::Optimizer) = isnothing(optimizer.model) function MOI.empty!(optimizer::Optimizer) + optimizer.solver = nothing optimizer.model = nothing optimizer.lin_cones = nothing return @@ -90,9 +91,6 @@ function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value) throw(MOI.UnsupportedAttribute(param)) end optimizer.options[param.name] = value - if !isnothing(optimizer.solver) - setproperty!(optimizer.solver, Symbol(param.name), value) - end return end @@ -264,11 +262,28 @@ function MOI.get(optimizer::Optimizer, attr::RawStatus) return getfield(optimizer.solver.stats, attr.name) end +# From the point of view of the solver, only a local solution is found. +# However, this this is a convex problem, this is actually a global minimum! +# We define this function instead of hard-coding `MOI.OPTIMAL` so that +# `BurerMonteiro` can override it since it is solving a non-convex formulation. +function termination_status(solver::SolverCore.AbstractOptimizationSolver) + if isnothing(solver.stats) + return MOI.OPTIMIZE_NOT_CALLED + end + status = NLPModelsJuMP.TERMINATION_STATUS[solver.stats.status] + if status == MOI.LOCALLY_SOLVED + status = MOI.OPTIMAL + elseif status == MOI.LOCALLY_INFEASIBLE + status = MOI.INFEASIBLE + end + return status +end + function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) - if isnothing(optimizer.solver.stats) - return MOI.OPTIMIZE_NOT_CALLED - end - return NLPModelsJuMP.TERMINATION_STATUS[optimizer.solver.stats.status] + if isnothing(optimizer.solver) + return MOI.OPTIMIZE_NOT_CALLED + end + return termination_status(optimizer.solver) end function MOI.get(model::Optimizer, ::MOI.ResultCount) @@ -282,7 +297,7 @@ end function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) return MOI.NO_SOLUTION - elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.LOCALLY_SOLVED + elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT else # TODO @@ -292,26 +307,34 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) - val = optimizer.solver.stats.objective + val = dual_obj(optimizer.model, optimizer.solver.stats.multipliers) return optimizer.max_sense ? -val : val end function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) - MOI.check_result_index_bounds(optimizer, attr) - return optimizer.solver.stats.solution[vi.value] + MOI.check_result_index_bounds(optimizer, attr) + return optimizer.solver.stats.multipliers[vi.value] end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) - val = Solvers.obj(optimizer.solver.model, optimizer.solver.X_lin, optimizer.solver.X)::T + val = optimizer.solver.stats.objective return optimizer.max_sense ? -val : val end -function MOI.get(::Optimizer, ::MOI.DualStatus) - # TODO - return MOI.NO_SOLUTION +function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) + if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) + return MOI.NO_SOLUTION + elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] + return MOI.FEASIBLE_POINT + else + # TODO + return MOI.UNKNOWN_RESULT_STATUS + end end +_solution(optimizer::Optimizer) = VectorizedSolution(optimizer.solver.stats.solution, optimizer.model.dim) + function MOI.get( optimizer::Optimizer{T}, attr::MOI.ConstraintDual, @@ -319,9 +342,7 @@ function MOI.get( ) where {T} MOI.check_result_index_bounds(optimizer, attr) lmi_id = optimizer.lmi_id[ci] - X = optimizer.solver.X[lmi_id] - n = optimizer.solver.model.msizes[lmi_id] - return [X[i, j] for j = 1:n for i = 1:j]::Vector{T} + return TriangleVectorization(_solution(optimizer)[MatrixIndex(lmi_id)]) end function MOI.get( @@ -331,5 +352,5 @@ function MOI.get( ) where {T} MOI.check_result_index_bounds(optimizer, attr) rows = MOI.Utilities.rows(optimizer.lin_cones, ci) - return optimizer.solver.X_lin[rows]::Vector{T} + return _solution(optimizer)[ScalarIndex][rows] end diff --git a/src/model.jl b/src/model.jl index 86b26e6..6deace8 100644 --- a/src/model.jl +++ b/src/model.jl @@ -6,26 +6,66 @@ import MutableArithmetics as MA import MathOptInterface as MOI import NLPModels -struct Solution{T} <: AbstractVector{T} - scalars::Vector{T} - matrices::Vector{LinearAlgebra.Symmetric{T,Matrix{T}}} +struct Dimensions + num_scalars::Int64 + side_dimensions::Vector{Int64} + offsets::Vector{Int64} end -function Base.zero(::Type{Solution{T}}, num_scalars::Integer, side_dimensions) where {T} - return Solution{T}( - zeros(T, num_scalars), - [LinearAlgebra.Symmetric(zeros(T, d, d)) for d in side_dimensions] - ) +Base.length(d::Dimensions) = d.offsets[end] + +struct ScalarIndex + value::Int64 +end + +struct MatrixIndex + value::Int64 +end + +abstract type AbstractSolution{T} <: AbstractVector{T} end + +Base.getindex(s::AbstractSolution, i::Union{Type{ScalarIndex},MatrixIndex}) = view(s, i) + +struct VectorizedSolution{T} <: AbstractSolution{T} + x::Vector{T} + dim::Dimensions +end + +Base.similar(s::VectorizedSolution) = VectorizedSolution(similar(s.x), s.dim) + +Base.size(s::VectorizedSolution) = (length(s.dim),) + +Base.to_index(s::VectorizedSolution, ::Type{ScalarIndex}) = Base.OneTo(s.dim.num_scalars) + +function Base.to_index(s::VectorizedSolution, mi::MatrixIndex) + i = mi.value + return (1 + s.dim.offsets[i]):s.dim.offsets[i+1] +end + +function Base.view(s::VectorizedSolution, ::Type{ScalarIndex}) + return view(s.x, Base.to_index(s, ScalarIndex)) end -struct Meta{T} <: NLPModels.AbstractNLPModelMeta{T,Solution{T}} - nvar::Int - x0::Solution{T} - ncon::Int - y0::Vector{T} - minimize::Bool +# `s[i] .= ...` calleds `copyto!(view(s, i), Broadcasted(...))` +function Base.view(s::VectorizedSolution, i::MatrixIndex) + v = view(s.x, Base.to_index(s, i)) + dim = s.dim.side_dimensions[i.value] + X = reshape(v, dim, dim) + return X end +Base.setindex!(s::VectorizedSolution, v, i::Integer) = setindex!(s.x, v, i) + +struct ShapedSolution{T,MT<:AbstractMatrix{T}} <: AbstractSolution{T} + scalars::Vector{T} + matrices::Vector{MT} +end + +Base.size(s::ShapedSolution) = (length(s.scalars) + sum(length, s.matrices, init = 0),) + +Base.view(s::ShapedSolution, ::Type{ScalarIndex}) = s.scalars +Base.view(s::ShapedSolution, i::MatrixIndex) = s.matrices[i.value] + """ Model @@ -48,7 +88,8 @@ The fields of the `struct` as related to the arrays of the above formulation as * The matrix ``A_{i,j}`` is given by `-A[i,j]`. """ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vector{T}} - meta::Meta{T} + meta::NLPModels.NLPModelMeta{T,Vector{T}} + dim::Dimensions C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} A::Matrix{A} b::Vector{T} @@ -74,17 +115,13 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vec model.d_lin = d_lin model.C_lin = C_lin model.msizes = msizes - model.meta = Meta{T}( - num_scalars(model) + sum( - Base.Fix1(side_dimension, model), - matrix_indices(model); - init = 0 - ), - zero(Solution{T}, num_scalars(model), msizes), - length(b), - zero(b), - true, + n = num_scalars(model) + model.meta = NLPModels.NLPModelMeta{T,Vector{T}}( + n + sum(abs2, msizes, init = 0), + ncon = length(b), ) + offsets = n .+ [0; cumsum(abs2.(msizes))] + model.dim = Dimensions(n, msizes, offsets) return model end end @@ -97,20 +134,12 @@ end # but these are not box constraints NLPModels.has_bounds(::Model) = false -struct ScalarIndex - value::Int64 -end - num_scalars(model::Model) = length(model.d_lin) function scalar_indices(model::Model) return MOI.Utilities.LazyMap{ScalarIndex}(ScalarIndex, Base.OneTo(num_scalars(model))) end -struct MatrixIndex - value::Int64 -end - num_matrices(model::Model) = length(model.C) function matrix_indices(model::Model) @@ -135,7 +164,7 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) return zero(T) end - return norm(model.A[i.value, :]) + return LinearAlgebra.norm(model.A[i.value, :]) end function NLPModels.obj(model::Model, X, i::MatrixIndex) @@ -166,7 +195,7 @@ function NLPModels.grad!(model::Model, _, g) return g end -dual_obj(model::Model, y) = -dot(model.b, y) + model.b_const +dual_obj(model::Model, y) = -LinearAlgebra.dot(model.b, y) + model.b_const function jtprod(model::Model, ::Type{ScalarIndex}, y) return -model.C_lin' * y @@ -232,7 +261,7 @@ NLPModels.grad(model::Model, i::MatrixIndex) = model.C[i.value] cons_constant(model::Model) = model.b -function NLPModels.cons!(model::Model, x, cx) +function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) NLPModels.jprod!(model, x, x, cx) cx .*= -1 cx .+= model.b @@ -245,7 +274,7 @@ function add_jprod!(model::Model, i::MatrixIndex, V, Jv) end end -function NLPModels.jprod!(model::Model, _, v, Jv) +function NLPModels.jprod!(model::Model, _::AbstractVector, v::AbstractVector, Jv::AbstractVector) LinearAlgebra.mul!(Jv, model.C_lin, v[ScalarIndex]) for i in matrix_indices(model) add_jprod!(model, i, v[i], Jv) diff --git a/src/schur.jl b/src/schur.jl index ed92f8b..ec25364 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -7,18 +7,18 @@ function _dot(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, # have the same sizes so we can safely use `@inbounds` result = zero(eltype(A)) @inbounds for i in axes(A, 2) - nzA = nzrange(A, i) + nzA = SparseArrays.nzrange(A, i) if !isempty(nzA) for j in axes(B, 2) - nzB = nzrange(B, j) + nzB = SparseArrays.nzrange(B, j) if !isempty(nzB) AW = zero(result) for k in nzA - AW += nonzeros(A)[k] * W[rowvals(A)[k], j] + AW += SparseArrays.nonzeros(A)[k] * W[SparseArrays.rowvals(A)[k], j] end WB = zero(result) for k in nzB - WB += W[i, rowvals(B)[k]] * nonzeros(B)[k] + WB += W[i, SparseArrays.rowvals(B)[k]] * SparseArrays.nonzeros(B)[k] end result += AW * WB end @@ -64,11 +64,11 @@ end function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) n = num_constraints(model) - BBBB = zeros(eltype(eltype(W)), n, n) - for mat_idx in matrix_indices(model) - BBBB += schur_complement(buffer, model, mat_idx, W[mat_idx.value]) + H = zeros(eltype(eltype(W)), n, n) + for i in matrix_indices(model) + H += schur_complement(buffer, model, i, W[i]) end - return BBBB + return H end ##### @@ -80,26 +80,28 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T dim = side_dimension(model, mat_idx) @assert dim == size(W, 1) == size(W, 2) tmp1 = Matrix{T}(undef, size(W, 2), dim) + tmp2 = Vector{T}(undef, num_constraints(model)) tmp = zeros(T, size(W, 2), dim) for ii = 1:n i = σ[ii,ilmi] Ai = model.A[ilmi, i] - if nnz(Ai) > 0 + if SparseArrays.nnz(Ai) > 0 if ii <= last_dense[ilmi] - mul!(tmp1, W, Ai) - mul!(tmp, tmp1, W) - tmp2 = jprod(model, mat_idx, tmp) + LinearAlgebra.mul!(tmp1, W, Ai) + LinearAlgebra.mul!(tmp, tmp1, W) + fill!(tmp2, zero(T)) + add_jprod!(model, mat_idx, tmp, tmp2) indi = σ[ii:end,ilmi] BBBB[indi,i] .= -tmp2[indi] BBBB[i,indi] .= -tmp2[indi] else - if !iszero(nnz(Ai)) - if nnz(Ai) > 1 + if !iszero(SparseArrays.nnz(Ai)) + if SparseArrays.nnz(Ai) > 1 @inbounds for jj = ii:n j = σ[jj,ilmi] Aj = model.A[ilmi, j] - if !iszero(nnz(Aj)) + if !iszero(SparseArrays.nnz(Aj)) ttt = _dot(Ai, Aj, W) if i >= j BBBB[i,j] = ttt @@ -110,17 +112,17 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T end else # A is symmetric - iiiiAi = jjjiAi = only(rowvals(Ai)) - vvvi = only(nonzeros(Ai)) + iiiiAi = jjjiAi = only(SparseArrays.rowvals(Ai)) + vvvi = only(SparseArrays.nonzeros(Ai)) @inbounds for jj = ii:n j = σ[jj,ilmi] Ajjj = model.A[ilmi, j] # As we sort the matrices in decreasing `nnz` order, # the rest of matrices is either zero or have only # one entry - if !iszero(nnz(Ajjj)) - iiijAj = jjjjAj = only(rowvals(Ajjj)) - vvvj = only(nonzeros(Ajjj)) + if !iszero(SparseArrays.nnz(Ajjj)) + iiijAj = jjjjAj = only(SparseArrays.rowvals(Ajjj)) + vvvj = only(SparseArrays.nonzeros(Ajjj)) ttt = vvvi * W[iiiiAi,iiijAj] * W[jjjiAi,jjjjAj] * vvvj if i >= j BBBB[i,j] = ttt @@ -140,34 +142,33 @@ end # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA_jW⟩ -function schur_complement(buffer, model::Model, w, W::AbstractVector) +function schur_complement(buffer, model::Model, W::AbstractVector) H = MA.Zero() if num_matrices(model) > 0 H = MA.add!!(H, schur_complement(buffer, model, W, MatrixIndex)) end if num_scalars(model) > 0 - H = MA.add!!(H, schur_complement(model, w, ScalarIndex)) + H = MA.add!!(H, schur_complement(model, W[ScalarIndex], ScalarIndex)) end if H isa MA.Zero n = num_constraints(model) - H = zeros(eltype(w), n, n) + H = zeros(eltype(W), n, n) end - return Hermitian(H, :L) + return LinearAlgebra.Hermitian(H, :L) end function schur_complement(model::Model, w, ::Type{ScalarIndex}) - return model.C_lin * spdiagm(w) * model.C_lin' + return model.C_lin * SparseArrays.spdiagm(w) * model.C_lin' end # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA(y)W⟩ -function eval_schur_complement!(buffer, result, model::Model, w, W, y) - result .= 0.0 - for mat_idx in matrix_indices(model) - i = mat_idx.value - result .-= jprod(model, mat_idx, W[i] * jtprod!(buffer[i], model, mat_idx, y) * W[i]) +function eval_schur_complement!(buffer, result, model::Model, W, y) + fill!(result, zero(eltype(result))) + for i in matrix_indices(model) + add_jprod!(model, i, -W[i] * jtprod!(buffer[i.value], model, i, y) * W[i], result) end - result .+= model.C_lin * (w .* (model.C_lin' * y)) + result .+= model.C_lin * (W[ScalarIndex] .* (model.C_lin' * y)) return result end From d7c9296953fe92920c128aabf7bcbde6ee483322 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 16:03:56 +0200 Subject: [PATCH 09/48] Remove objective constant --- src/MOI_wrapper.jl | 16 ++++++++-------- src/model.jl | 23 ++++++++++++++--------- 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 871942e..dd20682 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -75,6 +75,7 @@ function MOI.empty!(optimizer::Optimizer) optimizer.solver = nothing optimizer.model = nothing optimizer.lin_cones = nothing + optimizer.objective_constant = NaN return end @@ -193,7 +194,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} end for row in eachindex(back) lmi_id, i, j = back[row] - _add(lmi_id, 1, i, j, -psd_AC.constants[row]) + _add(lmi_id, 1, i, j, psd_AC.constants[row]) end for var = 1:n for k in SparseArrays.nzrange(psd_A, var) @@ -204,9 +205,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} end dest.max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE obj = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}()) - # objective_constant = MOI.constant(obj) # TODO # MK: done(?) - b_const = obj.constant - b_const = dest.max_sense ? -b_const : b_const + dest.objective_constant = MOI.constant(obj) b0 = zeros(T, n) for term in obj.terms b0[term.variable.value] += term.coefficient @@ -216,10 +215,9 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} AA = SparseArrays.SparseMatrixCSC{T,Int}[SparseArrays.sparse(IJV...) for IJV in A] dest.model = Model( - -AA[:,1], + AA[:,1], AA[:,2:end], b, - b_const, convert(SparseArrays.SparseVector{T,Int64}, SparseArrays.sparsevec(Cd_lin.constants)), C_lin, msizes, @@ -299,6 +297,8 @@ function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) return MOI.NO_SOLUTION elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT + elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE + return MOI.INFEASIBLE_POINT else # TODO return MOI.UNKNOWN_RESULT_STATUS @@ -308,7 +308,7 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = dual_obj(optimizer.model, optimizer.solver.stats.multipliers) - return optimizer.max_sense ? -val : val + return optimizer.objective_constant + (optimizer.max_sense ? -val : val) end function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) @@ -319,7 +319,7 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = optimizer.solver.stats.objective - return optimizer.max_sense ? -val : val + return optimizer.objective_constant + (optimizer.max_sense ? -val : val) end function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) diff --git a/src/model.jl b/src/model.jl index 6deace8..8376a1d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -69,15 +69,23 @@ Base.view(s::ShapedSolution, i::MatrixIndex) = s.matrices[i.value] """ Model -Model representing the problem: +Model representing the primal-dual pair of problems: ```math \\begin{aligned} -\\max {} & b^\\top y - b_\\text{const} +\\min {} & \\sum_{i=1}^\\text{nlmi} +\\langle C_i, X_i \\rangle + \\langle d_\\text{lin}, x \\rangle & +\\max {} & b^\\top y \\\\ -& \\sum_{j=1}^n y_j A_{i,j} \\preceq C_i +\\text{s.t. } & \\sum_{i=1}^\\text{nlmi} \\langle A_{i,j}, X_i \\rangle + (C_\\text{lin} x)_j = b_j +\\qquad +\\forall j \\in \\{1,\\ldots,m\\} & +\\text{s.t. } & \\sum_{j=1}^m y_j A_{i,j} \\preceq C_i \\qquad \\forall i \\in \\{1,\\ldots,\\text{nlmi}\\} \\\\ +& x \\ge 0, X_i \\succeq 0 +\\qquad +\\forall i \\in \\{1,\\ldots,\\text{nlmi}\\} & & C_\\text{lin}^\\top y \\le d_\\text{lin} \\end{aligned} ``` @@ -85,7 +93,7 @@ The fields of the `struct` as related to the arrays of the above formulation as * The ``i``th PSD constraint is of size `msize[i] × msisze[i]` * The matrix ``C_i`` is given by `C[i]`. -* The matrix ``A_{i,j}`` is given by `-A[i,j]`. +* The matrix ``A_{i,j}`` is given by `A[i,j]`. """ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vector{T}} meta::NLPModels.NLPModelMeta{T,Vector{T}} @@ -93,7 +101,6 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vec C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} A::Matrix{A} b::Vector{T} - b_const::T d_lin::SparseArrays.SparseVector{T,Int64} C_lin::SparseArrays.SparseMatrixCSC{T,Int64} msizes::Vector{Int64} @@ -102,7 +109,6 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vec C::Vector{SparseArrays.SparseMatrixCSC{T,Int}}, A::Matrix{AT}, b::Vector{T}, - b_const::T, d_lin::SparseArrays.SparseVector{T,Int64}, C_lin::SparseArrays.SparseMatrixCSC{T,Int64}, msizes::Vector{Int64}, @@ -111,7 +117,6 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vec model.C = C model.A = A model.b = b - model.b_const = b_const model.d_lin = d_lin model.C_lin = C_lin model.msizes = msizes @@ -184,7 +189,7 @@ function NLPModels.obj(model::Model, x, ::Type{ScalarIndex}) end function NLPModels.obj(model::Model, x) - return model.b_const + NLPModels.obj(model, x, MatrixIndex) + NLPModels.obj(model, x, ScalarIndex) + return NLPModels.obj(model, x, MatrixIndex) + NLPModels.obj(model, x, ScalarIndex) end function NLPModels.grad!(model::Model, _, g) @@ -195,7 +200,7 @@ function NLPModels.grad!(model::Model, _, g) return g end -dual_obj(model::Model, y) = -LinearAlgebra.dot(model.b, y) + model.b_const +dual_obj(model::Model, y) = -LinearAlgebra.dot(model.b, y) function jtprod(model::Model, ::Type{ScalarIndex}, y) return -model.C_lin' * y From bb932634dea4baa98ba567ba7e0f27294b910838 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 16:19:03 +0200 Subject: [PATCH 10/48] Flip sign of A and C --- src/MOI_wrapper.jl | 2 +- src/model.jl | 10 +++++----- src/schur.jl | 4 ++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index dd20682..97e4e18 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -182,7 +182,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} I, J, V, _, _ = A[lmi_id, k] push!(I, i) push!(J, j) - push!(V, v) + push!(V, -v) return end function _add(lmi_id, k, i, j, coef) diff --git a/src/model.jl b/src/model.jl index 8376a1d..3af62d0 100644 --- a/src/model.jl +++ b/src/model.jl @@ -163,7 +163,7 @@ end # Should be only used with `norm` NLPModels.jac(model::Model, ::Type{ScalarIndex}) = model.C_lin -NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = model.A[j.value, i.value] +NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = -model.A[j.value, i.value] NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) @@ -173,7 +173,7 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} end function NLPModels.obj(model::Model, X, i::MatrixIndex) - return -LinearAlgebra.dot(model.C[i.value], X) + return LinearAlgebra.dot(model.C[i.value], X) end function NLPModels.obj(model::Model, x, ::Type{MatrixIndex}) @@ -251,14 +251,14 @@ function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) end _zero!(buffer) for j in eachindex(y) - _add_mul!(buffer, model.A[mat_idx.value, j], y[j]) + _add_mul!(buffer, model.A[mat_idx.value, j], -y[j]) end return buffer end function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y, S) i = mat_idx.value - return jtprod!(buffer[i], model, mat_idx, y) + model.C[i] - S[i] + return jtprod!(buffer[i], model, mat_idx, y) - model.C[i] - S[i] end NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin @@ -275,7 +275,7 @@ end function add_jprod!(model::Model, i::MatrixIndex, V, Jv) for j in 1:num_constraints(model) - Jv[j] -= LinearAlgebra.dot(model.A[i.value, j], V) + Jv[j] += LinearAlgebra.dot(model.A[i.value, j], V) end end diff --git a/src/schur.jl b/src/schur.jl index ec25364..16e2135 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -93,8 +93,8 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T fill!(tmp2, zero(T)) add_jprod!(model, mat_idx, tmp, tmp2) indi = σ[ii:end,ilmi] - BBBB[indi,i] .= -tmp2[indi] - BBBB[i,indi] .= -tmp2[indi] + BBBB[indi,i] .= tmp2[indi] + BBBB[i,indi] .= tmp2[indi] else if !iszero(SparseArrays.nnz(Ai)) if SparseArrays.nnz(Ai) > 1 From c14e6e927f25766d08094c56fda1c746c6981a43 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 16:54:41 +0200 Subject: [PATCH 11/48] Reverse cons --- src/model.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/model.jl b/src/model.jl index 3af62d0..694372e 100644 --- a/src/model.jl +++ b/src/model.jl @@ -163,7 +163,7 @@ end # Should be only used with `norm` NLPModels.jac(model::Model, ::Type{ScalarIndex}) = model.C_lin -NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = -model.A[j.value, i.value] +NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = model.A[j.value, i.value] NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) @@ -206,10 +206,6 @@ function jtprod(model::Model, ::Type{ScalarIndex}, y) return -model.C_lin' * y end -function dual_cons(model::Model, ::Type{ScalarIndex}, y, S) - return model.d_lin - S + jtprod(model, ScalarIndex, y) -end - function buffer_for_jtprod(model::Model) if iszero(num_matrices(model)) return @@ -256,9 +252,13 @@ function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) return buffer end -function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y, S) +function dual_cons(model::Model, ::Type{ScalarIndex}, y) + return model.d_lin + jtprod(model, ScalarIndex, y) +end + +function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y) i = mat_idx.value - return jtprod!(buffer[i], model, mat_idx, y) - model.C[i] - S[i] + return jtprod!(buffer[i], model, mat_idx, y) - model.C[i] end NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin @@ -268,8 +268,7 @@ cons_constant(model::Model) = model.b function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) NLPModels.jprod!(model, x, x, cx) - cx .*= -1 - cx .+= model.b + cx .-= model.b return cx end From fdc0798a92d41e24dfddc2b40d22c9f472e30ec2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 17:10:06 +0200 Subject: [PATCH 12/48] Swap jtprod --- src/model.jl | 8 ++++---- src/schur.jl | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/model.jl b/src/model.jl index 694372e..97ab99e 100644 --- a/src/model.jl +++ b/src/model.jl @@ -203,7 +203,7 @@ end dual_obj(model::Model, y) = -LinearAlgebra.dot(model.b, y) function jtprod(model::Model, ::Type{ScalarIndex}, y) - return -model.C_lin' * y + return model.C_lin' * y end function buffer_for_jtprod(model::Model) @@ -247,18 +247,18 @@ function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) end _zero!(buffer) for j in eachindex(y) - _add_mul!(buffer, model.A[mat_idx.value, j], -y[j]) + _add_mul!(buffer, model.A[mat_idx.value, j], y[j]) end return buffer end function dual_cons(model::Model, ::Type{ScalarIndex}, y) - return model.d_lin + jtprod(model, ScalarIndex, y) + return model.d_lin - jtprod(model, ScalarIndex, y) end function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y) i = mat_idx.value - return jtprod!(buffer[i], model, mat_idx, y) - model.C[i] + return -jtprod!(buffer[i], model, mat_idx, y) - model.C[i] end NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin diff --git a/src/schur.jl b/src/schur.jl index 16e2135..7cd2738 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -167,7 +167,7 @@ end function eval_schur_complement!(buffer, result, model::Model, W, y) fill!(result, zero(eltype(result))) for i in matrix_indices(model) - add_jprod!(model, i, -W[i] * jtprod!(buffer[i.value], model, i, y) * W[i], result) + add_jprod!(model, i, W[i] * jtprod!(buffer[i.value], model, i, y) * W[i], result) end result .+= model.C_lin * (W[ScalarIndex] .* (model.C_lin' * y)) return result From 42e217ca4d29cd91dc26bbdf1d75efc6997bab37 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 17:36:24 +0200 Subject: [PATCH 13/48] Swap C --- src/MOI_wrapper.jl | 10 ++++------ src/model.jl | 8 ++++---- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 97e4e18..d1e23d7 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -141,8 +141,6 @@ function MOI.optimize!(model::Optimizer) ) if model.silent options[:verbose] = 0 - else - options[:verbose] = 1 end SolverCore.solve!(model.solver, model.model; options...) return @@ -182,7 +180,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} I, J, V, _, _ = A[lmi_id, k] push!(I, i) push!(J, j) - push!(V, -v) + push!(V, v) return end function _add(lmi_id, k, i, j, coef) @@ -200,7 +198,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} for k in SparseArrays.nzrange(psd_A, var) lmi_id, i, j = back[SparseArrays.rowvals(psd_A)[k]] col = 1 + var - _add(lmi_id, col, i, j, SparseArrays.nonzeros(psd_A)[k]) + _add(lmi_id, col, i, j, -SparseArrays.nonzeros(psd_A)[k]) end end dest.max_sense = MOI.get(src, MOI.ObjectiveSense()) == MOI.MAX_SENSE @@ -308,7 +306,7 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = dual_obj(optimizer.model, optimizer.solver.stats.multipliers) - return optimizer.objective_constant + (optimizer.max_sense ? -val : val) + return optimizer.objective_constant + (optimizer.max_sense ? val : -val) end function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) @@ -319,7 +317,7 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = optimizer.solver.stats.objective - return optimizer.objective_constant + (optimizer.max_sense ? -val : val) + return optimizer.objective_constant + (optimizer.max_sense ? val : -val) end function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) diff --git a/src/model.jl b/src/model.jl index 97ab99e..e942d53 100644 --- a/src/model.jl +++ b/src/model.jl @@ -185,7 +185,7 @@ function NLPModels.obj(model::Model, x, ::Type{MatrixIndex}) end function NLPModels.obj(model::Model, x, ::Type{ScalarIndex}) - return -LinearAlgebra.dot(model.d_lin, x[ScalarIndex]) + return LinearAlgebra.dot(model.d_lin, x[ScalarIndex]) end function NLPModels.obj(model::Model, x) @@ -200,7 +200,7 @@ function NLPModels.grad!(model::Model, _, g) return g end -dual_obj(model::Model, y) = -LinearAlgebra.dot(model.b, y) +dual_obj(model::Model, y) = LinearAlgebra.dot(model.b, y) function jtprod(model::Model, ::Type{ScalarIndex}, y) return model.C_lin' * y @@ -258,11 +258,11 @@ end function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y) i = mat_idx.value - return -jtprod!(buffer[i], model, mat_idx, y) - model.C[i] + return model.C[i] - jtprod!(buffer[i], model, mat_idx, y) end NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin -NLPModels.grad(model::Model, i::MatrixIndex) = model.C[i.value] +NLPModels.grad(model::Model, i::MatrixIndex) = -model.C[i.value] cons_constant(model::Model) = model.b From 138236021494793ee09251fca914b6f07b44a62d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 18:49:41 +0200 Subject: [PATCH 14/48] Add errors --- src/LowRankOpt.jl | 1 + src/errors.jl | 34 ++++++++++++++++++++++++++++++++++ src/model.jl | 13 ++++++++++--- 3 files changed, 45 insertions(+), 3 deletions(-) create mode 100644 src/errors.jl diff --git a/src/LowRankOpt.jl b/src/LowRankOpt.jl index ea44268..d8835a6 100644 --- a/src/LowRankOpt.jl +++ b/src/LowRankOpt.jl @@ -20,5 +20,6 @@ include("model.jl") include("schur.jl") include("MOI_wrapper.jl") include("BurerMonteiro.jl") +include("errors.jl") end # module LowRankOpt diff --git a/src/errors.jl b/src/errors.jl new file mode 100644 index 0000000..cb48027 --- /dev/null +++ b/src/errors.jl @@ -0,0 +1,34 @@ +""" + errors(model::Model, x, y) + +Return [the 6 standard DIMACS errors](https://plato.asu.edu/dimacs/node3.html). +""" +function errors(model::Model, x; + y = nothing, + primal_err = NLPModels.cons(model, x), + dual_slack = nothing, + dual_err = nothing, + pobj = NLPModels.obj(model, x), + dobj = dual_obj(model, y), +) + b_den = 1 + LinearAlgebra.norm(cons_constant(model), 1) + C_den = 1 + LinearAlgebra.norm(NLPModels.grad(model, ScalarIndex), 1) + sum(matrix_indices(model), init = zero(b_den)) do i + LinearAlgebra.norm(NLPModels.grad(model, i), 1) + end + obj_den = 1 + abs(pobj) + abs(dobj) + ( + LinearAlgebra.norm(primal_err) / b_den, + max(0, -LinearAlgebra.eigmin(x)) / b_den, + isnothing(dual_err) ? zero(b_den) : LinearAlgebra.norm(dual_err) / C_den, + isnothing(dual_slack) ? zero(b_den) : LinearAlgebra.norm(dual_err) / C_den, + (pobj - dobj) / obj_den, + LinearAlgebra.dot(x, dual_slack) / obj_den, + ) +end + +# As defined in https://plato.asu.edu/dimacs/node3.html +function LinearAlgebra.eigmin(x::AbstractSolution{T}) where {T} + return min(minimum(x[ScalarIndex], init = zero(T)) + minimum(matrix_indices(x), init = zero(T)) do i + LinearAlgebra.eigmin(LinearAlgebra.Symmetric(x[i])) + end) +end diff --git a/src/model.jl b/src/model.jl index e942d53..8f7d767 100644 --- a/src/model.jl +++ b/src/model.jl @@ -12,6 +12,7 @@ struct Dimensions offsets::Vector{Int64} end +num_matrices(d::Dimensions) = length(d.side_dimensions) Base.length(d::Dimensions) = d.offsets[end] struct ScalarIndex @@ -31,6 +32,10 @@ struct VectorizedSolution{T} <: AbstractSolution{T} dim::Dimensions end +LinearAlgebra.dot(x::VectorizedSolution, z::VectorizedSolution) = LinearAlgebra.dot(x.x, z.x) + +num_matrices(x::VectorizedSolution) = num_matrices(x.dim) + Base.similar(s::VectorizedSolution) = VectorizedSolution(similar(s.x), s.dim) Base.size(s::VectorizedSolution) = (length(s.dim),) @@ -55,6 +60,7 @@ function Base.view(s::VectorizedSolution, i::MatrixIndex) end Base.setindex!(s::VectorizedSolution, v, i::Integer) = setindex!(s.x, v, i) +Base.getindex(s::VectorizedSolution, i::Integer) = getindex(s.x, i) struct ShapedSolution{T,MT<:AbstractMatrix{T}} <: AbstractSolution{T} scalars::Vector{T} @@ -89,6 +95,7 @@ Model representing the primal-dual pair of problems: & C_\\text{lin}^\\top y \\le d_\\text{lin} \\end{aligned} ``` +This corresponds to [this primal-dual pair](https://plato.asu.edu/dimacs/node2.html). The fields of the `struct` as related to the arrays of the above formulation as follows: * The ``i``th PSD constraint is of size `msize[i] × msisze[i]` @@ -147,7 +154,7 @@ end num_matrices(model::Model) = length(model.C) -function matrix_indices(model::Model) +function matrix_indices(model::Union{Model,AbstractSolution}) return MOI.Utilities.LazyMap{MatrixIndex}(MatrixIndex, Base.OneTo(num_matrices(model))) end @@ -261,8 +268,8 @@ function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y) return model.C[i] - jtprod!(buffer[i], model, mat_idx, y) end -NLPModels.grad(model::Model, ::Type{ScalarIndex}) = -model.d_lin -NLPModels.grad(model::Model, i::MatrixIndex) = -model.C[i.value] +NLPModels.grad(model::Model, ::Type{ScalarIndex}) = model.d_lin +NLPModels.grad(model::Model, i::MatrixIndex) = model.C[i.value] cons_constant(model::Model) = model.b From 8615b951ea879143afe702f4f536a9d99414c02f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 22:23:45 +0200 Subject: [PATCH 15/48] Fix diff check tests --- src/BurerMonteiro.jl | 27 +++++++++---- src/MOI_wrapper.jl | 3 ++ src/factorization.jl | 2 +- test/BurerMonteiro.jl | 92 +++++++++++++++++++++++++++++++++++++------ test/runtests.jl | 1 + 5 files changed, 104 insertions(+), 21 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 04f80e4..de8ab7a 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -50,12 +50,12 @@ struct Model{T,AT} <: NLPModels.AbstractNLPModel{T,Vector{T}} end end -struct Solution{T,VT<:AbstractVector{T}} +struct Solution{T,VT<:AbstractVector{T}} <: AbstractVector{T} x::VT dim::Dimensions end -struct _OuterProduct{T,UT<:AbstractVector{T},VT<:AbstractVector{T}} +struct _OuterProduct{T,UT<:AbstractVector{T},VT<:AbstractVector{T}} <: AbstractVector{T} x::Solution{T,VT} v::Solution{T,UT} end @@ -63,6 +63,18 @@ end Base.eltype(::Type{<:Union{Solution{T},_OuterProduct{T}}}) where {T} = T Base.eltype(x::Union{Solution,_OuterProduct}) = eltype(typeof(x)) +Base.size(s::Solution) = size(s.x) +Base.getindex(s::Solution, i::Integer) = getindex(s.x, i) + +Base.size(s::_OuterProduct) = size(s.x) +function Base.show(io::IO, s::_OuterProduct) + print(io, "_OuterProduct(") + print(io, s.x) + print(io, ", ") + print(io, s.v) + print(io, ")") +end + function Base.getindex(s::Solution, ::Type{LRO.ScalarIndex}) return view(s.x, Base.OneTo(s.dim.num_scalars)) end @@ -82,8 +94,8 @@ function Base.getindex(s::Solution, mi::LRO.MatrixIndex) end function Base.getindex(s::_OuterProduct{T}, i::LRO.MatrixIndex) where {T} - U = s.x[i] - V = s.v[i] + U = s.x[i].factor + V = s.v[i].factor return LRO.AsymmetricFactorization(U, V, FillArrays.Fill(T(2), size(U, 2))) end @@ -125,10 +137,10 @@ function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, J ) for i in LRO.matrix_indices(model.model) U = JtV[i].factor + fill!(U, zero(eltype(U))) for j in eachindex(y) A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) - LinearAlgebra.mul!(U, A, X[i].factor) - U .*= (2y[j]) + U .+= A * X[i].factor .* (2y[j]) end end return Jtv @@ -146,9 +158,10 @@ function NLPModels.hprod!(model::Model{T}, ::AbstractVector, y, v::AbstractVecto Hvi .*= 2obj_weight for j in 1:model.meta.ncon A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) - Hvi .+= A * Vi .* (2y[j]) + Hvi .-= A * Vi .* (2y[j]) end end + Hv end struct Solver{T,ST} <: SolverCore.AbstractOptimizationSolver diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index d1e23d7..50b36a7 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -75,6 +75,7 @@ function MOI.empty!(optimizer::Optimizer) optimizer.solver = nothing optimizer.model = nothing optimizer.lin_cones = nothing + empty!(optimizer.lmi_id) optimizer.objective_constant = NaN return end @@ -271,6 +272,8 @@ function termination_status(solver::SolverCore.AbstractOptimizationSolver) status = MOI.OPTIMAL elseif status == MOI.LOCALLY_INFEASIBLE status = MOI.INFEASIBLE + elseif status == MOI.NORM_LIMIT + status = MOI.DUAL_INFEASIBLE end return status end diff --git a/src/factorization.jl b/src/factorization.jl index 74cdb1c..ebd01e4 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -12,7 +12,7 @@ end function Base.getindex(m::AbstractFactorization, i::Int, j::Int) left = left_factor(m) - right = left_factor(m) + right = right_factor(m) return sum( left[i, k] * m.scaling[k] * right[j, k]' for k in eachindex(m.scaling) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index fa135b9..f8df929 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -2,32 +2,98 @@ using Test using LowRankOpt import Percival -include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) -weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; -model = maxcut(weights, LowRankOpt.Optimizer) - set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) set_attribute(model, "ranks", [1]) optimize!(model) solution_summary(model) +function test_vecprod(f, len, J; tol = 1e-6) + v = ones(len) + @test f(v) ≈ J * v rtol = tol atol = tol + v = -ones(len) + @test f(v) ≈ J * v rtol = tol atol = tol + v = 2ones(len) + @test f(v) ≈ J * v rtol = tol atol = tol + v = -2ones(len) + @test f(v) ≈ J * v rtol = tol atol = tol + v = rand(len) + @test f(v) ≈ J * v rtol = tol atol = tol +end + import NLPModels, FiniteDiff -function jac_check(model, x; tol = 1e-6) +function jac_check(model, x; kws...) f(x) = NLPModels.cons(model, x) J = FiniteDiff.finite_difference_jacobian(f, x) - v = rand(model.meta.nvar) - @test NLPModels.jprod(model, x, v) ≈ J * v rtol = tol atol = tol - v = rand(model.meta.ncon) - @test NLPModels.jtprod(model, x, v) ≈ J' * v rtol = tol atol = tol + @testset "jprod" begin + test_vecprod(v -> NLPModels.jprod(model, x, v), model.meta.nvar, J; kws...) + end + @testset "jtprod" begin + test_vecprod(v -> NLPModels.jtprod(model, x, v), model.meta.ncon, J'; kws...) + end end +function hess_check(model, x; kws...) + obj_weight = rand() + y = rand(model.meta.ncon) + f(x) = obj_weight * NLPModels.obj(model, x) - dot(y, NLPModels.cons(model, x)) + J = FiniteDiff.finite_difference_hessian(f, x) + test_vecprod(v -> NLPModels.hprod(model, x, y, v; obj_weight), model.meta.nvar, J; kws...) +end -@testset "Diff check" begin +using NLPModelsTest +function diff_check(model) b = unsafe_backend(model) - using NLPModelsTest bm = b.solver.model x = rand(bm.meta.nvar) - @test isempty(NLPModelsTest.gradient_check(bm; x)) - jac_check(bm, x) + @testset "Gradient" begin + @test isempty(NLPModelsTest.gradient_check(bm; x)) + end + @testset "Jacobian" begin + jac_check(bm, x) + end + @testset "Hessian" begin + hess_check(bm, x) + end +end + +function full_check(model) + set_attribute(model, "solver", LRO.BurerMonteiro.Solver) + set_attribute(model, "sub_solver", Percival.PercivalSolver) + set_attribute(model, "max_iter", 0) + set_attribute(model, "max_eval", 0) + set_attribute(model, "ranks", [1]) + optimize!(model) + diff_check(model) +end + +@testset "Simple LP" begin + model = Model(LowRankOpt.Optimizer) + @variable(model, x) + @constraint(model, x + 1 >= 0) + @objective(model, Min, x) + full_check(model) +end + +@testset "Simple SDP" begin + model = Model(LowRankOpt.Optimizer) + @variable(model, x) + @constraint(model, x * ones(1, 1) in PSDCone()) + @objective(model, Min, x) + full_check(model) +end + +@testset "Simple SDP" begin + include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) + weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; + model = maxcut(weights, LowRankOpt.Optimizer) + set_attribute(model, "solver", LRO.BurerMonteiro.Solver) + set_attribute(model, "sub_solver", Percival.PercivalSolver) + set_attribute(model, "ranks", [1]) + set_attribute(model, "max_iter", 200) + set_attribute(model, "max_eval", 200) + set_attribute(model, "verbose", 2) + optimize!(model) + solution_summary(model) + diff_check(model) end diff --git a/test/runtests.jl b/test/runtests.jl index fe6daee..cda3aeb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,5 +6,6 @@ using Test include("sets.jl") +include("BurerMonteiro.jl") include("Bridges/runtests.jl") include("Test/runtests.jl") From ee9099a5877a67c2dc7364c3157d1aaed05fb2ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Jun 2025 23:37:36 +0200 Subject: [PATCH 16/48] Fix --- src/BurerMonteiro.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index de8ab7a..28f8424 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -110,7 +110,7 @@ function NLPModels.grad!(model::Model, x::AbstractVector, g::AbstractVector) for i in LRO.matrix_indices(model.model) C = NLPModels.grad(model.model, i) LinearAlgebra.mul!(G[i].factor, C, X[i].factor) - G[i].factor .*= -2 + G[i].factor .*= 2 end return g end From 24e482843182e41cf92736a077c2e610ac0bc449 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 09:28:58 +0200 Subject: [PATCH 17/48] Remove full_check --- test/BurerMonteiro.jl | 37 +++++++++++++++++++++++-------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index f8df929..c1c97f3 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -22,6 +22,13 @@ function test_vecprod(f, len, J; tol = 1e-6) end import NLPModels, FiniteDiff + +function grad_check(model, x; tol = 1e-6) + f(x) = NLPModels.obj(model, x) + g = FiniteDiff.finite_difference_gradient(f, x) + @test NLPModels.grad(model, x) ≈ g rtol = tol atol = tol +end + function jac_check(model, x; kws...) f(x) = NLPModels.cons(model, x) J = FiniteDiff.finite_difference_jacobian(f, x) @@ -47,6 +54,7 @@ function diff_check(model) bm = b.solver.model x = rand(bm.meta.nvar) @testset "Gradient" begin + grad_check(bm, x) @test isempty(NLPModelsTest.gradient_check(bm; x)) end @testset "Jacobian" begin @@ -57,30 +65,31 @@ function diff_check(model) end end -function full_check(model) - set_attribute(model, "solver", LRO.BurerMonteiro.Solver) - set_attribute(model, "sub_solver", Percival.PercivalSolver) - set_attribute(model, "max_iter", 0) - set_attribute(model, "max_eval", 0) - set_attribute(model, "ranks", [1]) - optimize!(model) - diff_check(model) -end - @testset "Simple LP" begin model = Model(LowRankOpt.Optimizer) @variable(model, x) @constraint(model, x + 1 >= 0) @objective(model, Min, x) - full_check(model) + set_attribute(model, "solver", LRO.BurerMonteiro.Solver) + set_attribute(model, "sub_solver", Percival.PercivalSolver) + optimize!(model) + diff_check(model) end @testset "Simple SDP" begin model = Model(LowRankOpt.Optimizer) @variable(model, x) - @constraint(model, x * ones(1, 1) in PSDCone()) - @objective(model, Min, x) - full_check(model) + @constraint(model, (1 - x) * ones(1, 1) in PSDCone()) + @objective(model, Max, x) + set_attribute(model, "solver", LRO.BurerMonteiro.Solver) + set_attribute(model, "sub_solver", Percival.PercivalSolver) + set_attribute(model, "max_iter", 10) + set_attribute(model, "ranks", [1]) + set_attribute(model, "verbose", 2) + optimize!(model) + solution_summary(model) + @test objective_value(model) ≈ 1 + diff_check(model) end @testset "Simple SDP" begin From cea54baaf2ba4814d79ea035e2e12628e927eca9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 12:16:01 +0200 Subject: [PATCH 18/48] Fixes --- test/BurerMonteiro.jl | 10 +++------- test/Project.toml | 1 + 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index c1c97f3..ebca505 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -1,13 +1,8 @@ using Test -using LowRankOpt +using JuMP +import LowRankOpt as LRO import Percival -set_attribute(model, "solver", LRO.BurerMonteiro.Solver) -set_attribute(model, "sub_solver", Percival.PercivalSolver) -set_attribute(model, "ranks", [1]) -optimize!(model) -solution_summary(model) - function test_vecprod(f, len, J; tol = 1e-6) v = ones(len) @test f(v) ≈ J * v rtol = tol atol = tol @@ -72,6 +67,7 @@ end @objective(model, Min, x) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) + set_attribute(model, "ranks", Int[]) optimize!(model) diff_check(model) end diff --git a/test/Project.toml b/test/Project.toml index 3ea19cf..85326a3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" From 95b5877033912c1a36c4fec250edc0095cb3f529 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 17:48:23 +0200 Subject: [PATCH 19/48] Fixes --- src/BurerMonteiro.jl | 19 +++++++++++++++++++ src/MOI_wrapper.jl | 35 ++++++++++++++++++++++++++++------- src/model.jl | 3 +++ test/BurerMonteiro.jl | 24 +++++++++++++++++++----- 4 files changed, 69 insertions(+), 12 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 28f8424..e2c00fc 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -4,6 +4,8 @@ import LinearAlgebra import FillArrays import SolverCore import NLPModels +import MathOptInterface as MOI +import NLPModelsJuMP import LowRankOpt as LRO struct Dimensions @@ -117,6 +119,8 @@ end function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) NLPModels.cons!(model.model, Solution(x, model.dim), cx) + @show cx + cx end function NLPModels.jprod!(model::Model, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) @@ -185,5 +189,20 @@ function SolverCore.solve!( SolverCore.solve!(solver.solver, solver.model, solver.stats; kws...) end +function MOI.get(solver::Solver, attr::MOI.SolverName) + return "BurerMonteiro with " * MOI.get(solver.solver, attr) +end + +function MOI.get(solver::Solver, ::MOI.TerminationStatus) + if isnothing(solver.stats) + return MOI.OPTIMIZE_NOT_CALLED + end + return NLPModelsJuMP.TERMINATION_STATUS[solver.stats.status] + # TODO if the dual is feasible, we can still claim that we found the optimal +end + +function MOI.get(solver::Solver, ::LRO.Solution) + return Solution(solver.stats.solution, solver.model.dim) +end end diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 50b36a7..07655ba 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -80,7 +80,18 @@ function MOI.empty!(optimizer::Optimizer) return end -MOI.get(::Optimizer, ::MOI.SolverName) = "LowRankOpt" +# /!\ FIXME type piracy +function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::MOI.SolverName) + return string(parentmodule(typeof(solver))) +end + +function MOI.get(optimizer::Optimizer, attr::MOI.SolverName) + if isnothing(optimizer.solver) + return "LowRankOpt with no solver loaded yet" + else + return MOI.get(optimizer.solver, attr) + end +end # MOI.RawOptimizerAttribute @@ -263,7 +274,8 @@ end # However, this this is a convex problem, this is actually a global minimum! # We define this function instead of hard-coding `MOI.OPTIMAL` so that # `BurerMonteiro` can override it since it is solving a non-convex formulation. -function termination_status(solver::SolverCore.AbstractOptimizationSolver) +# FIXME This is type piracy, this should be moved to an extension of SolverCore maybe ? +function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::MOI.TerminationStatus) if isnothing(solver.stats) return MOI.OPTIMIZE_NOT_CALLED end @@ -278,11 +290,11 @@ function termination_status(solver::SolverCore.AbstractOptimizationSolver) return status end -function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) +function MOI.get(optimizer::Optimizer, attr::MOI.TerminationStatus) if isnothing(optimizer.solver) return MOI.OPTIMIZE_NOT_CALLED end - return termination_status(optimizer.solver) + return MOI.get(optimizer.solver, attr) end function MOI.get(model::Optimizer, ::MOI.ResultCount) @@ -320,6 +332,7 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = optimizer.solver.stats.objective + @show val return optimizer.objective_constant + (optimizer.max_sense ? val : -val) end @@ -334,7 +347,13 @@ function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) end end -_solution(optimizer::Optimizer) = VectorizedSolution(optimizer.solver.stats.solution, optimizer.model.dim) +struct Solution <: MOI.AbstractModelAttribute end +MOI.is_set_by_optimize(::Solution) = true + +function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::Solution) + return VectorizedSolution(solver.solver.stats.solution, solver.model.dim) +end +MOI.get(optimizer::Optimizer, attr::Solution) = MOI.get(optimizer.solver, attr) function MOI.get( optimizer::Optimizer{T}, @@ -343,7 +362,8 @@ function MOI.get( ) where {T} MOI.check_result_index_bounds(optimizer, attr) lmi_id = optimizer.lmi_id[ci] - return TriangleVectorization(_solution(optimizer)[MatrixIndex(lmi_id)]) + sol = MOI.get(optimizer, Solution()) + return TriangleVectorization(sol[MatrixIndex(lmi_id)]) end function MOI.get( @@ -353,5 +373,6 @@ function MOI.get( ) where {T} MOI.check_result_index_bounds(optimizer, attr) rows = MOI.Utilities.rows(optimizer.lin_cones, ci) - return _solution(optimizer)[ScalarIndex][rows] + sol = MOI.get(optimizer, Solution()) + return sol[ScalarIndex][rows] end diff --git a/src/model.jl b/src/model.jl index 8f7d767..68b0b77 100644 --- a/src/model.jl +++ b/src/model.jl @@ -180,6 +180,9 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} end function NLPModels.obj(model::Model, X, i::MatrixIndex) + @show model.C[i.value] + @show X + @show LinearAlgebra.dot(model.C[i.value], X) return LinearAlgebra.dot(model.C[i.value], X) end diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index ebca505..bee5b61 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -1,4 +1,5 @@ using Test +using LinearAlgebra using JuMP import LowRankOpt as LRO import Percival @@ -61,7 +62,7 @@ function diff_check(model) end @testset "Simple LP" begin - model = Model(LowRankOpt.Optimizer) + model = Model(LRO.Optimizer) @variable(model, x) @constraint(model, x + 1 >= 0) @objective(model, Min, x) @@ -73,25 +74,38 @@ end end @testset "Simple SDP" begin - model = Model(LowRankOpt.Optimizer) + model = Model(LRO.Optimizer) @variable(model, x) - @constraint(model, (1 - x) * ones(1, 1) in PSDCone()) + @constraint(model, con_ref, Symmetric((1 - x) * ones(1, 1)) in PSDCone()) @objective(model, Max, x) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) - set_attribute(model, "max_iter", 10) + set_attribute(model, "max_iter", 20) set_attribute(model, "ranks", [1]) set_attribute(model, "verbose", 2) + @test solver_name(model) == "LowRankOpt with no solver loaded yet" optimize!(model) solution_summary(model) + @test solver_name(model) == "BurerMonteiro with Percival" + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test primal_status(model) == MOI.FEASIBLE_POINT + @test dual_status(model) == MOI.FEASIBLE_POINT + @test value(x) ≈ 1 + t = MOI.get(model, MOI.ConstraintDual(), con_ref) + @test t isa LRO.TriangleVectorization + @test t.matrix isa LRO.Factorization + @test_broken t.matrix ≈ ones(1, 1) + @test_broken only(dual(con_ref)) ≈ 1 + solution_summary(model) @test objective_value(model) ≈ 1 + @test_broken dualobjective_value(model) ≈ 1 diff_check(model) end @testset "Simple SDP" begin include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; - model = maxcut(weights, LowRankOpt.Optimizer) + model = maxcut(weights, LRO.Optimizer) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) set_attribute(model, "ranks", [1]) From c081befed8c624e8365ca3702634732c34dca9a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 17:49:47 +0200 Subject: [PATCH 20/48] SA v0.10 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 9bb671f..d30aa1e 100644 --- a/Project.toml +++ b/Project.toml @@ -24,6 +24,6 @@ MutableArithmetics = "1.6.4" NLPModels = "0.21.5" NLPModelsJuMP = "0.13.2" SolverCore = "0.3.8" -SparseArrays = "1.11.0" +SparseArrays = "1.10" Test = "1.10" julia = "1.10" From 0a2bbcc9f32afdf96466c17878ac2749f49010c9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 18:07:18 +0200 Subject: [PATCH 21/48] Fixes --- src/BurerMonteiro.jl | 7 ++++--- src/MOI_wrapper.jl | 1 - src/model.jl | 3 --- test/BurerMonteiro.jl | 33 +++++++++++++++++++++++++-------- 4 files changed, 29 insertions(+), 15 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index e2c00fc..4a6826e 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -118,9 +118,10 @@ function NLPModels.grad!(model::Model, x::AbstractVector, g::AbstractVector) end function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) - NLPModels.cons!(model.model, Solution(x, model.dim), cx) - @show cx - cx + X = Solution(x, model.dim) + # We don't call `cons!` as we don't want to include `-b` since the constraint + # is encoded as `b <= c(x) <= b` and we just need to specify `c(x)` here. + NLPModels.jprod!(model.model, X, X, cx) end function NLPModels.jprod!(model::Model, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 07655ba..2bfef13 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -332,7 +332,6 @@ end function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} MOI.check_result_index_bounds(optimizer, attr) val = optimizer.solver.stats.objective - @show val return optimizer.objective_constant + (optimizer.max_sense ? val : -val) end diff --git a/src/model.jl b/src/model.jl index 68b0b77..8f7d767 100644 --- a/src/model.jl +++ b/src/model.jl @@ -180,9 +180,6 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} end function NLPModels.obj(model::Model, X, i::MatrixIndex) - @show model.C[i.value] - @show X - @show LinearAlgebra.dot(model.C[i.value], X) return LinearAlgebra.dot(model.C[i.value], X) end diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index bee5b61..295eb57 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -64,12 +64,25 @@ end @testset "Simple LP" begin model = Model(LRO.Optimizer) @variable(model, x) - @constraint(model, x + 1 >= 0) - @objective(model, Min, x) + @constraint(model, con_ref, 1 - x in Nonnegatives()) + @objective(model, Max, x) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) + set_attribute(model, "verbose", 1) set_attribute(model, "ranks", Int[]) + + set_attribute(model, "max_iter", 0) + optimize!(model) + @test termination_status(model) == MOI.ITERATION_LIMIT + diff_check(model) + + set_attribute(model, "max_iter", 10) optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test value(x) ≈ 1 + @test dual(con_ref) ≈ 1 + @test objective_value(model) ≈ 1 + @test dual_objective_value(model) ≈ 1 diff_check(model) end @@ -80,13 +93,19 @@ end @objective(model, Max, x) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) - set_attribute(model, "max_iter", 20) set_attribute(model, "ranks", [1]) set_attribute(model, "verbose", 2) @test solver_name(model) == "LowRankOpt with no solver loaded yet" + + set_attribute(model, "max_iter", 0) optimize!(model) solution_summary(model) @test solver_name(model) == "BurerMonteiro with Percival" + @test termination_status(model) == MOI.ITERATION_LIMIT + diff_check(model) + + set_attribute(model, "max_iter", 10) + optimize!(model) @test termination_status(model) == MOI.LOCALLY_SOLVED @test primal_status(model) == MOI.FEASIBLE_POINT @test dual_status(model) == MOI.FEASIBLE_POINT @@ -94,11 +113,11 @@ end t = MOI.get(model, MOI.ConstraintDual(), con_ref) @test t isa LRO.TriangleVectorization @test t.matrix isa LRO.Factorization - @test_broken t.matrix ≈ ones(1, 1) - @test_broken only(dual(con_ref)) ≈ 1 + @test t.matrix ≈ ones(1, 1) + @test only(dual(con_ref)) ≈ 1 solution_summary(model) @test objective_value(model) ≈ 1 - @test_broken dualobjective_value(model) ≈ 1 + @test dual_objective_value(model) ≈ 1 diff_check(model) end @@ -109,8 +128,6 @@ end set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) set_attribute(model, "ranks", [1]) - set_attribute(model, "max_iter", 200) - set_attribute(model, "max_eval", 200) set_attribute(model, "verbose", 2) optimize!(model) solution_summary(model) From 7a48bbc102fba5bc2cfff23bb19cd0a6d82fdc36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 3 Jun 2025 22:48:41 +0200 Subject: [PATCH 22/48] Test with Dualizaton --- src/BurerMonteiro.jl | 4 ++-- test/BurerMonteiro.jl | 54 ++++++++++++++++++++++++++++++------------- test/Project.toml | 1 + 3 files changed, 41 insertions(+), 18 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 4a6826e..9284e29 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -41,8 +41,8 @@ struct Model{T,AT} <: NLPModels.AbstractNLPModel{T,Vector{T}} ncon = ncon, x0 = rand(n), y0 = rand(ncon), - lvar = fill(-Inf, n), - uvar = fill(Inf, n), + lvar = [fill(zero(T), dim.num_scalars); fill(typemin(T), n - dim.num_scalars)], + uvar = fill(typemax(T), n), lcon = LRO.cons_constant(model), ucon = LRO.cons_constant(model), minimize = true, diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 295eb57..69b37fc 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -2,6 +2,7 @@ using Test using LinearAlgebra using JuMP import LowRankOpt as LRO +using Dualization import Percival function test_vecprod(f, len, J; tol = 1e-6) @@ -47,6 +48,9 @@ end using NLPModelsTest function diff_check(model) b = unsafe_backend(model) + if b isa DualOptimizer + b = b.dual_problem.dual_model.model.optimizer + end bm = b.solver.model x = rand(bm.meta.nvar) @testset "Gradient" begin @@ -61,8 +65,8 @@ function diff_check(model) end end -@testset "Simple LP" begin - model = Model(LRO.Optimizer) +@testset "Simple LP $opt" for opt in [LRO.Optimizer, dual_optimizer(LRO.Optimizer)] + model = Model(dual_optimizer(LRO.Optimizer)) @variable(model, x) @constraint(model, con_ref, 1 - x in Nonnegatives()) @objective(model, Max, x) @@ -84,10 +88,10 @@ end @test objective_value(model) ≈ 1 @test dual_objective_value(model) ≈ 1 diff_check(model) -end +end; -@testset "Simple SDP" begin - model = Model(LRO.Optimizer) +@testset "Simple SDP $opt" for (is_dual, opt) in [(false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer))] + model = Model(opt) @variable(model, x) @constraint(model, con_ref, Symmetric((1 - x) * ones(1, 1)) in PSDCone()) @objective(model, Max, x) @@ -95,12 +99,20 @@ end set_attribute(model, "sub_solver", Percival.PercivalSolver) set_attribute(model, "ranks", [1]) set_attribute(model, "verbose", 2) - @test solver_name(model) == "LowRankOpt with no solver loaded yet" + if is_dual + @test solver_name(model) == "Dual model with LowRankOpt with no solver loaded yet attached" + else + @test solver_name(model) == "LowRankOpt with no solver loaded yet" + end set_attribute(model, "max_iter", 0) optimize!(model) solution_summary(model) - @test solver_name(model) == "BurerMonteiro with Percival" + if is_dual + @test solver_name(model) == "Dual model with BurerMonteiro with Percival attached" + else + @test solver_name(model) == "BurerMonteiro with Percival" + end @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) @@ -111,25 +123,35 @@ end @test dual_status(model) == MOI.FEASIBLE_POINT @test value(x) ≈ 1 t = MOI.get(model, MOI.ConstraintDual(), con_ref) - @test t isa LRO.TriangleVectorization - @test t.matrix isa LRO.Factorization - @test t.matrix ≈ ones(1, 1) + if !is_dual + @test t isa LRO.TriangleVectorization + @test t.matrix isa LRO.Factorization + @test t.matrix ≈ ones(1, 1) + end @test only(dual(con_ref)) ≈ 1 solution_summary(model) @test objective_value(model) ≈ 1 @test dual_objective_value(model) ≈ 1 diff_check(model) -end +end; -@testset "Simple SDP" begin +@testset "Simple SDP $opt" for (is_dual, opt) in [(false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer))] include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; - model = maxcut(weights, LRO.Optimizer) + model = maxcut(weights, opt) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) - set_attribute(model, "ranks", [1]) + set_attribute(model, "ranks", [is_dual ? 2 : 3]) set_attribute(model, "verbose", 2) + + set_attribute(model, "max_iter", 0) optimize!(model) - solution_summary(model) + @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) -end + + set_attribute(model, "max_iter", 10) + optimize!(model) + @test termination_status(model) == MOI.LOCALLY_SOLVED + @test objective_value(model) ≈ 18 + diff_check(model) +end; diff --git a/test/Project.toml b/test/Project.toml index 85326a3..0c5136f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,4 +1,5 @@ [deps] +Dualization = "191a621a-6537-11e9-281d-650236a99e60" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" From bf1f5dd9f28c87f477846dee373759d2c4d7be94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 10:11:19 +0200 Subject: [PATCH 23/48] Add ref to Kojima paper --- src/schur.jl | 30 ++++++++++++++++-------------- 1 file changed, 16 insertions(+), 14 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index 7cd2738..77c74d7 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -1,4 +1,6 @@ -# Adapted from Loraine.jl +# This code computes the Schur complement using the ideas detailed in [FKN97, Section 3] +# It was adapted from dapted from Michal Kocvara's code in +# https://github.com/kocvara/Loraine.jl/blob/bd2821ba830786a78f04081d7e8f5cac25e56cac/src/makeBBBB.jl # Computes `⟨A * W, W * B⟩` for symmetric sparse matrices `A` and `B` function _dot(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, W::AbstractMatrix) @@ -45,19 +47,19 @@ function buffer_for_schur_complement(model::Model, κ) return σ, last_dense end -function makeBBBB_rank1(n,nlmi,B,G) +function makeH_rank1(n,nlmi,B,G) tmp = zeros(Float64, n, n) - BBBB = zeros(Float64, n, n) + H = zeros(Float64, n, n) for ilmi = 1:nlmi BB = transpose(B[ilmi] * G[ilmi]) mul!(tmp,BB',BB) if ilmi == 1 - BBBB = tmp .^ 2 + H = tmp .^ 2 else - BBBB += tmp .^ 2 + H += tmp .^ 2 end end - return BBBB + return H end ######################### @@ -76,7 +78,7 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T σ, last_dense = buffer ilmi = mat_idx.value n = num_constraints(model) - BBBB = zeros(T, n, n) + H = zeros(T, n, n) dim = side_dimension(model, mat_idx) @assert dim == size(W, 1) == size(W, 2) tmp1 = Matrix{T}(undef, size(W, 2), dim) @@ -93,8 +95,8 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T fill!(tmp2, zero(T)) add_jprod!(model, mat_idx, tmp, tmp2) indi = σ[ii:end,ilmi] - BBBB[indi,i] .= tmp2[indi] - BBBB[i,indi] .= tmp2[indi] + H[indi,i] .= tmp2[indi] + H[i,indi] .= tmp2[indi] else if !iszero(SparseArrays.nnz(Ai)) if SparseArrays.nnz(Ai) > 1 @@ -104,9 +106,9 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T if !iszero(SparseArrays.nnz(Aj)) ttt = _dot(Ai, Aj, W) if i >= j - BBBB[i,j] = ttt + H[i,j] = ttt else - BBBB[j,i] = ttt + H[j,i] = ttt end end end @@ -125,9 +127,9 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T vvvj = only(SparseArrays.nonzeros(Ajjj)) ttt = vvvi * W[iiiiAi,iiijAj] * W[jjjiAi,jjjjAj] * vvvj if i >= j - BBBB[i,j] = ttt + H[i,j] = ttt else - BBBB[j,i] = ttt + H[j,i] = ttt end end end @@ -136,7 +138,7 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T end end end - return BBBB + return H end # [HKS24, (5b)] From 9651db2a8f4801441104a2bda438c21608a0150c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 10:11:33 +0200 Subject: [PATCH 24/48] Fix format --- src/BurerMonteiro.jl | 46 +++++++++++++++++------ src/MOI_wrapper.jl | 78 +++++++++++++++++++++++++++----------- src/errors.jl | 30 ++++++++++----- src/factorization.jl | 15 ++++++-- src/model.jl | 66 ++++++++++++++++++++++++-------- src/schur.jl | 87 ++++++++++++++++++++++++++++--------------- test/BurerMonteiro.jl | 43 ++++++++++++++++----- 7 files changed, 261 insertions(+), 104 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 9284e29..4f11d99 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -16,7 +16,8 @@ struct Dimensions end function Dimensions(model::LRO.Model, ranks) - side_dimensions = [LRO.side_dimension(model, i) for i in LRO.matrix_indices(model)] + side_dimensions = + [LRO.side_dimension(model, i) for i in LRO.matrix_indices(model)] num_scalars = LRO.num_scalars(model) offsets = num_scalars .+ [0; cumsum(side_dimensions .* ranks)] return Dimensions(num_scalars, side_dimensions, ranks, offsets) @@ -41,7 +42,10 @@ struct Model{T,AT} <: NLPModels.AbstractNLPModel{T,Vector{T}} ncon = ncon, x0 = rand(n), y0 = rand(ncon), - lvar = [fill(zero(T), dim.num_scalars); fill(typemin(T), n - dim.num_scalars)], + lvar = [ + fill(zero(T), dim.num_scalars); + fill(typemin(T), n - dim.num_scalars) + ], uvar = fill(typemax(T), n), lcon = LRO.cons_constant(model), ucon = LRO.cons_constant(model), @@ -57,7 +61,8 @@ struct Solution{T,VT<:AbstractVector{T}} <: AbstractVector{T} dim::Dimensions end -struct _OuterProduct{T,UT<:AbstractVector{T},VT<:AbstractVector{T}} <: AbstractVector{T} +struct _OuterProduct{T,UT<:AbstractVector{T},VT<:AbstractVector{T}} <: + AbstractVector{T} x::Solution{T,VT} v::Solution{T,UT} end @@ -74,7 +79,7 @@ function Base.show(io::IO, s::_OuterProduct) print(io, s.x) print(io, ", ") print(io, s.v) - print(io, ")") + return print(io, ")") end function Base.getindex(s::Solution, ::Type{LRO.ScalarIndex}) @@ -88,7 +93,7 @@ end function Base.getindex(s::Solution, mi::LRO.MatrixIndex) i = mi.value U = reshape( - view(s.x, (1 + s.dim.offsets[i]):s.dim.offsets[i+1]), + view(s.x, (1+s.dim.offsets[i]):s.dim.offsets[i+1]), s.dim.side_dimensions[i], s.dim.ranks[i], ) @@ -121,18 +126,28 @@ function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) X = Solution(x, model.dim) # We don't call `cons!` as we don't want to include `-b` since the constraint # is encoded as `b <= c(x) <= b` and we just need to specify `c(x)` here. - NLPModels.jprod!(model.model, X, X, cx) + return NLPModels.jprod!(model.model, X, X, cx) end -function NLPModels.jprod!(model::Model, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) +function NLPModels.jprod!( + model::Model, + x::AbstractVector, + v::AbstractVector, + Jv::AbstractVector, +) X = Solution(x, model.dim) V = Solution(v, model.dim) # The second argument is ignored as it is linear so it does # not matter that we give `x` - NLPModels.jprod!(model.model, X, _OuterProduct(X, V), Jv) + return NLPModels.jprod!(model.model, X, _OuterProduct(X, V), Jv) end -function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, Jtv::AbstractVector) +function NLPModels.jtprod!( + model::Model, + x::AbstractVector, + y::AbstractVector, + Jtv::AbstractVector, +) X = Solution(x, model.dim) JtV = Solution(Jtv, model.dim) LinearAlgebra.mul!( @@ -151,7 +166,14 @@ function NLPModels.jtprod!(model::Model, x::AbstractVector, y::AbstractVector, J return Jtv end -function NLPModels.hprod!(model::Model{T}, ::AbstractVector, y, v::AbstractVector, Hv::AbstractVector; obj_weight = one(T)) where {T} +function NLPModels.hprod!( + model::Model{T}, + ::AbstractVector, + y, + v::AbstractVector, + Hv::AbstractVector; + obj_weight = one(T), +) where {T} V = Solution(v, model.dim) HV = Solution(Hv, model.dim) fill!(Hv, zero(eltype(Hv))) @@ -166,7 +188,7 @@ function NLPModels.hprod!(model::Model{T}, ::AbstractVector, y, v::AbstractVecto Hvi .-= A * Vi .* (2y[j]) end end - Hv + return Hv end struct Solver{T,ST} <: SolverCore.AbstractOptimizationSolver @@ -187,7 +209,7 @@ function SolverCore.solve!( model::NLPModels.AbstractNLPModel; # Same as `solver.model.model` kws..., ) - SolverCore.solve!(solver.solver, solver.model, solver.stats; kws...) + return SolverCore.solve!(solver.solver, solver.model, solver.stats; kws...) end function MOI.get(solver::Solver, attr::MOI.SolverName) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 2bfef13..63f040a 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -81,7 +81,10 @@ function MOI.empty!(optimizer::Optimizer) end # /!\ FIXME type piracy -function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::MOI.SolverName) +function MOI.get( + solver::SolverCore.AbstractOptimizationSolver, + ::MOI.SolverName, +) return string(parentmodule(typeof(solver))) end @@ -134,22 +137,30 @@ end function MOI.supports( ::Optimizer, - ::Union{MOI.ObjectiveSense,MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}}, + ::Union{ + MOI.ObjectiveSense, + MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}, + }, ) where {T} return true end const SUPPORTED_CONES = Union{NNG,PSD} -function MOI.supports_constraint(::Optimizer{T}, ::Type{VAF{T}}, ::Type{<:SUPPORTED_CONES}) where {T} +function MOI.supports_constraint( + ::Optimizer{T}, + ::Type{VAF{T}}, + ::Type{<:SUPPORTED_CONES}, +) where {T} return true end SOLVER_OPTIONS = ["solver", "sub_solver", "ranks"] function MOI.optimize!(model::Optimizer) - options = Dict{Symbol, Any}( - Symbol(key) => model.options[key] for key in keys(model.options) if !(key in SOLVER_OPTIONS) + options = Dict{Symbol,Any}( + Symbol(key) => model.options[key] for + key in keys(model.options) if !(key in SOLVER_OPTIONS) ) if model.silent options[:verbose] = 0 @@ -168,21 +179,26 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} C_lin = -convert(SM, C_lin') n = MOI.get(src, MOI.NumberOfVariables()) nlmi = MOI.get(src, MOI.NumberOfConstraints{VAF{T},PSD}()) - A = Matrix{Tuple{Vector{Int64},Vector{Int64},Vector{T},Int64,Int64}}(undef, nlmi, n + 1) + A = Matrix{Tuple{Vector{Int64},Vector{Int64},Vector{T},Int64,Int64}}( + undef, + nlmi, + n + 1, + ) back = Vector{Tuple{Int64,Int64,Int64}}(undef, size(psd_A, 1)) empty!(dest.lmi_id) row = 0 msizes = Int64[] - for (lmi_id, ci) in enumerate(MOI.get(src, MOI.ListOfConstraintIndices{VAF{T},PSD}())) + for (lmi_id, ci) in + enumerate(MOI.get(src, MOI.ListOfConstraintIndices{VAF{T},PSD}())) dest.lmi_id[ci] = lmi_id set = MOI.get(src, MOI.ConstraintSet(), ci) d = set.side_dimension push!(msizes, d) - for k = 1:(n+1) + for k in 1:(n+1) A[lmi_id, k] = (Int64[], Int64[], T[], d, d) end - for j = 1:d - for i = 1:j + for j in 1:d + for i in 1:j row += 1 back[row] = (lmi_id, i, j) end @@ -206,7 +222,7 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} lmi_id, i, j = back[row] _add(lmi_id, 1, i, j, psd_AC.constants[row]) end - for var = 1:n + for var in 1:n for k in SparseArrays.nzrange(psd_A, var) lmi_id, i, j = back[SparseArrays.rowvals(psd_A)[k]] col = 1 + var @@ -223,12 +239,17 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} b = dest.max_sense ? b0 : -b0 # b = max_sense ? -b0 : b0 - AA = SparseArrays.SparseMatrixCSC{T,Int}[SparseArrays.sparse(IJV...) for IJV in A] + AA = SparseArrays.SparseMatrixCSC{T,Int}[ + SparseArrays.sparse(IJV...) for IJV in A + ] dest.model = Model( - AA[:,1], - AA[:,2:end], + AA[:, 1], + AA[:, 2:end], b, - convert(SparseArrays.SparseVector{T,Int64}, SparseArrays.sparsevec(Cd_lin.constants)), + convert( + SparseArrays.SparseVector{T,Int64}, + SparseArrays.sparsevec(Cd_lin.constants), + ), C_lin, msizes, ) @@ -238,8 +259,9 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} options["verb"] = 0 end dest.lin_cones = Cd_lin.sets - options = Dict{Symbol, Any}( - Symbol(key) => dest.options[key] for key in keys(dest.options) if key in SOLVER_OPTIONS && key != "solver" + options = Dict{Symbol,Any}( + Symbol(key) => dest.options[key] for key in keys(dest.options) if + key in SOLVER_OPTIONS && key != "solver" ) dest.solver = dest.options["solver"](dest.model; options...) return MOI.Utilities.identity_index_map(src) @@ -275,7 +297,10 @@ end # We define this function instead of hard-coding `MOI.OPTIMAL` so that # `BurerMonteiro` can override it since it is solving a non-convex formulation. # FIXME This is type piracy, this should be moved to an extension of SolverCore maybe ? -function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::MOI.TerminationStatus) +function MOI.get( + solver::SolverCore.AbstractOptimizationSolver, + ::MOI.TerminationStatus, +) if isnothing(solver.stats) return MOI.OPTIMIZE_NOT_CALLED end @@ -308,7 +333,8 @@ end function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) return MOI.NO_SOLUTION - elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] + elseif MOI.get(optimizer, MOI.TerminationStatus()) in + [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE return MOI.INFEASIBLE_POINT @@ -324,12 +350,19 @@ function MOI.get(optimizer::Optimizer{T}, attr::MOI.ObjectiveValue) where {T} return optimizer.objective_constant + (optimizer.max_sense ? val : -val) end -function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex) +function MOI.get( + optimizer::Optimizer, + attr::MOI.VariablePrimal, + vi::MOI.VariableIndex, +) MOI.check_result_index_bounds(optimizer, attr) return optimizer.solver.stats.multipliers[vi.value] end -function MOI.get(optimizer::Optimizer{T}, attr::MOI.DualObjectiveValue) where {T} +function MOI.get( + optimizer::Optimizer{T}, + attr::MOI.DualObjectiveValue, +) where {T} MOI.check_result_index_bounds(optimizer, attr) val = optimizer.solver.stats.objective return optimizer.objective_constant + (optimizer.max_sense ? val : -val) @@ -338,7 +371,8 @@ end function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) if attr.result_index > MOI.get(optimizer, MOI.ResultCount()) return MOI.NO_SOLUTION - elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] + elseif MOI.get(optimizer, MOI.TerminationStatus()) in + [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT else # TODO diff --git a/src/errors.jl b/src/errors.jl index cb48027..8f7349a 100644 --- a/src/errors.jl +++ b/src/errors.jl @@ -3,7 +3,9 @@ Return [the 6 standard DIMACS errors](https://plato.asu.edu/dimacs/node3.html). """ -function errors(model::Model, x; +function errors( + model::Model, + x; y = nothing, primal_err = NLPModels.cons(model, x), dual_slack = nothing, @@ -12,15 +14,20 @@ function errors(model::Model, x; dobj = dual_obj(model, y), ) b_den = 1 + LinearAlgebra.norm(cons_constant(model), 1) - C_den = 1 + LinearAlgebra.norm(NLPModels.grad(model, ScalarIndex), 1) + sum(matrix_indices(model), init = zero(b_den)) do i - LinearAlgebra.norm(NLPModels.grad(model, i), 1) - end + C_den = + 1 + + LinearAlgebra.norm(NLPModels.grad(model, ScalarIndex), 1) + + sum(matrix_indices(model), init = zero(b_den)) do i + return LinearAlgebra.norm(NLPModels.grad(model, i), 1) + end obj_den = 1 + abs(pobj) + abs(dobj) - ( + return ( LinearAlgebra.norm(primal_err) / b_den, max(0, -LinearAlgebra.eigmin(x)) / b_den, - isnothing(dual_err) ? zero(b_den) : LinearAlgebra.norm(dual_err) / C_den, - isnothing(dual_slack) ? zero(b_den) : LinearAlgebra.norm(dual_err) / C_den, + isnothing(dual_err) ? zero(b_den) : + LinearAlgebra.norm(dual_err) / C_den, + isnothing(dual_slack) ? zero(b_den) : + LinearAlgebra.norm(dual_err) / C_den, (pobj - dobj) / obj_den, LinearAlgebra.dot(x, dual_slack) / obj_den, ) @@ -28,7 +35,10 @@ end # As defined in https://plato.asu.edu/dimacs/node3.html function LinearAlgebra.eigmin(x::AbstractSolution{T}) where {T} - return min(minimum(x[ScalarIndex], init = zero(T)) + minimum(matrix_indices(x), init = zero(T)) do i - LinearAlgebra.eigmin(LinearAlgebra.Symmetric(x[i])) - end) + return min( + minimum(x[ScalarIndex], init = zero(T)) + + minimum(matrix_indices(x), init = zero(T)) do i + return LinearAlgebra.eigmin(LinearAlgebra.Symmetric(x[i])) + end, + ) end diff --git a/src/factorization.jl b/src/factorization.jl index ebd01e4..bd870c8 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -14,8 +14,7 @@ function Base.getindex(m::AbstractFactorization, i::Int, j::Int) left = left_factor(m) right = right_factor(m) return sum( - left[i, k] * m.scaling[k] * right[j, k]' for - k in eachindex(m.scaling) + left[i, k] * m.scaling[k] * right[j, k]' for k in eachindex(m.scaling) ) end @@ -186,7 +185,11 @@ function AsymmetricFactorization( right::AbstractMatrix{T}, scaling::AbstractVector{T}, ) where {T} - return AsymmetricFactorization{T,typeof(left),typeof(scaling)}(left, right, scaling) + return AsymmetricFactorization{T,typeof(left),typeof(scaling)}( + left, + right, + scaling, + ) end function AsymmetricFactorization( @@ -194,7 +197,11 @@ function AsymmetricFactorization( right::AbstractVector{T}, scaling::AbstractArray{T,0}, ) where {T} - return AsymmetricFactorization{T,typeof(left),typeof(scaling)}(left, right, scaling) + return AsymmetricFactorization{T,typeof(left),typeof(scaling)}( + left, + right, + scaling, + ) end left_factor(m::AsymmetricFactorization) = m.left diff --git a/src/model.jl b/src/model.jl index 8f7d767..485d50d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -25,14 +25,21 @@ end abstract type AbstractSolution{T} <: AbstractVector{T} end -Base.getindex(s::AbstractSolution, i::Union{Type{ScalarIndex},MatrixIndex}) = view(s, i) +function Base.getindex( + s::AbstractSolution, + i::Union{Type{ScalarIndex},MatrixIndex}, +) + return view(s, i) +end struct VectorizedSolution{T} <: AbstractSolution{T} x::Vector{T} dim::Dimensions end -LinearAlgebra.dot(x::VectorizedSolution, z::VectorizedSolution) = LinearAlgebra.dot(x.x, z.x) +function LinearAlgebra.dot(x::VectorizedSolution, z::VectorizedSolution) + return LinearAlgebra.dot(x.x, z.x) +end num_matrices(x::VectorizedSolution) = num_matrices(x.dim) @@ -40,11 +47,13 @@ Base.similar(s::VectorizedSolution) = VectorizedSolution(similar(s.x), s.dim) Base.size(s::VectorizedSolution) = (length(s.dim),) -Base.to_index(s::VectorizedSolution, ::Type{ScalarIndex}) = Base.OneTo(s.dim.num_scalars) +function Base.to_index(s::VectorizedSolution, ::Type{ScalarIndex}) + return Base.OneTo(s.dim.num_scalars) +end function Base.to_index(s::VectorizedSolution, mi::MatrixIndex) i = mi.value - return (1 + s.dim.offsets[i]):s.dim.offsets[i+1] + return (1+s.dim.offsets[i]):s.dim.offsets[i+1] end function Base.view(s::VectorizedSolution, ::Type{ScalarIndex}) @@ -67,7 +76,9 @@ struct ShapedSolution{T,MT<:AbstractMatrix{T}} <: AbstractSolution{T} matrices::Vector{MT} end -Base.size(s::ShapedSolution) = (length(s.scalars) + sum(length, s.matrices, init = 0),) +function Base.size(s::ShapedSolution) + return (length(s.scalars) + sum(length, s.matrices, init = 0),) +end Base.view(s::ShapedSolution, ::Type{ScalarIndex}) = s.scalars Base.view(s::ShapedSolution, i::MatrixIndex) = s.matrices[i.value] @@ -102,7 +113,8 @@ The fields of the `struct` as related to the arrays of the above formulation as * The matrix ``C_i`` is given by `C[i]`. * The matrix ``A_{i,j}`` is given by `A[i,j]`. """ -mutable struct Model{T,A<:AbstractMatrix{T}} <: NLPModels.AbstractNLPModel{T,Vector{T}} +mutable struct Model{T,A<:AbstractMatrix{T}} <: + NLPModels.AbstractNLPModel{T,Vector{T}} meta::NLPModels.NLPModelMeta{T,Vector{T}} dim::Dimensions C::Vector{SparseArrays.SparseMatrixCSC{T,Int}} @@ -149,13 +161,19 @@ NLPModels.has_bounds(::Model) = false num_scalars(model::Model) = length(model.d_lin) function scalar_indices(model::Model) - return MOI.Utilities.LazyMap{ScalarIndex}(ScalarIndex, Base.OneTo(num_scalars(model))) + return MOI.Utilities.LazyMap{ScalarIndex}( + ScalarIndex, + Base.OneTo(num_scalars(model)), + ) end num_matrices(model::Model) = length(model.C) function matrix_indices(model::Union{Model,AbstractSolution}) - return MOI.Utilities.LazyMap{MatrixIndex}(MatrixIndex, Base.OneTo(num_matrices(model))) + return MOI.Utilities.LazyMap{MatrixIndex}( + MatrixIndex, + Base.OneTo(num_matrices(model)), + ) end side_dimension(model::Model, i::MatrixIndex) = model.msizes[i.value] @@ -165,13 +183,20 @@ struct ConstraintIndex end num_constraints(model::Model) = length(model.b) function constraint_indices(model::Model) - return MOI.Utilities.LazyMap{ConstraintIndex}(ConstraintIndex, Base.OneTo(num_constraints(model))) + return MOI.Utilities.LazyMap{ConstraintIndex}( + ConstraintIndex, + Base.OneTo(num_constraints(model)), + ) end # Should be only used with `norm` NLPModels.jac(model::Model, ::Type{ScalarIndex}) = model.C_lin -NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) = model.A[j.value, i.value] -NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) = model.C_lin[i.value,:] +function NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) + return model.A[j.value, i.value] +end +function NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) + return model.C_lin[i.value, :] +end function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) return zero(T) @@ -196,7 +221,8 @@ function NLPModels.obj(model::Model, x, ::Type{ScalarIndex}) end function NLPModels.obj(model::Model, x) - return NLPModels.obj(model, x, MatrixIndex) + NLPModels.obj(model, x, ScalarIndex) + return NLPModels.obj(model, x, MatrixIndex) + + NLPModels.obj(model, x, ScalarIndex) end function NLPModels.grad!(model::Model, _, g) @@ -226,12 +252,15 @@ function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) end # FIXME: at some point, switch to dense return sum( - abs.(model.A[mat_idx.value, j]) - for j in 1:num_constraints(model) + abs.(model.A[mat_idx.value, j]) for j in 1:num_constraints(model) ) end -function _add_mul!(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, α) +function _add_mul!( + A::SparseArrays.SparseMatrixCSC, + B::SparseArrays.SparseMatrixCSC, + α, +) for col in axes(A, 2) range_A = SparseArrays.nzrange(A, col) it_A = iterate(range_A) @@ -285,7 +314,12 @@ function add_jprod!(model::Model, i::MatrixIndex, V, Jv) end end -function NLPModels.jprod!(model::Model, _::AbstractVector, v::AbstractVector, Jv::AbstractVector) +function NLPModels.jprod!( + model::Model, + _::AbstractVector, + v::AbstractVector, + Jv::AbstractVector, +) LinearAlgebra.mul!(Jv, model.C_lin, v[ScalarIndex]) for i in matrix_indices(model) add_jprod!(model, i, v[i], Jv) diff --git a/src/schur.jl b/src/schur.jl index 77c74d7..56f74d2 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -3,8 +3,14 @@ # https://github.com/kocvara/Loraine.jl/blob/bd2821ba830786a78f04081d7e8f5cac25e56cac/src/makeBBBB.jl # Computes `⟨A * W, W * B⟩` for symmetric sparse matrices `A` and `B` -function _dot(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, W::AbstractMatrix) - @assert LinearAlgebra.checksquare(W) == LinearAlgebra.checksquare(A) == LinearAlgebra.checksquare(B) +function _dot( + A::SparseArrays.SparseMatrixCSC, + B::SparseArrays.SparseMatrixCSC, + W::AbstractMatrix, +) + @assert LinearAlgebra.checksquare(W) == + LinearAlgebra.checksquare(A) == + LinearAlgebra.checksquare(B) # After these asserts, we know that `A`, `B` and `W` are square and # have the same sizes so we can safely use `@inbounds` result = zero(eltype(A)) @@ -16,11 +22,15 @@ function _dot(A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, if !isempty(nzB) AW = zero(result) for k in nzA - AW += SparseArrays.nonzeros(A)[k] * W[SparseArrays.rowvals(A)[k], j] + AW += + SparseArrays.nonzeros(A)[k] * + W[SparseArrays.rowvals(A)[k], j] end WB = zero(result) for k in nzB - WB += W[i, SparseArrays.rowvals(B)[k]] * SparseArrays.nonzeros(B)[k] + WB += + W[i, SparseArrays.rowvals(B)[k]] * + SparseArrays.nonzeros(B)[k] end result += AW * WB end @@ -38,8 +48,8 @@ function buffer_for_schur_complement(model::Model, κ) for mat_idx in matrix_indices(model) i = mat_idx.value nzA = [SparseArrays.nnz(model.A[i, j]) for j in 1:n] - σ[:,i] = sortperm(nzA, rev = true) - sorted = nzA[σ[:,i]] + σ[:, i] = sortperm(nzA, rev = true) + sorted = nzA[σ[:, i]] last_dense[i] = something(findlast(Base.Fix1(isless, κ), sorted), 0) end @@ -47,12 +57,12 @@ function buffer_for_schur_complement(model::Model, κ) return σ, last_dense end -function makeH_rank1(n,nlmi,B,G) +function makeH_rank1(n, nlmi, B, G) tmp = zeros(Float64, n, n) H = zeros(Float64, n, n) - for ilmi = 1:nlmi + for ilmi in 1:nlmi BB = transpose(B[ilmi] * G[ilmi]) - mul!(tmp,BB',BB) + mul!(tmp, BB', BB) if ilmi == 1 H = tmp .^ 2 else @@ -74,7 +84,12 @@ function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) end ##### -function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T} +function schur_complement( + buffer, + model, + mat_idx, + W::AbstractMatrix{T}, +) where {T} σ, last_dense = buffer ilmi = mat_idx.value n = num_constraints(model) @@ -83,10 +98,10 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T @assert dim == size(W, 1) == size(W, 2) tmp1 = Matrix{T}(undef, size(W, 2), dim) tmp2 = Vector{T}(undef, num_constraints(model)) - tmp = zeros(T, size(W, 2), dim) + tmp = zeros(T, size(W, 2), dim) - for ii = 1:n - i = σ[ii,ilmi] + for ii in 1:n + i = σ[ii, ilmi] Ai = model.A[ilmi, i] if SparseArrays.nnz(Ai) > 0 if ii <= last_dense[ilmi] @@ -94,46 +109,51 @@ function schur_complement(buffer, model, mat_idx, W::AbstractMatrix{T}) where {T LinearAlgebra.mul!(tmp, tmp1, W) fill!(tmp2, zero(T)) add_jprod!(model, mat_idx, tmp, tmp2) - indi = σ[ii:end,ilmi] - H[indi,i] .= tmp2[indi] - H[i,indi] .= tmp2[indi] + indi = σ[ii:end, ilmi] + H[indi, i] .= tmp2[indi] + H[i, indi] .= tmp2[indi] else if !iszero(SparseArrays.nnz(Ai)) if SparseArrays.nnz(Ai) > 1 - @inbounds for jj = ii:n - j = σ[jj,ilmi] + @inbounds for jj in ii:n + j = σ[jj, ilmi] Aj = model.A[ilmi, j] if !iszero(SparseArrays.nnz(Aj)) ttt = _dot(Ai, Aj, W) if i >= j - H[i,j] = ttt + H[i, j] = ttt else - H[j,i] = ttt + H[j, i] = ttt end - end - end + end + end else # A is symmetric iiiiAi = jjjiAi = only(SparseArrays.rowvals(Ai)) vvvi = only(SparseArrays.nonzeros(Ai)) - @inbounds for jj = ii:n - j = σ[jj,ilmi] + @inbounds for jj in ii:n + j = σ[jj, ilmi] Ajjj = model.A[ilmi, j] # As we sort the matrices in decreasing `nnz` order, # the rest of matrices is either zero or have only # one entry if !iszero(SparseArrays.nnz(Ajjj)) - iiijAj = jjjjAj = only(SparseArrays.rowvals(Ajjj)) + iiijAj = + jjjjAj = only(SparseArrays.rowvals(Ajjj)) vvvj = only(SparseArrays.nonzeros(Ajjj)) - ttt = vvvi * W[iiiiAi,iiijAj] * W[jjjiAi,jjjjAj] * vvvj + ttt = + vvvi * + W[iiiiAi, iiijAj] * + W[jjjiAi, jjjjAj] * + vvvj if i >= j - H[i,j] = ttt + H[i, j] = ttt else - H[j,i] = ttt + H[j, i] = ttt end end - end - end + end + end end end end @@ -169,7 +189,12 @@ end function eval_schur_complement!(buffer, result, model::Model, W, y) fill!(result, zero(eltype(result))) for i in matrix_indices(model) - add_jprod!(model, i, W[i] * jtprod!(buffer[i.value], model, i, y) * W[i], result) + add_jprod!( + model, + i, + W[i] * jtprod!(buffer[i.value], model, i, y) * W[i], + result, + ) end result .+= model.C_lin * (W[ScalarIndex] .* (model.C_lin' * y)) return result diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 69b37fc..080d2a7 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -30,19 +30,35 @@ function jac_check(model, x; kws...) f(x) = NLPModels.cons(model, x) J = FiniteDiff.finite_difference_jacobian(f, x) @testset "jprod" begin - test_vecprod(v -> NLPModels.jprod(model, x, v), model.meta.nvar, J; kws...) + test_vecprod( + v -> NLPModels.jprod(model, x, v), + model.meta.nvar, + J; + kws..., + ) end @testset "jtprod" begin - test_vecprod(v -> NLPModels.jtprod(model, x, v), model.meta.ncon, J'; kws...) + test_vecprod( + v -> NLPModels.jtprod(model, x, v), + model.meta.ncon, + J'; + kws..., + ) end end function hess_check(model, x; kws...) obj_weight = rand() y = rand(model.meta.ncon) - f(x) = obj_weight * NLPModels.obj(model, x) - dot(y, NLPModels.cons(model, x)) + f(x) = + obj_weight * NLPModels.obj(model, x) - dot(y, NLPModels.cons(model, x)) J = FiniteDiff.finite_difference_hessian(f, x) - test_vecprod(v -> NLPModels.hprod(model, x, y, v; obj_weight), model.meta.nvar, J; kws...) + return test_vecprod( + v -> NLPModels.hprod(model, x, y, v; obj_weight), + model.meta.nvar, + J; + kws..., + ) end using NLPModelsTest @@ -65,7 +81,8 @@ function diff_check(model) end end -@testset "Simple LP $opt" for opt in [LRO.Optimizer, dual_optimizer(LRO.Optimizer)] +@testset "Simple LP $opt" for opt in + [LRO.Optimizer, dual_optimizer(LRO.Optimizer)] model = Model(dual_optimizer(LRO.Optimizer)) @variable(model, x) @constraint(model, con_ref, 1 - x in Nonnegatives()) @@ -90,7 +107,10 @@ end diff_check(model) end; -@testset "Simple SDP $opt" for (is_dual, opt) in [(false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer))] +@testset "Simple SDP $opt" for (is_dual, opt) in [ + (false, LRO.Optimizer), + (true, dual_optimizer(LRO.Optimizer)), +] model = Model(opt) @variable(model, x) @constraint(model, con_ref, Symmetric((1 - x) * ones(1, 1)) in PSDCone()) @@ -100,7 +120,8 @@ end; set_attribute(model, "ranks", [1]) set_attribute(model, "verbose", 2) if is_dual - @test solver_name(model) == "Dual model with LowRankOpt with no solver loaded yet attached" + @test solver_name(model) == + "Dual model with LowRankOpt with no solver loaded yet attached" else @test solver_name(model) == "LowRankOpt with no solver loaded yet" end @@ -109,7 +130,8 @@ end; optimize!(model) solution_summary(model) if is_dual - @test solver_name(model) == "Dual model with BurerMonteiro with Percival attached" + @test solver_name(model) == + "Dual model with BurerMonteiro with Percival attached" else @test solver_name(model) == "BurerMonteiro with Percival" end @@ -135,7 +157,10 @@ end; diff_check(model) end; -@testset "Simple SDP $opt" for (is_dual, opt) in [(false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer))] +@testset "Simple SDP $opt" for (is_dual, opt) in [ + (false, LRO.Optimizer), + (true, dual_optimizer(LRO.Optimizer)), +] include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; model = maxcut(weights, opt) From 5501063857af3fd5e9e14155f8b386d227519884 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 10:31:18 +0200 Subject: [PATCH 25/48] Add types to arguments --- src/model.jl | 21 +++++++++++---------- src/schur.jl | 8 +++----- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/src/model.jl b/src/model.jl index 485d50d..074529a 100644 --- a/src/model.jl +++ b/src/model.jl @@ -204,11 +204,11 @@ function norm_jac(model::Model{T}, i::MatrixIndex) where {T} return LinearAlgebra.norm(model.A[i.value, :]) end -function NLPModels.obj(model::Model, X, i::MatrixIndex) +function NLPModels.obj(model::Model, X::AbstractMatrix, i::MatrixIndex) return LinearAlgebra.dot(model.C[i.value], X) end -function NLPModels.obj(model::Model, x, ::Type{MatrixIndex}) +function NLPModels.obj(model::Model, x::AbstractVector, ::Type{MatrixIndex}) result = zero(eltype(x)) for i in matrix_indices(model) result += NLPModels.obj(model, x[i], i) @@ -216,16 +216,16 @@ function NLPModels.obj(model::Model, x, ::Type{MatrixIndex}) return result end -function NLPModels.obj(model::Model, x, ::Type{ScalarIndex}) +function NLPModels.obj(model::Model, x::AbstractVector, ::Type{ScalarIndex}) return LinearAlgebra.dot(model.d_lin, x[ScalarIndex]) end -function NLPModels.obj(model::Model, x) +function NLPModels.obj(model::Model, x::AbstractVector) return NLPModels.obj(model, x, MatrixIndex) + NLPModels.obj(model, x, ScalarIndex) end -function NLPModels.grad!(model::Model, _, g) +function NLPModels.grad!(model::Model, _::AbstractVector, g::AbstractVector) copyto!(g[ScalarIndex], model.d_lin) for i in matrix_indices(model) copyto!(g[i], model.C[i.value]) @@ -233,9 +233,9 @@ function NLPModels.grad!(model::Model, _, g) return g end -dual_obj(model::Model, y) = LinearAlgebra.dot(model.b, y) +dual_obj(model::Model, y::AbstractVector) = LinearAlgebra.dot(model.b, y) -function jtprod(model::Model, ::Type{ScalarIndex}, y) +function jtprod(model::Model, ::Type{ScalarIndex}, y::AbstractVector) return model.C_lin' * y end @@ -256,6 +256,7 @@ function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) ) end +# Computes `A .+= B * α` function _add_mul!( A::SparseArrays.SparseMatrixCSC, B::SparseArrays.SparseMatrixCSC, @@ -288,11 +289,11 @@ function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) return buffer end -function dual_cons(model::Model, ::Type{ScalarIndex}, y) +function dual_cons(model::Model, ::Type{ScalarIndex}, y::AbstractVector) return model.d_lin - jtprod(model, ScalarIndex, y) end -function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y) +function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y::AbstractVector) i = mat_idx.value return model.C[i] - jtprod!(buffer[i], model, mat_idx, y) end @@ -308,7 +309,7 @@ function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) return cx end -function add_jprod!(model::Model, i::MatrixIndex, V, Jv) +function add_jprod!(model::Model, i::MatrixIndex, V::AbstractVector, Jv::AbstractVector) for j in 1:num_constraints(model) Jv[j] += LinearAlgebra.dot(model.A[i.value, j], V) end diff --git a/src/schur.jl b/src/schur.jl index 56f74d2..c63d735 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -1,4 +1,6 @@ # This code computes the Schur complement using the ideas detailed in [FKN97, Section 3] +# This is useful to compute search direction in primal-dual interior-point methods for semidefinite programs [FKN97] +# [FKN97] Fujisawa, Katsuki, Masakazu Kojima, and Kazuhide Nakata. "Exploiting sparsity in primal-dual interior-point methods for semidefinite programming." Mathematical Programming 79 (1997): 235-253. # It was adapted from dapted from Michal Kocvara's code in # https://github.com/kocvara/Loraine.jl/blob/bd2821ba830786a78f04081d7e8f5cac25e56cac/src/makeBBBB.jl @@ -63,11 +65,7 @@ function makeH_rank1(n, nlmi, B, G) for ilmi in 1:nlmi BB = transpose(B[ilmi] * G[ilmi]) mul!(tmp, BB', BB) - if ilmi == 1 - H = tmp .^ 2 - else - H += tmp .^ 2 - end + H .+= tmp .^ 2 end return H end From 022e85a3aefceb3cab226282677475c5c38a0ec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 10:33:30 +0200 Subject: [PATCH 26/48] fix format --- src/model.jl | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/model.jl b/src/model.jl index 074529a..5b66751 100644 --- a/src/model.jl +++ b/src/model.jl @@ -293,7 +293,12 @@ function dual_cons(model::Model, ::Type{ScalarIndex}, y::AbstractVector) return model.d_lin - jtprod(model, ScalarIndex, y) end -function dual_cons!(buffer, model::Model, mat_idx::MatrixIndex, y::AbstractVector) +function dual_cons!( + buffer, + model::Model, + mat_idx::MatrixIndex, + y::AbstractVector, +) i = mat_idx.value return model.C[i] - jtprod!(buffer[i], model, mat_idx, y) end @@ -309,7 +314,12 @@ function NLPModels.cons!(model::Model, x::AbstractVector, cx::AbstractVector) return cx end -function add_jprod!(model::Model, i::MatrixIndex, V::AbstractVector, Jv::AbstractVector) +function add_jprod!( + model::Model, + i::MatrixIndex, + V::AbstractVector, + Jv::AbstractVector, +) for j in 1:num_constraints(model) Jv[j] += LinearAlgebra.dot(model.A[i.value, j], V) end From dc329c7b319f8eddcacacb5039c14eda1b9e1874 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 11:39:17 +0200 Subject: [PATCH 27/48] Fix --- src/model.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/model.jl b/src/model.jl index 5b66751..c728bdb 100644 --- a/src/model.jl +++ b/src/model.jl @@ -317,7 +317,7 @@ end function add_jprod!( model::Model, i::MatrixIndex, - V::AbstractVector, + V::AbstractMatrix, Jv::AbstractVector, ) for j in 1:num_constraints(model) From c7cac6272dec9f3b5d7d31b6af5cb41fd9ee787a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 11:56:49 +0200 Subject: [PATCH 28/48] Alloc 20 iterations --- test/BurerMonteiro.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 080d2a7..10fcc2d 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -157,11 +157,11 @@ end; diff_check(model) end; -@testset "Simple SDP $opt" for (is_dual, opt) in [ +include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) +@testset "Max-CUT $opt" for (is_dual, opt) in [ (false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer)), ] - include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; model = maxcut(weights, opt) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) @@ -174,7 +174,7 @@ end; @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) - set_attribute(model, "max_iter", 10) + set_attribute(model, "max_iter", 20) optimize!(model) @test termination_status(model) == MOI.LOCALLY_SOLVED @test objective_value(model) ≈ 18 From e0cdf66d905a9024870651bc56e3e17613aa78e6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 21:24:00 +0200 Subject: [PATCH 29/48] Add tests --- src/BurerMonteiro.jl | 2 +- src/MOI_wrapper.jl | 22 +++++----------------- test/BurerMonteiro.jl | 16 +++++++++++++++- 3 files changed, 21 insertions(+), 19 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 4f11d99..675b300 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -216,7 +216,7 @@ function MOI.get(solver::Solver, attr::MOI.SolverName) return "BurerMonteiro with " * MOI.get(solver.solver, attr) end -function MOI.get(solver::Solver, ::MOI.TerminationStatus) +function MOI.get(solver::Solver, ::LRO.ConvexTerminationStatus) if isnothing(solver.stats) return MOI.OPTIMIZE_NOT_CALLED end diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 63f040a..80263f4 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -103,17 +103,11 @@ function MOI.supports(::Optimizer, param::MOI.RawOptimizerAttribute) end function MOI.set(optimizer::Optimizer, param::MOI.RawOptimizerAttribute, value) - if !MOI.supports(optimizer, param) - throw(MOI.UnsupportedAttribute(param)) - end optimizer.options[param.name] = value return end function MOI.get(optimizer::Optimizer, param::MOI.RawOptimizerAttribute) - if !MOI.supports(optimizer, param) - throw(MOI.UnsupportedAttribute(param)) - end return optimizer.options[param.name] end @@ -128,11 +122,6 @@ end MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent -function MOI.set(optimizer::Optimizer, ::MOI.ObjectiveSense, value::Bool) - optimizer.max_sense = value - return -end - # MOI.supports function MOI.supports( @@ -255,9 +244,6 @@ function MOI.copy_to(dest::Optimizer{T}, src::OptimizerCache{T}) where {T} ) # FIXME this does not work if an option is changed between `MOI.copy_to` and `MOI.optimize!` options = copy(dest.options) - if dest.silent - options["verb"] = 0 - end dest.lin_cones = Cd_lin.sets options = Dict{Symbol,Any}( Symbol(key) => dest.options[key] for key in keys(dest.options) if @@ -296,7 +282,9 @@ end # However, this this is a convex problem, this is actually a global minimum! # We define this function instead of hard-coding `MOI.OPTIMAL` so that # `BurerMonteiro` can override it since it is solving a non-convex formulation. -# FIXME This is type piracy, this should be moved to an extension of SolverCore maybe ? +struct ConvexTerminationStatus <: MOI.AbstractModelAttribute end +MOI.is_set_by_optimize(::ConvexTerminationStatus) = true + function MOI.get( solver::SolverCore.AbstractOptimizationSolver, ::MOI.TerminationStatus, @@ -315,11 +303,11 @@ function MOI.get( return status end -function MOI.get(optimizer::Optimizer, attr::MOI.TerminationStatus) +function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) if isnothing(optimizer.solver) return MOI.OPTIMIZE_NOT_CALLED end - return MOI.get(optimizer.solver, attr) + return MOI.get(optimizer.solver, ConvexTerminationStatus()) end function MOI.get(model::Optimizer, ::MOI.ResultCount) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 10fcc2d..811a339 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -89,7 +89,7 @@ end @objective(model, Max, x) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) - set_attribute(model, "verbose", 1) + set_silent(model) set_attribute(model, "ranks", Int[]) set_attribute(model, "max_iter", 0) @@ -104,6 +104,7 @@ end @test dual(con_ref) ≈ 1 @test objective_value(model) ≈ 1 @test dual_objective_value(model) ≈ 1 + @test abs(MOI.get(model, LRO.RawStatus(:solution))[1]) < 1e-6 diff_check(model) end; @@ -137,6 +138,7 @@ end; end @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) + @test MOI.get(backend(model), MOI.RawOptimizerAttribute("max_iter")) == 0 set_attribute(model, "max_iter", 10) optimize!(model) @@ -180,3 +182,15 @@ include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) @test objective_value(model) ≈ 18 diff_check(model) end; + +@testset "MOI runtests" begin + model = LRO.Optimizer() + MOI.set(model, MOI.RawOptimizerAttribute("solver"), LRO.BurerMonteiro.Solver) + MOI.set(model, MOI.RawOptimizerAttribute("sub_solver"), Percival.PercivalSolver) + config = MOI.Test.Config() + MOI.Test.runtests( + model, + config; + include = ["Silent"], + ) +end; \ No newline at end of file From 17349ecea000357caeb451738336dda2fb9c0b9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 4 Jun 2025 21:45:06 +0200 Subject: [PATCH 30/48] fix --- test/BurerMonteiro.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 811a339..7fdd5de 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -185,12 +185,16 @@ end; @testset "MOI runtests" begin model = LRO.Optimizer() - MOI.set(model, MOI.RawOptimizerAttribute("solver"), LRO.BurerMonteiro.Solver) - MOI.set(model, MOI.RawOptimizerAttribute("sub_solver"), Percival.PercivalSolver) - config = MOI.Test.Config() - MOI.Test.runtests( + MOI.set( + model, + MOI.RawOptimizerAttribute("solver"), + LRO.BurerMonteiro.Solver, + ) + MOI.set( model, - config; - include = ["Silent"], + MOI.RawOptimizerAttribute("sub_solver"), + Percival.PercivalSolver, ) -end; \ No newline at end of file + config = MOI.Test.Config() + MOI.Test.runtests(model, config; include = ["Silent"]) +end; From bd43c06fa92b4f5fb70d34a2f22210781a8b6660 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 10:14:53 +0200 Subject: [PATCH 31/48] Improve coverage of MOI wrapper --- src/MOI_wrapper.jl | 18 +++++++------- test/BurerMonteiro.jl | 55 ++++++++++++++++++++++++++++++++++++++++--- test/Project.toml | 1 + 3 files changed, 63 insertions(+), 11 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 80263f4..e469702 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -287,18 +287,17 @@ MOI.is_set_by_optimize(::ConvexTerminationStatus) = true function MOI.get( solver::SolverCore.AbstractOptimizationSolver, - ::MOI.TerminationStatus, + ::ConvexTerminationStatus, ) - if isnothing(solver.stats) - return MOI.OPTIMIZE_NOT_CALLED - end status = NLPModelsJuMP.TERMINATION_STATUS[solver.stats.status] if status == MOI.LOCALLY_SOLVED status = MOI.OPTIMAL elseif status == MOI.LOCALLY_INFEASIBLE - status = MOI.INFEASIBLE - elseif status == MOI.NORM_LIMIT + # Since we solve the dual, we need to dualize the status status = MOI.DUAL_INFEASIBLE + elseif status == MOI.NORM_LIMIT + # Since we solve the dual, we need to dualize the status + status = MOI.INFEASIBLE end return status end @@ -324,7 +323,8 @@ function MOI.get(optimizer::Optimizer, attr::MOI.PrimalStatus) elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT - elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE + elseif MOI.get(optimizer, MOI.TerminationStatus()) in + [MOI.INFEASIBLE, MOI.LOCALLY_INFEASIBLE] return MOI.INFEASIBLE_POINT else # TODO @@ -362,6 +362,8 @@ function MOI.get(optimizer::Optimizer, attr::MOI.DualStatus) elseif MOI.get(optimizer, MOI.TerminationStatus()) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED] return MOI.FEASIBLE_POINT + elseif MOI.get(optimizer, MOI.TerminationStatus()) == MOI.DUAL_INFEASIBLE + return MOI.INFEASIBLE_POINT else # TODO return MOI.UNKNOWN_RESULT_STATUS @@ -372,7 +374,7 @@ struct Solution <: MOI.AbstractModelAttribute end MOI.is_set_by_optimize(::Solution) = true function MOI.get(solver::SolverCore.AbstractOptimizationSolver, ::Solution) - return VectorizedSolution(solver.solver.stats.solution, solver.model.dim) + return VectorizedSolution(solver.stats.solution, solver.model.dim) end MOI.get(optimizer::Optimizer, attr::Solution) = MOI.get(optimizer.solver, attr) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 7fdd5de..09f8ce2 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -3,6 +3,7 @@ using LinearAlgebra using JuMP import LowRankOpt as LRO using Dualization +import SolverCore import Percival function test_vecprod(f, len, J; tol = 1e-6) @@ -104,7 +105,12 @@ end @test dual(con_ref) ≈ 1 @test objective_value(model) ≈ 1 @test dual_objective_value(model) ≈ 1 - @test abs(MOI.get(model, LRO.RawStatus(:solution))[1]) < 1e-6 + raw_sol = MOI.get(model, LRO.RawStatus(:solution)) + sol = MOI.get(model, LRO.Solution()) + @test raw_sol isa Vector{Float64} + @test sol isa LRO.BurerMonteiro.Solution{Float64,Vector{Float64}} + @test sol == raw_sol + @test abs(sol[1]) < 1e-6 diff_check(model) end; @@ -138,7 +144,10 @@ end; end @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) - @test MOI.get(backend(model), MOI.RawOptimizerAttribute("max_iter")) == 0 + if !is_dual # See https://github.com/jump-dev/Dualization.jl/issues/195 + @test MOI.supports(unsafe_backend(model), MOI.RawOptimizerAttribute("max_iter")) + end + @test MOI.get(unsafe_backend(model), MOI.RawOptimizerAttribute("max_iter")) == 0 set_attribute(model, "max_iter", 10) optimize!(model) @@ -160,11 +169,11 @@ end; end; include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) +weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; @testset "Max-CUT $opt" for (is_dual, opt) in [ (false, LRO.Optimizer), (true, dual_optimizer(LRO.Optimizer)), ] - weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; model = maxcut(weights, opt) set_attribute(model, "solver", LRO.BurerMonteiro.Solver) set_attribute(model, "sub_solver", Percival.PercivalSolver) @@ -183,6 +192,13 @@ include(joinpath(dirname(@__DIR__), "examples", "maxcut.jl")) diff_check(model) end; +@testset "ResultCount" begin + model = LRO.Optimizer() + @test MOI.get(model, MOI.ResultCount()) == 0 + @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION + @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION +end + @testset "MOI runtests" begin model = LRO.Optimizer() MOI.set( @@ -198,3 +214,36 @@ end; config = MOI.Test.Config() MOI.Test.runtests(model, config; include = ["Silent"]) end; + +struct ConvexSolver{T} <: SolverCore.AbstractOptimizationSolver + model::LRO.Model{T} + stats::SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any} +end + +function ConvexSolver(model::LRO.Model) + stats = SolverCore.GenericExecutionStats(model) + return ConvexSolver(model, stats) +end + +function SolverCore.solve!( + ::ConvexSolver, + ::LRO.Model, +) + return +end + +@testset "ConvexSolver" begin + model = maxcut(weights, LRO.Optimizer) + set_attribute(model, "solver", ConvexSolver) + b = unsafe_backend(model) + optimize!(model) + b.solver.stats.status = :first_order + @test termination_status(model) == MOI.OPTIMAL + @test MOI.get(model, LRO.Solution()) isa LRO.VectorizedSolution{Float64} + b.solver.stats.status = :infeasible + @test termination_status(model) == MOI.DUAL_INFEASIBLE + @test dual_status(model) == MOI.INFEASIBLE_POINT + b.solver.stats.status = :unbounded + @test termination_status(model) == MOI.INFEASIBLE + @test primal_status(model) == MOI.INFEASIBLE_POINT +end diff --git a/test/Project.toml b/test/Project.toml index 0c5136f..6bfa942 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -8,4 +8,5 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From af6d598d8287bb2dc009b42122bbbb2ae20f705f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 10:40:54 +0200 Subject: [PATCH 32/48] Fix format --- test/BurerMonteiro.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 09f8ce2..7d2bf9e 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -145,9 +145,15 @@ end; @test termination_status(model) == MOI.ITERATION_LIMIT diff_check(model) if !is_dual # See https://github.com/jump-dev/Dualization.jl/issues/195 - @test MOI.supports(unsafe_backend(model), MOI.RawOptimizerAttribute("max_iter")) + @test MOI.supports( + unsafe_backend(model), + MOI.RawOptimizerAttribute("max_iter"), + ) end - @test MOI.get(unsafe_backend(model), MOI.RawOptimizerAttribute("max_iter")) == 0 + @test MOI.get( + unsafe_backend(model), + MOI.RawOptimizerAttribute("max_iter"), + ) == 0 set_attribute(model, "max_iter", 10) optimize!(model) @@ -225,10 +231,7 @@ function ConvexSolver(model::LRO.Model) return ConvexSolver(model, stats) end -function SolverCore.solve!( - ::ConvexSolver, - ::LRO.Model, -) +function SolverCore.solve!(::ConvexSolver, ::LRO.Model) return end From d8d8c6f36ee0c1b4c92f850c54001927728212c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 11:12:13 +0200 Subject: [PATCH 33/48] Add coverage --- src/MOI_wrapper.jl | 8 ++++++-- test/BurerMonteiro.jl | 1 + 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index e469702..d429d07 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -302,11 +302,15 @@ function MOI.get( return status end -function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) +function MOI.get(optimizer::Optimizer, attr::ConvexTerminationStatus) if isnothing(optimizer.solver) return MOI.OPTIMIZE_NOT_CALLED end - return MOI.get(optimizer.solver, ConvexTerminationStatus()) + return MOI.get(optimizer.solver, attr) +end + +function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus) + return MOI.get(optimizer, ConvexTerminationStatus()) end function MOI.get(model::Optimizer, ::MOI.ResultCount) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 7d2bf9e..3c73137 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -241,6 +241,7 @@ end b = unsafe_backend(model) optimize!(model) b.solver.stats.status = :first_order + @test MOI.get(model, LRO.ConvexTerminationStatus()) == MOI.OPTIMAL @test termination_status(model) == MOI.OPTIMAL @test MOI.get(model, LRO.Solution()) isa LRO.VectorizedSolution{Float64} b.solver.stats.status = :infeasible From a37de470c4fe41eb1a6f223b9183d3d1b8eb3f3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 11:34:03 +0200 Subject: [PATCH 34/48] Complete coverage of src/BurerMonteiro --- src/BurerMonteiro.jl | 7 +++---- test/BurerMonteiro.jl | 3 +++ 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index 675b300..d2122a3 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -79,7 +79,8 @@ function Base.show(io::IO, s::_OuterProduct) print(io, s.x) print(io, ", ") print(io, s.v) - return print(io, ")") + print(io, ")") + return end function Base.getindex(s::Solution, ::Type{LRO.ScalarIndex}) @@ -217,11 +218,9 @@ function MOI.get(solver::Solver, attr::MOI.SolverName) end function MOI.get(solver::Solver, ::LRO.ConvexTerminationStatus) - if isnothing(solver.stats) - return MOI.OPTIMIZE_NOT_CALLED - end return NLPModelsJuMP.TERMINATION_STATUS[solver.stats.status] # TODO if the dual is feasible, we can still claim that we found the optimal + # and turn `LOCALLY_SOLVED` into `OPTIMAL` end function MOI.get(solver::Solver, ::LRO.Solution) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 3c73137..e3a403d 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -109,6 +109,9 @@ end sol = MOI.get(model, LRO.Solution()) @test raw_sol isa Vector{Float64} @test sol isa LRO.BurerMonteiro.Solution{Float64,Vector{Float64}} + outer = LRO.BurerMonteiro._OuterProduct(sol, sol) + @test length(outer) == length(sol) + @test sprint(show, outer) == "_OuterProduct($sol, $sol)" @test sol == raw_sol @test abs(sol[1]) < 1e-6 diff_check(model) From ce108921799db708f7c6ef51bafd057a78b6a003 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 11:40:51 +0200 Subject: [PATCH 35/48] Complete coverage of factorization --- src/factorization.jl | 7 ++++++- test/sets.jl | 14 ++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/src/factorization.jl b/src/factorization.jl index bd870c8..d3d858e 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -166,7 +166,7 @@ struct AsymmetricFactorization{ end if length(scaling) != size(left, 2) error( - "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", + "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(left, 2))` of factor", ) end return new{T,F,S}(left, right, scaling) @@ -176,6 +176,11 @@ struct AsymmetricFactorization{ right::AbstractVector{T}, scaling::AbstractArray{T,0}, ) where {T,F<:AbstractVector{T},S<:AbstractArray{T,0}} + if length(left) != length(right) + error( + "Length `$(length(left))` of left factor does not match the length `$(length(right))` of right factor", + ) + end return new{T,F,S}(left, right, scaling) end end diff --git a/test/sets.jl b/test/sets.jl index 5c7a4b8..d8d7ebf 100644 --- a/test/sets.jl +++ b/test/sets.jl @@ -46,11 +46,25 @@ function test_inconsistent_length() "Length `1` of diagonal does not match number of columns `2` of factor", ) @test_throws err LRO.Factorization(ones(1, 2), [1.0]) + err = ErrorException( + "Size `(2, 1)` of left factor does not match size `(2, 2)` of right factor", + ) + @test_throws err LRO.AsymmetricFactorization(ones(2, 1), ones(2, 2), [1.0]) + err = ErrorException( + "Length `1` of diagonal does not match number of columns `2` of factor", + ) + @test_throws err LRO.AsymmetricFactorization(ones(2, 2), ones(2, 2), [1.0]) + err = ErrorException( + "Length `2` of left factor does not match the length `1` of right factor", + ) + @test_throws err LRO.AsymmetricFactorization(ones(2), ones(1), ones(tuple())) end function test_factorizations() f = [1, 2] + g = [3, 4] _test_factorization(f * f', LRO.positive_semidefinite_factorization(f)) + _test_factorization(5 * f * g', LRO.AsymmetricFactorization(f, g, 5 * ones(Int, tuple()))) _test_factorization(2 * f * f', LRO.Factorization(f, 2)) F = [1 2; 3 4; 5 6] d = [7, 8] From a8e15e4e6c79891ccdaf15b9a25b2faae7d3c920 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 14:33:41 +0200 Subject: [PATCH 36/48] Fix format --- test/sets.jl | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/test/sets.jl b/test/sets.jl index d8d7ebf..e43c319 100644 --- a/test/sets.jl +++ b/test/sets.jl @@ -57,14 +57,21 @@ function test_inconsistent_length() err = ErrorException( "Length `2` of left factor does not match the length `1` of right factor", ) - @test_throws err LRO.AsymmetricFactorization(ones(2), ones(1), ones(tuple())) + @test_throws err LRO.AsymmetricFactorization( + ones(2), + ones(1), + ones(tuple()), + ) end function test_factorizations() f = [1, 2] g = [3, 4] _test_factorization(f * f', LRO.positive_semidefinite_factorization(f)) - _test_factorization(5 * f * g', LRO.AsymmetricFactorization(f, g, 5 * ones(Int, tuple()))) + _test_factorization( + 5 * f * g', + LRO.AsymmetricFactorization(f, g, 5 * ones(Int, tuple())), + ) _test_factorization(2 * f * f', LRO.Factorization(f, 2)) F = [1 2; 3 4; 5 6] d = [7, 8] From b72d43d597521a2476d7584e630975ac95c9136d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 14:54:07 +0200 Subject: [PATCH 37/48] Add tests for errors --- test/BurerMonteiro.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index e3a403d..ac4c713 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -253,4 +253,14 @@ end b.solver.stats.status = :unbounded @test termination_status(model) == MOI.INFEASIBLE @test primal_status(model) == MOI.INFEASIBLE_POINT + X = LRO.VectorizedSolution(collect(1:b.model.meta.nvar), b.model.dim) + y = collect(1:b.model.meta.ncon) + err = LRO.errors(b.solver.model, X; y, dual_slack = X, dual_err = X) + @test length(err) == 6 + @test err[1] ≈ 4.155017729878046 + @test err[2] ≈ 0.318237296391563 + @test err[3] ≈ 70/9 + @test err[4] ≈ 70/9 + @test err[5] ≈ 0.92 + @test err[6] ≈ 392.0 end From dac727921332934e401139049e7121f3cbc9e6d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 15:07:22 +0200 Subject: [PATCH 38/48] Relax tol --- test/BurerMonteiro.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index ac4c713..e490ab8 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -197,7 +197,7 @@ weights = [0 5 7 6; 5 0 0 1; 7 0 0 1; 6 1 1 0]; set_attribute(model, "max_iter", 20) optimize!(model) @test termination_status(model) == MOI.LOCALLY_SOLVED - @test objective_value(model) ≈ 18 + @test objective_value(model) ≈ 18 rtol = 1e-6 diff_check(model) end; From 849225e28f067f99d736b5801da487bcab0c2f5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 22:20:54 +0200 Subject: [PATCH 39/48] Add tests --- src/schur.jl | 109 ++++++++++++++++++------------------------ test/BurerMonteiro.jl | 18 +++++++ 2 files changed, 65 insertions(+), 62 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index c63d735..a7f96d4 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -53,25 +53,13 @@ function buffer_for_schur_complement(model::Model, κ) σ[:, i] = sortperm(nzA, rev = true) sorted = nzA[σ[:, i]] + # Last index for which nnz > κ last_dense[i] = something(findlast(Base.Fix1(isless, κ), sorted), 0) end return σ, last_dense end -function makeH_rank1(n, nlmi, B, G) - tmp = zeros(Float64, n, n) - H = zeros(Float64, n, n) - for ilmi in 1:nlmi - BB = transpose(B[ilmi] * G[ilmi]) - mul!(tmp, BB', BB) - H .+= tmp .^ 2 - end - return H -end - -######################### - function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) n = num_constraints(model) H = zeros(eltype(eltype(W)), n, n) @@ -81,24 +69,23 @@ function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) return H end -##### function schur_complement( buffer, - model, - mat_idx, + model::Model, + mat_idx::MatrixIndex, W::AbstractMatrix{T}, ) where {T} σ, last_dense = buffer ilmi = mat_idx.value - n = num_constraints(model) + n = model.meta.ncon H = zeros(T, n, n) dim = side_dimension(model, mat_idx) - @assert dim == size(W, 1) == size(W, 2) - tmp1 = Matrix{T}(undef, size(W, 2), dim) - tmp2 = Vector{T}(undef, num_constraints(model)) - tmp = zeros(T, size(W, 2), dim) + @assert dim == LinearAlgebra.checksquare(W) + tmp1 = Matrix{T}(undef, dim, dim) + tmp2 = Vector{T}(undef, n) + tmp = zeros(T, dim, dim) - for ii in 1:n + for ii in axes(H, 1) i = σ[ii, ilmi] Ai = model.A[ilmi, i] if SparseArrays.nnz(Ai) > 0 @@ -111,44 +98,42 @@ function schur_complement( H[indi, i] .= tmp2[indi] H[i, indi] .= tmp2[indi] else - if !iszero(SparseArrays.nnz(Ai)) - if SparseArrays.nnz(Ai) > 1 - @inbounds for jj in ii:n - j = σ[jj, ilmi] - Aj = model.A[ilmi, j] - if !iszero(SparseArrays.nnz(Aj)) - ttt = _dot(Ai, Aj, W) - if i >= j - H[i, j] = ttt - else - H[j, i] = ttt - end + if SparseArrays.nnz(Ai) > 1 + @inbounds for jj in ii:n + j = σ[jj, ilmi] + Aj = model.A[ilmi, j] + if !iszero(SparseArrays.nnz(Aj)) + ttt = _dot(Ai, Aj, W) + if i >= j + H[i, j] = ttt + else + H[j, i] = ttt end end - else - # A is symmetric - iiiiAi = jjjiAi = only(SparseArrays.rowvals(Ai)) - vvvi = only(SparseArrays.nonzeros(Ai)) - @inbounds for jj in ii:n - j = σ[jj, ilmi] - Ajjj = model.A[ilmi, j] - # As we sort the matrices in decreasing `nnz` order, - # the rest of matrices is either zero or have only - # one entry - if !iszero(SparseArrays.nnz(Ajjj)) - iiijAj = - jjjjAj = only(SparseArrays.rowvals(Ajjj)) - vvvj = only(SparseArrays.nonzeros(Ajjj)) - ttt = - vvvi * - W[iiiiAi, iiijAj] * - W[jjjiAi, jjjjAj] * - vvvj - if i >= j - H[i, j] = ttt - else - H[j, i] = ttt - end + end + elseif SparseArrays.nnz(Ai) == 1 + # A is symmetric + iiiiAi = jjjiAi = only(SparseArrays.rowvals(Ai)) + vvvi = only(SparseArrays.nonzeros(Ai)) + @inbounds for jj in ii:n + j = σ[jj, ilmi] + Ajjj = model.A[ilmi, j] + # As we sort the matrices in decreasing `nnz` order, + # the rest of matrices is either zero or have only + # one entry + if !iszero(SparseArrays.nnz(Ajjj)) + iiijAj = + jjjjAj = only(SparseArrays.rowvals(Ajjj)) + vvvj = only(SparseArrays.nonzeros(Ajjj)) + ttt = + vvvi * + W[iiiiAi, iiijAj] * + W[jjjiAi, jjjjAj] * + vvvj + if i >= j + H[i, j] = ttt + else + H[j, i] = ttt end end end @@ -159,6 +144,10 @@ function schur_complement( return H end +function schur_complement(model::Model, w, ::Type{ScalarIndex}) + return model.C_lin * SparseArrays.spdiagm(w) * model.C_lin' +end + # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA_jW⟩ @@ -177,10 +166,6 @@ function schur_complement(buffer, model::Model, W::AbstractVector) return LinearAlgebra.Hermitian(H, :L) end -function schur_complement(model::Model, w, ::Type{ScalarIndex}) - return model.C_lin * SparseArrays.spdiagm(w) * model.C_lin' -end - # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA(y)W⟩ diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index e490ab8..b4bc9a9 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -238,6 +238,21 @@ function SolverCore.solve!(::ConvexSolver, ::LRO.Model) return end +function schur_test(model, w, κ) + schur_buffer = LRO.buffer_for_schur_complement(model, κ) + jtprod_buffer = LRO.buffer_for_jtprod(model) + y = rand(model.meta.ncon) + H = LRO.schur_complement(schur_buffer, model, w) + Hy = similar(y) + LRO.eval_schur_complement!(jtprod_buffer, Hy, model, w, y) + @test Hy ≈ H * y +end + +function schur_test(model, κ) + w = rand(model.meta.nvar) + schur_test(model, LRO.VectorizedSolution(w, model.dim), κ) +end + @testset "ConvexSolver" begin model = maxcut(weights, LRO.Optimizer) set_attribute(model, "solver", ConvexSolver) @@ -263,4 +278,7 @@ end @test err[4] ≈ 70/9 @test err[5] ≈ 0.92 @test err[6] ≈ 392.0 + schur_test(b.model, 0) + #schur_test(b.model, 1) + #schur_test(b.model, 2) end From e2b217244cd17f6c2a356528873545bd57153520 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 22:30:20 +0200 Subject: [PATCH 40/48] Add more tests --- src/schur.jl | 1 + test/BurerMonteiro.jl | 10 +++++++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index a7f96d4..e5628cd 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -69,6 +69,7 @@ function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) return H end +# /!\ W needs to be symmetric function schur_complement( buffer, model::Model, diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index b4bc9a9..a5363f0 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -250,7 +250,11 @@ end function schur_test(model, κ) w = rand(model.meta.nvar) - schur_test(model, LRO.VectorizedSolution(w, model.dim), κ) + W = LRO.VectorizedSolution(w, model.dim) + for i in LRO.matrix_indices(model) + W[i] .= W[i] .+ W[i]' + end + schur_test(model, W, κ) end @testset "ConvexSolver" begin @@ -279,6 +283,6 @@ end @test err[5] ≈ 0.92 @test err[6] ≈ 392.0 schur_test(b.model, 0) - #schur_test(b.model, 1) - #schur_test(b.model, 2) + schur_test(b.model, 1) + schur_test(b.model, 2) end From f61638ed57f04672742058c759da41d0eb91c703 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 5 Jun 2025 22:53:25 +0200 Subject: [PATCH 41/48] Fix format --- src/schur.jl | 3 +-- test/BurerMonteiro.jl | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index e5628cd..3811e22 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -123,8 +123,7 @@ function schur_complement( # the rest of matrices is either zero or have only # one entry if !iszero(SparseArrays.nnz(Ajjj)) - iiijAj = - jjjjAj = only(SparseArrays.rowvals(Ajjj)) + iiijAj = jjjjAj = only(SparseArrays.rowvals(Ajjj)) vvvj = only(SparseArrays.nonzeros(Ajjj)) ttt = vvvi * diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index a5363f0..0a0b465 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -254,7 +254,7 @@ function schur_test(model, κ) for i in LRO.matrix_indices(model) W[i] .= W[i] .+ W[i]' end - schur_test(model, W, κ) + return schur_test(model, W, κ) end @testset "ConvexSolver" begin From d680192c96ff93e4dbecd7cfcbd4dc61222a01d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 11:05:56 +0200 Subject: [PATCH 42/48] Inplace schur complement --- src/schur.jl | 27 +++++++++++---------------- test/BurerMonteiro.jl | 6 ++++-- 2 files changed, 15 insertions(+), 18 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index 3811e22..f72731c 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -60,26 +60,24 @@ function buffer_for_schur_complement(model::Model, κ) return σ, last_dense end -function schur_complement(buffer, model::Model, W, ::Type{MatrixIndex}) - n = num_constraints(model) - H = zeros(eltype(eltype(W)), n, n) +function add_schur_complement!(buffer, model::Model, W, ::Type{MatrixIndex}, H) for i in matrix_indices(model) - H += schur_complement(buffer, model, i, W[i]) + add_schur_complement!(buffer, model, i, W[i], H) end return H end # /!\ W needs to be symmetric -function schur_complement( +function add_schur_complement!( buffer, model::Model, mat_idx::MatrixIndex, W::AbstractMatrix{T}, + H, ) where {T} σ, last_dense = buffer ilmi = mat_idx.value n = model.meta.ncon - H = zeros(T, n, n) dim = side_dimension(model, mat_idx) @assert dim == LinearAlgebra.checksquare(W) tmp1 = Matrix{T}(undef, dim, dim) @@ -144,24 +142,21 @@ function schur_complement( return H end -function schur_complement(model::Model, w, ::Type{ScalarIndex}) - return model.C_lin * SparseArrays.spdiagm(w) * model.C_lin' +function add_schur_complement!(model::Model, w, ::Type{ScalarIndex}, H) + H .+= model.C_lin * SparseArrays.spdiagm(w) * model.C_lin' + return H end # [HKS24, (5b)] # Returns the matrix equal to the sum, for each equation, of # ⟨A_i, WA_jW⟩ -function schur_complement(buffer, model::Model, W::AbstractVector) - H = MA.Zero() +function schur_complement!(buffer, model::Model, W::AbstractVector, H) + fill!(H, zero(eltype(H))) if num_matrices(model) > 0 - H = MA.add!!(H, schur_complement(buffer, model, W, MatrixIndex)) + add_schur_complement!(buffer, model, W, MatrixIndex, H) end if num_scalars(model) > 0 - H = MA.add!!(H, schur_complement(model, W[ScalarIndex], ScalarIndex)) - end - if H isa MA.Zero - n = num_constraints(model) - H = zeros(eltype(W), n, n) + add_schur_complement!(model, W[ScalarIndex], ScalarIndex, H) end return LinearAlgebra.Hermitian(H, :L) end diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 0a0b465..e32cf98 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -241,8 +241,10 @@ end function schur_test(model, w, κ) schur_buffer = LRO.buffer_for_schur_complement(model, κ) jtprod_buffer = LRO.buffer_for_jtprod(model) - y = rand(model.meta.ncon) - H = LRO.schur_complement(schur_buffer, model, w) + n = model.meta.ncon + y = rand(n) + H = zeros(n, n) + H = LRO.schur_complement(schur_buffer, model, w, H) Hy = similar(y) LRO.eval_schur_complement!(jtprod_buffer, Hy, model, w, y) @test Hy ≈ H * y From 0781cf93833110e099d88a8ad6dc982671d97cf7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 12:46:33 +0200 Subject: [PATCH 43/48] Fix --- test/BurerMonteiro.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index e32cf98..017f21a 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -244,7 +244,7 @@ function schur_test(model, w, κ) n = model.meta.ncon y = rand(n) H = zeros(n, n) - H = LRO.schur_complement(schur_buffer, model, w, H) + H = LRO.schur_complement!(schur_buffer, model, w, H) Hy = similar(y) LRO.eval_schur_complement!(jtprod_buffer, Hy, model, w, y) @test Hy ≈ H * y From e56f9e75b970cadb8f59ed74de96a7818126659c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 13:45:27 +0200 Subject: [PATCH 44/48] Fixes --- src/schur.jl | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index f72731c..beb23ae 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -93,9 +93,10 @@ function add_schur_complement!( LinearAlgebra.mul!(tmp, tmp1, W) fill!(tmp2, zero(T)) add_jprod!(model, mat_idx, tmp, tmp2) - indi = σ[ii:end, ilmi] - H[indi, i] .= tmp2[indi] - H[i, indi] .= tmp2[indi] + H[i, i] += tmp2[i] + indi = σ[(ii+1):end, ilmi] + H[indi, i] .+= tmp2[indi] + H[i, indi] .+= tmp2[indi] else if SparseArrays.nnz(Ai) > 1 @inbounds for jj in ii:n @@ -103,10 +104,9 @@ function add_schur_complement!( Aj = model.A[ilmi, j] if !iszero(SparseArrays.nnz(Aj)) ttt = _dot(Ai, Aj, W) - if i >= j - H[i, j] = ttt - else - H[j, i] = ttt + H[i, j] += ttt + if i != j + H[j, i] += ttt end end end @@ -128,10 +128,9 @@ function add_schur_complement!( W[iiiiAi, iiijAj] * W[jjjiAi, jjjjAj] * vvvj - if i >= j - H[i, j] = ttt - else - H[j, i] = ttt + H[i, j] += ttt + if i != j + H[j, i] += ttt end end end @@ -158,7 +157,7 @@ function schur_complement!(buffer, model::Model, W::AbstractVector, H) if num_scalars(model) > 0 add_schur_complement!(model, W[ScalarIndex], ScalarIndex, H) end - return LinearAlgebra.Hermitian(H, :L) + return H end # [HKS24, (5b)] From 7ad4403445b4e30fbf89a827ab423fe6e646f982 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 17:55:09 +0200 Subject: [PATCH 45/48] Add tests --- src/BurerMonteiro.jl | 6 ++--- src/model.jl | 51 ++++++++++++++++-------------------- src/schur.jl | 2 +- test/BurerMonteiro.jl | 60 +++++++++++++++++++++++++++++++++++-------- 4 files changed, 75 insertions(+), 44 deletions(-) diff --git a/src/BurerMonteiro.jl b/src/BurerMonteiro.jl index d2122a3..0dc6b36 100644 --- a/src/BurerMonteiro.jl +++ b/src/BurerMonteiro.jl @@ -33,7 +33,7 @@ struct Model{T,AT} <: NLPModels.AbstractNLPModel{T,Vector{T}} function Model(model::LRO.Model{T,AT}, ranks) where {T,AT} dim = Dimensions(model, ranks) n = length(dim) - ncon = LRO.num_constraints(model) + ncon = model.meta.ncon return new{T,AT}( model, dim, @@ -160,7 +160,7 @@ function NLPModels.jtprod!( U = JtV[i].factor fill!(U, zero(eltype(U))) for j in eachindex(y) - A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) + A = NLPModels.jac(model.model, j, i) U .+= A * X[i].factor .* (2y[j]) end end @@ -185,7 +185,7 @@ function NLPModels.hprod!( Hvi .+= C * Vi Hvi .*= 2obj_weight for j in 1:model.meta.ncon - A = NLPModels.jac(model.model, LRO.ConstraintIndex(j), i) + A = NLPModels.jac(model.model, j, i) Hvi .-= A * Vi .* (2y[j]) end end diff --git a/src/model.jl b/src/model.jl index c728bdb..c8d849d 100644 --- a/src/model.jl +++ b/src/model.jl @@ -79,6 +79,18 @@ end function Base.size(s::ShapedSolution) return (length(s.scalars) + sum(length, s.matrices, init = 0),) end +num_matrices(s::ShapedSolution) = length(s.matrices) + +function LinearAlgebra.norm2(s::ShapedSolution{T}) where {T} + # `LinearAlgebra.generic_norm2` starts by computing the ∞ norm and do a rescaling, we don't do that here + return √(LinearAlgebra.dot(s, s)) +end + +function LinearAlgebra.dot(a::ShapedSolution{T}, b::ShapedSolution{T}) where {T} + return LinearAlgebra.dot(a.scalars, b.scalars) + sum(eachindex(a.matrices); init = zero(T)) do i + LinearAlgebra.dot(a.matrices[i], b.matrices[i]) + end +end Base.view(s::ShapedSolution, ::Type{ScalarIndex}) = s.scalars Base.view(s::ShapedSolution, i::MatrixIndex) = s.matrices[i.value] @@ -151,7 +163,7 @@ mutable struct Model{T,A<:AbstractMatrix{T}} <: end function NLPModels.unconstrained(model::Model) - return iszero(num_constraints(model)) + return iszero(model.meta.ncon) end # TODO the scalar actually have lower bounds and the SDP variables too @@ -160,13 +172,6 @@ NLPModels.has_bounds(::Model) = false num_scalars(model::Model) = length(model.d_lin) -function scalar_indices(model::Model) - return MOI.Utilities.LazyMap{ScalarIndex}( - ScalarIndex, - Base.OneTo(num_scalars(model)), - ) -end - num_matrices(model::Model) = length(model.C) function matrix_indices(model::Union{Model,AbstractSolution}) @@ -178,24 +183,13 @@ end side_dimension(model::Model, i::MatrixIndex) = model.msizes[i.value] -struct ConstraintIndex - value::Int64 -end -num_constraints(model::Model) = length(model.b) -function constraint_indices(model::Model) - return MOI.Utilities.LazyMap{ConstraintIndex}( - ConstraintIndex, - Base.OneTo(num_constraints(model)), - ) -end - # Should be only used with `norm` NLPModels.jac(model::Model, ::Type{ScalarIndex}) = model.C_lin -function NLPModels.jac(model::Model, i::ConstraintIndex, j::MatrixIndex) - return model.A[j.value, i.value] +function NLPModels.jac(model::Model, j::Integer, i::MatrixIndex) + return model.A[i.value, j] end -function NLPModels.jac(model::Model, i::ConstraintIndex, ::Type{ScalarIndex}) - return model.C_lin[i.value, :] +function NLPModels.jac(model::Model, j::Integer, ::Type{ScalarIndex}) + return model.C_lin[j, :] end function norm_jac(model::Model{T}, i::MatrixIndex) where {T} if isempty(model.A) @@ -229,6 +223,7 @@ function NLPModels.grad!(model::Model, _::AbstractVector, g::AbstractVector) copyto!(g[ScalarIndex], model.d_lin) for i in matrix_indices(model) copyto!(g[i], model.C[i.value]) + @show @__LINE__ end return g end @@ -247,12 +242,13 @@ function buffer_for_jtprod(model::Model) end function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) - if iszero(num_constraints(model)) + if iszero(model.meta.ncon) + @show @__LINE__ return end # FIXME: at some point, switch to dense return sum( - abs.(model.A[mat_idx.value, j]) for j in 1:num_constraints(model) + abs.(model.A[mat_idx.value, j]) for j in 1:model.meta.ncon ) end @@ -279,9 +275,6 @@ end _zero!(A::SparseArrays.SparseMatrixCSC) = fill!(SparseArrays.nonzeros(A), 0.0) function jtprod!(buffer, model::Model, mat_idx::MatrixIndex, y) - if iszero(num_constraints(model)) - return MA.Zero() - end _zero!(buffer) for j in eachindex(y) _add_mul!(buffer, model.A[mat_idx.value, j], y[j]) @@ -320,7 +313,7 @@ function add_jprod!( V::AbstractMatrix, Jv::AbstractVector, ) - for j in 1:num_constraints(model) + for j in 1:model.meta.ncon Jv[j] += LinearAlgebra.dot(model.A[i.value, j], V) end end diff --git a/src/schur.jl b/src/schur.jl index beb23ae..16303d4 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -43,7 +43,7 @@ function _dot( end function buffer_for_schur_complement(model::Model, κ) - n = num_constraints(model) + n = model.meta.ncon σ = zeros(Int64, n, num_matrices(model)) last_dense = zeros(Int64, num_matrices(model)) diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 017f21a..52e3a75 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -1,5 +1,6 @@ using Test using LinearAlgebra +using SparseArrays using JuMP import LowRankOpt as LRO using Dualization @@ -84,7 +85,7 @@ end @testset "Simple LP $opt" for opt in [LRO.Optimizer, dual_optimizer(LRO.Optimizer)] - model = Model(dual_optimizer(LRO.Optimizer)) + model = Model(opt) @variable(model, x) @constraint(model, con_ref, 1 - x in Nonnegatives()) @objective(model, Max, x) @@ -105,6 +106,10 @@ end @test dual(con_ref) ≈ 1 @test objective_value(model) ≈ 1 @test dual_objective_value(model) ≈ 1 + b = unsafe_backend(model) + if !(b isa DualOptimizer) + @test isnothing(LRO.buffer_for_jtprod(b.model)) + end raw_sol = MOI.get(model, LRO.RawStatus(:solution)) sol = MOI.get(model, LRO.Solution()) @test raw_sol isa Vector{Float64} @@ -113,7 +118,9 @@ end @test length(outer) == length(sol) @test sprint(show, outer) == "_OuterProduct($sol, $sol)" @test sol == raw_sol - @test abs(sol[1]) < 1e-6 + if b isa DualOptimizer + @test abs(sol[1]) < 1e-6 + end diff_check(model) end; @@ -224,6 +231,19 @@ end MOI.Test.runtests(model, config; include = ["Silent"]) end; +@testset "No constraints" begin + model = LRO.Model( + [spzeros(1, 1)], + [ones(1, 1) for _ in 1:1, _ in 1:0], + zeros(0), + sparsevec(Int[], Float64[], 0), + sparse(Int[], Int[], Float64[], 0, 0), + [1], + ) + @test model.meta.ncon == 0 + @test LRO.norm_jac(model, LRO.MatrixIndex(1)) == 0 +end; + struct ConvexSolver{T} <: SolverCore.AbstractOptimizationSolver model::LRO.Model{T} stats::SolverCore.GenericExecutionStats{T,Vector{T},Vector{T},Any} @@ -248,6 +268,11 @@ function schur_test(model, w, κ) Hy = similar(y) LRO.eval_schur_complement!(jtprod_buffer, Hy, model, w, y) @test Hy ≈ H * y + for i in LRO.matrix_indices(model) + ret = LRO.dual_cons!(jtprod_buffer, model, i, y) + @test ret isa SparseMatrixCSC + end + @test LRO.dual_cons(model, LRO.ScalarIndex, y) isa SparseArrays.SparseVector end function schur_test(model, κ) @@ -274,16 +299,29 @@ end b.solver.stats.status = :unbounded @test termination_status(model) == MOI.INFEASIBLE @test primal_status(model) == MOI.INFEASIBLE_POINT - X = LRO.VectorizedSolution(collect(1:b.model.meta.nvar), b.model.dim) + x = LRO.VectorizedSolution(collect(1:b.model.meta.nvar), b.model.dim) + sim = similar(x) + @test sim isa typeof(x) + sim .= x + @test sim == x + X = LRO.ShapedSolution( + Vector(x[LRO.ScalarIndex]), + [Matrix(x[i]) for i in LRO.matrix_indices(b.model)], + ) y = collect(1:b.model.meta.ncon) - err = LRO.errors(b.solver.model, X; y, dual_slack = X, dual_err = X) - @test length(err) == 6 - @test err[1] ≈ 4.155017729878046 - @test err[2] ≈ 0.318237296391563 - @test err[3] ≈ 70/9 - @test err[4] ≈ 70/9 - @test err[5] ≈ 0.92 - @test err[6] ≈ 392.0 + @test NLPModels.jac(b.model, 1, LRO.MatrixIndex(1)) == sparse([1], [1], [-1], 4, 4) + @test NLPModels.jac(b.model, 1, LRO.ScalarIndex) == sparsevec([1, 2], [-1, 1], 8) + @test LRO.norm_jac(b.model, LRO.MatrixIndex(1)) == 4 + for xx in [x, X] + err = LRO.errors(b.solver.model, xx; y, dual_slack = xx, dual_err = xx) + @test length(err) == 6 + @test err[1] ≈ 4.155017729878046 + @test err[2] ≈ 0.318237296391563 + @test err[3] ≈ 70/9 + @test err[4] ≈ 70/9 + @test err[5] ≈ 0.92 + @test err[6] ≈ 392.0 + end schur_test(b.model, 0) schur_test(b.model, 1) schur_test(b.model, 2) From c7dc384fafe7af0ddb90b437b3541825d259a840 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 21:32:15 +0200 Subject: [PATCH 46/48] Cover all --- src/model.jl | 2 -- test/BurerMonteiro.jl | 7 +++++++ 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/model.jl b/src/model.jl index c8d849d..7f01f30 100644 --- a/src/model.jl +++ b/src/model.jl @@ -223,7 +223,6 @@ function NLPModels.grad!(model::Model, _::AbstractVector, g::AbstractVector) copyto!(g[ScalarIndex], model.d_lin) for i in matrix_indices(model) copyto!(g[i], model.C[i.value]) - @show @__LINE__ end return g end @@ -243,7 +242,6 @@ end function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) if iszero(model.meta.ncon) - @show @__LINE__ return end # FIXME: at some point, switch to dense diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 52e3a75..4519da3 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -242,6 +242,7 @@ end; ) @test model.meta.ncon == 0 @test LRO.norm_jac(model, LRO.MatrixIndex(1)) == 0 + @test isnothing(LRO.buffer_for_jtprod(model, LRO.MatrixIndex(1))) end; struct ConvexSolver{T} <: SolverCore.AbstractOptimizationSolver @@ -312,6 +313,12 @@ end @test NLPModels.jac(b.model, 1, LRO.MatrixIndex(1)) == sparse([1], [1], [-1], 4, 4) @test NLPModels.jac(b.model, 1, LRO.ScalarIndex) == sparsevec([1, 2], [-1, 1], 8) @test LRO.norm_jac(b.model, LRO.MatrixIndex(1)) == 4 + grad = similar(x) + NLPModels.grad!(b.model, X, grad) + @test Vector(grad) == [ + Vector(NLPModels.grad(b.model, LRO.ScalarIndex)); + NLPModels.grad(b.model, LRO.MatrixIndex(1))[:] + ] for xx in [x, X] err = LRO.errors(b.solver.model, xx; y, dual_slack = xx, dual_err = xx) @test length(err) == 6 From f32596e61fef6adcffdc7e7528348647b072202c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 21:32:23 +0200 Subject: [PATCH 47/48] Fix format --- src/model.jl | 9 ++++----- test/BurerMonteiro.jl | 6 ++++-- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/model.jl b/src/model.jl index 7f01f30..64c709e 100644 --- a/src/model.jl +++ b/src/model.jl @@ -87,8 +87,9 @@ function LinearAlgebra.norm2(s::ShapedSolution{T}) where {T} end function LinearAlgebra.dot(a::ShapedSolution{T}, b::ShapedSolution{T}) where {T} - return LinearAlgebra.dot(a.scalars, b.scalars) + sum(eachindex(a.matrices); init = zero(T)) do i - LinearAlgebra.dot(a.matrices[i], b.matrices[i]) + return LinearAlgebra.dot(a.scalars, b.scalars) + + sum(eachindex(a.matrices); init = zero(T)) do i + return LinearAlgebra.dot(a.matrices[i], b.matrices[i]) end end @@ -245,9 +246,7 @@ function buffer_for_jtprod(model::Model, mat_idx::MatrixIndex) return end # FIXME: at some point, switch to dense - return sum( - abs.(model.A[mat_idx.value, j]) for j in 1:model.meta.ncon - ) + return sum(abs.(model.A[mat_idx.value, j]) for j in 1:model.meta.ncon) end # Computes `A .+= B * α` diff --git a/test/BurerMonteiro.jl b/test/BurerMonteiro.jl index 4519da3..844c3f7 100644 --- a/test/BurerMonteiro.jl +++ b/test/BurerMonteiro.jl @@ -310,8 +310,10 @@ end [Matrix(x[i]) for i in LRO.matrix_indices(b.model)], ) y = collect(1:b.model.meta.ncon) - @test NLPModels.jac(b.model, 1, LRO.MatrixIndex(1)) == sparse([1], [1], [-1], 4, 4) - @test NLPModels.jac(b.model, 1, LRO.ScalarIndex) == sparsevec([1, 2], [-1, 1], 8) + @test NLPModels.jac(b.model, 1, LRO.MatrixIndex(1)) == + sparse([1], [1], [-1], 4, 4) + @test NLPModels.jac(b.model, 1, LRO.ScalarIndex) == + sparsevec([1, 2], [-1, 1], 8) @test LRO.norm_jac(b.model, LRO.MatrixIndex(1)) == 4 grad = similar(x) NLPModels.grad!(b.model, X, grad) From 3f2cbf7820602e5a18f6d1caed60f2708cfafc4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 6 Jun 2025 21:38:57 +0200 Subject: [PATCH 48/48] Add SparseArrays --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 6bfa942..9af24e3 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -9,4 +9,5 @@ NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"