Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Approximate PSD projections #101

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Approximate PSD projections #101

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
62af097
Initial commit with LOBPCG projection :tada:
nrontsis Dec 9, 2019
b9ca329
Add commented tests
nrontsis Dec 9, 2019
97e8463
Include custom lmul! and rmul! to support Julia 1.0
nrontsis Dec 9, 2019
d37fb16
Remove deprecated dependency
nrontsis Dec 9, 2019
0654276
Average lobpcg/exact iterations
nrontsis Dec 9, 2019
fb81907
Attempt to solve LinearAlgebra conflicts
nrontsis Dec 9, 2019
ece3780
Removing unnecessary dependency
nrontsis Dec 9, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 20 additions & 20 deletions src/MOIWrapper.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -359,7 +359,7 @@ function processconstraints!(optimizer::Optimizer, src::MOI.ModelLike, idxmap, r
set_constant = zeros(Float64, m)
# loop over constraints and modify A, l, u and constants
for (F, S) in MOI.get(src, MOI.ListOfConstraints())
processconstraints!((I, J, V), b, COSMOconvexSets, constant, set_constant, src, idxmap, rowranges, F, S)
processconstraints!(optimizer, (I, J, V), b, COSMOconvexSets, constant, set_constant, src, idxmap, rowranges, F, S)
end
optimizer.set_constant = set_constant
# subtract constant from right hand side
Expand All @@ -379,7 +379,7 @@ function processconstraints!(optimizer::Optimizer, src::MOI.ModelLike, idxmap, r
return nothing
end

function processconstraints!(triplets::SparseTriplets, b::Vector, COSMOconvexSets::Array{COSMO.AbstractConvexSet{Float64}}, constant::Vector{Float64}, set_constant::Vector{Float64},
function processconstraints!(optimizer, triplets::SparseTriplets, b::Vector, COSMOconvexSets::Array{COSMO.AbstractConvexSet{Float64}}, constant::Vector{Float64}, set_constant::Vector{Float64},
src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}},
F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
Expand All @@ -394,7 +394,7 @@ function processconstraints!(triplets::SparseTriplets, b::Vector, COSMOconvexSet
set_constant[rows] = MOI.constant(s)
end
processConstraint!(triplets, f, rows, idxmap, s)
processSet!(b, rows, COSMOconvexSets, s)
processSet!(b, rows, COSMOconvexSets, s, optimizer.inner.settings)
end
nothing
end
Expand Down Expand Up @@ -425,7 +425,7 @@ end


# process function like f(x) = coeff*x_1
function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunction, row::Int, idxmap::MOIU.IndexMap, s::MOI.AbstractSet)
function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunction, row::Int, idxmap, s::MOI.AbstractSet)
(I, J, V) = triplets
for term in f.terms
push!(I, row)
Expand All @@ -434,7 +434,7 @@ function processConstraint!(triplets::SparseTriplets, f::MOI.ScalarAffineFunctio
end
end

function processConstraint!(triplets::SparseTriplets, f::MOI.VectorOfVariables, rows::UnitRange{Int}, idxmap::MOIU.IndexMap, s::MOI.AbstractSet)
function processConstraint!(triplets::SparseTriplets, f::MOI.VectorOfVariables, rows::UnitRange{Int}, idxmap, s::MOI.AbstractSet)
(I, J, V) = triplets
cols = zeros(length(rows))
for (i, var) in enumerate(f.variables)
Expand Down Expand Up @@ -470,75 +470,75 @@ end
##############################

# process the following sets Union{LessThan, GreaterThan, EqualTo}
function processSet!(b::Vector, row::Int, cs, s::LessThan)
function processSet!(b::Vector, row::Int, cs, s::LessThan, settings::Settings)
b[row] += s.upper
push!(cs, COSMO.Nonnegatives{Float64}(1))
nothing
end
function processSet!(b::Vector, row::Int, cs, s::GreaterThan)
function processSet!(b::Vector, row::Int, cs, s::GreaterThan, settings::Settings)
b[row] -= s.lower
push!(cs, COSMO.Nonnegatives{Float64}(1))
nothing
end
function processSet!(b::Vector, row::Int, cs, s::EqualTo)
function processSet!(b::Vector, row::Int, cs, s::EqualTo, settings::Settings)
b[row] += s.value
push!(cs, COSMO.ZeroSet{Float64}(1))
nothing
end

function processSet!(b::Vector, row::Int, cs, s::MOI.Interval)
function processSet!(b::Vector, row::Int, cs, s::MOI.Interval, settings::Settings)
push!(cs, COSMO.Box{Float64}([s.lower], [s.upper]))
nothing
end

# process the following sets Zeros
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::Zeros)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::Zeros, settings::Settings)
push!(cs, COSMO.ZeroSet{Float64}(length(rows)))
nothing
end

