How to use max function within integral in the objective function

[christineymshen]

I'm trying to set up an optimization where my objective function is in the following form:

$$\text{min}_{\lambda, z} \ \int \left(q(x)-z\int p(x,\theta)\lambda(\theta)d\theta\right)^+ dx + \alpha z.$$

$$\text{s.t.} \ z \ge 0, \lambda(\theta) \ge 0, \int \lambda(\theta) d\theta =1.$$

i.e., I want to find $z \in \mathbb{R}$ and a function $\lambda(\theta)$ which minimizes the above function. I'm new to Julia, and I'm struggling to set this up. In particular, I'd like to ask:

1. is it a problem that InfiniteOpt package can solve?
2. how to implement the flooring at zero within the outer integral?

Any thoughts/ comments/ help would be greatly appreciated!!

The codes I've got so far as well as the error message are posted below. I got the error msg while executing the line for objective function. The long formula is just to the exact form of function p and q listed above. The errors are enclosed after the codes.

using InfiniteOpt, Distributions, Ipopt 

tau2 = 1;
mu = 1;
alpha = 0.05;

model = InfiniteModel(Ipopt.Optimizer)

@finite_parameter(model, z == 1)
@infinite_parameter(model, x in [-10000, 10000], num_supports=100, 
derivative_method = OrthogonalCollocation(2))

@infinite_parameter(model, delta in [-10000, 0], num_supports=100, 
derivative_method = OrthogonalCollocation(2))

@variable(model, lambda>=0, Infinite(delta))

@objective(model, Min, alpha * z + integral( max(0, 
    (4*pi*(1+tau2))^(-0.5) * exp(-(x-mu)^2/(4*(1+tau2))) - 
    z * integral( (4*pi)^(-0.5)*exp(-(x-delta)^2/2)*lambda, delta) ), x))

@constraint(model, z >= 0)

@constraint(model, lambda(delta) >= 0)

@constraint(model, expect(lambda(delta), delta) = 1)

optimize!(model)

