<p dir="auto"><em>Original ticket <a href="http://projects.scipy.org/scipy/ticket/491" rel="nofollow">http://projects.scipy.org/scipy/ticket/491</a> on 2007-08-31 by trac user Mark Poolman, assigned to unknown.</em></p>
<p dir="auto">There appears to be a problem in svd results for certain large matrices, as demonstrated by use of svd to calculate a kernel. Briefly, for some matrices if k = kernel(a) we find a.k != 0. This has been seen with some, but not all,<br>
combinations of various flavours of ubunti and x86_64 and i686 hardware.</p>
<p dir="auto">A more complete description is here:<br>
<a href="https://bugs.launchpad.net/ubuntu/+source/python-scipy/+bug/134277" rel="nofollow">https://bugs.launchpad.net/ubuntu/+source/python-scipy/+bug/134277</a></p>
<p dir="auto"><em>This is a follow up to <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="987205536" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14688" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/14688/hovercard" href="https://github.com/scipy/scipy/issues/14688">gh-14688</a>. That issue was originally about a kernel panic (fixed in macOS 12.0.1), and after that the same reproducer showed severe performance issues. This issue is about those performance issues. Note that while the reproducer is the same, it's not clear whether or not the kernel panic and the performance issues share a root cause or not.</em></p>
<h2 dir="auto">Issue reproducer</h2>
<p dir="auto">A reproducer (warning: do <em>NOT</em> run on macOS 11.x, it will crash the OS):</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="from time import perf_counter
import numpy as np
from scipy.sparse.linalg import eigsh


n_samples, n_features = 2000, 10
rng = np.random.default_rng(0)
X = rng.normal(size=(n_samples, n_features))
K = X @ X.T

for i in range(10):
    print(&quot;running eigsh...&quot;)
    tic = perf_counter()
    s, _ = eigsh(K, 3, which=&quot;LA&quot;, tol=0)
    toc = perf_counter()
    print(f&quot;computed {s} in {toc - tic:.3f} s&quot;)"><pre class="notranslate"><span class="pl-k">from</span> <span class="pl-s1">time</span> <span class="pl-k">import</span> <span class="pl-s1">perf_counter</span>
<span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
<span class="pl-k">from</span> <span class="pl-s1">scipy</span>.<span class="pl-s1">sparse</span>.<span class="pl-s1">linalg</span> <span class="pl-k">import</span> <span class="pl-s1">eigsh</span>


<span class="pl-s1">n_samples</span>, <span class="pl-s1">n_features</span> <span class="pl-c1">=</span> <span class="pl-c1">2000</span>, <span class="pl-c1">10</span>
<span class="pl-s1">rng</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">random</span>.<span class="pl-en">default_rng</span>(<span class="pl-c1">0</span>)
<span class="pl-v">X</span> <span class="pl-c1">=</span> <span class="pl-s1">rng</span>.<span class="pl-en">normal</span>(<span class="pl-s1">size</span><span class="pl-c1">=</span>(<span class="pl-s1">n_samples</span>, <span class="pl-s1">n_features</span>))
<span class="pl-v">K</span> <span class="pl-c1">=</span> <span class="pl-v">X</span> @ <span class="pl-v">X</span>.<span class="pl-v">T</span>

<span class="pl-k">for</span> <span class="pl-s1">i</span> <span class="pl-c1">in</span> <span class="pl-en">range</span>(<span class="pl-c1">10</span>):
    <span class="pl-en">print</span>(<span class="pl-s">"running eigsh..."</span>)
    <span class="pl-s1">tic</span> <span class="pl-c1">=</span> <span class="pl-en">perf_counter</span>()
    <span class="pl-s1">s</span>, <span class="pl-s1">_</span> <span class="pl-c1">=</span> <span class="pl-en">eigsh</span>(<span class="pl-v">K</span>, <span class="pl-c1">3</span>, <span class="pl-s1">which</span><span class="pl-c1">=</span><span class="pl-s">"LA"</span>, <span class="pl-s1">tol</span><span class="pl-c1">=</span><span class="pl-c1">0</span>)
    <span class="pl-s1">toc</span> <span class="pl-c1">=</span> <span class="pl-en">perf_counter</span>()
    <span class="pl-en">print</span>(<span class="pl-s">f"computed <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">s</span><span class="pl-kos">}</span></span> in <span class="pl-s1"><span class="pl-kos">{</span><span class="pl-s1">toc</span> <span class="pl-c1">-</span> <span class="pl-s1">tic</span>:.3f<span class="pl-kos">}</span></span> s"</span>)</pre></div>
