<p dir="auto">Run this script with <code class="notranslate">ulimit -v 1000000</code> (limit memory use to 1GB)</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="## dumb_convolve.py

import numpy as np
from numpy.lib.stride_tricks import as_strided

N = 16384
dt = np.dtype(float)
dtsz = dt.itemsize

# Assume these vectors contained actual data, not just ones
vect1 = np.ones((2*N-1,), dtype=dt)
vect2 = np.ones((N,), dtype=dt)

arr1 = as_strided(vect1, shape=(N, N), strides=(dtsz, dtsz))

vect3 = np.dot(arr1, vect2)  # MemoryError here"><pre class="notranslate"><span class="pl-c">## dumb_convolve.py</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">numpy</span>.<span class="pl-s1">lib</span>.<span class="pl-s1">stride_tricks</span> <span class="pl-k">import</span> <span class="pl-s1">as_strided</span>

<span class="pl-v">N</span> <span class="pl-c1">=</span> <span class="pl-c1">16384</span>
<span class="pl-s1">dt</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">dtype</span>(<span class="pl-s1">float</span>)
<span class="pl-s1">dtsz</span> <span class="pl-c1">=</span> <span class="pl-s1">dt</span>.<span class="pl-s1">itemsize</span>

<span class="pl-c"># Assume these vectors contained actual data, not just ones</span>
<span class="pl-s1">vect1</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">2</span><span class="pl-c1">*</span><span class="pl-v">N</span><span class="pl-c1">-</span><span class="pl-c1">1</span>,), <span class="pl-s1">dtype</span><span class="pl-c1">=</span><span class="pl-s1">dt</span>)
<span class="pl-s1">vect2</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-v">N</span>,), <span class="pl-s1">dtype</span><span class="pl-c1">=</span><span class="pl-s1">dt</span>)

<span class="pl-s1">arr1</span> <span class="pl-c1">=</span> <span class="pl-en">as_strided</span>(<span class="pl-s1">vect1</span>, <span class="pl-s1">shape</span><span class="pl-c1">=</span>(<span class="pl-v">N</span>, <span class="pl-v">N</span>), <span class="pl-s1">strides</span><span class="pl-c1">=</span>(<span class="pl-s1">dtsz</span>, <span class="pl-s1">dtsz</span>))

<span class="pl-s1">vect3</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">dot</span>(<span class="pl-s1">arr1</span>, <span class="pl-s1">vect2</span>)  <span class="pl-c"># MemoryError here</span></pre></div>
<p dir="auto"><code class="notranslate">np.dot</code> tries to allocate a 16384*16384 array and fails with MemoryError</p>
<p dir="auto">For cases like this, even when there is enough memory to allocate the big array, it would be much faster (re: cache etc) to do the calculation using the existing array and an unoptimised algorithm than copying to a big array then using an optimised BLAS.</p>
<p dir="auto">While trying to test how <code class="notranslate">dot</code> would behave on a very large array without actually constructing one, I found this discrepancy:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; x = as_strided(array([1]), shape=(2**33,), strides=(0,))
&gt;&gt;&gt; x.dtype
dtype('int64')
&gt;&gt;&gt; x = as_strided(x, shape=(2**33,), strides=(0,))
&gt;&gt;&gt; np.dot(x, x)
8589934592"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; x = as_strided(array([1]), shape=(2**33,), strides=(0,))
&gt;&gt;&gt; x.dtype
dtype('int64')
&gt;&gt;&gt; x = as_strided(x, shape=(2**33,), strides=(0,))
&gt;&gt;&gt; np.dot(x, x)
8589934592
</code></pre></div>
<p dir="auto">This is the expected result, but on floats it fails:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; x = as_strided(array([1.]), shape=(2**33,), strides=(0,))
&gt;&gt;&gt; np.dot(x, x)
Traceback (most recent call last):
  File &quot;&lt;ipython-input-37-7170710dec0d&gt;&quot;, line 1, in &lt;module&gt;
    np.dot(x, x)
ValueError: matrices are not aligned"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; x = as_strided(array([1.]), shape=(2**33,), strides=(0,))
&gt;&gt;&gt; np.dot(x, x)
Traceback (most recent call last):
  File "&lt;ipython-input-37-7170710dec0d&gt;", line 1, in &lt;module&gt;
    np.dot(x, x)
ValueError: matrices are not aligned
</code></pre></div>
<p dir="auto">This is unfortunate, since it makes it hard to test whether <code class="notranslate">dot</code> behaves correctly on arrays of size &gt; 2**31-1.</p>