Merge pull request #62 from chriscoey/newlinearpdform

WIP: use new primal-dual format (linear objectives only)

.gitignore

@@ -2,3 +2,4 @@
 *.jl.*.cov
 *.jl.mem
 deps/deps.jl
+scratch.jl

Manifest.toml

@@ -9,9 +9,9 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
 [[BinaryProvider]]
 deps = ["Libdl", "Pkg", "SHA", "Test"]
-git-tree-sha1 = "b530fbeb6f41ab5a83fbe3db1fcbe879334bcd2d"
+git-tree-sha1 = "48c147e63431adbcee69bc40b04c3f0fec0a4982"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.4.2"
+version = "0.5.0"
 
 [[Combinatorics]]
 deps = ["LinearAlgebra", "Polynomials", "Test"]
@@ -21,15 +21,15 @@ version = "0.7.0"
 
 [[Compat]]
 deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "ae262fa91da6a74e8937add6b613f58cd56cdad4"
+git-tree-sha1 = "ff2595695fc4f14427358ce2593f867085c45dcb"
 uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "1.1.0"
+version = "1.2.0"
 
 [[Conda]]
 deps = ["Compat", "JSON", "VersionParsing"]
-git-tree-sha1 = "a47f9a2c7b80095e6a935536795635522fe27f5d"
+git-tree-sha1 = "85b5bf3ffcf4f39abe019dab1dd00a0aead8d882"
 uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
-version = "1.0.1"
+version = "1.0.2"
 
 [[Dates]]
 deps = ["Printf"]
@@ -78,9 +78,9 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[MathOptInterface]]
 deps = ["Compat", "Unicode"]
-git-tree-sha1 = "f2791f990529b7aea1fcb9fd45b636cd3a9d1e79"
+git-tree-sha1 = "ba12e7ce825c1458c03f88aae84fa630d882d303"
 uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "0.6.0"
+version = "0.6.1"
 
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
@@ -91,9 +91,9 @@ uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 
 [[Polynomials]]
 deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "af883752c4935425a3ab30031a2069254c451b8b"
+git-tree-sha1 = "1a1eae52956658a6acae6fa1b6d6c3d488192895"
 uuid = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
-version = "0.5.0"
+version = "0.5.1"
 
 [[Printf]]
 deps = ["Unicode"]
@@ -139,7 +139,7 @@ deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[UUIDs]]
-deps = ["Random"]
+deps = ["Random", "SHA"]
 uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [[Unicode]]
@@ -149,4 +149,4 @@ uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 deps = ["Compat"]
 git-tree-sha1 = "c9d5aa108588b978bd859554660c8a5c4f2f7669"
 uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
-version = "1.1.2"
+version = "1.1.3"

README.md

