# How to multiply matrices with different structure?

**URL:** <https://discourse.julialang.org/t/how-to-multiply-matrices-with-different-structure/86428>\
**Category:** Optimization (Mathematical)\
**Created:** [August 27, 2022, 12:22pm UTC](https://discourse.julialang.org/t/how-to-multiply-matrices-with-different-structure/86428 "2022-08-27T12:22:32Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ivzktch99](https://avatars.discourse-cdn.com/v4/letter/i/b19c9b/32.png) [@ivzktch99](https://discourse.julialang.org/u/ivzktch99)\
**Post date:** [August 27, 2022, 12:22pm UTC](https://discourse.julialang.org/t/how-to-multiply-matrices-with-different-structure/86428/1 "2022-08-27T12:22:32Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![SteffenPL](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/steffenpl/32/206270_2.png) [@SteffenPL](https://discourse.julialang.org/u/SteffenPL)\
**Post date:** [August 27, 2022, 8:25pm UTC](https://discourse.julialang.org/t/how-to-multiply-matrices-with-different-structure/86428/2 "2022-08-27T20:25:37Z")

</div>

A manual way to do it would be

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

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 28, 2022, 9:42pm UTC](https://discourse.julialang.org/t/how-to-multiply-matrices-with-different-structure/86428/3 "2022-08-28T21:42:10Z")

</div>

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:

```nohighlight
@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:

```nohighlight
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]²

```