<p dir="auto">Running <code class="notranslate">scipy.test()</code> or <code class="notranslate">scipy.linalg.test()</code> will also show a significant performance impact.</p>
<h2 dir="auto">Performance impact</h2>
<p dir="auto">In situations where we hit the performance problem, the above code will show:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 1.062 s
running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.888 s
..."><pre class="notranslate"><code class="notranslate">running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 1.062 s
running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.888 s
...
</code></pre></div>
<p dir="auto">And if we don't hit that problem:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.018 s
running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.018 s
..."><pre class="notranslate"><code class="notranslate">running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.018 s
running eigsh...
computed [2096.59480134 2144.21662824 2188.11799492] in 0.018 s
...
</code></pre></div>
<p dir="auto">So a ~50x slowdown for this particular example.</p>
<p dir="auto">There is in general an impact on functions that use BLAS/LAPACK. The impact on the total time taken by <code class="notranslate">scipy.test()</code> was about 30% (311 sec. with default settings, 234 sec. when using <code class="notranslate">OPENBLAS_NUM_THREADS=1</code>) - note that this was just a single test on one build config, results may vary: <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="987205536" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14688" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/14688/hovercard?comment_id=968748706&amp;comment_type=issue_comment" href="https://github.com/scipy/scipy/issues/14688#issuecomment-968748706">#14688 (comment)</a>. The single-threaded case has similar timings as when running the test suite on a <code class="notranslate">scipy</code> install that doesn't show the problem at all (~240 sec. seems expected on arm64 macOS, and it doesn't depend on the threading setting (because test arrays are always small)).  Important: ensure <code class="notranslate">pytest-xdist</code> is not installed when looking at time taken by the test suite (see <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="946849269" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14425" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/14425/hovercard" href="https://github.com/scipy/scipy/issues/14425">gh-14425</a> for why).</p>
<h2 dir="auto">When the problem occurs</h2>
<p dir="auto">The discussion in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="987205536" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14688" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/14688/hovercard" href="https://github.com/scipy/scipy/issues/14688">gh-14688</a> showed that this problem gets hit when two copies of <code class="notranslate">libopenblas</code> get loaded. The following configurations showed a problem so far:</p>
<ul dir="auto">
<li>Installing both <code class="notranslate">numpy</code> and <code class="notranslate">scipy</code> from a wheel (e.g., <code class="notranslate">numpy</code> 1.21.4 from PyPI and the latest <code class="notranslate">1.8.0.dev0</code> wheel from <a href="https://anaconda.org/scipy-wheels-nightly/scipy/" rel="nofollow">https://anaconda.org/scipy-wheels-nightly/scipy/</a>)</li>
<li>Installing <code class="notranslate">numpy</code> 1.21.4 from PyPI and installing <code class="notranslate">scipy</code> locally when built against conda-forge's <code class="notranslate">openblas</code>.</li>
</ul>
<p dir="auto">These configurations did not show a problem:</p>
<ul dir="auto">
<li>Installing <code class="notranslate">numpy</code> 1.21.4 from PyPI and installing <code class="notranslate">scipy</code> locally when built against Homebrew's <code class="notranslate">openblas</code>.</li>
<li>Any situation where only a single <code class="notranslate">libopenblas</code> is loaded.</li>
</ul>
<p dir="auto">It is unclear right now what the exact root cause is. The situation when using conda-forge's <code class="notranslate">openblas</code> is very similar to that using Homebrew's <code class="notranslate">openblas</code>, but only one of those triggers the issue. The most important situation is installing both NumPy and SciPy from wheels though, that's what the vast majority of <code class="notranslate">pip</code>/PyPI users will get.</p>
<p dir="auto">A difference between conda-forge and Homebrew that may be relevant is that the former uses <code class="notranslate">@rpath</code> and the latter a hardcoded path to load <code class="notranslate">libopenblas</code>:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="% # conda-forge
% otool -L _fblas.cpython-39-darwin.so
_fblas.cpython-39-darwin.so:
	@rpath/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.100.5)

