Performance of numerical optimization with triple integral [closed]

Question

I'm trying to solve a numerical optimisation that looks something like this:

Joint[X_, Y_, Z_] := PDF[MultinormalDistribution[{0, 0, 0}, {{1, 1/2, 1/2}, {1/2, 1, 1/2}, {1/2, 1/2, 1}}], {X, Y, Z}]

Cond[z_] := ConditionalMultinormalDistribution[MultinormalDistribution[{0, 0}, {{1, 1/2}, {1/2, 1}}], {z}, {2}]

zstar[y_?NumberQ] := z /. FindRoot[1 - CDF[Cond[z], y] == .3, {z, y - 0.1}]

cap[y1_?NumberQ, y2_?NumberQ] := (1/2)*SurvivalFunction[BinormalDistribution[{0, 0}, {1, 1}, (1/2)], {y1, zstar[y1]}] + (1/2)*SurvivalFunction[BinormalDistribution[{0, 0}, {1, 1}, 1/2], {y2, zstar[y2]}]

q[y_?NumberQ] := NIntegrate[X*Joint[X, Y, Z] , {X, -∞, ∞}, {Y, y, ∞}, {Z, zstar[y], ∞}]

NMaximize[{(1/2)*q[y1] + (1/2)*q[y2], cap[y1, y2] == .1}, {t1, t2}]

The definition of ConditionalMultinormalDistribution follows this post

This numerical optmization problem is computationally very expensive and I would like to improve its performance

Answer 1

Note: In the updated question, there is a genuine optimization problem, which is still amenable to the techniques of the original answer.

Original answer

I worked on the integral before I realized the optimization problem was trivial. Yet I thought I might as well share what I got.

It's sometimes hard to figure out when pursuing a symbolic solution will be faster and more productive than pursuing a numerical one. Numerical solutions are usually fast and approximate with sufficient accuracy, but in this case they are not fast. Normal distributions have been a subject of intense study, and so there's a chance that Mathematica can produce an exact solution. In this case that turns out to be worth trying.

First, when working with numerics in which you might have to adjust the precision (through WorkingPrecision, say), it helps to use exact numbers for coefficients and parameters if available (or set the precision high enough with SetPrecision). In this case the distributions to use are

BinormalDistribution[{0, 0}, {1, 1}, 1/2]
MultinormalDistribution[{0, 0, 0}, {{1, 1/2, 1/2}, {1/2, 1, 1/2}, {1/2, 1/2, 1}}]

Second, cap[y, z] is the same as

CDF[BinormalDistribution[{0, 0}, {1, 1}, 1/2], -{y, z}]

and the built-in function is faster and more accurate.

Third, the triple integral can be solved exactly (with the indefinite limits y and z). The first (triple integral) took twice as long as the second (iterated integral):

Integrate[X*Joint2[X, Y, Z],
     {Y, y, Infinity}, {Z, z, Infinity}, {X, -Infinity, Infinity}]
Fold[Integrate[#1, #2] &, 
     X * Joint2[X, Y, Z],
     {{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, z, Infinity}}]
(* (E^(-(y^2/2)) (1 + Erf[(y - 2 z)/Sqrt[6]]) + E^(-(z^2/2)) Erfc[(2 y - z)/Sqrt[6]]) /
    (4 Sqrt[2 Pi]) *)

(The timings were 30.770272 and 16.735236 resp.)

---

So a maximization problem can be set up with the integration having having been completed beforehand. Most likely it would be quick. I can't really put this to the test on the problem in the OP, since there is only one number in the domain that satisfies the constraint cap[y, 1] == 0.1.

---

Update - Response to the additional problem

The new q can still be done symbolically as above, which will save a lot of time over a triple NIntegrate:

Fold[Integrate[#1, #2] &, 
     X*Joint2[X, Y, Z],
     {{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, zstar[y], Infinity}}]
(*(E^(-(y^2/2)) (1+Erf[(y-2 zstar[y])/Sqrt[6]])
     +E^(-(1/2) zstar[y]^2) \ Erfc[(2 y-zstar[y])/Sqrt[6]])/
  (4 Sqrt[2 \[Pi]])*)

Note: it helps to do the X integral first, probably in this case because it is a definite integral. This is true even for the triple Integrate in the original answer above.

The next most expensive part of the calculation is FindRoot. It can be sped up by about 50% by using InverseCDF instead of CDF:

z /. FindRoot[InverseCDF[Cond[z], 1 - .3] == y, {z, y - 0.1}]

But in this case FindRoot can be replaced by Solve, since Cond[z] == NormalDistribution[z/2, Sqrt[3]/2]. Both of the following return the same thing (although the first generates an error message about using inverse functions:

z /. First@Solve[1 - CDF[Cond[z], y] == .3, z]
z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z] // Simplify
(* -0.908288 + 2. y *)

So zstar can be defined by a linear equation (and could even be made exact, if desired).

With these changes NMaximize takes a little time:

q[y_?NumberQ] := (E^(-(y^2/2)) (1 + Erf[(y - 2 zstar[y])/Sqrt[6]]) + 
                  E^(-(1/2) zstar[y]^2) Erfc[(2 y - zstar[y])/Sqrt[6]]) / (4 Sqrt[2 Pi]);
zstar[y_?NumberQ] := Evaluate[z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z]];

NMaximize[{(1/2)*q[y1] + (1/2)*q[y2], cap[y1, y2] == .1}, {y1, y2}] // Timing
(* {2.683494, {0.0970647, {y1 -> 0.816697, y2 -> 0.816683}}} *)

Answer 2

Map to spherical coordinates

The canonical conversion from cartesian to spherical coordinates:

map = {a, b, c} -> {r Cos[ϕ] Sin[θ], r Sin[ϕ] Sin[θ], r Cos[θ]};

with the scale factor r^2 Sin[θ].

q2[y_?NumberQ, z_?NumberQ] :=Module[{map},
    map = {a, b, c} -> {r Cos[ϕ] Sin[θ], r Sin[ϕ] Sin[θ], r Cos[θ]}; 

    NIntegrate[X*Joint2[X, Y, Z] r^2 Sin[θ] /. Thread[map], 
        {X, -∞, ∞}, {Y, y, ∞}, {Z, z, ∞}, 
        Method -> {Automatic, "SymbolicProcessing" -> 0}, 
        WorkingPrecision -> 10]
]

Now let's call it:

NMaximize[{q2[y, 1], cap[y, 1] == .1}, {y}] // Timing

==> {43.089157, {0., {y -> 0.462727}}}