Does that mean I can’t define matrix product at all?
Correct. In the current version of JuMP, you cannot create expressions like X * Y * Z
where X
, Y
and Z
are arrays of JuMP variables.
One work-around is:
using JuMP, LinearAlgebra
n = 5
model = Model()
@variable(model, R[1:n, 1:n], Symmetric)
@variable(model, Q[1:n, 1:n], Symmetric)
@variable(model, RQ[1:n, 1:n], Symmetric)
@constraint(model, RQ .== R * Q)
@expression(model, P, Q' * QR)
This limitation is fixed in the work-in-progress nonlinear rewrite: https://github.com/jump-dev/JuMP.jl/pull/3106.