NumPy is a Python library that stores numerical data in contiguous, typed memory buffers (the ndarray) and dispatches arithmetic over those buffers to optimized C and Fortran kernels—BLAS for linear algebra, SIMD-vectorized loops for element-wise work—so that a single Python function call can process millions of values at hardware memory bandwidth rather than at Python interpreter speed.

<svg viewBox="0 0 760 470" xmlns="http://www.w3.org/2000/svg" font-family="ui-monospace,monospace">
  <rect width="760" height="470" fill="#0d1117"/>
  <text x="380" y="32" text-anchor="middle" fill="#e6edf3" font-size="21" font-weight="700">NumPy ndarray Memory Model</text>
  <text x="380" y="52" text-anchor="middle" fill="#8b98a5" font-size="13">Strides, views, and SIMD dispatch — why the same data can look like different arrays</text>

  <!-- Raw memory buffer -->
  <text x="30" y="88" fill="#8b98a5" font-size="12">Contiguous float64 buffer (8 bytes/element)</text>
  <g transform="translate(30,96)">
    <!-- 8 cells -->
    <rect x="0"   y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#1f6feb" stroke-width="1.5"/>
    <rect x="62"  y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#1f6feb" stroke-width="1.5"/>
    <rect x="124" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#1f6feb" stroke-width="1.5"/>
    <rect x="186" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#1f6feb" stroke-width="1.5"/>
    <rect x="248" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#3fb950" stroke-width="1.5"/>
    <rect x="310" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#3fb950" stroke-width="1.5"/>
    <rect x="372" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#3fb950" stroke-width="1.5"/>
    <rect x="434" y="0" width="60" height="36" rx="3" fill="#161b22" stroke="#3fb950" stroke-width="1.5"/>
    <text x="30"  y="23" text-anchor="middle" fill="#58a6ff" font-size="11">0</text>
    <text x="92"  y="23" text-anchor="middle" fill="#58a6ff" font-size="11">8</text>
    <text x="154" y="23" text-anchor="middle" fill="#58a6ff" font-size="11">16</text>
    <text x="216" y="23" text-anchor="middle" fill="#58a6ff" font-size="11">24</text>
    <text x="278" y="23" text-anchor="middle" fill="#3fb950" font-size="11">32</text>
    <text x="340" y="23" text-anchor="middle" fill="#3fb950" font-size="11">40</text>
    <text x="402" y="23" text-anchor="middle" fill="#3fb950" font-size="11">48</text>
    <text x="464" y="23" text-anchor="middle" fill="#3fb950" font-size="11">56</text>
    <text x="252" y="50" fill="#6e7681" font-size="10">byte offset</text>
  </g>

  <!-- View A: shape (2,4), strides (32,8) -->
  <text x="30" y="176" fill="#8b98a5" font-size="12">A = buf.reshape(2,4)  strides=(32, 8)  — zero copy</text>
  <g transform="translate(30,184)">
    <rect x="0"   y="0" width="120" height="30" rx="3" fill="#1f6feb" fill-opacity="0.15" stroke="#1f6feb" stroke-width="1"/>
    <rect x="122" y="0" width="120" height="30" rx="3" fill="#1f6feb" fill-opacity="0.15" stroke="#1f6feb" stroke-width="1"/>
    <rect x="244" y="0" width="120" height="30" rx="3" fill="#1f6feb" fill-opacity="0.15" stroke="#1f6feb" stroke-width="1"/>
    <rect x="366" y="0" width="120" height="30" rx="3" fill="#1f6feb" fill-opacity="0.15" stroke="#1f6feb" stroke-width="1"/>
    <rect x="0"   y="32" width="120" height="30" rx="3" fill="#3fb950" fill-opacity="0.15" stroke="#3fb950" stroke-width="1"/>
    <rect x="122" y="32" width="120" height="30" rx="3" fill="#3fb950" fill-opacity="0.15" stroke="#3fb950" stroke-width="1"/>
    <rect x="244" y="32" width="120" height="30" rx="3" fill="#3fb950" fill-opacity="0.15" stroke="#3fb950" stroke-width="1"/>
    <rect x="366" y="32" width="120" height="30" rx="3" fill="#3fb950" fill-opacity="0.15" stroke="#3fb950" stroke-width="1"/>
    <text x="60"  y="20" text-anchor="middle" fill="#58a6ff" font-size="11">A[0,0]</text>
    <text x="182" y="20" text-anchor="middle" fill="#58a6ff" font-size="11">A[0,1]</text>
    <text x="304" y="20" text-anchor="middle" fill="#58a6ff" font-size="11">A[0,2]</text>
    <text x="426" y="20" text-anchor="middle" fill="#58a6ff" font-size="11">A[0,3]</text>
    <text x="60"  y="52" text-anchor="middle" fill="#3fb950" font-size="11">A[1,0]</text>
    <text x="182" y="52" text-anchor="middle" fill="#3fb950" font-size="11">A[1,1]</text>
    <text x="304" y="52" text-anchor="middle" fill="#3fb950" font-size="11">A[1,2]</text>
    <text x="426" y="52" text-anchor="middle" fill="#3fb950" font-size="11">A[1,3]</text>
  </g>

  <!-- View B: A.T, strides swapped to (8,32) -->
  <text x="30" y="286" fill="#8b98a5" font-size="12">B = A.T  strides=(8, 32)  — zero copy, same buffer</text>
  <g transform="translate(30,294)">
    <rect x="0"  y="0"  width="56" height="64" rx="3" fill="#161b22" stroke="#f0883e" stroke-width="1.5"/>
    <rect x="58" y="0"  width="56" height="64" rx="3" fill="#161b22" stroke="#f0883e" stroke-width="1.5"/>
    <rect x="116" y="0" width="56" height="64" rx="3" fill="#161b22" stroke="#f0883e" stroke-width="1.5"/>
    <rect x="174" y="0" width="56" height="64" rx="3" fill="#161b22" stroke="#f0883e" stroke-width="1.5"/>
    <text x="28"  y="26" text-anchor="middle" fill="#f0883e" font-size="10">B[0,0]</text>
    <text x="28"  y="48" text-anchor="middle" fill="#f0883e" font-size="10">B[0,1]</text>
    <text x="86"  y="26" text-anchor="middle" fill="#f0883e" font-size="10">B[1,0]</text>
    <text x="86"  y="48" text-anchor="middle" fill="#f0883e" font-size="10">B[1,1]</text>
    <text x="144" y="26" text-anchor="middle" fill="#f0883e" font-size="10">B[2,0]</text>
    <text x="144" y="48" text-anchor="middle" fill="#f0883e" font-size="10">B[2,1]</text>
    <text x="202" y="26" text-anchor="middle" fill="#f0883e" font-size="10">B[3,0]</text>
    <text x="202" y="48" text-anchor="middle" fill="#f0883e" font-size="10">B[3,1]</text>
    <text x="310" y="36" fill="#8b98a5" font-size="11">B.shape = (4, 2)</text>
    <text x="310" y="54" fill="#8b98a5" font-size="11">B.strides = (8, 32) bytes</text>
    <text x="310" y="72" fill="#6e7681" font-size="10">0 bytes allocated</text>
  </g>

  <!-- SIMD dispatch bar -->
  <text x="30" y="384" fill="#8b98a5" font-size="12">Element-wise throughput — Python loop vs NumPy SIMD</text>
  <text x="30" y="400" fill="#6e7681" font-size="10">Python list sum (~100 ns/elem)</text>
  <rect x="200" y="390" width="480" height="14" rx="2" fill="#f85149" fill-opacity="0.7"/>
  <text x="688" y="401" fill="#f85149" font-size="10">200×</text>
  <text x="30" y="420" fill="#6e7681" font-size="10">NumPy AVX-512 (~0.5 ns/elem)</text>
  <rect x="200" y="410" width="2.4" height="14" rx="2" fill="#3fb950" fill-opacity="0.9"/>
  <text x="210" y="421" fill="#3fb950" font-size="10">1×</text>

  <text x="380" y="460" text-anchor="middle" fill="#6e7681" font-size="11">A.T, A[::2], A[1:] all share the same buffer — slicing allocates no memory</text>
