From 7e5eea6a700c22a6aaab077b3736fb1432c13c3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 19 Dec 2023 17:47:52 +0100 Subject: [PATCH] Add maxrank option with Barvinok-Pataki bound (#11) * Add maxrank option with Barvinok-Pataki bound * Fix format * Fix format * Clarify * Fix format * Use inverse_trimap * Fix format * Add bounds in runtests --- src/MOI_wrapper.jl | 1 + src/SDPLR.jl | 11 +++++- src/bounds.jl | 97 ++++++++++++++++++++++++++++++++++++++++++++++ test/bounds.jl | 78 +++++++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 5 files changed, 186 insertions(+), 2 deletions(-) create mode 100644 src/bounds.jl create mode 100644 test/bounds.jl diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 5a500bf..00152ae 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -263,6 +263,7 @@ function MOI.optimize!(model::Optimizer) CAinfo_type = reduce(vcat, vcat([model.Cinfo_type], model.Ainfo_type), init = Cchar[]) maxranks = default_maxranks( + model.params.maxrank, model.blktype, model.blksz, CAinfo_entptr, diff --git a/src/SDPLR.jl b/src/SDPLR.jl index b08285a..9a53d19 100644 --- a/src/SDPLR.jl +++ b/src/SDPLR.jl @@ -2,6 +2,8 @@ module SDPLR using SDPLR_jll +include("bounds.jl") + function solve_sdpa_file(file) return run(`$(SDPLR_jll.sdplr()) $file`) end @@ -22,6 +24,7 @@ Base.@kwdef mutable struct Parameters gaptol::Cdouble = 1.0e-3 checkbd::Cptrdiff_t = -1 typebd::Cptrdiff_t = 1 + maxrank::Function = default_maxrank end # See `macros.h` @@ -43,16 +46,19 @@ function default_R(blktype::Vector{Cchar}, blksz, maxranks) return rand(nr) - rand(nr) end -function default_maxranks(blktype, blksz, CAinfo_entptr, m) +function default_maxranks(maxrank, blktype, blksz, CAinfo_entptr, m) numblk = length(blktype) # See `getstorage` in `main.c` return map(eachindex(blktype)) do k if blktype[k] == Cchar('s') + # Because we do `1:m` and not `0:m`, we do not count the objective. + # `maxrank` will increment our count by `1` assuming that it is + # always part of the objective. cons = count(1:m) do i ind = datablockind(i, k, numblk) return CAinfo_entptr[ind+1] > CAinfo_entptr[ind] end - return Csize_t(min(isqrt(2cons) + 1, blksz[k])) + return Csize_t(maxrank(cons, blksz[k])) else @assert blktype[k] == Cchar('d') return Csize_t(1) @@ -85,6 +91,7 @@ function solve( CAinfo_type::Vector{Cchar}; params::Parameters = Parameters(), maxranks::Vector{Csize_t} = default_maxranks( + params.maxrank, blktype, blksz, CAinfo_entptr, diff --git a/src/bounds.jl b/src/bounds.jl new file mode 100644 index 0000000..c7aea11 --- /dev/null +++ b/src/bounds.jl @@ -0,0 +1,97 @@ +""" + pataki(m, n = 0) + +Return the upper bound of [B02; Proposition 13.1], [P12; Corollary 3.3.1(1)] +for the minimal rank of an optimal solution for an SDP of `m` constraints. + +!!! note + The argument `n` is ignore and is there for the sole purpose of allowing + the function to be used as `maxrank` optimizer attribute. + +It is the maximal `r` such that `div(r * (r + 1), 2) ≤ m`. + +!!! note + This bound can trivially be improved by computing `min(pataki(m, n), n)`. + This `min` with `n` is not done so that these two bounds can easily be + manipulated independently. This allows doing `min(pataki(m, n) + 1, n)`. + +[B02] Barvinok, "A Course in Convexity", 2002. +[P12] Pataki, "The geometry of semidefinite programming", 2012. +""" +function pataki(m, n = 0) + # Let `τ(r) = r * (r + 1) / 2` + # `inverse_trimap(m)[2]` is the smallest `r` such that + # `τ(r) ≥ m` while we want the largest `r` such that `τ(r) ≤ m`. + # `inverse_trimap(m + 1)[2]` is smallest `r` such that + # `τ(r) ≥ m + 1` hence `τ(r) > m` which is the negation of `τ(r) ≤ m`. + # Therefore, we can do: + return MOI.Utilities.inverse_trimap(m + 1)[2] - 1 +end + +""" + barvinok(m, n) + +Return the improved upper bound of [B02; Proposition 13.4] +for the minimal rank of an optimal solution for an SDP of `m` constraints. + +[B02; Proposition 13.4] states that if `div(r * (r + 1), 2) ≤ m` and `1 < r < n` +then there exists a solution of rank `≤ r - 1`. +So for this case, `barvinok(m, n)` returns `r - 1` while `pataki(m, n)` +returns `r`. + +!!! note + This bound can trivially be improved by computing `min(barvinok(m, n), n)`. + This `min` with `n` is not done so that these two bounds can easily be + manipulated independently. This allows doing `min(barvinok(m, n) + 1, n)`. + +[B02] Barvinok, "A Course in Convexity", 2002. +""" +function barvinok(m, n) + r = pataki(m, n) + if m == MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(r)) && 1 < r < n + r -= 1 + end + return r +end + +""" + default_maxrank(m, n) + +Return the value of `min(pataki(m + 1) + 1, n)` following the results of [BM05] +suggesting to use a `r > r_{m + 1}` where `r_{m}` is `pataki(m)`. So it means +`r > pataki(m + 1)` or equivalently `r ≥ pataki(m + 1) + 1`. + +> "It is interesting to note that the implementation of Burer and Monteiro [BM03] +> required only `r ≥ r_m`, but now with the insight provided by Theorem 3.4, our +> implementation in Section 6 implements r ≥ r_{m+1}." [BM05; p. 433] + +Although Theorem 3.4 requires `≥` [BM05; Theorem 5.3], the perturbed Augmented +Lagrangian implemented in SDPLR requires the strict version `>` [BM05; Theorem 5.3]. + +!!! note + For `m > 1` `min(barvinok(m, n) + 1, n)` is the same as `min(pataki(m - 1, n) + 1, n)`. + Without the `min(⋅ + 1, n)`, these can be different. For instance, with + `m = 3, n = 2`, `barvinok(3, 2) = 2` and `pataki(2, 2) = 1` but after passing + them through `min(⋅ + 1, n)`, they are both `2`. The constraint `m > 1` is also + important. For instance, with `m = 1` and `n = 2`, `SDPLR.barvinok(1, 2)` is `1` while + `SDPLR.pataki(0, 2)` is `0`. + +!!! note + When calling the SDPLR executable, it defaults to choosing + `min(isqrt(2m) + 1, n)` (see line 307 in the `getstorage` function in + `SDPLR-1.03-beta/source/main.c`). The value of `isqrt(2m)` is equal to + `MOI.Utilities.side_dimension_for_vectorized_dimension(m)`. The value + is equal to `pataki(m + 1)` for all `m` between `1` and `17` except for `8` and `13`. + However, starting from `18`, they are equal only around 50% of the times + but the error is at most one since `isqrt(2m) - 1 ≤ pataki(m + 1) ≤ isqrt(2m)` + always holds. The value `isqrt(2m)` is therefore always larger so it goes in + the direction of a more benign landscape which is a safe choice. + +[BM03] Burer, Samuel, and Monteiro, Renato DC. +"A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization." +Mathematical programming 95.2 (2003): 329-357. +[BM05] Burer, Samuel, and Monteiro, Renato DC. +"Local minima and convergence in low-rank semidefinite programming." +Mathematical programming 103.3 (2005): 427-444. +""" +default_maxrank(m, n) = min(isqrt(2m) + 1, n) diff --git a/test/bounds.jl b/test/bounds.jl new file mode 100644 index 0000000..48c6b07 --- /dev/null +++ b/test/bounds.jl @@ -0,0 +1,78 @@ +module TestBounds + +using Test +import SDPLR +import MathOptInterface as MOI + +const PATAKI = [ + 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 + 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 + 3 3 3 3 3 3 3 + 3 3 3 3 3 3 3 +] +const BARVINOK = [ + 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 + 2 2 1 1 1 1 1 + 2 2 2 2 2 2 2 + 2 2 2 2 2 2 2 + 3 3 3 2 2 2 2 + 3 3 3 3 3 3 3 +] +const DEFAULT_MAXRANK = [ + 1 2 2 2 2 2 2 + 1 2 3 3 3 3 3 + 1 2 3 3 3 3 3 + 1 2 3 3 3 3 3 + 1 2 3 4 4 4 4 + 1 2 3 4 4 4 4 + 1 2 3 4 4 4 4 +] + +τ(r) = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(r)) + +function test_bounds() + for m in 1:7 + for n in 1:7 + @test SDPLR.pataki(m, n) == PATAKI[m, n] + @test SDPLR.barvinok(m, n) == BARVINOK[m, n] + @test SDPLR.default_maxrank(m, n) == DEFAULT_MAXRANK[m, n] + end + end + for m in 1:100 + @test isqrt(2m) == + MOI.Utilities.side_dimension_for_vectorized_dimension(m) + @test isqrt(2m) - 1 <= SDPLR.pataki(m) + @test isqrt(2m) - 1 <= SDPLR.pataki(m + 1) + @test SDPLR.pataki(m) <= isqrt(2m) + @test SDPLR.pataki(m + 1) <= isqrt(2m) + r = SDPLR.pataki(m) + @test τ(r) ≤ m + @test τ(r + 1) > m + end + for m in 1:10 + for n in 1:10 + @test min(SDPLR.pataki(m, n) + 1, n) == + min(SDPLR.barvinok(m + 1, n) + 1, n) + end + end + return +end + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +end # module + +TestBounds.runtests() diff --git a/test/runtests.jl b/test/runtests.jl index 7b43729..be306d1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ include("test_vibra.jl") include("test_simple.jl") +include("bounds.jl") include("MOI_wrapper.jl")