From c9f1eb5f1565ab5a4b99c63d27ae5e81f001d317 Mon Sep 17 00:00:00 2001 From: Chris Coey Date: Sun, 19 Aug 2018 21:15:38 -0400 Subject: [PATCH] clean up cone functions --- examples/lp/lp.jl | 2 +- src/conedata.jl | 144 +++++++++++++++++++---------- src/nativeinterface.jl | 203 ++++++++--------------------------------- test/nativeexamples.jl | 16 ++-- 4 files changed, 146 insertions(+), 219 deletions(-) diff --git a/examples/lp/lp.jl b/examples/lp/lp.jl index d3a4f9ecc..527b09b9c 100644 --- a/examples/lp/lp.jl +++ b/examples/lp/lp.jl @@ -40,7 +40,7 @@ function build_lp!(alf::Alfonso.AlfonsoOpt, m::Int, n::Int; use_data=false, dens return Alfonso.load_data!(alf, A, b, c, cones, coneidxs) end -# alf = Alfonso.AlfonsoOpt(maxiter=100, verbose=false) +# alf = Alfonso.AlfonsoOpt(maxiter=100, verbose=true) # optionally use fixed data in folder # select the random matrix size, dense/sparse, sparsity fraction diff --git a/src/conedata.jl b/src/conedata.jl index f3f87272e..bf3709388 100644 --- a/src/conedata.jl +++ b/src/conedata.jl @@ -1,23 +1,89 @@ # export ConeData, NonnegData, SumOfSqrData - abstract type ConeData end - #= Nonnegative cone =# +# +# mutable struct NonnegData <: ConeData +# dim::Int +# pnt::Vector{Float64} +# has_g::Bool +# g::Vector{Float64} +# g_prev::Vector{Float64} +# has_Hi::Bool +# Hi::Diagonal{Float64} +# Hi_prev::Diagonal{Float64} +# +# function NonnegData(dim) +# k = new() +# k.dim = dim +# k.has_g = false +# k.g = Vector{Float64}(undef, dim) +# k.g_prev = copy(k.g) +# k.has_Hi = false +# k.Hi = Diagonal{Float64}(copy(k.g)) +# k.Hi_prev = copy(k.Hi) +# return k +# end +# end +# +# dimension(k::NonnegData) = k.dim +# +# barpar(k::NonnegData) = k.dim +# +# function load_txk(k::NonnegData, pnt::Vector{Float64}, save_prev::Bool) +# @assert length(pnt) == k.dim +# if save_prev # may want to use previously calculated values, so store +# k.g_prev .= k.g +# k.Hi_prev.diag .= k.Hi.diag +# end +# k.pnt = pnt +# k.has_g = false +# k.has_Hi = false +# return nothing +# end +# +# function use_prevk(k::NonnegData) +# k.g .= k.g_prev +# k.Hi.diag .= k.Hi_prev.diag +# return nothing +# end +# +# inconek(k::NonnegData) = all(x -> (x > 0.0), k.pnt) +# +# function calc_gk(k::NonnegData) +# if !k.has_g +# k.g .= -inv.(k.pnt) +# k.has_g = true +# end +# return k.g +# end +# +# function calc_Hik(k::NonnegData) +# if !k.has_Hi +# k.Hi.diag .= abs2.(k.pnt) +# k.has_Hi = true +# end +# return k.Hi +# end +# +# calc_Lk(k::NonnegData) = Diagonal(-k.g) +# + mutable struct NonnegData <: ConeData dim::Int - - pnt - pnt_prev - # invpnt # TODO store this and the square + pnt::Vector{Float64} + pnt_prev::Vector{Float64} function NonnegData(dim) k = new() k.dim = dim + k.pnt = Vector{Float64}(undef, dim) + k.pnt_prev = copy(k.pnt) + # TODO prealloc etc return k end end @@ -28,56 +94,48 @@ barpar(k::NonnegData) = k.dim function load_txk(k::NonnegData, pnt, save_prev) @assert length(pnt) == k.dim - if save_prev - k.pnt_prev = copy(k.pnt) + k.pnt_prev .= k.pnt end - - k.pnt = pnt + k.pnt .= pnt end -use_prevk(k::NonnegData) = (k.pnt = copy(k.pnt_prev)) +use_prevk(k::NonnegData) = (k.pnt .= k.pnt_prev) inconek(k::NonnegData) = all(x -> (x > 0.0), k.pnt) calc_gk(k::NonnegData) = -inv.(k.pnt) -calc_Hinvk(k::NonnegData) = Diagonal(k.pnt.^2) - -calc_HCholLk(k::NonnegData) = Diagonal(inv.(k.pnt)) # TODO just -g +calc_Hik(k::NonnegData) = Diagonal(abs2.(k.pnt)) +calc_Lk(k::NonnegData) = Diagonal(inv.(k.pnt)) # TODO just -g #= Sum of squares cone =# mutable struct SumOfSqrData <: ConeData dim::Int - ip - ipwt - - isnewpnt - pnt - - grad + ip::Matrix{Float64} + ipwt::Vector{Matrix{Float64}} + pnt::Vector{Float64} Hfact - hasgH - hasHi - Hi - hasL - L - - grad_prev - Hi_prev - L_prev - - # TODO prealloc etc + hasgH::Bool + hasHi::Bool + hasL::Bool + grad::Vector{Float64} + grad_prev::Vector{Float64} + Hi::Matrix{Float64} + Hi_prev::Matrix{Float64} + L::Matrix{Float64} + L_prev::Matrix{Float64} function SumOfSqrData(dim, ip, ipwt) k = new() k.dim = dim k.ip = ip k.ipwt = ipwt - k.isnewpnt = false + k.pnt = Vector{Float64}(undef, dim) + # TODO prealloc etc return k end end @@ -86,19 +144,16 @@ dimension(k::SumOfSqrData) = k.dim barpar(k::SumOfSqrData) = (size(k.ip, 2) + sum(size(k.ipwt[j], 2) for j in 1:length(k.ipwt))) -function load_txk(k::SumOfSqrData, pnt, save_prev) +function load_txk(k::SumOfSqrData, pnt::Vector{Float64}, save_prev::Bool) @assert length(pnt) == k.dim - if save_prev # may want to use previously calculated values, so store @assert k.hasgH k.grad_prev = calc_gk(k) - k.Hi_prev = calc_Hinvk(k) - k.L_prev = calc_HCholLk(k) + k.Hi_prev = calc_Hik(k) + k.L_prev = calc_Lk(k) end - - k.pnt = pnt - k.isnewpnt = true + k.pnt .= pnt k.hasgH = false k.hasHi = false k.hasL = false @@ -115,9 +170,6 @@ function use_prevk(k::SumOfSqrData) end function inconek(k::SumOfSqrData) - @assert k.isnewpnt - k.isnewpnt = false - F = cholesky!(Symmetric(k.ip'*Diagonal(k.pnt)*k.ip), check=false) # TODO do inplace. TODO could this cholesky of P'DP be faster? if !issuccess(F) return false @@ -154,16 +206,16 @@ function calc_gk(k::SumOfSqrData) return k.grad end -function calc_Hinvk(k::SumOfSqrData) +function calc_Hik(k::SumOfSqrData) @assert k.hasgH if !k.hasHi - k.Hi = inv(k.Hfact) + k.Hi = inv(k.Hfact) # TODO maybe should do L\L'\A k.hasHi = true end return k.Hi end -function calc_HCholLk(k::SumOfSqrData) +function calc_Lk(k::SumOfSqrData) @assert k.hasgH if !k.hasL k.L = k.Hfact.L diff --git a/src/nativeinterface.jl b/src/nativeinterface.jl index d078b86ad..8da76c8f9 100644 --- a/src/nativeinterface.jl +++ b/src/nativeinterface.jl @@ -218,33 +218,14 @@ function solve!(alf::AlfonsoOpt) dir_ty = similar(ty) dir_tx = similar(tx) dir_ts = similar(ts) - # Hic = similar(tx) - # HiAt = similar(A, n, m) - # # AHi = similar(A, m, n) - # Hirxrs = similar(tx) - # lhsdydtau = similar(A, m+1, m+1) - # rhsdydtau = similar(tx, m+1) - # dydtau = similar(rhsdydtau) sa_tx = similar(tx) sa_ts = similar(ts) - # G = [zeros(m,m) A -b; -A' zeros(n,n) c; b' -c' 0] - # Gp = [zeros(m,m) A; -A' zeros(n,n)] - # Gps = [zeros(m,m) -A; -A' zeros(n,n)] - # lhs = Symmetric([zeros(m,m) -A; -A' zeros(n,n)]) - # W = zeros(n, n) - Wi = similar(A, n, n) - # rhs1 = [-b; -c] - # rhs2 = similar(rhs1) - # yx1 = similar(rhs1) - # yx2 = similar(rhs1) + HiAt = similar(A, n, m) # TODO for very sparse LPs, this is good, but for sparse problems with dense hessians, unclear y1 = similar(b) x1 = similar(c) y2 = similar(b) x2 = similar(c) - AWi = similar(A, m, n) - Witc = similar(c) - Witxs = similar(c) #= main loop @@ -302,65 +283,15 @@ function solve!(alf::AlfonsoOpt) prediction phase =# # determine prediction direction - - # @time begin - # invmu = 1.0/mu - # calc_Hinv_vec!(cones, coneidxs, Hic, c) - # Hic .*= invmu - # calc_Hinv_At!(cones, coneidxs, HiAt, A) - # HiAt .*= invmu - # dir_ts .= invmu*(rhs_tx - ts) - # calc_Hinv_vec!(cones, coneidxs, Hirxrs, dir_ts) - # - # # TODO maybe can use special structure of lhsdydtau: top left mxm is symmetric (L*A)^2, then last row and col are skew-symmetric - # lhsdydtau .= [A*HiAt (-b - A*Hic); (b' - c'*HiAt) (mu/tau^2 + dot(c, Hic))] - # rhsdydtau .= [(rhs_ty - A*Hirxrs); (rhs_tau - kap + dot(c, Hirxrs))] - # dydtau .= lhsdydtau\rhsdydtau - # - # dir_ty .= dydtau[1:m] - # dir_tau = dydtau[m+1] - # dir_tx .= Hirxrs + HiAt*dir_ty - Hic*dir_tau - # dir_ts .= -rhs_tx - A'*dir_ty + c*dir_tau - # dir_kap = -rhs_tau + dot(b, dir_ty) - dot(c, dir_tx) - # end - - # calc_W!(cones, coneidxs, W) - # H = W'*W - # - # lhs = G + [zeros(m,m) zeros(m,n) zeros(m,1); zeros(n,m) mu*H zeros(n,1); zeros(1,m) zeros(1,n) mu/tau^2] - # rhs = -G*[ty; tx; tau] - # @time yxt = lhs\rhs - # - # dir_ty = yxt[1:m] - # dir_tx = yxt[m+1:m+n] - # dir_tau = yxt[end] - # dir_ts .= -rhs_tx - A'*dir_ty + c*dir_tau - # dir_kap = -rhs_tau + dot(b, dir_ty) - dot(c, dir_tx) - # - # @show dir_tau, dir_kap - - # lhs.data[m+1:m+n,m+1:m+n] .= W'*W - # rhs2[1:m] .= -rhs_ty - # rhs2[m+1:m+n] .= rhs_tx - ts - # - # @time begin - # F = factorize(lhs) - # yx1 .= F\rhs1 - # yx2 .= F\rhs2 - # end - - #10.1 - calc_Wi!(cones, coneidxs, Wi) # TODO maybe can be faster using factorizations in cones - Wi .*= inv(sqrt(mu)) - AWi .= A*Wi - FAW = cholesky(Symmetric(AWi*AWi')) # TODO use structure of W (upper triangular), and precomputed factorize(A) - Witc .= Wi'*c - Witxs .= Wi'*(ts - rhs_tx) - # TODO can parallelize 1 and 2 - y1 .= FAW\(b + AWi*Witc) - x1 .= Wi*(AWi'*y1 - Witc) - y2 .= FAW\(rhs_ty + AWi*Witxs) - x2 .= Wi*(AWi'*y2 - Witxs) + invmu = 1.0/mu + calc_Hinv_At!(cones, coneidxs, HiAt, A) + HiAt .*= invmu + FAW = cholesky(Symmetric(A*HiAt)) + # # TODO can parallelize 1 and 2 + y1 .= FAW\(b + HiAt'*c) + calc_Hinv_vec!(cones, coneidxs, x1, invmu*(A'*y1 - c)) + y2 .= FAW\(rhs_ty + HiAt'*(ts - rhs_tx)) + calc_Hinv_vec!(cones, coneidxs, x2, invmu*(A'*y2 - ts + rhs_tx)) dir_tau = (rhs_tau - kap - dot(b, y2) + dot(c, x2))/(mu/tau^2 + dot(b, y1) - dot(c, x1)) dir_ty .= y2 + dir_tau*y1 @@ -387,10 +318,7 @@ function solve!(alf::AlfonsoOpt) sa_ts .= ts + alpha*dir_ts sa_tk = (tau + alpha*dir_tau)*(kap + alpha*dir_kap) sa_mu = (dot(sa_tx, sa_ts) + sa_tk)/bnu - calc_g!(cones, coneidxs, g) - calc_Wi!(cones, coneidxs, Wi) - Witxs .= Wi'*(sa_ts + sa_mu*g) - nbhd_beta = sqrt(dot(Witxs, Witxs) + (sa_tk - sa_mu)^2)/sa_mu + nbhd_beta = calc_nbhd(cones, coneidxs, sa_ts, sa_mu, sa_tk) if nbhd_beta < beta # iterate is inside the beta-neighborhood @@ -462,37 +390,16 @@ function solve!(alf::AlfonsoOpt) while true ncorrsteps += 1 # calculate correction direction - # invmu = 1.0/mu - # calc_Hinv_vec!(cones, coneidxs, Hic, c) - # Hic .*= invmu - # calc_Hinv_At!(cones, coneidxs, HiAt, A) - # HiAt .*= invmu - # dir_ts .= -invmu*ts - calc_g!(cones, coneidxs, g) - # calc_Hinv_vec!(cones, coneidxs, Hirxrs, dir_ts) - # - # lhsdydtau .= [A*HiAt (-b - A*Hic); (b' - c'*HiAt) (mu/tau^2 + dot(c, Hic))] - # rhsdydtau .= [-A*Hirxrs; (-kap + mu/tau + dot(c, Hirxrs))] - # dydtau .= lhsdydtau\rhsdydtau - # - # dir_ty .= dydtau[1:m] - # dir_tau = dydtau[m+1] - # dir_tx .= Hirxrs + HiAt*dir_ty - Hic*dir_tau - # dir_ts .= -A'*dir_ty + c*dir_tau - # dir_kap = dot(b, dir_ty) - dot(c, dir_tx) - - #10.1 - calc_Wi!(cones, coneidxs, Wi) # TODO maybe can be faster using factorizations in cones - Wi .*= inv(sqrt(mu)) - AWi .= A*Wi - FAW = cholesky(Symmetric(AWi*AWi')) # TODO use structure of W (upper triangular), and precomputed factorize(A) - Witc .= Wi'*c calc_g!(cones, coneidxs, g) - Witxs .= Wi'*(ts + mu*g) + invmu = 1.0/mu + calc_Hinv_At!(cones, coneidxs, HiAt, A) + HiAt .*= invmu + FAW = cholesky(Symmetric(A*HiAt), check=false) # TODO can parallelize 1 and 2 - y1 .= FAW\(b + AWi*Witc) - x1 .= Wi*(AWi'*y1 - Witc) - y2 .= FAW\(AWi*Witxs) - x2 .= Wi*(AWi'*y2 - Witxs) + y1 .= FAW\(b + HiAt'*c) + calc_Hinv_vec!(cones, coneidxs, x1, invmu*(A'*y1 - c)) + y2 .= FAW\(HiAt'*(ts + mu*g)) + calc_Hinv_vec!(cones, coneidxs, x2, invmu*(A'*y2 - ts) - g) dir_tau = (mu/tau - kap - dot(b, y2) + dot(c, x2))/(mu/tau^2 + dot(b, y1) - dot(c, x1)) dir_ty .= y2 + dir_tau*y1 @@ -540,10 +447,7 @@ function solve!(alf::AlfonsoOpt) # finish if allowed and current iterate is in the eta-neighborhood, or if taken max steps if (ncorrsteps == alf.maxcorrsteps) || alf.corrcheck - calc_g!(cones, coneidxs, g) - calc_Wi!(cones, coneidxs, Wi) - Witxs .= Wi'*(ts + mu*g) - if sqrt(dot(Witxs, Witxs) + (tau*kap - mu)^2) <= mu*eta + if calc_nbhd(cones, coneidxs, ts, mu, tau*kap) <= eta break elseif ncorrsteps == alf.maxcorrsteps # nbhd_eta > eta, so corrector failed @@ -623,58 +527,29 @@ function calc_g!(cones, coneidxs, g) return g end -# TODO store inv(W) inside the cones -function calc_Wi!(cones, coneidxs, Wi) +function calc_Hinv_vec!(cones, coneidxs, Hi_vec, v) + for k in eachindex(cones) + Hi_vec[coneidxs[k]] .= calc_Hik(cones[k])*v[coneidxs[k]] + end + return Hi_vec +end + +# TODO could save At submatrix inside cone objects +function calc_Hinv_At!(cones, coneidxs, Hi_At, A) for k in eachindex(cones) - Wi[coneidxs[k],coneidxs[k]] .= inv(calc_HCholLk(cones[k]))' + Hi_At[coneidxs[k],:] .= calc_Hik(cones[k])*A[:,coneidxs[k]]' end - return Wi + return Hi_At end -# function calc_Hinv_vec!(cones, coneidxs, Hi_vec, v) -# for k in eachindex(cones) -# Hi_vec[coneidxs[k]] .= calc_Hinvk(cones[k])*v[coneidxs[k]] -# end -# return Hi_vec -# end - -# # TODO could save At submatrix inside cone objects -# function calc_Hinv_At!(cones, coneidxs, Hi_At, A) -# for k in eachindex(cones) -# Hi_At[coneidxs[k],:] .= calc_Hinvk(cones[k])*A[:,coneidxs[k]]' -# end -# return Hi_At -# end -# -# function calc_A_Hinv!(cones, coneidxs, A_Hi, A) -# for k in eachindex(cones) -# A_Hi[:,coneidxs[k]] .= A[:,coneidxs[k]]*calc_Hinvk(cones[k]) -# end -# return A_Hi -# end - -# function calc_nbhd(cones, coneidxs, ts, mu, tk) -# # sqrt(sum(abs2, L\(ts + mu*g)) + (tau*kap - mu)^2)/mu -# sumsqr = (tk - mu)^2 -# for k in eachindex(cones) -# sumsqr += sum(abs2, calc_HCholLk(cones[k])\(ts[coneidxs[k]] + mu*calc_gk(cones[k]))) -# end -# return sqrt(sumsqr)/mu -# end - -# function calc_L!(cones, coneidxs, L) -# for k in eachindex(cones) -# L[coneidxs[k],coneidxs[k]] .= calc_HCholLk(cones[k]) -# end -# return L -# end -# -# function calc_W!(cones, coneidxs, W) -# for k in eachindex(cones) -# W[coneidxs[k],coneidxs[k]] .= calc_HCholLk(cones[k])' -# end -# return W -# end +function calc_nbhd(cones, coneidxs, ts, mu, tk) + # sqrt(sum(abs2, L\(ts + mu*g)) + (tau*kap - mu)^2)/mu + sumsqr = (tk - mu)^2 + for k in eachindex(cones) + sumsqr += sum(abs2, calc_Lk(cones[k])\(ts[coneidxs[k]] + mu*calc_gk(cones[k]))) + end + return sqrt(sumsqr)/mu +end function getbetaeta(maxcorrsteps, bnu) if maxcorrsteps <= 2 diff --git a/test/nativeexamples.jl b/test/nativeexamples.jl index a55b4f0b5..92dcf514e 100644 --- a/test/nativeexamples.jl +++ b/test/nativeexamples.jl @@ -99,14 +99,14 @@ end # @test Alfonso.get_dobj(alf) ≈ -0.25 atol=1e-4 rtol=1e-4 # end -@testset "Motzkin" begin - alf = Alfonso.AlfonsoOpt(verbose=verbflag) - build_namedpoly!(alf, :motzkin, 7) - @time Alfonso.solve!(alf) - @test Alfonso.get_status(alf) == :Optimal - @test Alfonso.get_pobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 - @test Alfonso.get_dobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 -end +# @testset "Motzkin" begin +# alf = Alfonso.AlfonsoOpt(verbose=verbflag) +# build_namedpoly!(alf, :motzkin, 7) +# @time Alfonso.solve!(alf) +# @test Alfonso.get_status(alf) == :Optimal +# @test Alfonso.get_pobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 +# @test Alfonso.get_dobj(alf) ≈ 0 atol=1e-4 rtol=1e-4 +# end @testset "Reaction-diffusion" begin alf = Alfonso.AlfonsoOpt(verbose=verbflag)