How to define a lower triangular decision matrix

Is there an efficient way to define a lower triangular matrix (other than defining n by n matrix and adding a constraint that sets every upper triangular element to be zero).

I tried the following but it resulted in an error when I multiply it by another matrix B

using JuMP
using LinearAlgebra
model = Model()
n = 5
@variable(model, A[i=1:n, j=1:i])
@variable(model, B[1:n, 1:n], Symmetric)
C = A * B’

You were pretty close:

using JuMP
import LinearAlgebra
model = Model()
@variable(model, A[1:3,1:3], Symmetric)
Al = LinearAlgebra.LowerTriangular(A)

julia> @expression(model, Al * [1, 2, 3])
3-element Vector{AffExpr}:
 A[1,1]
 A[1,2] + 2 A[2,2]
 A[1,3] + 2 A[2,3] + 3 A[3,3]
1 Like

Great Thanks. I don’t see a need for A to be symmetric. Right?

If you define A to be symmetric, JuMP will only add 6 variables, because it knows that A[i, j] is the same as A[j, i]. If you omit Symmetric, it will add 9 variables.

2 Likes

@odow I am still struggling to use this lower triangular matrix in my model. Here is what I’m trying to set up. I tried wrapping them into expression but that didn’t help. Any idea?

using JuMP
using LinearAlgebra
model = Model()
n=5
@variable(model, A[1:n,1:n], Symmetric)
Al = LowerTriangular(A)
@variable(model, B[1:n,1:n], Symmetric)
@variable(model, c[1:n])
X = Al * B + c * c'
# @expression(model, X, Al * B + c * c')

@variable(model, Q[1:n,1:n])
@variable(model, P[1:n,1:n])
Y = Q*P'
# @expression(model, Y, Q*P)
@constraint(model, X .== Y)

What’s the problem?

Do you mean in the code above? It gives me the following error

ArgumentError: Cannot set a non-diagonal index in a symmetric matrix

Apparently It doesn’t like multiplying the lower triangular matrix Al by another matrix B. I also tried to put them in expressions like what you did in your answer but that didn’t help.

Ah. Okay the issue is LinearAlgebra matrix types and expressions · Issue #66 · jump-dev/MutableArithmetics.jl · GitHub.

You can work around the problem as Al = collect(LowerTriangular(A)). But I should take another look at what would be required to fix this properly.

The underlying problem is that LinearAlgebra makes some assumptions about JuMP variables that aren’t true. For example, VariableRef * VariableRef is a QuadExpr, whereas Float64 * Float64 is still a Float64.

1 Like