LinearAlgebra matrix types and expressions

[odow]

As reported on Discourse, affine expressions don't play well with LinearAlgebra.LowerTriangular:

julia> using JuMP

julia> model = Model()
A JuMP Model
Feasibility problem with:
Variables: 0
ariableModel mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.

julia> @variable(model, x[1:3])
3-element Array{VariableRef,1}:
 x[1]
 x[2]
 x[3]

julia> using LinearAlgebra

julia> A = LowerTriangular(rand(3, 3))
3×3 LowerTriangular{Float64,Array{Float64,2}}:
 0.246249   ⋅         ⋅      
 0.925955  0.871257   ⋅      
 0.388344  0.724686  0.435545

julia> b = x .+ rand(3)
3-element Array{GenericAffExpr{Float64,VariableRef},1}:
 x[1] + 0.41008732270098425
 x[2] + 0.6152593235089165 
 x[3] + 0.7871105854964822 

julia> c = rand(3)'
1×3 Adjoint{Float64,Array{Float64,1}}:
 0.0817674  0.517406  0.660573

julia> A * b
ERROR: MethodError: Cannot `convert` an object of type GenericQuadExpr{Float64,VariableRef} to an object of type GenericAffExpr{Float64,VariableRef}
Closest candidates are:
  convert(::Type{GenericAffExpr{T,V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:314
  convert(::Type{GenericAffExpr{T,V}}, ::V) where {T, V} at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:311
  convert(::Type{T}, ::GenericAffExpr{T,VarType} where VarType) where T at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:318
  ...
Stacktrace:
 [1] setindex! at ./array.jl:769 [inlined]
 [2] lmul!(::LowerTriangular{GenericAffExpr{Float64,VariableRef},Array{GenericAffExpr{Float64,VariableRef},2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:760
 [3] *(::LowerTriangular{Float64,Array{Float64,2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:1802
 [4] top-level scope at none:0

julia> @expression(model, A * b)
3-element Array{GenericAffExpr{Float64,VariableRef},1}:
 0.24624865480766212 x[1] + 0.10098345156879301                                                 
 0.9259545359165902 x[1] + 0.8712570245873634 x[2] + 0.9157712241268796                         
 0.3883444595770791 x[1] + 0.724685976949409 x[2] + 0.4355450311946607 x[3] + 0.9479470481617296

julia> @expression(model, (A * b))
3-element Array{GenericAffExpr{Float64,VariableRef},1}:
 0.24624865480766212 x[1] + 0.10098345156879301                                                 
 0.9259545359165902 x[1] + 0.8712570245873634 x[2] + 0.9157712241268796                         
 0.3883444595770791 x[1] + 0.724685976949409 x[2] + 0.4355450311946607 x[3] + 0.9479470481617296

julia> @expression(model, (A * b) .* c)
ERROR: MethodError: Cannot `convert` an object of type GenericQuadExpr{Float64,VariableRef} to an object of type GenericAffExpr{Float64,VariableRef}
Closest candidates are:
  convert(::Type{GenericAffExpr{T,V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:314
  convert(::Type{GenericAffExpr{T,V}}, ::V) where {T, V} at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:311
  convert(::Type{T}, ::GenericAffExpr{T,VarType} where VarType) where T at /Users/oscar/.julia/dev/JuMP/src/aff_expr.jl:318
  ...
Stacktrace:
 [1] setindex! at ./array.jl:769 [inlined]
 [2] lmul!(::LowerTriangular{GenericAffExpr{Float64,VariableRef},Array{GenericAffExpr{Float64,VariableRef},2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:760
 [3] *(::LowerTriangular{Float64,Array{Float64,2}}, ::Array{GenericAffExpr{Float64,VariableRef},1}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:1802
 [4] top-level scope at /Users/oscar/.julia/packages/MutableArithmetics/hS4h3/src/rewrite.jl:224
 [5] top-level scope at /Users/oscar/.julia/dev/JuMP/src/macros.jl:45

[blegat]

This issue is that in the product of a LowerTriangular of real with an array of variables, LinearAlgebra decides to promote both to AffExpr. The same issue is all over LinearAlgebra and SparseArrays. The short term solution is to overwrite these methods in
https://github.com/JuliaOpt/MutableArithmetics.jl/blob/master/src/dispatch.jl
The long term solution is to fix this in LinearAlgebra

[odow]

I should have added, the even simpler work-around is to force people to call collect(A) on a LowerTriangular matrix. Given the number of different matrix types in LinearAlgebra, I'm hesitant to have to add them all to MutableArithmetics.

Not really sure what a good answer is.

[odow]

Another case from discourse (https://discourse.julialang.org/t/constraint-seems-linear-but-the-compiler-determined-it-was-quadratic-not-sure-why-or-how-to-fix/84426/6) is

julia> using JuMP, LinearAlgebra

julia> model = Model();

julia> @variable(model, an[1:2], Bin);

julia> @variable(model, ACR[1:2, 1:2], Int);

julia> flows = ones(2, 2);

julia> @constraint(model, ACR .<= Diagonal(1 .- an) * flows)
ERROR: MethodError: Cannot `convert` an object of type QuadExpr to an object of type AffExpr
Closest candidates are:
  convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{T, V}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:521
  convert(::Type{GenericAffExpr{T, V}}, ::GenericAffExpr{S, V}) where {S, T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:528
  convert(::Type{GenericAffExpr{T, V}}, ::Union{Number, UniformScaling}) where {T, V} at /Users/oscar/.julia/packages/JuMP/Y4piv/src/aff_expr.jl:517
  ...
Stacktrace:
  [1] setindex!
    @ ./array.jl:845 [inlined]
  [2] setindex!
    @ ./multidimensional.jl:645 [inlined]
  [3] macro expansion
    @ ./broadcast.jl:984 [inlined]
  [4] macro expansion
    @ ./simdloop.jl:77 [inlined]
  [5] copyto!
    @ ./broadcast.jl:983 [inlined]
  [6] copyto!
    @ ./broadcast.jl:936 [inlined]
  [7] materialize!
    @ ./broadcast.jl:894 [inlined]
  [8] materialize!
    @ ./broadcast.jl:891 [inlined]
  [9] lmul!(D::Diagonal{AffExpr, Vector{AffExpr}}, B::Matrix{AffExpr})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:212
 [10] *(D::Diagonal{AffExpr, Vector{AffExpr}}, A::Matrix{Float64})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/diagonal.jl:201
 [11] macro expansion
    @ ~/.julia/packages/MutableArithmetics/Lnlkl/src/rewrite.jl:294 [inlined]
 [12] macro expansion
    @ ~/.julia/packages/JuMP/Y4piv/src/macros.jl:823 [inlined]
 [13] top-level scope
    @ REPL[20]:1

[odow]

So I don't know if we really need to support things like LowerTriangular(A) * x in general.

A solution is for the user to use MA.operate(*, LowerTriangular(A), x) outside the macros and to use a macro with LowerTriangular(A) * x. The overloads in https://github.com/jump-dev/MutableArithmetics.jl/blob/master/src/dispatch.jl are pretty scary, and rely on intercepting the eltype of various matrices. That's going to be impossible to get right in the long-run.

The macros do work in almost all circumstances, especially now that #170 doesn't do additional rewriting.

[blegat]

Yes, the overload would have been easier if LinearAlgebra was implementing their default dispatch methods like we do in MOI and MA with foo_fallback so thatfoo does not have too many methods and is easy to overload

[odow]

Even Base gets some of these wrong: #176 (comment)