Question on programmatically defining a nonlinear expression involving product and sum

Specific Domains Optimization (Mathematical)
1

I have a n-dimensional vector m, a set of JuMP variables defined by

\{x_{i,k_{i}}\}_{i\in\{1,2,\ldots,n\},k_{i}\in\{1,2,\ldots,m_{i}\}}=\{x_{1,1},x_{1,2},\ldots,x_{1,m_{1}},x_{2,1},x_{2,2},\ldots,x_{2,m_{2}},\ldots,x_{n,1},x_{n,2},\ldots x_{n,m_{n}}\}

(where m_i means m[i]), a set of n multidimensional array defined by \{A[i]\}_{i\in\{1,2,\ldots,n\}} with each A[i]\in[0,1]^{m_{1}\times m_{2}\times\ldots\times m_{n}}. I want to define the following n nonlinear expressions in JuMP

\sum_{k_{1}\in[1:m_{1}]}\cdots\sum_{k_{n}\in[1:m_{n}]}\left(A[i]\right)[k_{1},k_{2},\ldots,k_{n}]\left(x_{1,k_{1}}\times x_{2,k_{2}}\times\ldots\times x_{n,k_{n}}\right), \quad i=1,2,\ldots,n \quad \text{(EXPR)}(i)

I am defining the set of variables in JuMP as follows:

# Sample data 

using JuMP

n = 3

m = [3, 1, 2]

function generate_A(n::Int, m::Vector{Int})

    dims = Tuple(m)

    # A[i] has dimensions m[1] x m[2] x ... x m[n].
    A = Vector{Array{Float64}}(undef, n)
    for i in 1:n
        A[i] = rand(dims...) 
    end

    return A
    
end

A = generate_A(n, m) # so A[i][k_1, k_2, ..., k_m] gives $\sum_{k_{1}\in[1:m_{1}]}\cdots\sum_{k_{n}\in[1:m_{n}]}A[i][k_{1},k_{2},\ldots,k_{n}]\left(x_{1,k_{1}}\times x_{2,k_{2}}\times\ldots\times x_{n,k_{n}}\right)$

# define JuMP model

myModel = Model(Gurobi.Optimizer)

# define {x_{i,k_i}}_{i ∈ 1:n, k_i ∈ 1:m[i]}

@variable(myModel, x[i=1:n, k_i=1:m[i]])

When n and entries of m are small, I can define the expressions (\text{EXPR})(i) for i \in \{1,\ldots,n\} in JuMP manually without issue. However, it gets cumbersome for larger values of n or when m has large entries.

What is the best way to define the expressions (\text{EXPR})(i) for i\in \{1,\ldots,n\} in JuMP programmatically, i.e., one inputs n, m, A and the goal is to generate the nonlinear expressions (\text{EXPR})(i) for i\in \{1,\ldots,n\} in JuMP? Any tips/suggestions will be much appreciated!

2

You could do something like this:

using JuMP
m, n = [3, 1, 2], 3
generate_A(n::Int, m::Vector{Int}) = [rand(m...) for i in 1:n]
A = generate_A(n, m)
model = Model()
@variable(model, x[i in 1:n, k_i in 1:m[i]])
@expression(
    model,
    [i in 1:n],
    sum(
        A[i][k...] * prod(x[j,k[j]] for j in 1:n)
        for k in Iterators.product(map(i->1:i, m)...)
    )
)
2 Likes
3

Hi @odow, thanks for the quick response! Sorry for the confusion, I am actually trying to define n total nonlinear expressions, with i in A[i] corresponding to the index i \in \{1,2,\ldots,n \}, i.e., I am trying to define the nonlinear expressions:

\sum_{k_{1}\in[1:m_{1}]}\cdots\sum_{k_{n}\in[1:m_{n}]}\left(A[i]\right)[k_{1},k_{2},\ldots,k_{n}]\left(x_{1,k_{1}}\times x_{2,k_{2}}\times\ldots\times x_{n,k_{n}}\right)

for i = 1,2,\ldots,n.

I have also corrected my original question accordingly. Please let me know your thoughts.

4

I think my code example above does that? I do’t see anywhere else that i appears.

1 Like
5

You are correct, your code already accounts for that. Perhaps remove the question What is i in A[i] in your first post? So I can make it a solution to the question, because this will come handy to other people with a similar question.

1 Like
6

@odow, another side question, if I want to drop the sum over k_{i} and product involving x_{i,k_{i}} from the expression, i.e., if I want to model the nonlinear expression

\sum_{k_{1}=1}^{m_{1}}\sum_{k_{2}=1}^{m_{2}}\ldots\sum_{k_{i-1}=1}^{m_{i-1}}\sum_{k_{i+1}=1}^{m_{i+1}}\ldots\sum_{k_{n}=1}^{m_{n}}(A[i])[k_{1},k_{2},\ldots,k_{n}]\prod_{j=1,j\neq i}^{n}x_{j,k_{j}}

is it possible to model that in JuMP?

7

What is k_i in your A[i](k...) expression?

is it possible to model that in JuMP ?

Pretty much any algebraic expression is possible.

Perhaps something like:

function iterators(m, i)
    ret = [1:m_i for m_i in m]
    ret[i] = 1:1 # ???
    return ret
end
@expression(
    model,
    [i in 1:n],
    sum(
        A[i][k...] * prod(x[j,k[j]] for j in 1:n if j != i)
        for k in Iterators.product(iterators(m, i)...)
    )
)
1 Like
8

Hi @odow, thanks for your response! The term k_{i} in A[i](k_{1},k_{2},\ldots,k_{i-1},k_{i},k_{i+1},\ldots,k_{n}) corresponds to any integer in \{1,2,\ldots,m_{i}\}. Recall that in the initial sum k_{i}\in\{1,2,\ldots,m_{i}\} and in this version k_{i} wil take one of those m_{i} values. In other words, in this case, we are modeling m_{i}\times n nonlinear expressions:

\sum_{k_{1}=1}^{m_{1}}\sum_{k_{2}=1}^{m_{2}}\ldots\sum_{k_{i-1}=1}^{m_{i-1}}\sum_{k_{i+1}=1}^{m_{i+1}}\ldots\sum_{k_{n}=1}^{m_{n}}(A[i])[k_{1},k_{2},\ldots,k_{n}]\prod_{j=1,j\neq i}^{n}x_{j,k_{j}},

for i=1,2,\ldots,n and k_{i}=1,2,\ldots,m_{i}. I think your code can be modified as follows to fit this as follows:

function iterators(m, i, k_i)
    ret = [1:m_i for m_i in m]
    ret[i] = k_i:k_i
    return ret
end
@expression(
    model,
    [i in 1:n],
    sum(
        A[i][k...] * prod(x[j,k[j]] for j in 1:n if j != i)
        for k in Iterators.product(iterators(m, i, k_i)...)
    )
)

Please let me know if the modification is correct. Thanks again @odow!

9

It’s not correct. Your k_i is not defined. I think you need:

function iterators(m, i, k_i)
    ret = [1:m_i for m_i in m]
    ret[i] = k_i:k_i
    return ret
end
@expression(
    model,
    [i in 1:n, k_i in 1:m[i]],
    sum(
        A[i][k...] * prod(x[j,k[j]] for j in 1:n if j != i)
        for k in Iterators.product(iterators(m, i, k_i)...)
    )
)
1 Like
10

Yes, I ran into this issue while running the code, so was able to fix it basically the same way suggested by you earlier today. Thanks very much @odow !

1 Like