Obscure error in mvncdf()

This error is kind of reproducible but kind of not reproducible. The following is a copy paste that of the output from matlab. Essentially, for a given value of the first two inputs to mvncdf() it throws an error

mvncdf(z_gridvals(1976,:)-z_gridspacing_down(1976,:),z_gridvals(1976,:)+z_gridspacing_up(1976,:),Mew,Sigma)
Error using gpuArray/log
LOG: needs to return a complex result, but this is not supported for real input X on the GPU. Use LOG(COMPLEX(X))
instead.

Error in internal.stats.bvncdf (line 70)
log(t1.*exp(-bsq./(2*asq)) - t2.*sqrt(2*pi).*Phi(-b./a)));

Error in mvncdf (line 249)
y = y + (-1)^i * internal.stats.bvncdf(X, rho, tol/4);

That is the error, but something as minor as the next lines of code being

aaa1=z_gridvals(1976,:)-z_gridspacing_down(1976,:);
aaa2=z_gridvals(1976,:)-z_gridspacing_up(1976,:);
mvncdf(aaa1,aaa2,Mew,Sigma)

ans =

0

and suddenly the error is gone.

This is really weird, because in principle the inputs are the same in both cases.

The values of the four inputs to mvncdf() are as follows (but just putting them in this way does not cause error)

aaa1=[4.546541072119, -19.52465262]
aaa2=[4.548657475153, -17.361228432903]
Mew = [3.6604, -5.0249]
Sigma =[3.1531,  -4.7301;  -4.7301, 7.1776]

As an attachment there is a mat file containing z_gridvals, etc.

If I just run

clear all
load mvnerror.mat
mvncdf(z_gridvals(1976,:)-z_gridspacing_down(1976,:),z_gridvals(1976,:)+z_gridspacing_up(1976,:),Mew,Sigma)

it is enough to generate the error message described at the beginning of this post.

Side note: these are mostly gpuArrays, not sure if that is in any way relevant (so are aaa1 and aaa2 when created as above, so seems unimportant).

[Matlab Release is R2023a]

Answer 1

[Paul]

Walking through the debugger for this code:

%{
aaa1=[4.546541072119, -19.52465262];
aaa2=[4.548657475153, -17.361228432903];
Mew = [3.6604, -5.0249];
Sigma =[3.1531,  -4.7301;  -4.7301, 7.1776];
mvncdf(aaa1,aaa2,Mew,Sigma)

The line 70 in internal.stats.bvncdf is actually a continuation line. The full line is

            p2 = exp(-shk/2 + ...
                log(t1.*exp(-bsq./(2*asq)) - t2.*sqrt(2*pi).*Phi(-b./a)));
%}

The argument to the log function is small but negative, presumably rounding error (copied from my command window):

arglog = -2.964393875047479e-323;

The log of a negative has an imaginary component equal to pi

log(arglog)

ans = -7.4265e+02 + 3.1416e+00i

and the argument to the exp() is (again, copied from my command window)

argexp = -7.438000043200474e+02 + 3.141592653589793e+00i

argexp = -7.4380e+02 + 3.1416e+00i

Normally, I would expect the exp(argexp) to have a small imaginary component, because mathematically, exp(argexp) can be expressed as exp(real(argexp))*exp(1i*imag(argexp))

format long e
exp(real(argexp)), exp(1i*imag(argexp))

ans = 
    9.881312916824931e-324

ans = 
     -1.000000000000000e+00 + 1.224646799147353e-16i

But, presumably the product of 1e-324 and 1e-16 rounds down to zero, so we get

exp(real(argexp))*exp(1i*imag(argexp))

ans = 
   -9.881312916824931e-324

whicih is the same as what exp() calculates on its own

exp(argexp)

ans = 
   -9.881312916824931e-324

So, mvncdf doesn't suffer with complex numbers for these double inputs because the imaginary part magically disappears (which isn't to say that bvncdf shouldn't be more careful about the argument to that log, maybe there are some inputs where the imaginary component will still stick around, unless an imaginary piece is dealt with later in the code)

From the error message, it looks like that on GPU code needs to be more careful to account for the possibility of a negative input to log.

