Merge pull request #74 from chriscoey/linsyscaches

use abstraction layer for linear system solvers

.gitignore

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

Manifest.toml

@@ -15,9 +15,9 @@ version = "0.8.10"
 
 [[BinaryProvider]]
 deps = ["Libdl", "Pkg", "SHA", "Test"]
-git-tree-sha1 = "48c147e63431adbcee69bc40b04c3f0fec0a4982"
+git-tree-sha1 = "38be61810f280c3c478f5c38aaf387f8f9199275"
 uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.0"
+version = "0.5.1"
 
 [[Combinatorics]]
 deps = ["LinearAlgebra", "Polynomials", "Test"]
@@ -83,6 +83,12 @@ version = "0.9.0"
 deps = ["LinearAlgebra", "Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[IterativeSolvers]]
+deps = ["LinearAlgebra", "Pkg", "Printf", "Random", "RecipesBase", "SparseArrays", "Test"]
+git-tree-sha1 = "a7b8056ecd43aa0d112d4b5d04180194d12c3c62"
+uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
+version = "0.7.1"
+
 [[JSON]]
 deps = ["Dates", "Distributed", "Mmap", "Pkg", "Sockets", "Test", "Unicode"]
 git-tree-sha1 = "fec8e4d433072731466d37ed0061b3ba7f70eeb9"
@@ -108,7 +114,7 @@ uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
 [[MathOptInterface]]
 deps = ["Compat", "Unicode"]
-git-tree-sha1 = "ba12e7ce825c1458c03f88aae84fa630d882d303"
+git-tree-sha1 = "14bbbdbc746e70bbae6cc96d1fa58f14e02fae4d"
 repo-rev = "master"
 repo-url = "https://github.com/JuliaOpt/MathOptInterface.jl.git"
 uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
@@ -145,6 +151,12 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[RecipesBase]]
+deps = ["Pkg", "Random", "Test"]
+git-tree-sha1 = "0b3cb370ee4dc00f47f1193101600949f3dcf884"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "0.6.0"
+
 [[Reexport]]
 deps = ["Pkg"]
 git-tree-sha1 = "7b1d07f411bc8ddb7977ec7f377b97b158514fe0"

Project.toml

@@ -9,6 +9,7 @@ DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

examples/envelope/envelope.jl

@@ -9,15 +9,23 @@ available at https://arxiv.org/abs/1712.01792
 =#
 
 using Hypatia
-using SparseArrays
 using LinearAlgebra
+using SparseArrays
 using DelimitedFiles
 using Random
 
-function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int, d::Int; use_data::Bool=false, dense::Bool=false, rseed::Int=1)
-    @assert deg <= d
+function build_envelope!(
+    npoly::Int,
+    deg::Int,
+    n::Int,
+    d::Int;
+    use_data::Bool = false,
+    dense::Bool = false,
+    rseed::Int = 1,
+    )
 
     # generate interpolation
+    @assert deg <= d
     (L, U, pts, P0, P, w) = Hypatia.interpolate(n, d, calc_w=true)
     LWts = fill(binomial(n+d-1, n), n)
     wtVals = 1.0 .- pts.^2
@@ -27,7 +35,7 @@ function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int,
     if dense
         A = repeat(Array(1.0I, U, U), outer=(1, npoly))
     else
-        A = repeat(sparse(1.0I, U, U), outer=(1, npoly)) # TODO maybe construct without repeat
+        A = repeat(sparse(1.0I, U, U), outer=(1, npoly))
     end
     G = Diagonal(-1.0I, npoly*U) # TODO uniformscaling
     b = w
