BUG: huge speed and quality regression for stats.multivariate_normal.cdf

[HDembinski]

Describe your issue.

There is a huge speed regression (180x slower) in computing the cdf for the multivariate normal since scipy-1.16.0. The numerical precision is also worse and that lets unit tests in the iminuit package fail that used to pass since the Python-3.9 days.

Tested on MacOS, Linux.

Reproducing Code Example

import numpy as np
from scipy import stats


def mvnorm(mux, muy, sx, sy, rho):
    C = np.empty((2, 2))
    C[0, 0] = sx**2
    C[0, 1] = C[1, 0] = sx * sy * rho
    C[1, 1] = sy**2
    m = [mux, muy]
    return stats.multivariate_normal(m, C)

def model(xy, mux, muy, sx, sy, rho):
    return mvnorm(mux, muy, sx, sy, rho).cdf(xy.T)

rng = np.random.default_rng(42)

%%timeit

model(rng.normal(size=(2, 1000)), 0, 0, 1, 1, 0)

v1.15.3: 1.63 ms ± 7.17 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
v1.16.0: 286 ms ± 21.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
v1.16.1: 265 ms ± 5.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Error message

None

SciPy/NumPy/Python version and system information

1.15.3 2.2.5 sys.version_info(major=3, minor=12, micro=8, releaselevel='final', serial=0)
C:\Users\hansd.Powercat\AppData\Roaming\Python\Python312\site-packages\scipy\__config__.py:154: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "3.0.12",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "10.3.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.17.0",
      "include directory": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\pip-build-env-dihuxz0e\\overlay\\Lib\\site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "windows"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.28",
      "detection method": "pkgconfig",
      "include directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-3dk08vlg/cp312-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include",
      "lib directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-3dk08vlg/cp312-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib",
      "openblas configuration": "OpenBLAS 0.3.28 DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=24",
      "pc file directory": "D:/a/scipy/scipy"
    },
    "lapack": {
      "name": "scipy-openblas",
      "found": true,
      "version": "0.3.28",
      "detection method": "pkgconfig",
      "include directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-3dk08vlg/cp312-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include",
      "lib directory": "C:/Users/runneradmin/AppData/Local/Temp/cibw-run-3dk08vlg/cp312-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib",
      "openblas configuration": "OpenBLAS 0.3.28 DYNAMIC_ARCH NO_AFFINITY Haswell MAX_THREADS=24",
      "pc file directory": "D:/a/scipy/scipy"
    },
    "pybind11": {
      "name": "pybind11",
      "version": "2.13.6",
      "detection method": "config-tool",
      "include directory": "unknown"
    }
  },
  "Python Information": {
    "path": "C:\\Users\\runneradmin\\AppData\\Local\\Temp\\cibw-run-3dk08vlg\\cp312-win_amd64\\build\\venv\\Scripts\\python.exe",
    "version": "3.12"
  }
}

[lucascolley]

from gh-22298:

1. the replacement code is obscenely slow. Probably needs some cythonizing or pythranizing, as discussed in ENH: stats.multivariate_t: add cdf method #17410 and elsewhere

And see also #22298 (comment).

Is this the same as #23413 (comment)?

[rkern]

I think it has to do with the bivariate special case in the old version of the code (and the lack of one, currently), yes.

[mdhaber]

Yes, copying from there:

That would make sense. That PR did not replace the bivariate-specific version of the code. I had an LLM translate the bivariate code - partly to spark the discussion of whether that would be an acceptable use of LLMs - but I believe the conclusion was that it was not. So if someone wants to manually translate bvnu from Matlab to Python (to be completed with Pythran) or C/C++, that could help here. (Should be quite easy - the Matlab is basically Python.)

[ev-br]

Here's a manual, no-LLM translation of the bvnu routine to python:
https://github.com/scipy/scipy/compare/main...ev-br:scipy:mvn_2D?expand=1

[mdhaber]

Do you want to open a PR with it?

[ev-br]

My bandwidth is limited ATM, so if someone beats me to it, great.

[mdhaber]

I wasn't getting good agreement with SciPy 1.15 mvnun. I'm sure there's just a simple bug somewhere, but I thought it would be easiest to just try a line by line myself. Here it is along with a script to test it against _mvn from SciPy 1.15:

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.stats import _mvn

rng = np.random.default_rng(8348324598925153)