When this code is executed

%{
aaa1=z_gridvals(1976,:)-z_gridspacing_down(1976,:);
aaa2=z_gridvals(1976,:)-z_gridspacing_up(1976,:);
%}

are aaa1 and aaa2 ordinary Matlab doubles (i.e., not gpuarray)?

[Robert Kirkby]

No, aaa1 and aaa2 are just gpuarray (I created them exactly using the code described in original post, so using "aaa1=z_gridvals(1976,:)-z_gridspacing_down(1976,:);". whos reveals them to be gpuArrays (as expected given z_gridvals and z_gridspacing_down are both gpuArrays)

That is part of what confused me, it seems to happen or not happen, but in both cases it is gpuArray.

[Paul]

Copy/paste from the Question, here's the line that causes an error

mvncdf(z_gridvals(1976,:)-z_gridspacing_down(1976,:),z_gridvals(1976,:)+z_gridspacing_up(1976,:),Mew,Sigma)
%                                                                      ^

Here are tthe lines that don't cause a problem, also copy/pasted from the question

aaa1=z_gridvals(1976,:)-z_gridspacing_down(1976,:);
aaa2=z_gridvals(1976,:)-z_gridspacing_up(1976,:);
%                      ^
mvncdf(aaa1,aaa2,Mew,Sigma)

The second argument to the first call to mvncdf uses a +, but aaa2 is defined with a - , so not the same call to mvncdf.

[Robert Kirkby]

Oops, quite right, sorry. My mistake.

[Robert Kirkby]

So then issue was because on gpu, not cpu. That makes more sense now.

I am guessing the different behaviour of mvncdf() in this (admittedly unusual) instance on cpu vs gpu is not intended? Which is considered the correct behaviour? And is it likely that one will be cleaned up in a future release?

[Paul]

Right, the issue appears to be that on the GPU the log() method for the gpuArray class has a different behavior than the log() method for the double class. I searched the doc but couldn't find anything specific to the gpuArray log() method. Based on the error message in the question, it appears that the former requires the input have the complex attribute if the output will be complex, while the latter has no constraint.

Is the different behaviour between GPU and CPU intended? Only TMW knows, but I doubt it.

Which is correct? For this particular case, I suspect that the CPU/double implementation is giving the correct numerical result, but I have no idea whether or not that's true for all inputs, or whether or not it should work the way it does, i.e., don't check that arglog is valid. It seems like it would be safer for mvncdf to detect when arglog is invalid and take appropriate action, whatever that may be.

Is it likely that mvncdf will be cleaned up? Probably not, at least not until someone brings the issue to the attention of TMW via a tech support request. I've had success in the recent past getting TMW to address another issue with mvncdf where it was returning NaN when it should have returned 0 (see this Answers thread) and suggest you submit this issue to tech support if you feel strongly about it. If you do submit to tech support, please update this thread with their response if you don't mind.

Finally, all of the above is based on my assumption that mvncdf works the same way for double and gpuArray inputs, with the only difference being that the functions called from inside mvncdf, like log and exp, are dispatched to the corresponding methods of the the input class.

[Robert Kirkby]

Thanks!!! I will raise to tech support and reference it here once I do.

[Robert Kirkby]

Has been submitted to tech support as

Case Number 06658477 Subject mvncdf(): different behaviour on GPU vs CPU

I will report back on the outcome once there is one.

[Walter Roberson]

-2.964393875047479e-323 / eps(0)

ans = -6

Almost but not exactly the smallest floating point number you can get.

[Robert Kirkby]

MathWorks support said

"Upon further investigation, I found that during the computation the value slightly falls below zero and then we take the log of it, which results in the error. The only fix to this issue is to ensure that the value does not go below zero.

The workaround to this issue is to run on the CPU and avoid using GPU Arrays.

Please know that we have identified this issue, and it will be addressed in future releases."

So hopefully it will be fixed soon in a future release :)

[Paul]

That's the workaround for this particular case. But I'd still be worried that there might be some combination of inputs where, on the CPU, the result for p2 on line 69-70 does have a small imaginary part. Not sure what would happen after line 70 should that occur.