diff --git a/README.md b/README.md index bb70d3ec..398bcc12 100644 --- a/README.md +++ b/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 diff --git a/benchmarks/README.md b/benchmarks/README.md index b72896e7..57229d1f 100644 --- a/benchmarks/README.md +++ b/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 | diff --git a/mdopt/backend/array.py b/mdopt/backend/array.py index 2a13512a..1a9e944f 100644 --- a/mdopt/backend/array.py +++ b/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== --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 # ---------------------------------------------------------------------- diff --git a/mdopt/contractor/contractor.py b/mdopt/contractor/contractor.py index 25f297a9..79d23de8 100644 --- a/mdopt/contractor/contractor.py +++ b/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) diff --git a/mdopt/utils/utils.py b/mdopt/utils/utils.py index 667aaa90..7dde4e06 100644 --- a/mdopt/utils/utils.py +++ b/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 diff --git a/poetry.lock b/poetry.lock index fa7bc7a9..3d82c815 100644 --- a/poetry.lock +++ b/poetry.lock @@ -1,4 +1,4 @@ -# This file is automatically @generated by Poetry 2.4.1 and should not be changed by hand. +# This file is automatically @generated by Poetry 2.4.3 and should not be changed by hand. [[package]] name = "accessible-pygments" @@ -4627,4 +4627,4 @@ files = [ [metadata] lock-version = "2.1" python-versions = ">=3.11,<4.0" -content-hash = "dd4a5e7d24bc35695c8f673dd13db5bd2ad2b19cd8c638a876466d4bcbcc89f2" +content-hash = "011aaabbb7188b795e0fbb3cf9c3b0bcc63ee44aaf486d05b3365665c405ed19" diff --git a/pyproject.toml b/pyproject.toml index 79c240c9..5f513554 100644 --- a/pyproject.toml +++ b/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", diff --git a/tests/backend/test_array.py b/tests/backend/test_array.py index 472ba3d8..bfd13d25 100644 --- a/tests/backend/test_array.py +++ b/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 diff --git a/tests/decoding/test_convergence.py b/tests/decoding/test_convergence.py index 25c5ed98..7d72bc6d 100644 --- a/tests/decoding/test_convergence.py +++ b/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