Unexpected array VariableRef[c[1] c[2] c[3]] in nonlinear expression

Specific Domains Optimization (Mathematical)
1

Code:

function b̂ₜfunc(
  x̂ₜ₊₁::AbstractVector,
  Σ̂ₜ₊₁::AbstractMatrix,
  γ::AbstractFloat,
  ξ::AbstractFloat,
  b̂ₜ::AbstractVector
)
  n_assets = length(x̂ₜ₊₁)
  model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
  @variable(model, c[1:n_assets])
  @constraint(model, sum(c)==0)
  @NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*(c')' + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs(item) for item in c))
  optimize!(model)
  return value.(c)
end

Reproducible example:

x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
gam = 0.02;
kasi = 0.03;
b = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, gam, kasi, b)

Throws:

ERROR: Unexpected array VariableRef[c[1] c[2] c[3]] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.

Other variants that I’ve tried:

@NLobjective(model, Min, sum((γ*c'*Σ̂ₜ₊₁)[item]*c[item] for item=1:n_assets) + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs(item) for item in c))

# And

@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*(c')' + sum(c[item]*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁)[item] for item=1:n_assets) + ξ*sum(abs(item) for item in c))

# And
@NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*c + only(c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁)) + ξ*sum(abs(item) for item in c))

Throws the same error.

When I try @NLobjective(model, Min, γ*c'*Σ̂ₜ₊₁*c) I get the same error. Hence, the problem is with the matrix multiplication. I do not understand. This is the same as this example. Why it does not work in my case? Got it. Because I used @NLobjective.

2

The problem was with @NLobjective. I replaced it with @objective. Also, there is no need to use only. Furthermore, sum(abs(item) for item in c) can be replaced with sum(abs.(c)). Final functions:

function b̂ₜfunc(
  x̂ₜ₊₁::AbstractVector,
  Σ̂ₜ₊₁::AbstractMatrix,
  γ::AbstractFloat,
  ξ::AbstractFloat,
  b̂ₜ::AbstractVector
)
  n_assets = length(x̂ₜ₊₁)
  model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
  @variable(model, c[1:n_assets])
  @constraint(model, sum(c)==0)
  @objective(model, Min, γ*c'*Σ̂ₜ₊₁*c + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(abs.(c)))
  optimize!(model)
  return value.(c)
end

x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
gam = 0.02;
kasi = 0.03;
b = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, gam, kasi, b)
# 3-element Vector{Float64}:
#   1.4843778267918864e20
#  -6.863558855303889e18
#  -1.4157422382388475e20

Wow :thinking: Such big values.

1 Like
3

You should always check the termination_status(model) before querying values. Did the problem solve correctly?

And yes, the @objective macro now supports nonlinear and it fixes a number of limitations from the @NL macros.

1 Like
4

If this can be implied by the function that you mentioned:

julia> termination_status(model)
NORM_LIMIT::TerminationStatusCode = 17

It’s not clear. I do not know what’s the meaning of NORM_LIMIT::TerminationStatusCode = 17. I guess it’s not converged. I don’t know.

Ah, Grrreat!! Thank you so much. :heart:

5

NORM_LIMIT means that the solver did not terminate with an optimal solution. See JuMP · JuMP

If you turn on printing, you’ll see Ipopt print:

EXIT: Iterates diverging; problem might be unbounded.

The issue is that Ipopt assumes problems are smooth and twice differentiable. abs(c) violates that assumption, so Ipopt is not guaranteed to find a solution.

One option is the replace the L1 penalty sum(abs.(c)) with the L2 sum(c.^2).

Another is to use a reformulation (see documentation Tips and tricks · JuMP):

using JuMP, Ipopt
function b̂ₜfunc(
    x̂ₜ₊₁::AbstractVector,
    Σ̂ₜ₊₁::AbstractMatrix,
    γ::AbstractFloat,
    ξ::AbstractFloat,
    b̂ₜ::AbstractVector
)
    n_assets = length(x̂ₜ₊₁)
    model = Model(Ipopt.Optimizer)
    @variable(model, c[1:n_assets])
    @constraint(model, sum(c) == 0)
    @variable(model, t[1:n_assets] >= 0)
    @constraint(model, t .>= c)
    @constraint(model, t .>= -c)
    @objective(model, Min, γ*c'*Σ̂ₜ₊₁*c + c'*((2γ*b̂ₜ'*Σ̂ₜ₊₁)'-x̂ₜ₊₁) + ξ*sum(t))
    optimize!(model)
    return value.(c)
end

x̂ₜ₊₁ = [0.9884, 1.02564, 1.01561];
Σ̂ₜ₊₁ = rand(3, 3);
γ = 0.02;
ξ = 0.03;
b̂ₜ = [0.5, 0.1, 0.4];
b̂ₜfunc(x̂ₜ₊₁, Σ̂ₜ₊₁, γ, ξ, b̂ₜ)
1 Like
6

Thank you so much. I shall rely much more on the provided documentation. Thanks.

Since I’m trying to implement a proposed algorithm, I should follow the same signature (Although I’m aware that it does not change the result of the model). However, I got your point on this. Thank you so much. :heart:

This is the most interesting part! However, it’s a little bit ambiguous for me. For example, you’ve set t[1:n_assets] >= 0, and later you noted t .>= -c. Isn’t this redundant? Or I’m missing something here. Oh, I got it. Because we do not know whether c_i would be a positive or negative value. Hence, we need both t .>= c and t .>= -c to ensure that t_i would result in the absolute value of c_i. Thank you so much! It’s a beautiful trick.

1 Like
7

Yes, technically the t >= 0 is redundant. But it doesn’t hurt to have it.

1 Like