% # Homebrew
% otool -L  /opt/homebrew/lib/python3.9/site-packages/scipy/linalg/_fblas.cpython-*-darwin.so
/opt/homebrew/lib/python3.9/site-packages/scipy/linalg/_fblas.cpython-39-darwin.so:
	/opt/homebrew/opt/openblas/lib/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/gcc/lib/gcc/11/libgfortran.5.dylib (compatibility version 6.0.0, current version 6.0.0)
	/opt/homebrew/opt/gcc/lib/gcc/11/libgcc_s.1.1.dylib (compatibility version 1.0.0, current version 1.1.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.100.5)"><pre class="notranslate"><code class="notranslate">% # conda-forge
% otool -L _fblas.cpython-39-darwin.so
_fblas.cpython-39-darwin.so:
	@rpath/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.100.5)

% # Homebrew
% otool -L  /opt/homebrew/lib/python3.9/site-packages/scipy/linalg/_fblas.cpython-*-darwin.so
/opt/homebrew/lib/python3.9/site-packages/scipy/linalg/_fblas.cpython-39-darwin.so:
	/opt/homebrew/opt/openblas/lib/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/opt/homebrew/opt/gcc/lib/gcc/11/libgfortran.5.dylib (compatibility version 6.0.0, current version 6.0.0)
	/opt/homebrew/opt/gcc/lib/gcc/11/libgcc_s.1.1.dylib (compatibility version 1.0.0, current version 1.1.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.100.5)
</code></pre></div>
<p dir="auto">That may not be the only difference, e.g. compilers used to build <code class="notranslate">libopenblas</code> and <code class="notranslate">scipy</code> were not the same. Also <code class="notranslate">libopenblas</code> can be built with either <code class="notranslate">pthreads</code> or <code class="notranslate">openmp</code> - numpy and scipy wheels use <code class="notranslate">pthreads</code>, while conda-forge and Homebrew both use <code class="notranslate">openmp</code>.</p>
<p dir="auto">To check if two <code class="notranslate">libopenblas</code> libraries get loaded, use:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="❯ python -m threadpoolctl -i scipy.linalg
[
  {
    &quot;user_api&quot;: &quot;blas&quot;,
    &quot;internal_api&quot;: &quot;openblas&quot;,
    &quot;prefix&quot;: &quot;libopenblas&quot;,
    &quot;filepath&quot;: &quot;/Users/ogrisel/mambaforge/envs/tmp/lib/python3.9/site-packages/numpy/.dylibs/libopenblas64_.0.dylib&quot;,
    &quot;version&quot;: &quot;0.3.18&quot;,
    &quot;threading_layer&quot;: &quot;pthreads&quot;,
    &quot;architecture&quot;: &quot;armv8&quot;,
    &quot;num_threads&quot;: 8
  },
  {
    &quot;user_api&quot;: &quot;blas&quot;,
    &quot;internal_api&quot;: &quot;openblas&quot;,
    &quot;prefix&quot;: &quot;libopenblas&quot;,
    &quot;filepath&quot;: &quot;/Users/ogrisel/mambaforge/envs/tmp/lib/python3.9/site-packages/scipy/.dylibs/libopenblas.0.dylib&quot;,
    &quot;version&quot;: &quot;0.3.17&quot;,
    &quot;threading_layer&quot;: &quot;pthreads&quot;,
    &quot;architecture&quot;: &quot;armv8&quot;,
    &quot;num_threads&quot;: 8
  }
]"><pre class="notranslate"><code class="notranslate">❯ python -m threadpoolctl -i scipy.linalg
[
  {
    "user_api": "blas",
    "internal_api": "openblas",
    "prefix": "libopenblas",
    "filepath": "/Users/ogrisel/mambaforge/envs/tmp/lib/python3.9/site-packages/numpy/.dylibs/libopenblas64_.0.dylib",
    "version": "0.3.18",
    "threading_layer": "pthreads",
    "architecture": "armv8",
    "num_threads": 8
  },
  {
    "user_api": "blas",
    "internal_api": "openblas",
    "prefix": "libopenblas",
    "filepath": "/Users/ogrisel/mambaforge/envs/tmp/lib/python3.9/site-packages/scipy/.dylibs/libopenblas.0.dylib",
    "version": "0.3.17",
    "threading_layer": "pthreads",
    "architecture": "armv8",
    "num_threads": 8
  }
]
</code></pre></div>
<h2 dir="auto">Context: why do 2 libopenblas copies get loaded</h2>
<p dir="auto">The reason is that the NumPy and SciPy wheels both vendor a copy of <code class="notranslate">libopenblas</code> within them, and extension modules that need <code class="notranslate">libopenblas</code> are depending directly on that vendored copy:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="% cd /path/to/site-packages/scipy/linalg
% otool -L _fblas.cpython-39-darwin.so 
_fblas.cpython-39-darwin.so:
	@loader_path/../.dylibs/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	@loader_path/../.dylibs/libgfortran.5.dylib (compatibility version 6.0.0, current version 6.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.60.1)