ERROR: MethodError: no method matching isless(::NLPExpr, ::Int64)
Closest candidates are:
  isless(::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at ~/.julia/packages/StatsBase/XgjIN/src/statmodels.jl:90
  isless(::ForwardDiff.Dual{Tx}, ::Integer) where Tx at ~/.julia/packages/ForwardDiff/QdStj/src/dual.jl:144
  isless(::ForwardDiff.Dual{Tx}, ::Real) where Tx at ~/.julia/packages/ForwardDiff/QdStj/src/dual.jl:144
  ...
Stacktrace:
 [1] max(x::Int64, y::NLPExpr)
   @ Base ./operators.jl:480
 [2] macro expansion
   @ ~/.julia/packages/MutableArithmetics/9dpep/src/rewrite.jl:322 [inlined]
 [3] macro expansion
   @ ~/.julia/packages/JuMP/yYfHy/src/macros.jl:1293 [inlined]
 [4] top-level scope
   @ REPL[13]:1

[pulsipher]

Hi @christineymshen and welcome.

Yes, we can set this up in InfiniteOpt without too much issue. First, we can reformulate using a common linear programming trick shown here to address the max function.

In InfiniteOpt, I get:

using InfiniteOpt, Ipopt

# Initialize the model
model = InfiniteModel(Ipopt.Optimizer)

# Define the infinite parameters
@infinite_parameter(model, x ∈ [-10000, 10000], num_supports  = 100)
@infinite_parameter(model, θ ∈ [-10000, 0], num_supports = 100)

# Define the decision variables and their bounds
@variable(model, λ ≥ 0, Infinite(θ)) # creates λ(θ)
@variable(model, z ≥ 0)  # creates z
@variable(model, aux ≥ 0,  Infinite(x)) # creates auxiliary variable for handling the max operation

# Define the known functions
q = # put function here
p = # put function here

# Define the objective function
α = 0.05
@objective(model, Min, ∫(aux, x) + α * z)
@constraint(model, aux ≥ q - z * ∫(p * λ, θ))

# Add the constraints
@constraint(model, ∫(λ, θ) == 1)

# Optimize and get the results
optimize!(model)
if has_values(model)
   opt_z = value(z)
   opt_λ = value(λ)
end

I should note that the use of x seems unnecessary here since it is not used by the decision variables. Depending on the choices of q and p you can probably simplify the expression by integrating over x analytically first.

The code you provided had a number of errors, I would recommend working through a few of the tutorials provided by JuMP.jl here and then review the InfiniteOpt tutorials. These will help you get a much better feel for the syntax and workflows.

[christineymshen]

Thanks a lot for the explanations! I think I'll go with the Monte Carlo approach. It'd be great if you could walk me through how to do it in InfiniteOpt!!

[christineymshen]

Um I just tried to do a 2-dimension run with similar setup. However, the results are weird. I wonder whether you might have time to have a quick look at my codes below to see whether anything's obviously wrong? In particular, one thing I don't understand is ... say my parameter $\theta$ is 2-dimensional. And for each dimension, I've used 59 support points. The decision variable $\lambda$ is a function of $\theta$. But the results I got for $\lambda$ has a length of $4761 = 69^2$. I wonder whether I'm missing something here ...

Many thanks again!!

using InfiniteOpt, Ipopt
using JLD2

# Initialize the model
model = InfiniteModel(Ipopt.Optimizer)

θ_supp = [-200;-100;-50;-10:1:-6;-5:0.1:0;]
x_supp = [-200;-100;-50;-8:0.2:8;50;100;200;]

# Define the infinite parameters
@infinite_parameter(model, x[1:2] ∈ [-200, 200],independent=true, supports = x_supp)
@infinite_parameter(model, θ[1:2] ∈ [-200, 0], independent=true, supports = θ_supp)

# Define the decision variables and their bounds
@variable(model, λ ≥ 0, Infinite(θ)) # creates λ(θ)
@variable(model, z ≥ 0)  # creates z
@variable(model, aux ≥ 0,  Infinite(x)) # creates auxiliary variable for handling the max operation

τ2 = 1
μ = [2;2]
α = 0.05

# Define the known functions
q = (4π*(1+τ2))^(-1) * exp(-sum((x-μ) .* (x-μ))/4/(1+τ2))
p = (4π)^(-1) * exp(-sum((x-θ) .* (x-θ))/4)

# Define the objective function
@objective(model, Min, ∫(aux, x)+ α*z )
@constraint(model, aux ≥ q - z * ∫(p * λ, θ))

# Add the constraints
@constraint(model, ∫(λ, θ) == 1)

# Optimize and get the results
optimize!(model)
if has_values(model)
   opt_z = value(z)
   opt_λ = value(λ)
end

# save results
jldsave(string("out/res_k3_model.jld2"); model=model)
jldsave(string("out/res_k3_res.jld2"); μ,τ2,opt_z,opt_λ,θ_supp,x_supp)

[pulsipher]

This starts to get into some interesting behavior with multi-dimensional integrals. It would be worth reading https://infiniteopt.github.io/InfiniteOpt.jl/stable/guide/measure/ to about working with these. A key takeaway is that currently, InfiniteOpt only supports MC integration for multi-dimensional integrals and the integral syntax only allows for independent interval domains.

In the above code, the integrals will generate 10 uniform MC samples (the default amount) for the integrals and will ignore the points you provide since they are not known to be samples from a uniform distribution. Hence, the discretization will then have additional points for the discretization (going from 59 to 69). There are a few ways we can deal with this.

Option 1: Only use the generated MC points
Here will let InfiniteOpt generate the points. This only will work for square independent domains currently.

using InfiniteOpt, Ipopt

# Initialize the model
model = InfiniteModel(Ipopt.Optimizer)

# Define the infinite parameters
@infinite_parameter(model, x[1:2] ∈ [-200, 200], independent=false)
@infinite_parameter(model, θ[1:2] ∈ [-200, 0], independent=false)

# Define the decision variables and their bounds
@variable(model, λ ≥ 0, Infinite(θ)) # creates λ(θ)
@variable(model, z ≥ 0)  # creates z
@variable(model, aux ≥ 0,  Infinite(x)) # creates auxiliary variable for handling the max operation

τ2 = 1
μ = [2;2]
α = 0.05

# Define the known functions
q = (4π*(1+τ2))^(-1) * exp(-sum((x-μ) .* (x-μ))/4/(1+τ2))
p = (4π)^(-1) * exp(-sum((x-θ) .* (x-θ))/4)

# Define the objective function
set_multi_integral_defaults(num_supports = 50) # use 50 MC supports for each integral
@objective(model, Min, ∫(aux, x)+ α*z )
@constraint(model, aux ≥ q - z * ∫(p * λ, θ))

# Add the constraints
@constraint(model, ∫(λ, θ) == 1)

# Optimize and get the results
optimize!(model)
if has_values(model)
   opt_z = value(z)
   opt_λ = value(λ)
end

Option 2: Use nested 1D integrals (the combinatorics here make this quickly intractable)
We can simply nest our independent dimensions of the integral such that the default trapezoid rule can be used.

using InfiniteOpt, Ipopt

# Initialize the model
model = InfiniteModel(Ipopt.Optimizer)

θ_supp = [-200;-100;-50;-10:1:-6;-5:0.1:0;]
x_supp = [-200;-100;-50;-8:0.2:8;50;100;200;]

# Define the infinite parameters
@infinite_parameter(model, x[1:2] ∈ [-200, 200],independent=true, supports = x_supp)
@infinite_parameter(model, θ[1:2] ∈ [-200, 0], independent=true, supports = θ_supp)

# Define the decision variables and their bounds
@variable(model, λ ≥ 0, Infinite(θ)) # creates λ(θ)
@variable(model, z ≥ 0)  # creates z
@variable(model, aux ≥ 0,  Infinite(x)) # creates auxiliary variable for handling the max operation

τ2 = 1
μ = [2;2]
α = 0.05

# Define the known functions
q = (4π*(1+τ2))^(-1) * exp(-sum((x-μ) .* (x-μ))/4/(1+τ2))
p = (4π)^(-1) * exp(-sum((x-θ) .* (x-θ))/4)

# Define the objective function
@objective(model, Min, ∫(∫(aux, x[1]), x[2])+ α*z )
@constraint(model, aux ≥ q - z * ∫(∫(p * λ, θ[1]), θ[2]))

# Add the constraints
@constraint(model, ∫(∫(λ, θ[1]), θ[2]) == 1)

# Optimize and get the results
optimize!(model)
if has_values(model)
   opt_z = value(z)
   opt_λ = value(λ)
end

Option 3: Manually handle the MC integration (enables arbitrary domains of integration)
We can directly enforce a MC integration rule by manually controlling what the supports do. Here is what the model would look like using a uniform MC integration scheme. This will enable you to consider nonsquare domains (e.g., halfspaces).

using InfiniteOpt, Ipopt

# Initialize the model
model = InfiniteModel(Ipopt.Optimizer)

num_samples = # TODO put the desired number of MC samples here
θ_supp = # TODO generate 2D uniform random samples from your domain here using num_samples
x_supp = # TODO generate 2D uniform random samples from your domain here using num_samples
V = # TODO put the volume of the domain here

# Define the infinite parameters
@infinite_parameter(model, x[1:2] ∈ [-200, 200], independent=false, supports = x_supp)
@infinite_parameter(model, θ[1:2] ∈ [-200, 0], independent=false, supports = θ_supp)

# Define the decision variables and their bounds
@variable(model, λ ≥ 0, Infinite(θ)) # creates λ(θ)
@variable(model, z ≥ 0)  # creates z
@variable(model, aux ≥ 0,  Infinite(x)) # creates auxiliary variable for handling the max operation

τ2 = 1
μ = [2;2]
α = 0.05

# Define the known functions
q = (4π*(1+τ2))^(-1) * exp(-sum((x-μ) .* (x-μ))/4/(1+τ2))
p = (4π)^(-1) * exp(-sum((x-θ) .* (x-θ))/4)

# Define the objective function
@objective(model, Min,  V / num_samples * support_sum(aux, x)+ α*z )
@constraint(model, aux ≥ q - z * V / num_samples * support_sum(p * λ, θ))

# Add the constraints
@constraint(model, V / num_samples * support_sum(λ, θ) == 1)

# Optimize and get the results
optimize!(model)
if has_values(model)
   opt_z = value(z)
   opt_λ = value(λ)
end

You might get better results using importance sampling.

Option 4 Extend InfiniteOpt to do what you want
InfiniteOpt is extendable, so we can define methods to have it do what you want without hacking so much. The easiest way to handle this would be to extend the measure evaluation techniques: https://infiniteopt.github.io/InfiniteOpt.jl/stable/develop/extensions/#Measure-Evaluation-Techniques. And the most complete way would be to define a new domain type (e.g., for half domains) and then extend the measures to handle it: https://infiniteopt.github.io/InfiniteOpt.jl/stable/develop/extensions/#Infinite-Domains.

Option 5: Discretize manually and use JuMP.jl (most flexible approach)
Instead of using InfiniteOpt, we can manually discretize the problem and pose it directly in JuMP. This is the most flexible approach, but requires you to do all the bookkeeping. JuMP has tutorials/examples to help you get started. Eventually, you could even translate this work to PLASMO to start getting into decomposition approaches to speed things up (hopefully).

Final thoughts
The above approaches each have their pros and cons. Moreover, they each will have different accuracy.

Admittedly, our current multi-dimensional integral support has room for improvement and many of the current limitations will be addressed with #258.

[christineymshen]

Thanks so much for your detailed explanations!! This is very helpful. I'm starting with option 1 and 2, as for now I want to get a initial understanding of my results in 2 dimension first before diving in deeper, and these two seems to require minimal changes to what I already have. Though I can see they are the least scalable ... I've tried a little bit and have a few questions to follow up:

1. Not sure why, but option 1 doesn't seem to work for me. I tried a few times, every time the process is killed automatically after running for, say 5-10 mins. I'm using a Macbook, and running it using VScode.
2. For option 2, I have to significantly cut down the number of support points I use in order for the program to successfully run. E.g., I was using about 35 points for $\theta$ and 55 point for $x$. I had it run overnight but it was stuck in the initialization phase and wasn't able to start the optimization (i.e., I only got a usual warning msg which pops up at the beginning of every run, but didn't get any other printouts regarding the optimization process). So I had to kill it manually. Now I'm running with just 12 points for $\theta$ and 10 points for $x$, and able to get some results.
3. I wonder whether you have any suggestions in terms of how I make use of the optimization results. The results I get back from the optimizer are the vector opt_\lambda and a real number opt_z. To make use of these results, I need to be able to calculate $p(x)$ which requires integrating over $\lambda(\theta)$. I was calculating a Riemann sum manually when it was just 1-dimensional. But now that it's 2-dimensional, I wonder whether there are smarter ways to do this. For example, here are the super naive codes I'm using currently to check that $\lambda(\theta)$ indeed integrates to 1.