# process the following sets Union{Zeros, Nonnegatives, Nonpositives}
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonnegatives)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonnegatives, settings::Settings)
push!(cs, COSMO.Nonnegatives{Float64}(length(rows)))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonpositives)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.Nonpositives, settings::Settings)
push!(cs, COSMO.Nonnegatives{Float64}(length(rows)))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::SOC)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::SOC, settings::Settings)
push!(cs, COSMO.SecondOrderCone{Float64}(length(rows)))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeSquare)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeSquare, settings::Settings)
push!(cs, COSMO.PsdCone{Float64}(length(rows)))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeTriangle)
push!(cs, COSMO.PsdConeTriangle{Float64}(length(rows)))
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PositiveSemidefiniteConeTriangle, settings::Settings)
push!(cs, settings.psd_projector(length(rows)))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.ExponentialCone)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.ExponentialCone, settings::Settings)
push!(cs, COSMO.ExponentialCone{Float64}())
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualExponentialCone)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualExponentialCone, settings::Settings)
push!(cs, COSMO.DualExponentialCone{Float64}())
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PowerCone)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.PowerCone, settings::Settings)
push!(cs, COSMO.PowerCone{Float64}(s.exponent))
nothing
end

function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualPowerCone)
function processSet!(b::Vector, rows::UnitRange{Int}, cs, s::MOI.DualPowerCone, settings::Settings)
push!(cs, COSMO.DualPowerCone{Float64}(s.exponent))
nothing
end
Expand Down
58 changes: 52 additions & 6 deletions src/algebra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
using LinearAlgebra
import LinearAlgebra.lmul!
import LinearAlgebra.rmul!
const IdentityMatrix = UniformScaling{Bool}


Expand Down Expand Up @@ -111,11 +113,11 @@ end

lmul!(L::IdentityMatrix, M::AbstractMatrix) = L.λ ? M : M .= zero(eltype(M))

function lmul!(L::Diagonal, x::AbstractVector)
(length(L.diag) == length(x)) || throw(DimensionMismatch())
@. x = x * L.diag
return nothing
end
function lmul!(L::Diagonal, x::AbstractVector)
(length(L.diag) == length(x)) || throw(DimensionMismatch())
@. x = x * L.diag
return nothing
end

lmul!(L::IdentityMatrix, x::AbstractVector) = L.λ ? x : x .= zero(eltype(x))

Expand Down Expand Up @@ -217,4 +219,48 @@ function is_neg_def!(X::AbstractMatrix{T}, tol::T=zero(T)) where T
end

is_pos_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(copy(X), tol)
is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)
is_neg_def(X::AbstractMatrix{T}, tol::T=zero(T)) where T = is_pos_def!(-X, tol)

function populate_upper_triangle!(A::AbstractMatrix{T}, x::AbstractVector{T}, scaling_factor::T=T(1/sqrt(2))) where T
k = 0
for j in 1:size(A, 2)
for i in 1:j
k += 1
if i != j
A[i, j] = scaling_factor * x[k]
else
A[i, j] = x[k]
end
end
end
nothing
end

function extract_upper_triangle!(A::AbstractMatrix{T}, x::AbstractVector{T}, scaling_factor::T=T(sqrt(2))) where T
k = 0
for j in 1:size(A, 2)
for i in 1:j
k += 1
if i != j
x[k] = scaling_factor * A[i, j]
else
x[k] = A[i, j]
end
end
end
nothing
end

function populate_upper_triangle(x::AbstractVector{T}, scaling_factor::T=T(1/sqrt(2))) where T
n = Int(1/2*(sqrt(8*length(x) + 1) - 1)) # Solution of (n^2 + n)/2 = length(x) obtained by WolframAlpha
A = zeros(n, n)
populate_upper_triangle!(A, x, scaling_factor)
return Symmetric(A)
end

