← Back
Home › Systems › Numerical and mathematical computing

array-computing

Lets scientists write whole-array expressions and get near hand-tuned speed by splitting work between a high-level language and compiled numerical kernels.

NumPy, Julia, MATLAB

Open the array-computing pagearray-computing Auto vectorizationcompilers: Loop and SLP vectorization: dependence tests, cost model, turning scalar code into SIMDBytecode interpretercompilers: Register/stack VM dispatch loop for executing compiled plans or rules (SQLite VDBE, nftables)Dataflow analysiscompilers: Compute a fact per program point (live variables, reaching definitions, available expressions) by iterating transfer functions over the CFG to a fixed point.›FFI bindings generationcompilers: Describe a library's C-ABI surface once (header, attributes or IDL) and generate the declarations, converters and wrappers for each caller language.›Generational garbage collectioncompilers: Split the heap by age, collect the small young region often and cheaply, and rarely touch old data; wins whenever most objects die young.›JIT compilationcompilers: Run cold code in a cheap tier, count executions, and compile only what gets hot, trading warm-up time and memory for peak speed.›Reference countingcompilers: Keep a count of owners in each object, free it when the count hits zero; deterministic and local, but pays on every copy and leaks cycles.›mmap and the page cachehardware-os: How read() and mmap reach file data through the kernel page cache, what readahead and madvise do, and why databases use their own buffer pool instead.›SIMD programminghardware-os: Process 4-16 lanes per instruction with AVX2, AVX-512 or NEON via Rust std::arch, with runtime dispatch and a portable fallback.›Pythonlanguages: The fastest language to write and the home of ML and data science; slow as an interpreter, fast when C/CUDA libraries do the work.›PCA and SVDmachine-learning: Find orthogonal directions of maximum variance by an SVD of the centred data; randomized SVD gets the top k components in a few passes.›Kernel fusionml-training: Merge chains of memory-bound GPU ops into one kernel, so intermediates stay in registers instead of round-tripping through HBM.›Tiled matrix multiplicationml-training: Block GEMM into tiles held in shared memory and registers, so each byte loaded from DRAM feeds many FMAs; the core of fast GEMM.›Fast Fourier transformnumerics: Computes the DFT in O(n log n) by recursively splitting into even and odd halves; the basis of fast convolution and spectral methods.›Floating-point errornumerics: How IEEE 754 rounding error accumulates and how to contain it: compensated (Kahan/Neumaier) and pairwise summation, error-free transformations with FMA, ULP comparisons.›LU and Cholesky factorizationnumerics: Factor A = PLU (any square matrix, partial pivoting) or A = LLᵀ (SPD) once in O(n³), then solve each right-hand side in O(n²).›PRNGs (PCG, xoshiro)numerics: Small, fast, statistically strong non-cryptographic generators: PCG32 and xoshiro256**, with seeding, parallel streams and unbiased ranges.›Sparse matrices (CSR)numerics: Compressed sparse row stores only nonzeros as row pointers, column indices and values; SpMV streams them once, so it is memory-bandwidth bound.›Alias analysiscompilers: Decide whether two pointers may refer to the same memory (Andersen, Steensgaard, type-based) to enable optimizationsLoop optimizationscompilers: LICM, induction variables, strength reduction, unrolling, loop interchange and tilingControl-flow graphcompilers: Split a linear instruction list into basic blocks and connect them by jumps; the graph every optimizer, SSA builder and dataflow pass walks.›Topological sortgraphs: Order the vertices of a DAG so every edge points forward; the core of build systems, schedulers and dependency resolution, with cycle detection for free.›SSA formcompilers: Give every variable exactly one definition, merging at joins with phi functions; def-use becomes explicit and optimizations turn sparse and cheap.›Calling conventions abicompilers: System V, Windows x64 and AAPCS64 ABIs: argument registers, stack frames, struct passing, FFIIdl code generationcompilers: Generate typed client/server stubs and codecs from an interface definitionCopying gccompilers: Cheney semispace copying collection: bump allocation, forwarding pointers, compaction for freeMark-sweep garbage collectioncompilers: Find everything reachable from the roots, free the rest in place; handles cycles, never moves objects, and is the base of most production collectors.›Dirty page trackinghardware-os: Write-protection or PML logs of modified guest/process pagesDeoptimizationcompilers: Speculate with guards, and on failure rebuild the slow-tier frames from compile-time metadata and resume there, so optimized code may assume the common case.›Executable memory under W^Xcompilers: How JITs and code generators emit machine code under W^X: RW to RX flips, dual mappings, per-thread write toggles, and the instruction-cache flush.›+ 8 moreUnion-findgraphs: Disjoint-set forest with union by size and path compression; merges sets and answers same-set queries in effectively constant time.›Dominator treecompilers: Dominators via Lengauer-Tarjan or Cooper-Harvey-Kennedy; dominance frontiers for SSABFS and DFSgraphs: The two basic graph traversals; BFS gives hop distances and levels, DFS gives discovery/finish order for cycles, components and topological order.›Recursive-descent parsingcompilers: Hand-written top-down parser with one function per grammar rule; simple, fast, with full control over error messages and recovery.›Template enginescompilers: Compile templates with loops, partials and auto-escaping into render functionsProtobuf wire formatnetworking: Protocol Buffers wire format: tags, varints, length-delimited fields and schema evolutionArena allocatorsystems: Allocate by bumping a pointer through big chunks and free everything at once: a few instructions per allocation for objects sharing a lifetime.›Slab allocatorsystems: Pre-carve pages into equal-size slots per size class with an in-place free list: O(1) alloc and free, no external fragmentation, cache-warm reuse.›Stack unwindingcompilers: Frame-pointer, DWARF CFI and compact unwinding of call stacks+ 4 moreName resolutioncompilers: Bind every use of a name to its declaration with a scope stack, handling shadowing, hoisting and closure capture in one AST walk.›Html sanitizationsecurity: Allowlist-based HTML sanitizing of untrusted markupDebug info (DWARF)compilers: DWARF stores a tree of typed entries (DIEs), a line-number bytecode, unwind tables and location expressions that map machine state back to source; pick it for ELF/Mach-O debug builds.›Executable formats elf pecompilers: ELF, PE and Mach-O: sections, segments, symbol tables, relocations, dynamic loading+ 3 moreVarint (LEB128)compression: Variable-length integers in 7-bit groups with a continuation bit; small numbers take 1 byte, the wire format of Protobuf, WebAssembly and DWARF.›