tmp = reshape(opt_λ, (12,12))
values = (tmp[1:(end-1), 1:(end-1)] + tmp[2:end, 1:(end-1)] +
tmp[1:(end-1), 2:end] + tmp[2:end, 2:end])/4

tmp = θ_supp[2:end]-θ_supp[1:(end-1)]
areas = [[i*j] for i in tmp, j in tmp]

sum(values .* areas)

4. Related to the last question. If I were to use option 3, I wonder how I could do integration myself after I get the optimization results ...

Again thanks a lot!!

[pulsipher]

Happy to help. Before I dive into your questions, I should reemphasize that high-dimensional integrals are computationally very challenging. With that in mind, I'll answer your questions.

1. I had a typo before with independent = true which should have been independent=false. I have edited my response and you'll find that this now computes rapidly without issue. The previous code did a lot of combinatorics (leading to 50^4 supports) which is why it gave the warning and then took a very long time or simply crashed by using too much memory. More on this below.
2. Nesting quadrature schemes is known to not be a scalable approach for high-dimensional integral evaluation. In this case, the warning correctly indicates this leads to over 28 million support points. Part of the problem here (for building the model) is our temporary nonlinear expression interface patch is not very efficient for evaluating integrals whose approximations have millions of terms. This should be resolved by [Nonlinear] add new nonlinear interface jump-dev/JuMP.jl#3106; however, even with it built, I am not hopeful Ipopt would solve it.
3. Most quadrature schemes (like Riemann sums or Gauss quadrature) generally don't handle multi-dimensional integrals well since they don't scale well. However, some colleagues of mine have told me that there are some new approaches that look promising, but I haven't looked into it yet. The easiest approach would be to use Monte Carlo integration which corresponds to option 3 above.
4. I highly recommend you use option 3. You basically already are doing the approximation rule manually with my code above. Thus, doing a post calculation integration via MC integration should be straightforward. You would just use the same sample points and follow the approximation scheme described here.