Skip to content

Commit

Permalink
Add maxrank option with Barvinok-Pataki bound (#11)
Browse files Browse the repository at this point in the history
* Add maxrank option with Barvinok-Pataki bound

* Fix format

* Fix format

* Clarify

* Fix format

* Use inverse_trimap

* Fix format

* Add bounds in runtests
  • Loading branch information
blegat authored Dec 19, 2023
1 parent 572e5c4 commit 7e5eea6
Show file tree
Hide file tree
Showing 5 changed files with 186 additions and 2 deletions.
1 change: 1 addition & 0 deletions src/MOI_wrapper.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
11 changes: 9 additions & 2 deletions src/SDPLR.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,8 @@ module SDPLR

using SDPLR_jll

include("bounds.jl")

function solve_sdpa_file(file)
return run(`$(SDPLR_jll.sdplr()) $file`)
end
Expand All @@ -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`
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand Down
97 changes: 97 additions & 0 deletions src/bounds.jl
Original file line number Diff line number Diff line change
@@ -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)
78 changes: 78 additions & 0 deletions test/bounds.jl
Original file line number Diff line number Diff line change
@@ -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()
1 change: 1 addition & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
include("test_vibra.jl")
include("test_simple.jl")
include("bounds.jl")
include("MOI_wrapper.jl")

0 comments on commit 7e5eea6

Please sign in to comment.