Skip to content

Commit cc15af2

Browse files
committed
Add README.
1 parent 84249b2 commit cc15af2

File tree

7 files changed

+86
-16
lines changed

7 files changed

+86
-16
lines changed

.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
*.swp
2+
Manifest.toml

README.md

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
# MathOptChordalDecomposition.jl
2+
3+
MathOptChordalDecomposition.jl is a [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) layer that implements chordal decomposition of
4+
sparse semidefinite constraints.
5+
6+
## Basic Usage
7+
8+
The `sdplib` directory contains three semidefinite programming problems from the [SDPLIB library](https://github.com/vsdp/SDPLIB). The function `construct_model`, defined below, reads one of the problems and constructs a [JuMP.jl](https://github.com/jump-dev/JuMP.jl) model.
9+
10+
```julia-repl
11+
julia> using FileIO, LinearAlgebra, JuMP
12+
13+
julia> import MathOptChordalDecomposition, MosekTools, Mosek
14+
15+
julia> function construct_model(f, name::String)
16+
# load data
17+
data = load("./sdplib/$name.jld2");
18+
F = data["F"]
19+
c = data["c"]
20+
m = data["m"]
21+
n = data["n"]
22+
23+
# construct model
24+
model = JuMP.Model(f)
25+
set_silent(model)
26+
@variable(model, x[1:m])
27+
@objective(model, Min, c' * x)
28+
@constraint(model, con1, Symmetric(-Matrix(F[1]) + sum(Matrix(F[k + 1]) .* x[k] for k in 1:m)) in JuMP.PSDCone())
29+
return model
30+
end
31+
construct_model (generic function with 1 method)
32+
```
33+
34+
Solve the problem using the [Mosek.jl](https://github.com/MOSEK/Mosek.jl) optimizer.
35+
36+
```julia-repl
37+
julia> model = construct_model(Mosek.Optimizer, "mcp124-1");
38+
39+
julia> @time JuMP.optimize!(model)
40+
6.005076 seconds (2.53 k allocations: 515.859 KiB)
41+
42+
julia> MOI.get(model, MOI.ObjectiveValue())
43+
141.9904770422396
44+
```
45+
46+
Solve the problem using [Mosek.jl](https://github.com/MOSEK/Mosek.jl) and MathOptChordalDecomposition.jl.
47+
48+
```julia-repl
49+
julia> model = construct_model("mcp124-1") do
50+
MathOptChordalDecomposition.Optimizer(Mosek.Optimizer)
51+
end;
52+
53+
julia> @time JuMP.optimize!(model)
54+
0.041175 seconds (230.72 k allocations: 11.800 MiB)
55+
56+
julia> MOI.get(model.moi_backend.optimizer.model.inner, MOI.ObjectiveValue())
57+
141.99047611570586
58+
```

sdplib/mcp124-1.jld2

189 KB
Binary file not shown.

sdplib/mcp124-2.jld2

194 KB
Binary file not shown.

sdplib/mcp124-3.jld2

204 KB
Binary file not shown.

sdplib/mcp124-4.jld2

224 KB
Binary file not shown.

src/MathOptChordalDecomposition.jl

Lines changed: 26 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -298,23 +298,11 @@ function MOI.get(
298298
return # ???
299299
end
300300

301-
function idx(i::Int, j::Int)
302-
@assert i <= j
303-
x = i + j * (j - 1) ÷ 2
304-
return x
305-
end
306-
307-
function col(x::Int)
308-
j = ceil(Int, (sqrt(1 + 8x) - 1.0) / 2)
309-
return j
310-
end
311-
312-
function row(x::Int)
313-
j = col(x)
314-
i = x - j * (j - 1) ÷ 2
315-
return i
316-
end
301+
# Utilities
317302

303+
# `(V, A, b) = decode(f, S)` satisfies
304+
# f(Vᵢ) = Aᵢ + b
305+
# for all 1 ≤ i ≤ n.
318306
function decode(f::MOI.VectorAffineFunction{T}, S::MOI.PositiveSemidefiniteConeTriangle) where {T}
319307
n = S.side_dimension
320308
index = Dict{MOI.VariableIndex, Int}()
@@ -365,10 +353,32 @@ function decode(f::MOI.VectorAffineFunction{T}, S::MOI.PositiveSemidefiniteConeT
365353
return V, A, b
366354
end
367355

356+
# `S = sparsitypattern(A)` is a binary matrix with the same
357+
# sparsity pattern as A.
368358
function sparsitypattern(A::AbstractMatrix{T}) where {T}
369359
S = sparse(Symmetric(A, :U))
370360
nonzeros(S) .= one(T)
371361
return S
372362
end
373363

364+
# upper triangular indices
365+
# idx ∘ (row, col) = id
366+
# (row, col) ∘ idx = id
367+
function idx(i::Int, j::Int)
368+
@assert i <= j
369+
x = i + j * (j - 1) ÷ 2
370+
return x
371+
end
372+
373+
function col(x::Int)
374+
j = ceil(Int, (sqrt(1 + 8x) - 1.0) / 2)
375+
return j
376+
end
377+
378+
function row(x::Int)
379+
j = col(x)
380+
i = x - j * (j - 1) ÷ 2
381+
return i
382+
end
383+
374384
end

0 commit comments

Comments
 (0)