Updates

Author: blegat

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
-

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`.

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}

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_")