</svg>

NumPy's performance advantage over pure Python is not primarily algorithmic—it is an interpreter dispatch cost measured at approximately 100 ns per element in a Python for-loop versus 0.5 ns per element in a vectorized NumPy kernel, a 200× ratio that compounds across any loop over a large array. The ndarray stores elements in a single contiguous block of typed memory: a 1,000,000-element float64 array occupies exactly 7.6 MB (8 bytes × 10^6) with no Python object overhead per element, while an equivalent Python list holds 1,000,000 pointer-sized references (8 bytes each) plus a separate float object (24 bytes each) on the heap—28× more memory and zero cache locality.

The stride mechanism is what makes reshape, transpose, and slicing allocation-free. Every ndarray carries a strides tuple—one integer per dimension—specifying how many bytes to advance in the buffer to step one index forward. Calling .T on a (10,000 × 10,000) float64 array, which holds 0.75 GB of data, costs zero bytes of allocation: NumPy swaps the stride tuple from (80000, 8) to (8, 80000) and returns a view pointing at the same buffer. Similarly, arr[::2] sets stride to 2 × itemsize (16 bytes for float64) with no copy; only operations like np.ascontiguousarray() or boolean fancy indexing force a materialization.

Broadcasting eliminates explicit loops over shape-mismatched arrays by stretching dimensions of size 1 virtually rather than physically. Adding a (1000, 1) column vector to a (1, 1000) row vector produces a (1000, 1000) result (7.6 MB) without allocating any intermediate array for the broadcasts; NumPy's inner loop strides over the virtual dimension with a step of zero bytes, reusing the same 8-byte value 1,000 times per row. The three broadcasting rules—pad shapes left with ones, treat size-1 dimensions as repeating, require all other dimensions to match—make this zero-copy expansion predictable and general across any number of dimensions.