How it is built

Architecture overview

An array-computing environment lets a scientist write C = A*B + 1 on whole matrices and get speed close to hand-tuned C. Its central object is the n-dimensional array: one flat buffer of numbers plus metadata (element type, shape, strides) that says how to read it. The one idea that makes it work is a split of labour: a slow, flexible language handles control flow and orchestration, and a small set of heavily optimized kernels (elementwise loops, matrix multiplication, factorizations, transforms, random numbers) handle the O(n) to O(n³) numeric work. Environments differ in how they bridge the two. NumPy and MATLAB keep the language interpreted and call precompiled kernels; Julia compiles the user's own loops to machine code, specialized on the argument types, so the user's code and the library code are equally fast.

The cost of the split is memory traffic: each whole-array operation reads and writes whole buffers, so chains of operations pay for temporaries. That makes fusion (one loop for a whole expression) the central optimization, done by the user (NumExpr-style or Julia dot syntax), by a compiler (JIT, GPU graph compilers), or by a library (lazy expressions in C++).

Parts:

 user program  (A, B, C = A*B .+ 1; loops; indexing)
      │
      ▼
 ┌─────────────────┐ interpreted      ┌──────────────────────────────────┐
 │ language front  │ calls ─────────► │ kernel library                   │
 │ end + runtime   │                  │ BLAS/LAPACK, FFT, sparse, RNG,   │
 │ (types, GC/RC)  │ compiled         │ ufunc loops (SIMD, dispatch)     │
 │                 │ loops ────────┐  └──────────────┬───────────────────┘
 └───────┬─────────┘               ▼                 │ strided buffers
         │                ┌────────────────┐         ▼
         │  array API     │ JIT: infer →   │   ┌───────────┐   ┌─────────┐
         └───────────────►│ SSA → opt →    │──►│ array     │◄──│ FFI:    │
                          │ vectorize →    │   │ buffer +  │   │ C, GPU  │
                          │ machine code   │   │ metadata  │   │ libs    │
                          └────────────────┘   └───────────┘   └─────────┘