function extract_upper_triangle(A::AbstractMatrix{T}, scaling_factor::T=T(sqrt(2))) where T
n = size(A, 2)
x = zeros(Int(n*(n + 1)/2))
extract_upper_triangle!(A, x, scaling_factor)
return x
end
48 changes: 11 additions & 37 deletions src/convexset.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,11 +317,11 @@ mutable struct PsdConeTriangle{T} <: AbstractConvexCone{T}
clique_ind::Int64

function PsdConeTriangle{T}(dim::Int, tree_ind::Int64, clique_ind::Int64) where{T}
dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative"))
side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension), tree_ind, clique_ind)
dim >= 0 || throw(DomainError(dim, "dimension must be nonnegative"))
side_dimension = Int(sqrt(0.25 + 2 * dim) - 0.5);
new(dim, side_dimension, zeros(side_dimension, side_dimension), PsdBlasWorkspace{T}(side_dimension), tree_ind, clique_ind)

end
end
end
PsdConeTriangle(dim) = PsdConeTriangle{DefaultFloat}(dim)
PsdConeTriangle{T}(dim::Int64) where{T} = PsdConeTriangle{T}(dim, 0, 0)
Expand Down Expand Up @@ -357,54 +357,28 @@ function project!(x::AbstractArray, cone::Union{PsdConeTriangle{T}, DensePsdCone
end

# Notice that we are using a (faster) in-place version that modifies the input
function in_dual!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
function in_dual!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where{T}
n = cone.sqrt_dim
populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
return COSMO.is_pos_def!(cone.X, tol)
end
in_dual(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_dual!(x, cone, tol)
in_dual(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where {T} = in_dual!(x, cone, tol)

function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where{T}
function in_pol_recc!(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where{T}
n = cone.sqrt_dim
populate_upper_triangle!(cone.X, x, 1 / sqrt(2))
Xs = Symmetric(cone.X)
return COSMO.is_neg_def!(cone.X, tol)
end
in_pol_recc(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}, tol::T) where {T} = in_pol_recc!(x, cone, tol)
in_pol_recc(x::AbstractVector{T}, cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}}, tol::T) where {T} = in_pol_recc!(x, cone, tol)


function allocate_memory!(cone::Union{PsdConeTriangle{T}, DensePsdConeTriangle{T}}) where {T}
cone.X = zeros(cone.sqrt_dim, cone.sqrt_dim)
end

function populate_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
k = 0
for j in 1:size(A, 2)
for i in 1:j
k += 1
if i != j
A[i, j] = scaling_factor * x[k]
else
A[i, j] = x[k]
end
end
end
nothing
end

function extract_upper_triangle!(A::AbstractMatrix, x::AbstractVector, scaling_factor::Float64)
k = 0
for j in 1:size(A, 2)
for i in 1:j
k += 1
if i != j
x[k] = scaling_factor * A[i, j]
else
x[k] = A[i, j]
end
end
end
nothing
function reset_iteration_counters!(cone::AbstractConvexSet)
nothing
end

"""
Expand Down Expand Up @@ -848,7 +822,7 @@ function scale!(cone::AbstractConvexCone{T}, ::AbstractVector{T}) where{T}
return nothing
end

function rectify_scaling!(E, work, set::Union{SecondOrderCone{T}, PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, DensePsdConeTriangle{T}, PowerCone{T}, DualPowerCone{T}, ExponentialCone{T}, DualExponentialCone{T}}) where{T}
function rectify_scaling!(E, work, set::Union{SecondOrderCone{T}, PsdCone{T}, DensePsdCone{T}, PsdConeTriangle{T}, PsdConeTriangleLOBPCG{T}, DensePsdConeTriangle{T}, PowerCone{T}, DualPowerCone{T}, ExponentialCone{T}, DualExponentialCone{T}}) where{T}
return rectify_scalar_scaling!(E, work)
end
rectify_scaling!(E, work, set::Union{ZeroSet{<:Real}, Nonnegatives{<:Real}, Box{<:Real}}) = false
Expand Down
Loading