@@ -3,3 +3,29 @@
 Under construction. Only works on Julia v0.7+.
 
 An interior point solver for general convex conic optimization problems. An extension of methods in [CVXOPT](https://github.com/cvxopt/cvxopt/blob/master/src/python/coneprog.py) and [Alfonso](https://github.com/dpapp-github/alfonso). 
+
+Solves a pair of primal and dual cone programs:
+
+primal (over x,s):
+```
+  min  c'x :          duals
+    b - Ax == 0       (y)
+    h - Gx == s in K  (z)
+```
+dual (over z,y):
+```
+  max  -b'y - h'z :      duals
+    c + A'y + G'z == 0   (x)
+                z in K*  (s)
+```
+where K is a convex cone defined as a Cartesian product of recognized primitive cones, and K* is its dual cone.
+
+The primal-dual optimality conditions are:
+```
+         b - Ax == 0
+         h - Gx == s
+  c + A'y + G'z == 0
+            s'z == 0
+              s in K
+              z in K*
+```

examples/envelope/envelope.jl

@@ -21,28 +21,30 @@ function build_envelope!(alf::Alfonso.AlfonsoOpt, npoly::Int, deg::Int, n::Int,
     (L, U, pts, P0, P, w) = Alfonso.interpolate(n, d, calc_w=true)
     LWts = fill(binomial(n+d-1, n), n)
     wtVals = 1.0 .- pts.^2
-    PWts = [Array((qr(Diagonal(sqrt.(wtVals[:,j]))*P[:,1:LWts[j]])).Q) for j in 1:n]
+    PWts = [Array((qr(Diagonal(sqrt.(wtVals[:, j])) * P[:, 1:LWts[j]])).Q) for j in 1:n]
 
     # set up problem data
     if dense
-        A = repeat(Array(1.0I, U, U), outer=(1, npoly)) # dense A
+        A = repeat(Array(1.0I, U, U), outer=(1, npoly))
     else
-        A = repeat(sparse(1.0I, U, U), outer=(1, npoly)) # sparse A
+        A = repeat(sparse(1.0I, U, U), outer=(1, npoly)) # TODO maybe construct without repeat
     end
+    G = Diagonal(-1.0I, npoly*U) # TODO uniformscaling
     b = w
+    h = zeros(npoly*U)
     if use_data
         # use provided data in data folder
         c = vec(readdlm(joinpath(@__DIR__, "data/c$(size(A,2)).txt"), ',', Float64))
     else
         # generate random data
         Random.seed!(rseed)
         LDegs = binomial(n+deg, n)
-        c = vec(P0[:,1:LDegs]*rand(-9:9, LDegs, npoly))
+        c = vec(P0[:, 1:LDegs]*rand(-9:9, LDegs, npoly))
     end
 
     cone = Alfonso.Cone([Alfonso.SumOfSquaresCone(U, [P, PWts...]) for k in 1:npoly], [1+(k-1)*U:k*U for k in 1:npoly])
 
-    return Alfonso.load_data!(alf, A, b, c, cone)
+    return Alfonso.load_data!(alf, c, A, b, G, h, cone)
 end
 
 # alf = Alfonso.AlfonsoOpt(maxiter=100, verbose=true)

examples/lp/lp.jl

@@ -32,9 +32,11 @@ function build_lp!(alf::Alfonso.AlfonsoOpt, m::Int, n::Int; use_data::Bool=false
     if tosparse && !issparse(A)
         A = sparse(A)
     end
+    G = Diagonal(-1.0I, n) # TODO uniformscaling
+    h = zeros(n)
     cone = Alfonso.Cone([Alfonso.NonnegativeCone(n)], [1:n])
 
-    return Alfonso.load_data!(alf, A, b, c, cone)
+    return Alfonso.load_data!(alf, c, A, b, G, h, cone)
 end
 
 # alf = Alfonso.AlfonsoOpt(maxiter=100, verbose=true)

examples/namedpoly/namedpoly.jl

@@ -68,9 +68,12 @@ function build_namedpoly!(alf::Alfonso.AlfonsoOpt, polyname::Symbol, d::Int)
     A = ones(1, U)
     b = [1.0,]
     c = [fn(pts[j,:]...) for j in 1:U]
+    G = Diagonal(-1.0I, U) # TODO uniformscaling?
+    h = zeros(U)
+
     cone = Alfonso.Cone([Alfonso.SumOfSquaresCone(U, [P0, PWts...])], [1:U])
 
-    return Alfonso.load_data!(alf, A, b, c, cone)
+    return Alfonso.load_data!(alf, c, A, b, G, h, cone)
 end
 
 # alf = Alfonso.AlfonsoOpt(maxiter=100, verbose=false)

src/Alfonso.jl

@@ -6,7 +6,15 @@ module Alfonso
 
     include("interpolation.jl")
     include("cone.jl")
-    for primcone in ["nonnegative", "sumofsquares", "secondorder", "exponential", "power", "rotatedsecondorder", "positivesemidefinite"]
+    for primcone in [
+        "nonnegative",
+        "sumofsquares",
+        "secondorder",
+        "exponential",
+        # "power",
+        "rotatedsecondorder",
+        "positivesemidefinite",
+        ]
         include(joinpath(@__DIR__, "primitivecones", primcone * ".jl"))
     end
     include("nativeinterface.jl")

src/cone.jl

@@ -14,11 +14,11 @@ end
 # calculate complexity parameter of the barrier (sum of the primitive cone barrier parameters)
 barrierpar(cone::Cone) = sum(barrierpar_prm(prm) for prm in cone.prms)
 
-function getintdir!(arr::Vector{Float64}, cone::Cone)
+function getintdir!(dir::Vector{Float64}, cone::Cone)
     for k in eachindex(cone.prms)
-        getintdir_prm!(view(arr, cone.idxs[k]), cone.prms[k])
+        getintdir_prm!(view(dir, cone.idxs[k]), cone.prms[k])
     end
-    return arr
+    return dir
 end
 
 # TODO can parallelize the functions acting on Cone
@@ -38,18 +38,32 @@ function calcg!(g::Vector{Float64}, cone::Cone)
     return g
 end
 
-function calcHiarr!(Hi_mat::AbstractMatrix{Float64}, mat::AbstractMatrix{Float64}, cone::Cone)
+function calcHarr!(prod::AbstractMatrix{Float64}, arr::AbstractMatrix{Float64}, cone::Cone)
     for k in eachindex(cone.prms)
-        calcHiarr_prm!(view(Hi_mat, cone.idxs[k], :), view(mat, cone.idxs[k], :), cone.prms[k])
+        calcHarr_prm!(view(prod, cone.idxs[k], :), view(arr, cone.idxs[k], :), cone.prms[k])
     end
-    return Hi_mat
+    return prod
 end
 
-function calcHiarr!(Hi_vec::AbstractVector{Float64}, vec::AbstractVector{Float64}, cone::Cone)
+function calcHarr!(prod::AbstractVector{Float64}, arr::AbstractVector{Float64}, cone::Cone)
     for k in eachindex(cone.prms)
-        calcHiarr_prm!(view(Hi_vec, cone.idxs[k]), view(vec, cone.idxs[k]), cone.prms[k])
+        calcHarr_prm!(view(prod, cone.idxs[k]), view(arr, cone.idxs[k]), cone.prms[k])
     end
-    return Hi_vec
+    return prod
+end
+
+function calcHiarr!(prod::AbstractMatrix{Float64}, arr::AbstractMatrix{Float64}, cone::Cone)
+    for k in eachindex(cone.prms)
+        calcHiarr_prm!(view(prod, cone.idxs[k], :), view(arr, cone.idxs[k], :), cone.prms[k])
+    end
+    return prod
+end
+
+function calcHiarr!(prod::AbstractVector{Float64}, arr::AbstractVector{Float64}, cone::Cone)
+    for k in eachindex(cone.prms)
+        calcHiarr_prm!(view(prod, cone.idxs[k]), view(arr, cone.idxs[k]), cone.prms[k])
+    end
+    return prod
 end
 
 # utilities for converting between smat and svec forms (lower triangle) for symmetric matrices