@@ -44,17 +52,42 @@ function build_envelope!(alf::Hypatia.HypatiaOpt, npoly::Int, deg::Int, n::Int,
 
     cone = Hypatia.Cone([Hypatia.DualSumOfSquaresCone(U, [P, PWts...]) for k in 1:npoly], [1+(k-1)*U:k*U for k in 1:npoly])
 
-    return Hypatia.load_data!(alf, c, A, b, G, h, cone)
+    return (c, A, b, G, h, cone)
 end
 
-# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=true)
+function run_envelope()
+    # optionally use fixed data in folder
+    # select number of polynomials and degrees for the envelope
+    # select dimension and SOS degree (to be squared)
+    (c, A, b, G, h, cone) =
+        # build_envelope!(2, 5, 1, 5, use_data=true)
+        # build_envelope!(2, 5, 2, 8)
+        # build_envelope!(3, 5, 3, 5)
+        build_envelope!(2, 3, 1, 4, dense=false)
 
-# optionally use fixed data in folder
-# select number of polynomials and degrees for the envelope
-# select dimension and SOS degree (to be squared)
-# build_envelope!(alf, 2, 5, 1, 5, use_data=true)
-# build_envelope!(alf, 2, 5, 2, 8)
-# build_envelope!(alf, 3, 5, 3, 5)
-# build_envelope!(alf, 2, 3, 3, 5, dense=false)
+    Hypatia.check_data(c, A, b, G, h, cone)
+    (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
 
-# @time Hypatia.solve!(alf)
+    opt = Hypatia.Optimizer(maxiter=100, verbose=false)
+    Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
+    Hypatia.solve!(opt)
+
+    x = zeros(length(c))
+    x[dukeep] = Hypatia.get_x(opt)
+    y = zeros(length(b))
+    y[prkeep] = Hypatia.get_y(opt)
+    s = Hypatia.get_s(opt)
+    z = Hypatia.get_z(opt)
+
+    status = Hypatia.get_status(opt)
+    solvetime = Hypatia.get_solvetime(opt)
+    pobj = Hypatia.get_pobj(opt)
+    dobj = Hypatia.get_dobj(opt)
+
+    # @show status
+    # @show x
+    # @show pobj
+    # @show dobj
+    return nothing
+end

examples/lp/lp.jl

@@ -11,7 +11,16 @@ using SparseArrays
 using DelimitedFiles
 using Random
 
-function build_lp!(alf::Hypatia.HypatiaOpt, m::Int, n::Int; use_data::Bool=false, dense::Bool=false, nzfrac::Float64=1/sqrt(n), tosparse::Bool=false, rseed::Int=1)
+function build_lp!(
+    m::Int,
+    n::Int;
+    use_data::Bool = false,
+    dense::Bool = false,
+    nzfrac::Float64 = 1/sqrt(n),
+    tosparse::Bool = false,
+    rseed::Int = 1,
+    )
+
     # set up problem data
     if use_data
         # use provided data in data folder
@@ -35,17 +44,43 @@ function build_lp!(alf::Hypatia.HypatiaOpt, m::Int, n::Int; use_data::Bool=false
     end
     G = Diagonal(-1.0I, n) # TODO uniformscaling
     h = zeros(n)
+
     cone = Hypatia.Cone([Hypatia.NonnegativeCone(n)], [1:n])
 
-    return Hypatia.load_data!(alf, c, A, b, G, h, cone)
+    return (c, A, b, G, h, cone)
 end
 
-# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=true)
+function run_lp()
+    # optionally use fixed data in folder
+    # select the random matrix size, dense/sparse, sparsity fraction
+    (c, A, b, G, h, cone) =
+        # build_lp!(500, 1000, use_data=true)
+        # build_lp!(500, 1000)
+        build_lp!(15, 20)
+
+    Hypatia.check_data(c, A, b, G, h, cone)
+    (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
 
-# optionally use fixed data in folder
-# select the random matrix size, dense/sparse, sparsity fraction
-# build_lp!(alf, 500, 1000)
-# build_lp!(alf, 500, 1000, use_data=true)
+    opt = Hypatia.Optimizer(maxiter=100, verbose=false)
+    Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
+    Hypatia.solve!(opt)
 
-# solve it
-# @time Hypatia.solve!(alf)
+    x = zeros(length(c))
+    x[dukeep] = Hypatia.get_x(opt)
+    y = zeros(length(b))
+    y[prkeep] = Hypatia.get_y(opt)
+    s = Hypatia.get_s(opt)
+    z = Hypatia.get_z(opt)
+
+    status = Hypatia.get_status(opt)
+    solvetime = Hypatia.get_solvetime(opt)
+    pobj = Hypatia.get_pobj(opt)
+    dobj = Hypatia.get_dobj(opt)
+
+    # @show status
+    # @show x
+    # @show pobj
+    # @show dobj
+    return nothing
+end

examples/namedpoly/namedpoly.jl

@@ -12,7 +12,7 @@ using Hypatia
 using LinearAlgebra
 
 # list of currently available named polynomials, see https://people.sc.fsu.edu/~jburkardt/py_src/polynomials/polynomials.html
-const polys = Dict{Symbol,NamedTuple}(
+polys = Dict{Symbol,NamedTuple}(
     :butcher => (n=6, lbs=[-1,-0.1,-0.1,-1,-0.1,-0.1], ubs=[0,0.9,0.5,-0.1,-0.05,-0.03], deg=3,
         fn=((u,v,w,x,y,z) -> z*v^2+y*w^2-u*x^2+x^3+x^2-(1/3)*u+(4/3)*x)
         ),
@@ -48,7 +48,11 @@ const polys = Dict{Symbol,NamedTuple}(
         ),
 )
 
-function build_namedpoly!(alf::Hypatia.HypatiaOpt, polyname::Symbol, d::Int)
+function build_namedpoly!(
+    polyname::Symbol,
+    d::Int,
+    )
+
     # get data for named polynomial
     (n, lbs, ubs, deg, fn) = polys[polyname]
     if d < ceil(Int, deg/2)
@@ -63,34 +67,58 @@ function build_namedpoly!(alf::Hypatia.HypatiaOpt, polyname::Symbol, d::Int)
     pts .+= (ubs + lbs)'/2
     wtVals = (pts .- lbs') .* (ubs' .- pts)
     LWts = fill(binomial(n+d-1, n), n)
-    PWts = [Diagonal(sqrt.(wtVals[:,j]))*P0[:,1:LWts[j]] for j in 1:n]
+    PWts = [Diagonal(sqrt.(wtVals[:, j])) * P0[:, 1:LWts[j]] for j in 1:n]
 
     # set up problem data
     A = ones(1, U)
     b = [1.0,]
-    c = [fn(pts[j,:]...) for j in 1:U]
+    c = [fn(pts[j, :]...) for j in 1:U]
     G = Diagonal(-1.0I, U) # TODO uniformscaling?
     h = zeros(U)
 
     cone = Hypatia.Cone([Hypatia.DualSumOfSquaresCone(U, [P0, PWts...])], [1:U])
 
-    return Hypatia.load_data!(alf, c, A, b, G, h, cone)
+    return (c, A, b, G, h, cone)
 end
 
-# alf = Hypatia.HypatiaOpt(maxiter=100, verbose=false)
-
-# select the named polynomial to minimize and the SOS degree (to be squared)
-# build_namedpoly!(alf, :butcher, 2)
-# build_namedpoly!(alf, :caprasse, 4)
-# build_namedpoly!(alf, :goldsteinprice, 7)
-# build_namedpoly!(alf, :heart, 2)
-# build_namedpoly!(alf, :lotkavolterra, 3)
-# build_namedpoly!(alf, :magnetism7, 2)
-# build_namedpoly!(alf, :motzkin, 7)
-# build_namedpoly!(alf, :reactiondiffusion, 4)
-# build_namedpoly!(alf, :robinson, 8)
-# build_namedpoly!(alf, :rosenbrock, 4)
-# build_namedpoly!(alf, :schwefel, 3)
-
-# solve it
-# @time Hypatia.solve!(alf)
+function run_namedpoly()
+    # select the named polynomial to minimize and the SOS degree (to be squared)
+    (c, A, b, G, h, cone) =
+        # build_namedpoly!(:butcher, 2)
+        # build_namedpoly!(:caprasse, 4)
+        # build_namedpoly!(:goldsteinprice, 7)
+        # build_namedpoly!(:heart, 2)
+        # build_namedpoly!(:lotkavolterra, 3)
+        # build_namedpoly!(:magnetism7, 2)
+        # build_namedpoly!(:motzkin, 7)
+        build_namedpoly!(:reactiondiffusion, 3)
+        # build_namedpoly!(:robinson, 8)
+        # build_namedpoly!(:rosenbrock, 4)
+        # build_namedpoly!(:schwefel, 3)
+
+    Hypatia.check_data(c, A, b, G, h, cone)
+    (c1, A1, b1, G1, prkeep, dukeep, Q2, RiQ1) = Hypatia.preprocess_data(c, A, b, G, useQR=true)
+    L = Hypatia.QRSymmCache(c1, A1, b1, G1, h, cone, Q2, RiQ1)
+
+    opt = Hypatia.Optimizer(maxiter=100, verbose=false)
+    Hypatia.load_data!(opt, c1, A1, b1, G1, h, cone, L)
+    Hypatia.solve!(opt)
+
+    x = zeros(length(c))
+    x[dukeep] = Hypatia.get_x(opt)
+    y = zeros(length(b))
+    y[prkeep] = Hypatia.get_y(opt)
+    s = Hypatia.get_s(opt)
+    z = Hypatia.get_z(opt)
+
+    status = Hypatia.get_status(opt)
+    solvetime = Hypatia.get_solvetime(opt)
+    pobj = Hypatia.get_pobj(opt)
+    dobj = Hypatia.get_dobj(opt)
+
+    # @show status
+    # @show x
+    # @show pobj
+    # @show dobj
+    return nothing
+end

src/Hypatia.jl

@@ -4,14 +4,17 @@ Copyright 2018, Chris Coey and contributors
 
 module Hypatia
     using Printf
-    using SparseArrays
     using LinearAlgebra
-    using ForwardDiff
-    using DiffResults
+    using SparseArrays
 
+    import FFTW
+    import Combinatorics
     include("interpolation.jl")
+
+    import ForwardDiff
+    import DiffResults
     include("cone.jl")
-    for primcone in [
+    for primitivecone in [
         "orthant",
         "dualsumofsquares",
         "secondorder",
@@ -21,11 +24,21 @@ module Hypatia
         "positivesemidefinite",
         "ellinfinity",
         ]
-        include(joinpath(@__DIR__, "primitivecones", primcone * ".jl"))
+        include(joinpath(@__DIR__, "primitivecones", primitivecone * ".jl"))
     end
+
+    import IterativeSolvers
+    include("linearsystem.jl")
+    for linsyssolver in [
+        "qrsymm",
+        "naive",
+        ]
+        include(joinpath(@__DIR__, "linsyssolvers", linsyssolver * ".jl"))
+    end
+
     include("nativeinterface.jl")
 
     import MathOptInterface
-    MOI = MathOptInterface
+    const MOI = MathOptInterface
     include("mathoptinterface.jl")
 end

src/cone.jl

@@ -1,15 +1,18 @@
+#=
+Copyright 2018, Chris Coey and contributors
+=#
 
 # cone object
 abstract type PrimitiveCone end
 
 # TODO reorder primitive cones so easiest ones to check incone are first
 mutable struct Cone
     prms::Vector{PrimitiveCone}
-    idxs::Vector{AbstractVector{Int}}
+    idxs::Vector{UnitRange{Int}}
 end
-Cone() = Cone(PrimitiveCone[], AbstractVector{Int}[])
+Cone() = Cone(PrimitiveCone[], UnitRange{Int}[])
 
-function addprimitivecone!(cone::Cone, prm::PrimitiveCone, idx::AbstractVector{Int})
+function addprimitivecone!(cone::Cone, prm::PrimitiveCone, idx::UnitRange{Int})
     @assert dimension(prm) == length(idx)
     push!(cone.prms, prm)
     push!(cone.idxs, idx)