Retaining the property of Symmetric?

General Usage
, , ,
1

It doesn’t keep the Symmetric property

import JuMP, Gurobi, LinearAlgebra
model = JuMP.Model(Gurobi.Optimizer)
JuMP.@variable(model, X[1:2, 1:2] >= 0, Symmetric)
JuMP.optimize!(model)
Xt = JuMP.value.(X) # ⚠️ Can this behavior be improved?

# julia> typeof(X)
# LinearAlgebra.Symmetric{JuMP.VariableRef, Matrix{JuMP.VariableRef}}

# julia> typeof(Xt)
# Matrix{Float64} (alias for Array{Float64, 2})

Xt = LinearAlgebra.Symmetric(Xt) # Since JuMP do not support, I have to do this myself

Here is another example

import JuMP, Gurobi, LinearAlgebra
model = JuMP.Model(Gurobi.Optimizer)
JuMP.@variable(model, x); JuMP.@variable(model, y)
trivial_value = 999
JuMP.@expression(model, X, LinearAlgebra.Symmetric([0 x; trivial_value y]))
JuMP.optimize!(model)
Xt = JuMP.value.(X) # ⚠️ Can this behavior be improved?

julia> typeof(X)
LinearAlgebra.Symmetric{JuMP.AffExpr, Matrix{JuMP.AffExpr}}

julia> typeof(Xt)
Matrix{Float64} (alias for Array{Float64, 2})
2

This is a property of Julia, not JuMP:

julia> using LinearAlgebra

julia> x = LinearAlgebra.Symmetric([1 2; 2 3])
2×2 Symmetric{Int64, Matrix{Int64}}:
 1  2
 2  3

julia> x.^2
2×2 Matrix{Int64}:
 1  4
 4  9
2 Likes
3

Related issue on redundant operations in elementwise operations: map and broadcast on symmetric matrices is inefficient · Issue #643 · JuliaLang/LinearAlgebra.jl. However, would it be breaking to go farther and change the return type to Symmetric, even if most people would welcome it?

1 Like
4

The current behavior is

julia> S
3×3 Symmetric{Int64, Matrix{Int64}}:
  0  -8  -2
 -8  -2   3
 -2   3   8

julia> S^2
3×3 Symmetric{Int64, Matrix{Int64}}:
  68  10  -40
  10  77   34
 -40  34   77

julia> S * S
3×3 Matrix{Int64}:
  68  10  -40
  10  77   34
 -40  34   77

julia> transpose(S) * S
3×3 Matrix{Int64}:
  68  10  -40
  10  77   34
 -40  34   77

julia> S .^ 2
3×3 Matrix{Int64}:
  0  64   4
 64   4   9
  4   9  64
5

These two probably need to stay the way they are. The product of two symmetric matrices is not necessarily symmetric. In your specific examples, the result is symmetric because the factors are identical, but that’s a property of the values, not their types. It’s generally desirable for the return type to only depend on the argument types, not the argument values.

For S.^2 it would be reasonable to return Symmetric, but changing that now might be considered too breaking, that is, not backward compatible.

2 Likes
6

I think it’s fine to change return types across versions, given that types are tightly coupled with performance. This is somewhat breaking if people are explicitly dispatching on a Matrix, but the return type is not documented to be stable across versions.

3 Likes
7

I have a novel idea: Why should we write S .^ 2?
A more thoughtful method could have been

julia> S
3×3 Symmetric{Int64, Matrix{Int64}}:
 -5  -1   3
 -1  -5   1
  3   1  -5

julia> U = UpperTriangular(S.data)
3×3 UpperTriangular{Int64, Matrix{Int64}}:
 -5  -1   3
  ⋅  -5   1
  ⋅   ⋅  -5

julia> U = U .^ 2
3×3 UpperTriangular{Int64, Matrix{Int64}}:
 25   1   9
  ⋅  25   1
  ⋅   ⋅  25

julia> ret = Symmetric(U.data)
3×3 Symmetric{Int64, Matrix{Int64}}:
 25   1   9
  1  25   1
  9   1  25

julia> ret == S .^ 2
true

I don’t know whether this is appropriate.

8

This is precisely how many algebraic operations are carried out for Symmetric, and seems appropriate once we specialize broadcasting for the type.

1 Like