How the parts fit together

Array object and views

The array is a data pointer, an element type, a shape and one stride per dimension (bytes to step along that axis). With strides, slicing, transposing, reversing, broadcasting (a stride of 0 repeats a value) and reshaping of compatible layouts all create a new metadata object over the same buffer; nothing is copied (NumPy documentation: "Copies and views"; the van der Walt, Colbert and Varoquaux paper). The cost is that every kernel must handle arbitrary strides, or first copy to a contiguous layout. Memory order is a language choice: NumPy defaults to row-major (C order), Julia and MATLAB use column-major (Fortran order, matching BLAS and LAPACK). Large arrays can be backed by a file (mmap and the page cache), giving arrays bigger than RAM; Virtual memory and the TLB explains why huge pages help large strided scans.

Language front end and execution tiers

Python is slow as an interpreter, so an array library designs for few calls with a lot of work each. The interpreter loop (Bytecode interpreter) dispatches A*B to the array type's operator, which finds a loop for the dtypes and runs it. Inline caches and, in JIT-ed Python implementations, Deoptimization speed up the interpretive overhead around the calls, but the arithmetic never runs in the interpreter. Julia reverses the design: a function is compiled the first time it is called with a given tuple of argument types (JIT compilation); type inference determines concrete types, and the compiler emits specialized machine code through LLVM. Writing a loop is then as fast as calling a library, so the library itself is written in the language. Code that defeats inference (values whose type depends on runtime data) falls back to boxed values and dynamic dispatch, which the Julia performance tips treat as the main thing to avoid.

Numba is a third route: compile selected Python numeric functions with LLVM on first call, specialized on argument types.

Compiler middle end and code generation

For compiled array code the compiler's job is to make the loop nest over strided data look like the loop a human would write. It starts from SSA form in a Control-flow graph, uses Dataflow analysis (type inference is itself a dataflow problem) and then:

  • Alias analysis (no page yet) must prove that the output array does not overlap the inputs; without it, the compiler cannot reorder loads and stores or vectorize. Languages give hints (Fortran forbids aliasing of arguments; restrict in C; Julia's @inbounds and @simd).
  • Loop optimizations (no page yet): invariant code motion, unrolling, interchange so the inner loop runs along the contiguous axis, tiling for caches.
  • Auto-vectorization (no page yet): turn the inner loop into SIMD instructions after dependence tests and a cost model.

Back end: Instruction selection, Register allocation, Machine code emission (no page yet) and, for a JIT, Executable memory under W^X for installing the result. All of these exist so that y .= a .* x .+ b becomes one tight loop with no temporaries.

Elementwise and reduction engine

Where there is no compiler in the loop (NumPy), the same effect is obtained by hand-written kernels: a ufunc is an inner loop over strided data, with versions for each dtype. NumPy builds several kernel variants per function from "universal intrinsics" and selects the best one the running CPU supports at start-up (NumPy documentation, "CPU/SIMD optimizations"); see SIMD programming for the dispatch pattern. Reductions deserve care: summing n numbers left to right accumulates rounding error that grows with n, so NumPy sums with partial pairwise summation along the fast axis (the numpy.sum notes say this is not guaranteed along other axes) (Floating-point error). The order of a reduction changes the result in the last bits, which matters when comparing runs. Expression evaluation order is a Topological sort of the expression graph when a library or compiler schedules fused expressions.

Fusion

Interpreter-driven environments allocate a temporary per operation: a*x + b creates a*x, then adds. Julia's dot syntax makes fusion part of the language: nested dotted calls such as sin.(cos.(X)) become a single loop and a single result allocation, and the manual states that this is a syntactic guarantee and not an optional optimization. On GPUs Kernel fusion does the same job for memory-bound elementwise chains and is the first thing a GPU array library such as CuPy or a graph compiler does. CUDA graphs (no page yet) record the resulting launch sequence once so repeated iterations pay one launch.

Dense linear algebra

A*B on matrices calls BLAS, not a loop in the host language: NumPy's matmul uses an optimized BLAS library when possible, and Julia does the same for * on dense matrices. The kernels follow the structure of Tiled matrix multiplication: block into panels that fit the cache levels, pack, and run a small register-blocked micro-kernel written for each CPU (Goto and van de Geijn 2008). OpenBLAS selects its kernels at run time for the CPU it finds when built with dynamic architecture support, and BLIS makes the micro-kernel the only architecture-specific piece. Higher routines such as A\b, inv, eigen and svd call LAPACK: LU and Cholesky factorization for solves and PCA and SVD for low-rank structure. The environment's contribution is the dispatch (symmetric positive definite? triangular? pick Cholesky or a triangular solve) and using a factorization instead of an explicit inverse.

Special-purpose libraries

Fast Fourier transform routines come from a planned library (see the FFT library system); Sparse matrices (CSR) store only nonzeros and have their own multiply and solve paths, with a different performance profile (memory-bound, not compute-bound); PRNGs (PCG, xoshiro) supply random numbers with seeding rules that allow independent parallel streams. NumPy's modern Generator uses PCG64 by default; the choice of generator is an explicit object, not a global state.

GPU path

The array API is mirrored on the device (CuPy mirrors NumPy; JAX and PyTorch do their own variants). Operations launch kernels written in CUDA and GPU languages, where performance depends on coalesced access and bank-conflict-free shared memory (GPU memory coalescing and banking, no page yet). Host-device transfers are slower than any single kernel, so the main rule is to keep arrays on the device between operations.

Memory management

Arrays are large, so when a buffer is freed matters. Reference counting (CPython, so NumPy) frees a buffer the instant the last reference disappears: deterministic and local, which keeps peak memory low, with cycles the unsolved case. Tracing collectors (Mark-sweep garbage collection, Generational garbage collection) free in batches. Julia's collector is non-moving, generational and mark-sweep (Julia developer documentation), so array buffers are never moved and can be handed to C directly. A moving collector (Copying garbage collection, no page yet) would force pinning for FFI. Copy-on-write semantics (MATLAB's value semantics, shared buffers that copy at first write) need write detection; Dirty-page tracking (no page yet) is the OS-level mechanism for it.