# function p = bvnu( dh, dk, r )
def bvnu(dh, dk, r):
    # if dh ==  inf | dk ==  inf:p = 0;
    if (dh == math.inf) or (dk == math.inf):
        p = 0
    # elseif dh == -inf, if dk == -inf, p = 1; else p = phid(-dk); end
    elif dh == -math.inf:
        if dk == -math.inf:
            p = 1
        else:
            p = phid(-dk)
    # elseif dk == -inf, p = phid(-dh);
    elif dk == -math.inf:
        p = phid(-dh)
    # elseif r == 0, p = phid(-dh)*phid(-dk);
    elif r == 0:
        p = phid(-dh) * phid(-dk)
    # else, tp = 2*pi; h = dh; k = dk; hk = h*k; bvn = 0;
    else:
        tp = 2*math.pi
        h = dh
        k = dk
        hk = h*k
        bvn = 0
        # if abs(r) < 0.3      % Gauss Legendre points and weights, n =  6
        #     w(1:3) = [0.1713244923791705 0.3607615730481384 0.4679139345726904];
        #     x(1:3) = [0.9324695142031522 0.6612093864662647 0.2386191860831970];
        if abs(r) < 0.3:
            w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904]
            x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970]
        # elseif abs(r) < 0.75 % Gauss Legendre points and weights, n = 12
        #     w(1:3) = [.04717533638651177 0.1069393259953183 0.1600783285433464];
        #     w(4:6) = [0.2031674267230659 0.2334925365383547 0.2491470458134029];
        #     x(1:3) = [0.9815606342467191 0.9041172563704750 0.7699026741943050];
        #     x(4:6) = [0.5873179542866171 0.3678314989981802 0.1252334085114692];
        elif abs(r) < 0.75:
            w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464,
                 0.2031674267230659, 0.2334925365383547, 0.2491470458134029]
            x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
                 0.5873179542866171, 0.3678314989981802, 0.1252334085114692]
        # else,                % Gauss Legendre points and weights, n = 20
        #     w(1:3) = [.01761400713915212 .04060142980038694 .06267204833410906];
        #     w(4:6) = [.08327674157670475 0.1019301198172404 0.1181945319615184];
        #     w(7:9) = [0.1316886384491766 0.1420961093183821 0.1491729864726037];
        #     w(10) =   0.1527533871307259;
        #     x(1:3) = [0.9931285991850949 0.9639719272779138 0.9122344282513259];
        #     x(4:6) = [0.8391169718222188 0.7463319064601508 0.6360536807265150];
        #     x(7:9) = [0.5108670019508271 0.3737060887154196 0.2277858511416451];
        #     x(10) =   0.07652652113349733;
        else:
            w = [.01761400713915212, .04060142980038694, .06267204833410906,
                 .08327674157670475, 0.1019301198172404, 0.1181945319615184,
                 0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
                 0.1527533871307259]
            x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
                 0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
                 0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
                 0.07652652113349733]
        # end, w = [w  w]; x = [1-x 1+x];
        w, x = np.asarray(w), np.asarray(x)
        w = np.concatenate((w, w))
        x = np.concatenate((1-x, 1+x))
        # if abs(r) < 0.925, hs = ( h*h + k*k )/2; asr = asin(r)/2;
        if abs(r) < 0.925:
            hs = ( h*h + k*k )/2
            asr = np.asin(r)/2
        #     sn = sin(asr*x); bvn = exp((sn*hk-hs)./(1-sn.^2))*w';
            sn = np.sin(asr*x)
            bvn = np.exp((sn*hk-hs) / (1-sn**2)) @ w
        #     bvn = bvn*asr/tp + phid(-h)*phid(-k);
            bvn = bvn*asr/tp + phid(-h)*phid(-k)
        # else, if r < 0, k = -k; hk = -hk; end
        else:
            if r < 0:
                k = -k
                hk = -hk
            # if abs(r) < 1, as = 1-r^2; a = sqrt(as); bs = (h-k)^2;
            if abs(r) < 1:
                as_ = 1-r**2
                a = math.sqrt(as_)
                bs = (h-k)**2
                # asr = -( bs/as + hk )/2; c = (4-hk)/8 ; d = (12-hk)/80;
                asr = -( bs/as_ + hk )/2
                c = (4-hk)/8
                d = (12-hk)/80
                # if asr > -100, bvn = a*exp(asr)*(1-c*(bs-as)*(1-d*bs)/3+c*d*as^2); end
                if asr > -100:
                    bvn = a*math.exp(asr)*(1-c*(bs-as_)*(1-d*bs)/3+c*d*as_**2)
                # if hk  > -100, b = sqrt(bs); sp = sqrt(tp)*phid(-b/a);
                if hk  > -100:
                    b = math.sqrt(bs)
                    sp = math.sqrt(tp)*phid(-b/a)
                    # bvn = bvn - exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 );
                    bvn = bvn - math.exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 )

                # end, a = a/2; xs = (a*x).^2; asr = -( bs./xs + hk )/2;
                a = a/2
                xs = (a*x)**2
                asr = -( bs / xs + hk )/2;
                # ix = find( asr > -100 ); xs = xs(ix); sp = ( 1 + c*xs.*(1+5*d*xs) );
                ix = asr > -100
                xs = xs[ix]
                sp = 1 + c*xs * (1+5*d*xs)
                # rs = sqrt(1-xs); ep = exp( -(hk/2)*xs./(1+rs).^2 )./rs;
                rs = np.sqrt(1-xs)
                ep = np.exp( -(hk/2)*xs / (1+rs)**2 )/rs
                # bvn = ( a*( (exp(asr(ix)).*(sp-ep))*w(ix)' ) - bvn )/tp;
                bvn = ( a*( (np.exp(asr[ix]) * (sp-ep)) @ w[ix] ) - bvn )/tp
            # end
            # if r > 0, bvn =  bvn + phid( -max( h, k ) );
            if r > 0:
                bvn =  bvn + phid( -max( h, k ) )
            # elseif h >= k, bvn = -bvn;
            elif h >= k:
                bvn = -bvn
            # else, if h < 0, L = phid(k)-phid(h); else, L = phid(-h)-phid(-k); end
            else:
                if h < 0:
                    L = phid(k)-phid(h)
                else:
                    L = phid(-h)-phid(-k)
                # bvn =  L - bvn;
                bvn =  L - bvn
            # end
        # end, p = max( 0, min( 1, bvn ) );
        p = max( 0, min( 1, bvn ) )
    # end
    return p

