# Registering the user-defined function "multivariate normal cdf" for nonlinear modeling in JuMP

Hello Julia-community,

I would like to use the multivariate normal cumulative distribution function within a constraint in JuMP.

For the function, I’m using the `mvnormcdf` from this package MvNormalCDF.jl. Since this requires nonlinear modeling, I have to define the gradient and register the function as per the documentation here.

I have the mathematical formulations on how to define the gradient but I’m having issues to write this in Julia/JuMP. The point where I’m stuck now, is that the 'mvnormcdf` returns two values (probability and error), for my nonlinear constraint however I only need the probability value.

Here is a simpflied version of my code:

``````using JuMP, LinearAlgebra
using Distributions, NLopt
using Random
using MvNormalCDF

begin
ex2 = Model(NLopt.Optimizer)
set_optimizer_attribute(ex2, "algorithm", :LD_MMA)

@variable(ex2, x[1:2] >= 0)

@objective(ex2, Min, x + x)

set_start_value(x, 5)
set_start_value(x, 5)

p = 0.9 # probability

μ = [5; 5] # mean
vari = [4; 4] # variance
ρ = 0.5 # correlation
Σ = [vari ρ; ρ vari] # covariance matrix
a = [-Inf; -Inf]

norm_cdf(x) = cdf(Normal(), x)
norm_pdf(x) = exp(-(x^2) / 2) / sqrt(2 * pi)
register(ex2, :cdf, 1, norm_cdf, norm_pdf; autodiff = true)

# Define the mvnormcdf and its gradient, then register the function
prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)
function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
g = norm_pdf((x - μ)/sqrt(vari)) * cdf((y-(μ+ρ/sqrt(vari)*(x-μ)))/sqrt(vari^2*(1- ρ^2)))
g = norm_pdf((y - μ)/sqrt(vari))* cdf((x-(μ+ρ/sqrt(vari)*(y-μ)))/sqrt(vari^2*(1- ρ^2)))
return
end
register(ex2, :mvncdf, 2, prob, ∇prob)

# Use the mvncdf in the nonlinear constraint
# This is where the type error occurs
@NLconstraint(ex2, mvncdf(x, x) >= p)

JuMP.optimize!(ex2)

println(value.(x))

end
``````

The error I get is that the optimizer is getting a tuple instead of a float

``````TypeError: in typeassert, expected Float64, got a value of type Tuple{Float64,Float64}
in top-level scope at example.jl:42
in optimize! at JuMP/e0Uc2/src/optimizer_interface.jl:106
in optimize! at JuMP/e0Uc2/src/optimizer_interface.jl:106
in #optimize!#100 at JuMP/e0Uc2/src/optimizer_interface.jl:130
in optimize! at MathOptInterface/YDdD3/src/Utilities/cachingoptimizer.jl:252
in optimize! at MathOptInterface/YDdD3/src/Bridges/bridge_optimizer.jl:319
in optimize! at NLopt/9omMe/src/MOI_wrapper.jl:942
in  at NLopt/9omMe/src/MOI_wrapper.jl:907
in eval_constraint_jacobian at JuMP/e0Uc2/src/nlp.jl:583
in macro expansion at base/timing.jl:233
in macro expansion at JuMP/e0Uc2/src/nlp.jl:585
in _forward_eval_all at JuMP/e0Uc2/src/nlp.jl:503
in forward_eval at JuMP/e0Uc2/src/_Derivatives/forward.jl:163
in eval_and_check_return_type at JuMP/e0Uc2/src/_Derivatives/forward.jl:5
in eval_objective at JuMP/e0Uc2/src/nlp.jl:1169

``````

Any hints on how to solve this?

Welcome. Moving this over to the Optimization topic area.

1 Like

You can set up `prob(x,y)` as the first value returned from `mvnormcdf`:

``````prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)
``````

Also, your `cdf` in `∇prob` need the first argument to work, e.g.:

``````cdf(Normal(),(y-(μ+ρ/sqrt(vari)*(x-μ)))/sqrt(vari^2*(1- ρ^2))) #I added Normal(), not sure if this is correct for you
``````

Then this works:

``````
using JuMP, LinearAlgebra
using Distributions, NLopt
using Random
using MvNormalCDF

begin
ex2 = Model(NLopt.Optimizer)
set_optimizer_attribute(ex2, "algorithm", :LD_MMA)

@variable(ex2, x[1:2] >= 0)

@objective(ex2, Min, x + x)

set_start_value(x, 5)
set_start_value(x, 5)

p = 0.9 # probability

μ = [5; 5] # mean
vari = [4; 4] # variance
ρ = 0.5 # correlation
Σ = [vari ρ; ρ vari] # covariance matrix
a = [-Inf; -Inf]

norm_cdf(x) = cdf(Normal(), x)
norm_pdf(x) = exp(-(x^2) / 2) / sqrt(2 * pi)
register(ex2, :cdf, 1, norm_cdf, norm_pdf; autodiff = true)

# Define the mvnormcdf and its gradient, then register the function
prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)
function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
g = norm_pdf((x - μ)/sqrt(vari)) * cdf(Normal(),(y-(μ+ρ/sqrt(vari)*(x-μ)))/sqrt(vari^2*(1- ρ^2)))
g = norm_pdf((y - μ)/sqrt(vari))* cdf(Normal(),(x-(μ+ρ/sqrt(vari)*(y-μ)))/sqrt(vari^2*(1- ρ^2)))
return
end

register(ex2, :mvncdf, 2, prob, ∇prob)

# Use the mvncdf in the nonlinear constraint
# This is where the type error occurs
@NLconstraint(ex2,mvncdf(x, x) >= p)

JuMP.optimize!(ex2)

println(value.(x))

end
``````

And I got:

``````[8.247724657052249, 8.247724657052249]
``````
2 Likes

Thank you! This worked.

1 Like