Foreign-function layer

Most kernels are C, C++ or Fortran. FFI bindings generation describes the call surface once and generates the glue; Calling conventions and ABI (no page yet) fix how arguments reach the callee, and IDL code generation (no page yet) covers generators for the declarations. The data contract is the array's pointer, shape, strides and dtype: Python's buffer protocol (PEP 3118) and the DLPack exchange format pass such metadata between libraries without copying, so a NumPy array can become a GPU library's tensor on the same memory.

Data and control flow

C = A @ B + 1 in NumPy (interpreted environment):

  1. The interpreter evaluates A @ B: dispatches to matmul, checks dtypes and shapes, allocates the result.
  2. matmul calls the BLAS gemm for the dtype; the library packs panels and runs the micro-kernel on all cores. Time goes to flops.
  3. + 1 creates a second temporary: a ufunc loop reads the whole matrix and writes the whole matrix. Time goes to memory bandwidth; for large arrays this pass can cost a noticeable fraction of the matrix product.
  4. The old temporary is released by reference counting as soon as nothing refers to it.

y .= a .* x .+ b in Julia (compiled environment):

  1. The dot syntax is lowered to one broadcast! call that fuses all three operations.
  2. The first call with these argument types triggers type inference, SSA optimization, LLVM and machine-code generation; later calls reuse the compiled method. First-call latency is seconds for big code, microseconds afterwards.
  3. The inner loop reads a, x, b, writes y, vectorized across elements where the compiler can prove independence.
  4. No temporary is allocated.