# function p = phid(z), p = erfc( -z/sqrt(2) )/2; % Normal cdf
def phid(z):
    return math.erfc(-z / math.sqrt(2)) / 2

def bvn(a, b, A):
    # covariance matrix is written [[s1**2, rho*s1*s2], [rho*s1*s2, s2**2]]
    # e.g. https://en.wikipedia.org/wiki/Multivariate_normal_distribution
    # therefore, s12 = rho*s1*s2 -> rho = s12/(s1*s2)
    s1 = np.sqrt(A[0, 0])
    s2 = np.sqrt(A[1, 1])
    s12 = A[0, 1]
    r = s12 / (s1 * s2)
    # the x and y coordinates seem to be normalized by the standard devs
    xl, xu = a[0] / s1, b[0] / s1
    yl, yu = a[1] / s2, b[1] / s2
    p = bvnu(xl, yl, r) - bvnu(xu, yl, r) - bvnu(xl, yu, r) + bvnu(xu, yu, r) 
    p = max( 0, min( p, 1 ) )
    return p

errs = []
rerrs = []
for i in range(100000):
    r = 10**rng.uniform(-2, 2) * np.sign(rng.normal())
    w = rng.random(2)
    w = w / np.linalg.norm(w)
    v = stats.special_ortho_group.rvs(2, random_state=rng)
    A = np.asarray([[1, r], [r, 1]])
    A = v @ np.diag(w) @ v.T

    dist = stats.multivariate_normal(cov=A)
    a = np.asarray([-np.inf, -np.inf])
    b = dist.rvs()

    res = bvn(a, b, A)
    ref = _mvn.mvnun(a, b, np.zeros(2), A, 100000, 1e-5, 1e-5)[0]
    err = np.abs(res - ref)
    errs.append(err)
    rerrs.append(err/ref)
    # np.testing.assert_allclose(res, ref, atol=2e-16, rtol=1e-15)

log_errs = np.log10(np.maximum(errs, 1e-18))
log_rerrs = np.log10(np.maximum(rerrs, 1e-18))

plt.hist(log_errs, bins=30)
plt.xlabel("log10 absolute error")

plt.hist(log_rerrs, bins=30)
plt.xlabel("log10 relative error")

I think the test is reasonably robust. The correlation coefficient values are loguniformly distributed over a pretty wide range, and the points to check the CDF are sampled from the resulting distribution. You can tell that there is some testing in the tails since the absolute errors are all small but the relative errors sometimes get larger - and I think that's to be expected, since there is at least one place in the code where subtractive cancellation can occur.

If that looks good, I can drop it (or your version after adjustment) into _quato and open a PR.

[ev-br]

LGTM!
The only minor tweak is there's no need to hardcode the roots and wieghts for the quadrature, can reuse special.roots_legendre instead: x, w = roots_legendre(n); x = 1 + x.

[mdhaber]

Yeah, I checked that out. I was going for as direct of a translation as possible to make it easy to spot the source of differences, if there were any. I will think about switching when I make the PR.