Warn about Accelerate-linked NumPy/SciPy, and correct figures added by #543

README.md

@@ -32,6 +32,39 @@ Otherwise, you can clone the repository and use [poetry](https://python-poetry.o
 poetry install
 ```
 
+### A note on NumPy's BLAS on Apple silicon
+
+On Apple silicon running macOS 14 or later, pip installs the NumPy and
+SciPy wheels tagged `macosx_14_0_arm64`, which link Apple's Accelerate
+framework. On mdopt's matrices (rank-deficient, with singular values
+spanning many orders of magnitude) Accelerate's LAPACK corrupted memory:
+`numpy.linalg.qr` died with SIGBUS on a 636x304 matrix, and
+`numpy.linalg.svd` tripped malloc's heap check. How much this showed in the
+decoding depended on the bond dimension. On the [[72,12,6]]
+bivariate-bicycle code in the natural qubit order, at `chi_max=400`, where
+the SVDs act on matrices of about 800x1600, 21 of 22 shots with a
+non-trivial error decoded wrongly while mdopt still reduced its SVDs by QR
+first, against 3 of 24 on OpenBLAS, and every unit test still passed. The
+same test shot decoded correctly at `chi_max` 64, 128 and 256. Not seeing
+the fault at those smaller bond dimensions is no guarantee, because it
+depends on the state of the heap. mdopt's SVD helpers call both libraries,
+and mdopt warns at import when it finds Accelerate behind either one. The
+same versions are also published as OpenBLAS builds, the `macosx_11_0_arm64`
+NumPy wheel and the `macosx_12_0_arm64` SciPy wheel, which install on the
+same machines:
+
+```bash
+pip download "numpy==$(python -c 'import numpy; print(numpy.__version__)')" \
+    --platform macosx_11_0_arm64 --only-binary=:all: --no-deps -d /tmp/openblas-wheels
+pip download "scipy==$(python -c 'import scipy; print(scipy.__version__)')" \
+    --platform macosx_12_0_arm64 --only-binary=:all: --no-deps -d /tmp/openblas-wheels
+pip install --force-reinstall --no-deps /tmp/openblas-wheels/*.whl
+```
+
+Set `MDOPT_ALLOW_ACCELERATE=1` to silence the warning if you must keep
+Accelerate. Use `OMP_NUM_THREADS=1` (or `OPENBLAS_NUM_THREADS=1`) for the
+per-process BLAS of worker pools.
+
 ## Minimal example
 
 ```python

benchmarks/README.md

@@ -10,7 +10,7 @@ and a change is only accepted if the fingerprints still match.
 | `surface_bitflip` | 5x5 surface code, bit-flip noise, 6 shots, χ = 64, natural qubit order | verdicts exact, posterior entries within 1e-2 |
 | `css_optimised` | the same code under `qubit_order_strategy="Optimised"` (reverse Cuthill-McKee) | as above |
 | `shor_depolarising` | Shor code, depolarising noise, 40 shots, χ = 128: the dense-readout path | as above |
-| `classical_ldpc` | random (3,4) LDPC code, XOR constraints + dephasing-DMRG readout | overlaps exact to 1e-10 |
+| `classical_ldpc` | three random (3,4) LDPC codes, XOR constraints + dephasing-DMRG readout | overlaps exact to 1e-10 |
 | `dmrg_ground_state` | DMRG on a 24-site transverse-field Ising chain | energy exact to 1e-10 |
 | `dem_d3` | circuit-level detector error model, d=3 r=3 p=0.8%, eight busiest syndromes, χ = 32 | verdicts exact, class masses within 1e-2 |
 | `dem_d5` | the same at d=5 r=5 p=0.5%, busiest syndrome, χ = 32 | as above |

mdopt/backend/array.py

@@ -60,6 +60,80 @@ def _load_backend():
 GPU = _xp.__name__ == "cupy"
 
 
+def _lapack_vendor(module) -> str:
+    """The LAPACK library a NumPy or SciPy module was built against, lower-cased.
+
+    Read from ``show_config(mode="dicts")``, which every NumPy and SciPy
+    release mdopt supports provides (NumPy 2.3+, SciPy 1.15+). Returns ``""``
+    when the configuration cannot be read.
+    """
+    try:
+        deps = module.show_config(mode="dicts")["Build Dependencies"]
+        return str(deps["lapack"]["name"]).lower()
+    except Exception:  # pylint: disable=broad-except
+        return ""
+
+
+def _warn_if_accelerate(numpy_module) -> None:
+    """Warn once when NumPy's LAPACK is Apple's Accelerate framework.
+
+    On the NumPy 2.x wheels for macOS 14+ on Apple silicon
+    (``macosx_14_0_arm64``, which link Accelerate) the decoders'
+    rank-deficient, wide-spectrum matrices made ``linalg.qr`` die with SIGBUS
+    and ``linalg.svd`` trip malloc's heap-corruption check inside dgesdd.
+    While mdopt still reduced its SVDs by QR first, this also turned a
+    [[72,12,6]] decode at chi_max=400 into wrong verdicts (21 of 22 shots with
+    a non-trivial error) while every unit test passed. The OpenBLAS build of
+    the same NumPy version has none of this; it is the ``macosx_11_0_arm64``
+    wheel::
+
+        pip download numpy==<version> --platform macosx_11_0_arm64 \\
+            --only-binary=:all: --no-deps -d /tmp/numpy-openblas
+        pip install --force-reinstall --no-deps /tmp/numpy-openblas/*.whl
+
+    Set MDOPT_ALLOW_ACCELERATE=1 to silence the warning.
+    """
+    if os.getenv("MDOPT_ALLOW_ACCELERATE") == "1":
+        return
+    if "accelerate" in _lapack_vendor(numpy_module):
+        warnings.warn(
+            "NumPy is built against Apple's Accelerate LAPACK, which corrupted "
+            "memory on mdopt's matrices and, with an earlier SVD code path, led "
+            "to wrong decoding verdicts (see "
+            "mdopt.backend.array._warn_if_accelerate). Install the "
+            "OpenBLAS build of NumPy (the macosx_11_0_arm64 wheel) or set "
+            "MDOPT_ALLOW_ACCELERATE=1 to silence this warning.",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+
+
+def _warn_if_scipy_accelerate() -> None:
+    """The same warning for SciPy, whose LAPACK the SVD helpers also call."""
+    if os.getenv("MDOPT_ALLOW_ACCELERATE") == "1":
+        return
+    try:
+        scipy = importlib.import_module("scipy")
+    except ImportError:
+        return
+    if "accelerate" in _lapack_vendor(scipy):
+        warnings.warn(
+            "SciPy is built against Apple's Accelerate LAPACK (see "
+            "mdopt.backend.array._warn_if_accelerate); install the OpenBLAS "
+            "build (the macosx_12_0_arm64 wheel) or set "
+            "MDOPT_ALLOW_ACCELERATE=1 to silence this warning.",
+            RuntimeWarning,
+            stacklevel=2,
+        )
+
+
+# Both checks run on every backend: with CuPy selected the orthogonality-centre
+# moves and host-side contractions still call NumPy's linalg and BLAS, and the
+# SVD fallbacks and qr call SciPy's LAPACK.
+_warn_if_accelerate(importlib.import_module("numpy"))
+_warn_if_scipy_accelerate()
+
+
 # ----------------------------------------------------------------------
 # Introspection helpers
 # ----------------------------------------------------------------------

mdopt/contractor/contractor.py

@@ -23,8 +23,9 @@ def _contract_cached(subscripts, path, backend, *tensors):
     rebuilds path metadata on every call even when ``optimize`` is explicit.
     An expression with an explicit path is shape-independent, so it is
     cached per (subscripts, path) only -- keying on operand shapes made
-    truncation's data-dependent bond dimensions miss 7-14% of calls on
-    large codes -- and built from whatever shapes the first call carries.
+    truncation's data-dependent bond dimensions miss the cache (1-2% of the
+    zip-up calls on the circuit-level DEM workloads, 23% on the RCM
+    surface-code one) -- and built from whatever shapes the first call carries.
     """
     key = (subscripts, path)
     expression = _EXPRESSIONS.get(key)

mdopt/utils/utils.py

@@ -16,7 +16,8 @@ def _to_numpy(a):
 
     Resolved through the backend's own transfer, which honours MDOPT_BACKEND
     and the CUDA device probe; a per-call ``import cupy`` here once cost
-    ~18% of a decoding run.
+    about 11,000 failed imports per decode: from about 4% of the decode on
+    surface_bitflip to about a third on dem_d3.
     """
     host = xp.to_host(a)
     try:
@@ -93,12 +94,14 @@ def svd(
                 #   (the macOS 14+ arm64 wheels) it corrupted memory on the
                 #   decoders' rank-deficient matrices: numpy.linalg.qr died
                 #   with SIGBUS, and a [[72,12,6]] decode at chi_max=400
-                #   returned wrong verdicts on 24 of 24 shots while every
-                #   unit test and benchmark fingerprint still passed;
+                #   returned wrong verdicts on 21 of 22 shots with a
+                #   non-trivial error, while every unit test and benchmark
+                #   fingerprint still passed;
                 # - with OpenBLAS it is memory-safe but not equivalent under
                 #   truncation: it moved a chi_max=64 surface-code posterior
                 #   entry by 0.017, beyond the benchmark suite's tolerance;
-                # - it saved 0-7% on the benchmark workloads.
+                # - it saved at most 7% on the benchmark workloads, and
+                #   nothing measurable on most of them.
                 # Reintroducing it needs `benchmarks/bench_suite.py --check`
                 # and tests/decoding/test_convergence.py to pass.
                 # No finiteness pre-scan: a non-finite input makes the SVD

pyproject.toml

@@ -23,7 +23,7 @@ classifiers = [
 ]
 
 dependencies = [
-    "scipy>=1.9.2,<2.0.0",
+    "scipy>=1.15.0,<2.0.0",
     "opt-einsum>=3.3.0,<4.0.0",
     "more-itertools>=8.12,<11.0",
     "threadpoolctl>=3.6.0,<4.0.0",

tests/backend/test_array.py

@@ -45,11 +45,23 @@ def get_device_count():
 
 
 def test_defaults_to_numpy(monkeypatch):
-    """With no MDOPT_BACKEND set, NumPy is used and nothing is warned about."""
+    """With no MDOPT_BACKEND set, NumPy is used and nothing unexpected is warned about.
+
+    The import deliberately warns when NumPy or SciPy links Apple's Accelerate
+    framework, which the stock macOS 14+ wheels do (the macOS CI runners among
+    them); those warnings are allowed, any other warning is not.
+    """
     module, caught = _reload_backend(monkeypatch, None, None)
     assert module.GPU is False
     assert module.backend_name() == "numpy"
-    assert not caught
+    unexpected = [
+        w
+        for w in caught
+        if not (
+            issubclass(w.category, RuntimeWarning) and "Accelerate" in str(w.message)
+        )
+    ]
+    assert not unexpected, [str(w.message) for w in unexpected]
 
 
 def test_falls_back_when_cupy_missing(monkeypatch):
@@ -79,7 +91,14 @@ def test_uses_cupy_when_device_present(monkeypatch):
     module, caught = _reload_backend(monkeypatch, "cupy", _fake_cupy(num_devices=1))
     assert module.GPU is True
     assert module.backend_name() == "cupy"
-    assert not caught
+    unexpected = [
+        w
+        for w in caught
+        if not (
+            issubclass(w.category, RuntimeWarning) and "Accelerate" in str(w.message)
+        )
+    ]
+    assert not unexpected, [str(w.message) for w in unexpected]
 
 
 def test_stream_is_a_context_manager_on_numpy(monkeypatch):
@@ -88,3 +107,86 @@ def test_stream_is_a_context_manager_on_numpy(monkeypatch):
     with module.stream():
         pass
     module.synchronize()
+
+
+def test_accelerate_lapack_warns_unless_allowed(monkeypatch):
+    """A NumPy built against Accelerate triggers a RuntimeWarning at import
+    time; MDOPT_ALLOW_ACCELERATE=1 silences it; other vendors are silent."""
+    from types import SimpleNamespace
+
+    from mdopt.backend import array as backend
+
+    def fake_numpy(vendor):
+        return SimpleNamespace(
+            show_config=lambda mode="dicts": {
+                "Build Dependencies": {
+                    "lapack": {"name": vendor},
+                    "blas": {"name": vendor},
+                }
+            }
+        )
+
+    monkeypatch.delenv("MDOPT_ALLOW_ACCELERATE", raising=False)
+    with pytest.warns(RuntimeWarning, match="Accelerate"):
+        backend._warn_if_accelerate(fake_numpy("accelerate"))
+    with warnings.catch_warnings():
+        warnings.simplefilter("error")
+        backend._warn_if_accelerate(fake_numpy("scipy-openblas"))
+        monkeypatch.setenv("MDOPT_ALLOW_ACCELERATE", "1")
+        backend._warn_if_accelerate(fake_numpy("accelerate"))
+    # An unreadable configuration (no show_config, or one without the dicts
+    # mode) reports an unknown vendor instead of raising.
+    assert backend._lapack_vendor(SimpleNamespace()) == ""
+    assert backend._lapack_vendor(SimpleNamespace(show_config=lambda: None)) == ""
+
+
+def test_scipy_shaped_config_is_detected(monkeypatch):
+    """SciPy reports its LAPACK through show_config(mode="dicts"), with
+    the vendor capitalised ("Accelerate"); the SciPy check must warn on it."""
+    from types import SimpleNamespace
+
+    from mdopt.backend import array as backend
+
+    def show_config(mode="stdout"):
+        if mode == "stdout":
+            return None
+        return {
+            "Build Dependencies": {
+                "blas": {"name": "Accelerate"},
+                "lapack": {"name": "Accelerate"},
+            }
+        }
+
+    monkeypatch.delenv("MDOPT_ALLOW_ACCELERATE", raising=False)
+    monkeypatch.setitem(sys.modules, "scipy", SimpleNamespace(show_config=show_config))
+    with pytest.warns(RuntimeWarning, match="SciPy"):
+        backend._warn_if_scipy_accelerate()
+
+
+def test_accelerate_warnings_fire_on_the_gpu_backend(monkeypatch):
+    """With CuPy selected, NumPy and SciPy LAPACK still run on the host (centre
+    moves, host-side contractions, the SVD fallbacks, qr), so an
+    Accelerate-linked build must still be warned about."""
+    import numpy
+    import scipy
+
+    def show_config(mode="stdout"):
+        # A synthetic config, independent of what the installed NumPy and
+        # SciPy report.
+        if mode == "stdout":
+            return None
+        return {
+            "Build Dependencies": {
+                "blas": {"name": "accelerate"},
+                "lapack": {"name": "accelerate"},
+            }
+        }
+
+    monkeypatch.delenv("MDOPT_ALLOW_ACCELERATE", raising=False)
+    monkeypatch.setattr(numpy, "show_config", show_config)
+    monkeypatch.setattr(scipy, "show_config", show_config)
+    module, caught = _reload_backend(monkeypatch, "cupy", _fake_cupy(num_devices=1))
+    assert module.GPU is True
+    messages = [str(w.message) for w in caught]
+    assert any("NumPy is built against" in m for m in messages), messages
+    assert any("SciPy is built against" in m for m in messages), messages

tests/decoding/test_convergence.py

@@ -29,9 +29,10 @@ def test_bb_72_12_6_natural_order_single_error_converges_in_chi():
 
     This instance exposed the QR pre-reduction of ``svd`` (since removed):
     with it, the decode returned a flat or wrongly peaked
-    posterior at chi_max=400 on NumPy/Accelerate builds -- 24 of 24 shots
-    -- while every unit test and benchmark fingerprint still passed. Each
-    chi_max=400 decode takes 30-60 minutes on a laptop.
+    posterior at chi_max=400 on NumPy/Accelerate builds, as did 21 of the 22
+    shots with a non-trivial error in the natural-order experiment, while
+    every unit test and benchmark fingerprint still passed. The test takes
+    about 20 minutes on an idle laptop.
     """
     code = create_bb_code(6, 6, "x**3 + y + y**2", "y**3 + x + x**2")
     error = "I" * 58 + "Z" + "I" * 13