% cd ../../numpy/linalg
% otool -L _umath_linalg.cpython-39-darwin.so 
_umath_linalg.cpython-39-darwin.so:
	@loader_path/../.dylibs/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.60.1)"><pre class="notranslate"><code class="notranslate">% cd /path/to/site-packages/scipy/linalg
% otool -L _fblas.cpython-39-darwin.so 
_fblas.cpython-39-darwin.so:
	@loader_path/../.dylibs/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	@loader_path/../.dylibs/libgfortran.5.dylib (compatibility version 6.0.0, current version 6.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.60.1)

% cd ../../numpy/linalg
% otool -L _umath_linalg.cpython-39-darwin.so 
_umath_linalg.cpython-39-darwin.so:
	@loader_path/../.dylibs/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
	/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1292.60.1)
</code></pre></div>
<p dir="auto">This is how we have been shipping wheels for years, and it works fine across Windows, Linux and macOS. It seems like a weird thing to do of course (if you know how package managers work but are new to PyPI/wheels) - it's a long story, but the tl;dr is that PyPI wasn't designed with non-Python dependencies in mind, so the usual approach is to bundle those all into a wheel (it tends to work, unless you have complex non-Python dependencies). It'd be very much nontrivial to do any kind of unbundling here, and doing so would break situations where <code class="notranslate">numpy</code> and <code class="notranslate">scipy</code> are not installed in the same way (e.g., the former from conda-forge/Homebrew, the latter from PyPI).</p>
<h2 dir="auto">Possible root causes</h2>
<p dir="auto">The kernel panic had to do with spin locks apparently. It is not clear if the performance issues are also due to that, or have a completely different root cause. It does seem to be the case that two copies of the same shared library with the same version (all are <code class="notranslate">libopenblas.0.dylib</code>) cause a conflict at the OS level somehow. Anything beyond that is speculation at this point.</p>
<h2 dir="auto">Can we work around the problem?</h2>
<p dir="auto">If we release wheels for macOS 12, many people are going to hit this problem. A 50x slowdown for some code using linalg functionality for the default install configuration of <code class="notranslate">pip install numpy scipy</code> does not seem acceptable - that will lead too many users on wild goose chases. On the other hand it should be pointed out that if users build SciPy 1.7.2 from source on a native arm64 Python install, they will anyway hit the same problem. So not releasing any wheels isn't much better; at best it signals to users that they shouldn't use <code class="notranslate">arm64</code> just yet but stick with <code class="notranslate">x86_64</code> (but that does have some performance implications as well).</p>
<p dir="auto">At this point it looks like controlling the number of threads that OpenBLAS uses is the way we can work around this problem (or let users do so). Ways to control threading:</p>
<ul dir="auto">
<li>Use <code class="notranslate">threadpoolctl</code> (see the README at <a href="https://github.com/joblib/threadpoolctl">https://github.com/joblib/threadpoolctl</a> for how)</li>
<li>Set an environment variable to control the behavior, e.g. <code class="notranslate">OPENBLAS_NUM_THREADS</code></li>
<li>Rebuild the <code class="notranslate">libopenblas</code> we bundle in the wheel to have a max number of threads of 1, 2, or 4.</li>
</ul>
<p dir="auto">SciPy doesn't have a <code class="notranslate">threadpoolctl</code> runtime dependency, and it doesn't seem desirable to add one just for this issue. Note though that <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="948989375" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14441" data-hovercard-type="pull_request" data-hovercard-url="/scipy/scipy/pull/14441/hovercard" href="https://github.com/scipy/scipy/pull/14441">gh-14441</a> aims to add it as an <em>optional</em> dependency to improve test suite parallelism, and longer term we perhaps do want that dependency. Also, scikit-learn has a hard dependency on it, so many users will already have it installed.</p>
<p dir="auto">Rebuilding <code class="notranslate">libopenblas</code> with a low max number of threads does not allow users who know what they are doing or don't suffer from the problem to optimize threading behavior for their own code. It was pointed out in <a class="issue-link js-issue-link" data-error-text="Failed to load title" data-id="987205536" data-permission-text="Title is private" data-url="https://github.com/scipy/scipy/issues/14688" data-hovercard-type="issue" data-hovercard-url="/scipy/scipy/issues/14688/hovercard?comment_id=969143657&amp;comment_type=issue_comment" href="https://github.com/scipy/scipy/issues/14688#issuecomment-969143657">#14688 (comment)</a> that this is undesirable.</p>
<p dir="auto">Setting an environment variable is also not a great thing to do (a library should normally never ever do this), but if it works to do so in <code class="notranslate">scipy/__init__.py</code> then that may be the most pragmatic solution right now. However, this must be done before <code class="notranslate">libopenblas</code> is first loaded or it won't take effect. So if users import numpy first, then setting an env var will already have no effect on that copy of <code class="notranslate">libopenblas</code>. It needs testing whether this then still works around the problem or not.</p>
<p dir="auto"><em>Note: I wanted to have everything in one place, but let's discuss the release strategy on the mailing list (<a href="https://mail.python.org/archives/list/scipy-dev@python.org/thread/LLN2O4G2XI2MPILRW2XRRVCUK336WGKF/" rel="nofollow">link to thread</a>), and the actual performance issue here.</em></p>
<h2 dir="auto">Testing on other macOS arm64 build/install configurations</h2>
<p dir="auto"><em><strong>Request</strong>: if you have a build config on macOS arm64 that is not covered by the above summary yet, please run the following and reply on this issue with the results:</em></p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="% python -m threadpoolctl -i scipy.linalg

% cd /PATH/TO/scipy/linalg
% otool -L _fblas.cpython-*-darwin.so

% cd /PATH/TO/numpy/linalg
% otool -L _umath_linalg.cpython-*-darwin.so

% # Run the reproducer (again, only on macOS 12 - you will trigger an OS
% # crash on macOS 11.x!) and report if the time per `eigsh` call is ~0.02 sec. or ~1 sec.

% pip list    # if using pip for everything
% conda list  # if using conda"><pre class="notranslate"><code class="notranslate">% python -m threadpoolctl -i scipy.linalg

% cd /PATH/TO/scipy/linalg
% otool -L _fblas.cpython-*-darwin.so

% cd /PATH/TO/numpy/linalg
% otool -L _umath_linalg.cpython-*-darwin.so

% # Run the reproducer (again, only on macOS 12 - you will trigger an OS
% # crash on macOS 11.x!) and report if the time per `eigsh` call is ~0.02 sec. or ~1 sec.

% pip list    # if using pip for everything
% conda list  # if using conda
</code></pre></div>