How to multiply matrices with different structure?

Specific Domains Optimization (Mathematical)
1

I am defining a variable X as a lower trianglar matrix (to reduce the number of variables). How can I multiply it with another dense matrix Y? The following gives me an error

using JuMP
using LinearAlgebra

model = Model()
n = 4
@variable(model, X[i=1:n, j=1:i])  
@variable(model, Y[1:n, 1:n], Symmetric)

Z = X*Y
2

A manual way to do it would be

Z = [ sum( X[i,j] * Y[j,k] for j in 1:i) for i in 1:n, k in 1:n  ]

Alternatively, you could also define

X_ = [ (i>=j) ? X[i,j] : 0 for i in 1:n, j in 1:n ]
Z = X_ * Y    
# or Z = Matrix{QuadExpr}( X_ * Y ) if the eltype matters.

Of course, the most elegant way would be to use LowerTriagonal. But this seems
impossible since the LinearAlgebra package’s multiplication functions promote to a common type for the output before they perform the multiplication. This is a problem in this setup, since *(::VariableRef, ::VariableRef)::QuadExpr is defined, but *(::QuadExpr, ::QuadExpr)::QuadExpr is not (since that would lead to nonlinear terms.)

1 Like
3

Note that X isn’t really a matrix. It’s a JuMP.Containers.SparseAxisArray which in some cases acts like a matrix, but is more like a dictionary with keys.

One way to get a lower triangular matrix is:

@variable(model, X[1:n, 1:n], Symmetric)
Xl = LinearAlgebra.LowerTriangular(X)

But you’ll still run into the problem that @SteffenPL reports where JuMP cannot always multiply two matrices outside of the macros.

Here’s an example:

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

julia> n = 3
3

julia> @variable(model, X[1:n, 1:n], Symmetric)
3×3 Symmetric{VariableRef, Matrix{VariableRef}}:
 X[1,1]  X[1,2]  X[1,3]
 X[1,2]  X[2,2]  X[2,3]
 X[1,3]  X[2,3]  X[3,3]

julia> Xl = LinearAlgebra.LowerTriangular(X)
3×3 LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}}:
 X[1,1]  ⋅       ⋅
 X[1,2]  X[2,2]  ⋅
 X[1,3]  X[2,3]  X[3,3]

julia> X * Xl
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
  [1] setindex!(A::Symmetric{QuadExpr, Matrix{QuadExpr}}, v::VariableRef, i::Int64, j::Int64)
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:225
  [2] setindex!
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:257 [inlined]
  [3] _setindex!
    @ ./abstractarray.jl:1297 [inlined]
  [4] setindex!
    @ ./abstractarray.jl:1267 [inlined]
  [5] copyto_unaliased!(deststyle::IndexCartesian, dest::LowerTriangular{QuadExpr, Symmetric{QuadExpr, Matrix{QuadExpr}}}, srcstyle::IndexCartesian, src::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./abstractarray.jl:984
  [6] copyto!
    @ ./abstractarray.jl:950 [inlined]
  [7] copyto_axcheck!(dest::LowerTriangular{QuadExpr, Symmetric{QuadExpr, Matrix{QuadExpr}}}, src::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./abstractarray.jl:1056
  [8] AbstractMatrix{QuadExpr}(A::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./array.jl:541
  [9] AbstractArray
    @ ./boot.jl:475 [inlined]
 [10] convert
    @ ./abstractarray.jl:15 [inlined]
 [11] *(A::Symmetric{VariableRef, Matrix{VariableRef}}, B::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1699
 [12] top-level scope
    @ REPL[12]:1

julia> Xl * X
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
  [1] setindex!(A::Symmetric{QuadExpr, Matrix{QuadExpr}}, v::VariableRef, i::Int64, j::Int64)
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/symmetric.jl:225
  [2] setindex!
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:257 [inlined]
  [3] _setindex!
    @ ./abstractarray.jl:1297 [inlined]
  [4] setindex!
    @ ./abstractarray.jl:1267 [inlined]
  [5] copyto_unaliased!(deststyle::IndexCartesian, dest::LowerTriangular{QuadExpr, Symmetric{QuadExpr, Matrix{QuadExpr}}}, srcstyle::IndexCartesian, src::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./abstractarray.jl:984
  [6] copyto!
    @ ./abstractarray.jl:950 [inlined]
  [7] copyto_axcheck!(dest::LowerTriangular{QuadExpr, Symmetric{QuadExpr, Matrix{QuadExpr}}}, src::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./abstractarray.jl:1056
  [8] AbstractMatrix{QuadExpr}(A::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}})
    @ Base ./array.jl:541
  [9] AbstractArray
    @ ./boot.jl:475 [inlined]
 [10] convert
    @ ./abstractarray.jl:15 [inlined]
 [11] *(A::LowerTriangular{VariableRef, Symmetric{VariableRef, Matrix{VariableRef}}}, B::Symmetric{VariableRef, Matrix{VariableRef}})
    @ LinearAlgebra /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/LinearAlgebra/src/triangular.jl:1657
 [12] top-level scope
    @ REPL[13]:1

julia> @expression(model, Xl * X)
3×3 Matrix{QuadExpr}:
 X[1,1]²                                        …  X[1,1]*X[1,3]
 X[1,2]*X[1,1] + X[2,2]*X[1,2]                     X[1,2]*X[1,3] + X[2,2]*X[2,3]
 X[1,3]*X[1,1] + X[2,3]*X[1,2] + X[3,3]*X[1,3]     X[1,3]² + X[2,3]² + X[3,3]²
2 Likes