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

Specific Domains Optimization (Mathematical)
,
1

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[1] + x[2])

    set_start_value(x[1], 5)
    set_start_value(x[2], 5)

    p = 0.9 # probability

    μ = [5; 5] # mean
    vari = [4; 4] # variance
    ρ = 0.5 # correlation
    Σ = [vari[1] ρ; ρ vari[2]] # 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[1] = norm_pdf((x - μ[1])/sqrt(vari[1])) * cdf((y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2)))
        g[2] = norm_pdf((y - μ[2])/sqrt(vari[2]))* cdf((x-(μ[1]+ρ/sqrt(vari[2])*(y-μ[2])))/sqrt(vari[1]^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[1], x[2]) >= 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?

2

Welcome. Moving this over to the Optimization topic area.

1 Like
3

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

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

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

cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^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[1] + x[2])

    set_start_value(x[1], 5)
    set_start_value(x[2], 5)

    p = 0.9 # probability

    μ = [5; 5] # mean
    vari = [4; 4] # variance
    ρ = 0.5 # correlation
    Σ = [vari[1] ρ; ρ vari[2]] # 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]-μ)[1]
    function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
        g[1] = norm_pdf((x - μ[1])/sqrt(vari[1])) * cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2)))
        g[2] = norm_pdf((y - μ[2])/sqrt(vari[2]))* cdf(Normal(),(x-(μ[1]+ρ/sqrt(vari[2])*(y-μ[2])))/sqrt(vari[1]^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[1], x[2]) >= p)

    JuMP.optimize!(ex2)

    println(value.(x))

end

And I got:

[8.247724657052249, 8.247724657052249]
3 Likes
4

Thank you! This worked.

1 Like