Add maxrank option with Barvinok-Pataki bound

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,

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,

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)

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()

test/runtests.jl

@@ -1,3 +1,4 @@
 include("test_vibra.jl")
 include("test_simple.jl")
+include("bounds.jl")
 include("MOI_wrapper.jl")