From ea4c989b8a4705d236bdc2b2df069588f6f0776d Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 17 Apr 2025 12:38:23 +1200 Subject: [PATCH 01/11] [Utilities] add distance_to_set for PositiveSemidefiniteCone --- src/Utilities/distance_to_set.jl | 36 +++++++++++++++++++++++++++++++ test/Utilities/distance_to_set.jl | 22 +++++++++++++++++++ 2 files changed, 58 insertions(+) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index 3dc1464d4d..78a25bfc1f 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -489,3 +489,39 @@ function distance_to_set( elements = [x[i] for i in eachindex(x) if !(i in pairs[k])] return LinearAlgebra.norm(elements, 2) end + +function _reshape(x::AbstractVector, set::MOI.PositiveSemidefiniteConeSquare) + n = MOI.side_dimension(set) + return reshape(x, (n, n)) +end + +function _reshape(x::AbstractVector, set::MOI.PositiveSemidefiniteConeTriangle) + n = MOI.side_dimension(set) + X = zeros(eltype(x), n, n) + k = 1 + for i in 1:n + for j in 1:i + X[j,i] = X[i,j] = x[k] + k += 1 + end + end + return LinearAlgebra.Symmetric(X) +end + +function distance_to_set( + ::ProjectionUpperBoundDistance, + x::AbstractVector{T}, + set::Union{ + MOI.PositiveSemidefiniteConeSquare, + MOI.PositiveSemidefiniteConeTriangle, + }, +) where {T<:Real} + _check_dimension(x, set) + λ, U = LinearAlgebra.eigen(_reshape(x, set)) + if minimum(λ) >= 0 + return 0.0 + end + λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ)) + A = LinearAlgebra.Symmetric(U * λ_negative * U') + return LinearAlgebra.norm(A, 2) +end diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 558080efaf..3e4fa4d16a 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -308,6 +308,28 @@ function test_sos2() return end +function test_positivesemidefiniteconesquare() + _test_set( + MOI.PositiveSemidefiniteConeSquare(2), + [1.0, 0.0, 0.0, 1.0] => 0.0, + [1.0, -1.0, -1.0, 1.0] => 0.0, + [1.0, -2.0, -2.0, 1.0] => 1.0; + mismatch = [1.0], + ) + return +end + +function test_positivesemidefiniteconetriangle() + _test_set( + MOI.PositiveSemidefiniteConeTriangle(2), + [1.0, 0.0, 1.0] => 0.0, + [1.0, -1.0, 1.0] => 0.0, + [1.0, -2.0, 1.0] => 1.0; + mismatch = [1.0], + ) + return +end + end TestFeasibilityChecker.runtests() From 67480b7509da6c2426ab7c0775484d45f8fd68f4 Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 17 Apr 2025 12:58:21 +1200 Subject: [PATCH 02/11] Fix formatting --- src/Utilities/distance_to_set.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index 78a25bfc1f..160f087677 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -501,7 +501,7 @@ function _reshape(x::AbstractVector, set::MOI.PositiveSemidefiniteConeTriangle) k = 1 for i in 1:n for j in 1:i - X[j,i] = X[i,j] = x[k] + X[j, i] = X[i, j] = x[k] k += 1 end end From bce7f9f77e91c5490cd85b7eceb6fd9958c478d2 Mon Sep 17 00:00:00 2001 From: odow Date: Thu, 17 Apr 2025 14:50:34 +1200 Subject: [PATCH 03/11] Update --- src/Utilities/distance_to_set.jl | 19 +++++++++++++++++-- test/Utilities/distance_to_set.jl | 8 +++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index 160f087677..ec23e31eb3 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -508,6 +508,20 @@ function _reshape(x::AbstractVector, set::MOI.PositiveSemidefiniteConeTriangle) return LinearAlgebra.Symmetric(X) end +""" + distance_to_set( + ::ProjectionUpperBoundDistance, + x::AbstractVector, + set::Union{ + MOI.PositiveSemidefiniteConeSquare, + MOI.PositiveSemidefiniteConeTriangle, + }, + ) + +Let ``X`` be `x` reshaped into the appropriate matrix. The returned distance is +``||X - Y||_2^2`` where ``Y`` is the eigen decomposition of ``X`` with negative +eigen values removed. +""" function distance_to_set( ::ProjectionUpperBoundDistance, x::AbstractVector{T}, @@ -518,10 +532,11 @@ function distance_to_set( ) where {T<:Real} _check_dimension(x, set) λ, U = LinearAlgebra.eigen(_reshape(x, set)) - if minimum(λ) >= 0 + if minimum(λ) >= zero(T) return 0.0 end λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ)) A = LinearAlgebra.Symmetric(U * λ_negative * U') - return LinearAlgebra.norm(A, 2) + # Σ A^2 is needed to match the definition of Utilities.set_dot + return sum(Base.Fix2(^, 2), A) end diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 3e4fa4d16a..78374eca1d 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -30,7 +30,7 @@ function _test_set(set, pairs...; mismatch = nothing) ) end for (x, d) in pairs - @test MOI.Utilities.distance_to_set(x, set) ≈ d + @test ≈(MOI.Utilities.distance_to_set(x, set), d; atol = 1e-12) end return end @@ -313,7 +313,8 @@ function test_positivesemidefiniteconesquare() MOI.PositiveSemidefiniteConeSquare(2), [1.0, 0.0, 0.0, 1.0] => 0.0, [1.0, -1.0, -1.0, 1.0] => 0.0, - [1.0, -2.0, -2.0, 1.0] => 1.0; + [1.0, -2.0, -2.0, 1.0] => 1.0 + [1.0 1.1; 1.1 -2.3] => 2.6330532015051946; mismatch = [1.0], ) return @@ -324,7 +325,8 @@ function test_positivesemidefiniteconetriangle() MOI.PositiveSemidefiniteConeTriangle(2), [1.0, 0.0, 1.0] => 0.0, [1.0, -1.0, 1.0] => 0.0, - [1.0, -2.0, 1.0] => 1.0; + [1.0, -2.0, 1.0] => 1.0 + [1.0, 1.1, -2.3] => 2.6330532015051946; mismatch = [1.0], ) return From 3f918fc253f0b722a9755030b2ec426d7208bf14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 17 Apr 2025 09:51:23 +0200 Subject: [PATCH 04/11] Fix --- test/Utilities/distance_to_set.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 78374eca1d..28c175e509 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -313,7 +313,7 @@ function test_positivesemidefiniteconesquare() MOI.PositiveSemidefiniteConeSquare(2), [1.0, 0.0, 0.0, 1.0] => 0.0, [1.0, -1.0, -1.0, 1.0] => 0.0, - [1.0, -2.0, -2.0, 1.0] => 1.0 + [1.0, -2.0, -2.0, 1.0] => 1.0, [1.0 1.1; 1.1 -2.3] => 2.6330532015051946; mismatch = [1.0], ) From 9bd6ee5fe03f3e388ef99871a5d1d245780810c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 17 Apr 2025 09:56:30 +0200 Subject: [PATCH 05/11] Missing comma --- test/Utilities/distance_to_set.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 28c175e509..497e175a3b 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -325,7 +325,7 @@ function test_positivesemidefiniteconetriangle() MOI.PositiveSemidefiniteConeTriangle(2), [1.0, 0.0, 1.0] => 0.0, [1.0, -1.0, 1.0] => 0.0, - [1.0, -2.0, 1.0] => 1.0 + [1.0, -2.0, 1.0] => 1.0, [1.0, 1.1, -2.3] => 2.6330532015051946; mismatch = [1.0], ) From 89471bb90492805df3b2ee8800a5b10bb29d1b79 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 17 Apr 2025 11:01:30 +0200 Subject: [PATCH 06/11] Replace eigenvalue by its square --- test/Utilities/distance_to_set.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 497e175a3b..82013c974e 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -314,7 +314,7 @@ function test_positivesemidefiniteconesquare() [1.0, 0.0, 0.0, 1.0] => 0.0, [1.0, -1.0, -1.0, 1.0] => 0.0, [1.0, -2.0, -2.0, 1.0] => 1.0, - [1.0 1.1; 1.1 -2.3] => 2.6330532015051946; + [1.0, 1.1, 1.1, -2.3] => 6.9329523025; mismatch = [1.0], ) return @@ -326,7 +326,7 @@ function test_positivesemidefiniteconetriangle() [1.0, 0.0, 1.0] => 0.0, [1.0, -1.0, 1.0] => 0.0, [1.0, -2.0, 1.0] => 1.0, - [1.0, 1.1, -2.3] => 2.6330532015051946; + [1.0, 1.1, -2.3] => 6.9329523025; mismatch = [1.0], ) return From 0799b1283898a13e17c45357159d2de62a4ac546 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 17 Apr 2025 14:41:21 +0200 Subject: [PATCH 07/11] Fix --- src/Utilities/distance_to_set.jl | 22 +++++++++++++++------- test/Utilities/distance_to_set.jl | 4 ++-- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index ec23e31eb3..a16b924b52 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -531,12 +531,20 @@ function distance_to_set( }, ) where {T<:Real} _check_dimension(x, set) - λ, U = LinearAlgebra.eigen(_reshape(x, set)) - if minimum(λ) >= zero(T) - return 0.0 + # Let + # 1) `U` be the matrix of eigenvectors. + # 2) `λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ))` + # 3) `A = LinearAlgebra.Symmetric(U * λ_negative * U')` + # The upper bound distance is the Frobenius norm of `A` which is the square + # root of the sum of the squares of the negative eigenvalues. + # Indeed: + # `⟨U * λ_negative * U', U * λ_negative * U'⟩` + # which is equal to + # `⟨λ_negative * U' * U, U' * U * λ_negative⟩` + # Since `U'U = I` this is equal to + # `⟨λ_negative, λ_negative⟩` + # So we need the return: + return √sum(LinearAlgebra.eigvals(_reshape(x, set))) do λ + min(λ, zero(T))^2 end - λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ)) - A = LinearAlgebra.Symmetric(U * λ_negative * U') - # Σ A^2 is needed to match the definition of Utilities.set_dot - return sum(Base.Fix2(^, 2), A) end diff --git a/test/Utilities/distance_to_set.jl b/test/Utilities/distance_to_set.jl index 82013c974e..75f8e357eb 100644 --- a/test/Utilities/distance_to_set.jl +++ b/test/Utilities/distance_to_set.jl @@ -314,7 +314,7 @@ function test_positivesemidefiniteconesquare() [1.0, 0.0, 0.0, 1.0] => 0.0, [1.0, -1.0, -1.0, 1.0] => 0.0, [1.0, -2.0, -2.0, 1.0] => 1.0, - [1.0, 1.1, 1.1, -2.3] => 6.9329523025; + [1.0, 1.1, 1.1, -2.3] => 2.633053201505194; mismatch = [1.0], ) return @@ -326,7 +326,7 @@ function test_positivesemidefiniteconetriangle() [1.0, 0.0, 1.0] => 0.0, [1.0, -1.0, 1.0] => 0.0, [1.0, -2.0, 1.0] => 1.0, - [1.0, 1.1, -2.3] => 6.9329523025; + [1.0, 1.1, -2.3] => 2.633053201505194; mismatch = [1.0], ) return From 420cf5362e89a30a33bc7ebd881dffd849cfde53 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 17 Apr 2025 14:44:03 +0200 Subject: [PATCH 08/11] simpler explanation --- src/Utilities/distance_to_set.jl | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index a16b924b52..e3ad955a26 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -531,19 +531,16 @@ function distance_to_set( }, ) where {T<:Real} _check_dimension(x, set) - # Let - # 1) `U` be the matrix of eigenvectors. - # 2) `λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ))` - # 3) `A = LinearAlgebra.Symmetric(U * λ_negative * U')` - # The upper bound distance is the Frobenius norm of `A` which is the square - # root of the sum of the squares of the negative eigenvalues. - # Indeed: - # `⟨U * λ_negative * U', U * λ_negative * U'⟩` - # which is equal to - # `⟨λ_negative * U' * U, U' * U * λ_negative⟩` - # Since `U'U = I` this is equal to - # `⟨λ_negative, λ_negative⟩` - # So we need the return: + # We should return the norm of `A` defined by: + # ```julia + # λ, U = LinearAlgebra.eigvals(_reshape(x, set)) + # λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ)) + # A = LinearAlgebra.Symmetric(U * λ_negative * U') + # ``` + # The norm should correspond to `MOI.Utilities.set_dot` + # so it's the Frobenius norm. + # The Frobenius norm of `A` which is Euclidean + # norm of the vector of eigenvalues. return √sum(LinearAlgebra.eigvals(_reshape(x, set))) do λ min(λ, zero(T))^2 end From b7b88ccbe0be2ca5c56f880511bf0dad25ebb615 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 19 Apr 2025 08:38:26 +1200 Subject: [PATCH 09/11] Update distance_to_set.jl --- src/Utilities/distance_to_set.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index e3ad955a26..039cc28876 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -541,7 +541,6 @@ function distance_to_set( # so it's the Frobenius norm. # The Frobenius norm of `A` which is Euclidean # norm of the vector of eigenvalues. - return √sum(LinearAlgebra.eigvals(_reshape(x, set))) do λ - min(λ, zero(T))^2 - end + eigvals = LinearAlgebra.eigvals(_reshape(x, set)) + return sqrt(sum(min.(zero(T), eigvals).^2)) end From 109ab3dd0ca5928121160d15478233485aadfc40 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 19 Apr 2025 10:31:56 +1200 Subject: [PATCH 10/11] Update distance_to_set.jl --- src/Utilities/distance_to_set.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index 039cc28876..f770e69854 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -542,5 +542,6 @@ function distance_to_set( # The Frobenius norm of `A` which is Euclidean # norm of the vector of eigenvalues. eigvals = LinearAlgebra.eigvals(_reshape(x, set)) - return sqrt(sum(min.(zero(T), eigvals).^2)) + neg_eigvals = min.(zero(T), eigvals) + return LinearAlgebra.norm(neg_eigvals, 2) end From c4cb4b69404c974dcc6bf2502c13dabc0493ce74 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Sat, 19 Apr 2025 10:59:44 +1200 Subject: [PATCH 11/11] Update distance_to_set.jl --- src/Utilities/distance_to_set.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/Utilities/distance_to_set.jl b/src/Utilities/distance_to_set.jl index f770e69854..12655c4222 100644 --- a/src/Utilities/distance_to_set.jl +++ b/src/Utilities/distance_to_set.jl @@ -533,15 +533,14 @@ function distance_to_set( _check_dimension(x, set) # We should return the norm of `A` defined by: # ```julia - # λ, U = LinearAlgebra.eigvals(_reshape(x, set)) + # λ, U = LinearAlgebra.eigen(_reshape(x, set)) # λ_negative = LinearAlgebra.Diagonal(min.(zero(T), λ)) # A = LinearAlgebra.Symmetric(U * λ_negative * U') + # LinearAlgebra.norm(A, 2) # ``` - # The norm should correspond to `MOI.Utilities.set_dot` - # so it's the Frobenius norm. - # The Frobenius norm of `A` which is Euclidean - # norm of the vector of eigenvalues. + # The norm should correspond to `MOI.Utilities.set_dot` so it's the + # Frobenius norm, which is the Euclidean norm of the vector of eigenvalues. eigvals = LinearAlgebra.eigvals(_reshape(x, set)) - neg_eigvals = min.(zero(T), eigvals) - return LinearAlgebra.norm(neg_eigvals, 2) + eigvals .= min.(zero(T), eigvals) + return LinearAlgebra.norm(eigvals, 2) end