Calling a C library with an array:

  1. The caller passes the data pointer, length, strides and dtype (buffer protocol, DLPack, or a pointer plus size in Julia's ccall).
  2. The binding converts types, checks contiguity (copies if the library needs a different layout) and calls through the ABI.
  3. Results are written into a caller-provided buffer or a new array that adopts the library's allocation.

Scaling

In the order you hit the limits:

  1. Temporaries and passes over memory. Fuse expressions, use in-place operations, preallocate outputs, process in blocks that fit the cache.
  2. Kernel quality. Reach a tuned BLAS and LAPACK, match layouts to the kernels (contiguous, column-major for Fortran libraries), use the CPU's vector instructions.
  3. Threads. Parallel BLAS uses all cores by default; running several independent jobs on top of it oversubscribes the machine, so set the library's thread count or run one level of parallelism.
  4. Interpreter overhead. Move loops into array operations or into a compiled function; the signal is a profile dominated by the interpreter and not the kernels.
  5. Compilation latency. In JIT environments, cache compiled code, avoid recompiling because of changing types, and keep the hot kernels small and type-stable.
  6. Single-node memory and compute. Move to a GPU (keep data resident), memory-map files, process chunked arrays out of core, then distribute across nodes with a library that splits arrays into chunks and schedules tasks.

Failure modes

FailureSymptomCauseFix
Hidden temporariesmemory spikes, slow chains of operationseach operator allocates a full arrayfuse (dot syntax, compiled kernels), in-place operators, blocking
View/copy confusionwrites change another array, or changes are silently lostslices are views, fancy indexing copiesknow which operations return views; copy explicitly where needed
Silent dtype surprisesinteger overflow wraps, float32 results change, unexpected upcastfixed-width integer arithmetic, mixed-type promotionchoose dtypes explicitly, check ranges, test with small and large values
Broadcasting mistakewrong shape of result with no error (e.g. a column plus a row gives a matrix)rules extend size-1 and missing axesassert shapes at function boundaries
Type instability (JIT)code 10 to 100 times slower than expected (order of)inference cannot determine a concrete type, so values are boxedfind with the language's inference reports, make return types depend only on argument types
Wrong memory orderloops 5 to 10 times slower than the library routine (order of)inner loop runs against the contiguous axisiterate along the contiguous axis, transpose once, or call a library kernel
Thread oversubscriptionhigh CPU, low throughput, noisy timingslibrary threads times application threads exceed the coresone level of parallelism; set thread counts explicitly
Result differs between machineslast-digit differences, rarely a failed testdifferent BLAS kernels, reduction order, fused multiply-addcompare with tolerances; pin the BLAS if bitwise equality is a requirement
First-call latencyfirst request takes secondsJIT compilationwarm up at start, cache compiled code, ahead-of-time compile the hot kernels
GPU transfer dominatesdevice idle, copies in the profilearrays moved between host and device every operationkeep data on the device, fuse kernels

Real implementations

namewhat it isrepolicencenotes
NumPythe array type and ufunc/kernel library for Pythonnumpy/numpyBSD-3-Clausestrided views, SIMD dispatch, BLAS/LAPACK through a linked library
Julialanguage with type-specialized JIT; arrays and BLAS integration in the standard libraryJuliaLang/juliaMITdot-syntax fusion, LLVM code generation, non-moving generational GC
GNU Octaveinterpreted, MATLAB-compatible array languagegnu-octave/octaveGPL-3.0project site octave.org
MATLABproprietary array language with an execution engine and bundled math librariesmathworks.comproprietarycolumn-major arrays with value semantics
OpenBLASBLAS and LAPACK implementation with per-CPU kernels, from GotoBLAS2OpenMathLib/OpenBLASBSD-3-Clausethe default BLAS in many NumPy distributions
BLISframework for building BLAS-like libraries from one micro-kernel per architectureflame/blisBSD-3-Clause2020 SIAM Best Paper Prize for its design, per its README
NumbaLLVM JIT for numeric Python functionsnumba/numbaBSD-2-Clausetype-specialized compilation of NumPy-style loops
CuPyNumPy-compatible arrays on NVIDIA and AMD GPUscupy/cupyMITkernel fusion, vendor libraries
ndarrayn-dimensional arrays in Rustrust-ndarray/ndarrayMIT or Apache-2.0strided views, BLAS optional
xtensorC++ arrays with lazy expressions and NumPy-style semanticsxtensor-stack/xtensorBSD-3-Clauseexpression templates fuse loops at compile time

Sources

Similar systems

Share the most building blocks with this one.

eigensolver 24 %ml-compiler 17 %tensor-runtime 17 %statistics-environment 17 %circuit-sim 17 %