The BLAS and LAPACK libraries, not NumPy itself, supply the peak floating-point throughput for matrix operations. When np.dot() or the @ operator detects a 2D float array, it dispatches to the linked BLAS DGEMM routine—OpenBLAS, MKL, or BLIS depending on how NumPy was compiled. MKL's DGEMM on a modern CPU reaches ~82 GFLOPs for a (1024 × 1024) matrix multiply (2.15 GFLOPs of arithmetic) using AVX-512 FMA instructions that process 16 float32 or 8 float64 values per clock cycle. NumPy's own element-wise kernels also use AVX-512 where available, packing 16 float32 values into a single 512-bit register to achieve peak SIMD utilization.

Fancy indexing and boolean masking are the two cases where NumPy must allocate a new buffer. Integer array indexing—arr[[0, 5, 3, 2]]—gathers non-contiguous elements into a fresh contiguous allocation because the result cannot be expressed as a stride pattern. Boolean indexing—arr[arr > 0]—similarly materializes the selected elements since their count is data-dependent. Knowing this boundary prevents the common mistake of applying fancy indexing inside a tight loop, where the allocation cost accumulates; a single boolean mask followed by a scalar assignment (arr[mask] = 0) is an in-place operation that avoids the gather allocation entirely.

Numba and Cython extend NumPy's reach to cases the vectorized API cannot express. Numba's @jit(nopython=True) compiles Python loops over ndarrays to LLVM IR at first call, typically matching or exceeding hand-written C for stencil operations, custom reductions, and recurrences where each output depends on the previous one—patterns that have no NumPy primitive. Cython typed memoryviews (double[::1]) provide the same access with static dispatch and no JIT latency. For GPU execution, CuPy mirrors the NumPy API but dispatches to CUDA kernels; the bottleneck then shifts to PCIe Gen4 ×16 transfer bandwidth (~32 GB/s), so a 1 GB array transfer takes ~31 ms—enough to erase any GPU speedup unless the kernel runtime dominates.

NUMPY DISPATCH FLOWCHART

  Python call: result = A + B
         │
         ▼
  ┌─────────────────────┐
  │  __add__ ufunc      │  checks dtype, shape, broadcast rules
  │  np.add(A, B)       │
  └────────┬────────────┘
           │
     ┌─────┴─────┐
     │           │
  contiguous?  strided/
     │YES      broadcast?
     ▼           │YES
  ┌──────────┐   ▼
  │ AVX-512  │ ┌─────────────────┐
  │ kernel   │ │ inner-loop with │
  │ 16 fp32  │ │ stride=0 reuse  │
  │ per instr│ └────────┬────────┘
  └────┬─────┘          │
       └────────┬────────┘
                ▼
  ┌─────────────────────┐
  │  Output ndarray     │  contiguous C-order buffer
  │  (new allocation)   │  shape = broadcast(A,B).shape
  └─────────────────────┘

Read NumPy through a stride arithmetic lens rather than a fast math library lens. The entire API—reshape, slice, transpose, broadcast, fancy-index—is a set of operations on the (data pointer, shape, strides, dtype) quadruple, and every performance characteristic follows from whether the requested operation can be expressed by changing that quadruple without touching the underlying bytes. The moment an operation cannot—gather indexing, boolean selection, contiguity coercion—NumPy allocates, and that allocation cost is the right place to focus when profiling a NumPy-heavy pipeline.