diff --git a/src/Bridges/Variable/bridges/set_dot.jl b/src/Bridges/Variable/bridges/set_dot.jl index b2643321f3..ba593a717a 100644 --- a/src/Bridges/Variable/bridges/set_dot.jl +++ b/src/Bridges/Variable/bridges/set_dot.jl @@ -1,6 +1,6 @@ struct DotProductsBridge{T,S,V} <: SetMapBridge{T,S,MOI.SetWithDotProducts{S,V}} variables::Vector{MOI.VariableIndex} - constraint::MOI.ConstraintIndex{MOI.VectorOfVariables, S} + constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,S} set::MOI.SetWithDotProducts{S,V} end @@ -28,10 +28,7 @@ function bridge_constrained_variable( return BT(variables, constraint, set) end -function MOI.Bridges.map_set( - bridge::DotProductsBridge{T,S}, - set::S, -) where {T,S} +function MOI.Bridges.map_set(bridge::DotProductsBridge{T,S}, set::S) where {T,S} return MOI.SetWithDotProducts(set, bridge.vectors) end @@ -49,14 +46,27 @@ function MOI.Bridges.map_function( ) where {T} scalars = MOI.Utilities.eachscalar(func) if i.value in eachindex(bridge.set.vectors) - return MOI.Utilities.set_dot(bridge.set.vectors[i.value], scalars, bridge.set.set) + return MOI.Utilities.set_dot( + bridge.set.vectors[i.value], + scalars, + bridge.set.set, + ) else - return convert(MOI.ScalarAffineFunction{T}, scalars[i.value - length(bridge.vectors)]) + return convert( + MOI.ScalarAffineFunction{T}, + scalars[i.value-length(bridge.vectors)], + ) end end -function MOI.Bridges.inverse_map_function(bridge::DotProductsBridge{T}, func) where {T} +function MOI.Bridges.inverse_map_function( + bridge::DotProductsBridge{T}, + func, +) where {T} m = length(bridge.set.vectors) - return MOI.Utilities.operate(vcat, T, MOI.Utilities.eachscalar(func)[(m+1):end]) + return MOI.Utilities.operate( + vcat, + T, + MOI.Utilities.eachscalar(func)[(m+1):end], + ) end - diff --git a/src/Utilities/functions.jl b/src/Utilities/functions.jl index b2579159f0..2d558f92e4 100644 --- a/src/Utilities/functions.jl +++ b/src/Utilities/functions.jl @@ -638,7 +638,8 @@ end A type that allows iterating over the scalar-functions that comprise an `AbstractVectorFunction`. """ -struct ScalarFunctionIterator{F<:MOI.AbstractVectorFunction,C,S} <: AbstractVector{S} +struct ScalarFunctionIterator{F<:MOI.AbstractVectorFunction,C,S} <: + AbstractVector{S} f::F # Cache that can be used to store a precomputed datastructure that allows # an efficient implementation of `getindex`. diff --git a/src/sets.jl b/src/sets.jl index d7b8607ecf..2af3b3e4cd 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -1804,73 +1804,129 @@ function Base.getproperty( end """ - SetWithDotProducts(set, vectors::Vector{V}) + SetWithDotProducts(set::MOI.AbstractSet, vectors::AbstractVector) -Given a set `set` of dimension `d` and `m` vectors `v_1`, ..., `v_m` given in `vectors`, this is the set: -``\\{ ((\\langle v_1, x \\rangle, ..., \\langle v_m, x \\rangle, x) \\in \\mathbb{R}^{m + d} : x \\in \\text{set} \\}.`` +Given a set `set` of dimension `d` and `m` vectors `a_1`, ..., `a_m` given in `vectors`, this is the set: +``\\{ ((\\langle a_1, x \\rangle, ..., \\langle a_m, x \\rangle) \\in \\mathbb{R}^{m} : x \\in \\text{set} \\}.`` """ -struct SetWithDotProducts{S,V} <: AbstractVectorSet +struct SetWithDotProducts{S,A,V<:AbstractVector{A}} <: AbstractVectorSet set::S - vectors::Vector{V} + vectors::V end function Base.copy(s::SetWithDotProducts) return SetWithDotProducts(copy(s.set), copy(s.vectors)) end -function dimension(s::SetWithDotProducts) - return length(s.vectors) + dimension(s.set) -end +dimension(s::SetWithDotProducts) = length(s.vectors) function dual_set(s::SetWithDotProducts) return LinearCombinationInSet(s.set, s.vectors) end -function dual_set_type(::Type{SetWithDotProducts{S,V}}) where {S,V} - return LinearCombinationInSet{S,V} +function dual_set_type(::Type{SetWithDotProducts{S,A,V}}) where {S,A,V} + return LinearCombinationInSet{S,A,V} end """ - LinearCombinationInSet{S,V}(set::S, matrices::Vector{V}) + LinearCombinationInSet(set::MOI.AbstractSet, matrices::AbstractVector) -Given a set `set` of dimension `d` and `m` vectors `v_1`, ..., `v_m` given in `vectors`, this is the set: -``\\{ ((y, x) \\in \\mathbb{R}^{m + d} : \\sum_{i=1}^m y_i v_i + x \\in \\text{set} \\}.`` +Given a set `set` of dimension `d` and `m` vectors `a_1`, ..., `a_m` given in `vectors`, this is the set: +``\\{ (y \\in \\mathbb{R}^{m} : \\sum_{i=1}^m y_i a_i \\in \\text{set} \\}.`` """ -struct LinearCombinationInSet{S,V} <: AbstractVectorSet +struct LinearCombinationInSet{S,A,V<:AbstractVector{A}} <: AbstractVectorSet set::S - vectors::Vector{V} + vectors::V end -dimension(s::LinearCombinationInSet) = length(s.vectors) + simension(s.set) +dimension(s::LinearCombinationInSet) = length(s.vectors) function dual_set(s::LinearCombinationInSet) - return SetWithDotProducts(s.side_dimension, s.matrices) + return SetWithDotProducts(s.side_dimension, s.vectors) +end + +function dual_set_type(::Type{LinearCombinationInSet{S,A,V}}) where {S,A,V} + return SetWithDotProducts{S,A,V} end -function dual_set_type(::Type{LinearCombinationInSet{S,V}}) where {S,V} - return SetWithDotProducts{S,V} +abstract type AbstractFactorization{T,F} <: AbstractMatrix{T} end + +function Base.size(m::AbstractFactorization) + n = size(m.factor, 1) + return (n, n) end """ - struct LowRankMatrix{T} - diagonal::Vector{T} - factor::Matrix{T} + struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, + } <: AbstractMatrix{T} + factor::F + scaling::D end -`factor * Diagonal(diagonal) * factor'`. -""" -struct LowRankMatrix{T} <: AbstractMatrix{T} - diagonal::Vector{T} - factor::Matrix{T} +Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. +If `factor` is a vector and `diagonal` is a scalar, this corresponds to +the matrix `diagonal * factor * factor'`. +If `factor` is a matrix and `diagonal` is a vector, this corresponds to +the matrix `factor * Diagonal(scaling) * factor'`. +""" +struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, +} <: AbstractFactorization{T,F} + factor::F + scaling::D + function Factorization( + factor::AbstractMatrix{T}, + scaling::AbstractVector{T}, + ) where {T} + if length(scaling) != size(factor, 2) + error( + "Length `$(length(scaling))` of diagonal does not match number of columns `$(size(factor, 2))` of factor", + ) + end + return new{T,typeof(factor),typeof(scaling)}(factor, scaling) + end + function Factorization(factor::AbstractVector{T}, scaling::T) where {T} + return new{T,typeof(factor),typeof(scaling)}(factor, scaling) + end end -function Base.size(m::LowRankMatrix) - n = size(m.factor, 1) - return (n, n) +function Base.getindex(m::Factorization, i::Int, j::Int) + return sum( + m.factor[i, k] * m.scaling[k] * m.factor[j, k]' for + k in eachindex(m.scaling) + ) +end + +""" + struct Factorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, + D<:Union{T,AbstractVector{T}}, + } <: AbstractMatrix{T} + factor::F + scaling::D + end + +Matrix corresponding to `factor * Diagonal(diagonal) * factor'`. +If `factor` is a vector and `diagonal` is a scalar, this corresponds to +the matrix `diagonal * factor * factor'`. +If `factor` is a matrix and `diagonal` is a vector, this corresponds to +the matrix `factor * Diagonal(scaling) * factor'`. +""" +struct PositiveSemidefiniteFactorization{ + T, + F<:Union{AbstractVector{T},AbstractMatrix{T}}, +} <: AbstractFactorization{T,F} + factor::F end -function Base.getindex(m::LowRankMatrix, i::Int, j::Int) - return sum(m.factor[i, k] * m.diagonal[k] * m.factor[j, k]' for k in eachindex(m.diagonal)) +function Base.getindex(m::PositiveSemidefiniteFactorization, i::Int, j::Int) + return sum(m.factor[i, k] * m.factor[j, k]' for k in axes(m.factor, 2)) end struct TriangleVectorization{T,M<:AbstractMatrix{T}} <: AbstractVector{T} diff --git a/test/sets.jl b/test/sets.jl index bf52326974..122a562f7b 100644 --- a/test/sets.jl +++ b/test/sets.jl @@ -8,6 +8,7 @@ module TestSets using Test import MathOptInterface as MOI +import LinearAlgebra include("dummy.jl") @@ -353,6 +354,41 @@ function test_sets_reified() return end +function _test_factorization(A, B) + @test size(A) == size(B) + @test A ≈ B + d = LinearAlgebra.checksquare(A) + n = div(d * (d + 1), 2) + vA = MOI.TriangleVectorization(A) + @test length(vA) == n + @test eachindex(vA) == Base.OneTo(n) + vB = MOI.TriangleVectorization(B) + @test length(vB) == n + @test eachindex(vA) == Base.OneTo(n) + k = 0 + for j in 1:d + for i in 1:j + k += 1 + @test vA[k] == vB[k] + @test vA[k] == A[i, j] + end + end + return +end + +function test_factorizations() + f = [1, 2] + _test_factorization(f * f', MOI.PositiveSemidefiniteFactorization(f)) + _test_factorization(2 * f * f', MOI.Factorization(f, 2)) + F = [1 2; 3 4; 5 6] + d = [7, 8] + _test_factorization(F * F', MOI.PositiveSemidefiniteFactorization(F)) + return _test_factorization( + F * LinearAlgebra.Diagonal(d) * F', + MOI.Factorization(F, d), + ) +end + function runtests() for name in names(@__MODULE__; all = true) if startswith("$name", "test_")