Array Morphisms

A unified backend for numerical computing combining fusion-based lazy evaluation, Mathematics of Arrays (MoA) index transformations, and automatic memory reuse.

Rationale

Just as the Scheme programming language demonstrates that a small set of orthogonal primitives suffices for all of computation, array morphisms demonstrate that index transformations suffice for most of what we call "matrix operations."

The proliferation of reshape, transpose, slice, and broadcast in numerical multi-dimensional array libraries such as NumPy is an indication that the underlying abstraction, the coordinate-free linear map, is not a first-class citizen. A matrix is a linear mapping plus an arbitrary choice of basis; by making the mapping itself first-class and the coordinate representation derived, operations that appear distinct collapse into a single compositional abstraction over index spaces.

The Mathematics of Arrays (MoA), developed by Lenore Mullin, provides the formal foundation for this approach. MoA treats arrays not as containers of data but as functions from index spaces to values, and array operations as transformations of those index spaces. Libraries such as NumPy already exploit this insight as an internal optimization -- transpose, slice, and reshape can return views backed by the same storage rather than copies. What MoA makes possible, and what array morphisms realize, is promoting these optimizations from private implementation details to first-class abstractions. Index transformation become objects that can be inspected, composed, and extended. Users are no longer dependent on library authors having anticipated a particular combination of operations. Any composition of affine index functions is itself an affine index function, and the index algebra guarantees this closure. MoA's reduction theory (the Psi Calculus) then provides the rules for simplifying chains of such transformations before any data is ever accessed.

Array morphisms are the realization of MoA's index-transformation algebra as first-class Scheme objects. An array morphism pairs a shape, which serves as a description of the index space, with an index function that maps coordinates in that space to coordinates in some underlying storage. Morphisms compose: applying a transpose to a slice to a reshape produces a single new morphism whose index function is the composition of three affine maps, with no intermediate allocation. The representation is stratified: at the bottom lies a flat, contiguous buffer; above it, any number of morphisms may be stacked, each a pure description of how to interpret that buffer. Morphism evaluation, the act of reading or writing data, occurs only when explicitly requested, and only at the coordinates actually needed.

Key Features

• Index Functions as First-Class Objects - Inspectable, composable, and extensible algebraic data
• Unified Abstraction - No split between "generalized" and "specialized" arrays; all operations compose uniformly
• Zero-Copy Views - Structural operations (reshape, transpose, slice) via affine index functions
• Memory Reuse - Two-phase trace/replay execution with buffer planning via graph coloring
• BLAS Integration - Transparent dispatch to optimized linear algebra kernels
• Element-wise Kernel Backends - Pluggable whole-array kernels for activations, binary arithmetic, reductions and copies, with results identical to the Scheme code they replace
• Type Safety - Multiple element types (f64, f32, s64, s32, u32, u64) with automatic promotion
• Category-Theoretic Foundation - Morphisms as structure-preserving transformations between arrays

Array Morphisms vs. SRFI-231

Index Functions as First-Class Objects

In SRFI-231, the indexer of a specialized array is available through `array-indexer` as a Scheme procedure; the library tracks its affine structure internally, but that structure is not exposed to the programmer as data. Array morphisms make index functions inspectable algebraic data -- tagged records with known structure (affine, window, compute, composed). This makes the composition law programmable.

Unified Abstraction Rather Than Two-Tier Hierarchy

SRFI-231 distinguishes between generalized arrays (domain plus opaque getter/setter) and specialized arrays (domain plus flat buffer plus affine indexer). Array morphisms replace this with a single stratified structure where structural and computational morphisms compose uniformly.

Algebraic Simplification

SRFI-231 already composes the affine indexers of specialized arrays, so a chain of operations such as `array-extract`, `array-permute`, `array-reverse`, and `specialized-array-reshape` yields a single affine indexer. Array morphisms extend this idea to symbolic rewriting over all morphism kinds, including non-affine window and compute morphisms: `reshape(reshape(x, s1), s2)` rewrites to `reshape(x, s2)`; `transpose(transpose(x, p1), p2)` rewrites to the composed permutation, and if the result is identity, the operations cancel entirely.

Zero-Copy Structural Operations

SRFI-231's `specialized-array-reshape` succeeds without copying whenever there is an affine indexer on the new domain that visits the elements of the existing array in the same order; the array need not be contiguous (see the commentary at https://srfi.schemers.org/srfi-231/commentary.html). When no such indexer exists, it signals an error or, if requested, copies the array. Array morphisms treat reshape as an index transformation that is always valid, possibly non-affine, and leave the decision of whether to materialize a copy to the execution strategy.

Deferred Evaluation with Buffer Reuse

SRFI-231 generalized arrays compute elements on demand, and the programmer decides when to materialize them, for example with `array-copy`. Array morphisms add a two-phase execution model with trace/replay and automatic buffer planning via liveness analysis.

Batch Awareness

SRFI-231 has no notion of batch dimension. Array morphisms track the batch axis as a first-class attribute, enabling batch-aware broadcasting, convolution, and reduction.

Type System

Array morphisms support multiple element types with automatic promotion:

'f64 - Double precision floating point (64-bit)

'f32 - Single precision floating point (32-bit)

's64 - 64-bit signed integer

's32 - 32-bit signed integer

'u64 - 64-bit unsigned integer

'u32 - 32-bit unsigned integer

Type promotion rules:

• Mixed operations promote to higher precision type
• Transcendental functions promote integers to floating point
• Reductions preserve dtype (mean promotes to float)

Core Concepts

Morphisms vs Arrays

Array morphisms represent computation as structure-preserving transformations. There are two types:

Concrete Arrays - Materialized data with shape, dtype, and strides

Abstract Morphisms - Deferred computations represented as expression trees

;; Concrete array - data is stored
(define concrete (morph-from-list '(1.0 2.0 3.0) '(3) 'f64))

;; Abstract morphism - represents computation
(define abstract (morph+ concrete concrete))

;; Realization materializes the morphism
(define result (realize abstract))  ; Now concrete

Index Functions

Index functions describe transformations algebraically:

Affine Index Functions - Pure transformations (reshape, transpose, slice)

Compute Index Functions - Element-wise arithmetic operations

Window Index Functions - Convolution and pooling operations

Reduction Index Functions - Aggregate operations (sum, mean, max)

Zero-Copy Views

Structural operations manipulate array views without copying:

(define x (morph-from-list '(0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0) '(8) 'f64))

;; Non-contiguous slice (stride 2)
(define strided (morph-slice x '(0) '(8) 2))  ; (0.0 2.0 4.0 6.0)

;; Reshape works even on non-contiguous views
(define as-2x2 (morph-reshape strided #(2 2)))  ; ((0.0 2.0) (4.0 6.0))

Installation

chicken-install array-morphisms

Module Imports

(import array-morphisms-core)       ; Core data types and utilities
(import array-morphisms-basic-ops)  ; Arithmetic operations
(import array-morphisms-structural-ops)  ; Reshape, transpose, slice
(import array-morphisms-realization)     ; Materialization
(import array-morphisms-context)         ; Memory reuse contexts
(import array-morphisms-batch-ops)       ; Batch operations
(import array-morphisms-blas-exec)       ; BLAS backend records and dispatch
(import array-morphisms-activation-exec) ; Element-wise kernel backends
(import array-morphisms-grad)            ; Reverse-mode automatic differentiation
(import array-morphisms-grad-check)      ; Numerical gradient checking
(import array-morphisms-ssa)             ; SSA compilation, VJP and replay

Core API

Morphism Creation

Create a concrete morphism from a nested list.

Create a concrete morphism from a typed vector.

Morphism Properties

Get the shape of morphism m.

Get the data type of morphism m.

Get the total number of elements in morphism m.

Get the number of dimensions in morphism m.

Check if morphism is concrete (materialized).

Check if morphism is abstract (deferred).

Check if morphism carries a batch dimension.

Get the batch size (error if not batched).

Realization

Materialize morphism to concrete array. Uses context if active.

Alias for realize.

Convert concrete morphism to nested Scheme list.

Arithmetic Operations

All operations are lazy by default and create morphism expressions.

Comparison Operations

Functional Operations

Apply function element-wise to morphism.

Reduce morphism along specified axes. Operations: 'sum, 'mean, 'max, 'min, 'prod.

Structural Operations

Reshape morphism to new shape (zero-copy). Supports -1 for automatic inference.

Transpose morphism by permuting dimensions. Default reverses all axes.

Extract subarray via affine transformation. Supports negative indices and step.

Stack morphisms along new dimension.

Concatenate morphisms along existing dimension.

Remove dimensions of size 1.

Add dimension of size 1 at specified axis.

Pad morphism with constant, edge, or reflect mode.

Matrix Operations (BLAS)

Matrix-matrix multiply morphism: (M,K) x (K,N) -> (M,N).

Matrix-vector multiply morphism: (M,N) x (N,) -> (M,).

Vector dot product morphism: (N,) . (N,) -> scalar.

AXPY morphism: result = alpha * x + y.

Convolution Operations

Image-to-column morphism for convolution.

Column-to-image morphism (adjoint of im2col).

Batch Operations

Add a size-1 batch dimension to non-batched morphism.

Remove size-1 batch dimension.

Extract a single element from the batch dimension (zero-copy).

Stack list of morphisms into a batch.

Apply function to each element along batch dimension.

Reduce across batch dimension.

Combine two batched morphisms element-wise.

Batched matrix multiply with automatic shape handling.

Memory Reuse Context

Create execution context for buffer reuse.

Realize morphism using context for buffer planning.

Complete trace phase and compute buffer assignment.

Reset context for replay phase.

Get statistics about allocations in context.

Display buffer assignment plan.

BLAS Backend Configuration

Matrix multiply, matrix-vector multiply, dot, axpy and conv2d are dispatched through a pluggable BLAS backend. Four tiers are available:

Pure Scheme

built in; used when no backend is registered.

microBLAS

built in (array-morphisms-micro-blas-backend), with no external dependencies. It is registered automatically at load time if no other backend has been registered first.

System BLAS

the separate egg array-morphisms-blas, which needs the blas egg and a system BLAS library. Opt-in.

crunch

the separate egg array-morphisms-crunch, whose GEMM and convolution kernels are compiled from Scheme to C by the crunch compiler and split across threads. Opt-in.

Register a BLAS backend for accelerated operations, replacing any previously registered backend.

Return the registered BLAS backend, or #f.

Check if BLAS optimization is enabled.

Enable BLAS optimization.

Disable BLAS optimization (use pure Scheme).

Check if a BLAS backend is registered.

Element-wise Kernel Backends

Module array-morphisms-activation-exec.

Element-wise operations, reductions and strided copies can be handed to a pluggable activation backend, registered separately from the BLAS backend. A backend holds whole-array kernels, one per op and dtype ('f32 or 'f64), for:

• the activations relu, sigmoid and tanh, and the derivative ops relu-deriv, sigmoid-deriv and tanh-deriv emitted by the backward pass;
• the binary ops add, sub, mul and div;
• reductions of row-major 2-D arrays over axis 0 or 1;
• copies of strided views into row-major order.

When an SSA replay plan is compiled, bindings with the unary and binary ops become ri-activation-unary or ri-activation-binary instructions, which use the backend's kernel for their op and dtype at execution time. Reductions and strided copies use their kernels wherever arrays are realized. When no backend is registered, no kernel exists for the op and dtype, or the operand and output dtypes differ, the same Scheme code runs as without a backend. Kernels must compute exactly what that Scheme code computes, so a backend changes speed, not results. Bindings produced by the element-wise fusion pass carry a composed function and never use a kernel.

No activation backend is registered by default. The array-morphisms-crunch egg provides an unthreaded and a threaded backend; see the example below.

Create an activation backend called name (a symbol) with no kernels.

Install kernel, a procedure (size in out), as the implementation of unary op on dtype arrays. in and out are SRFI-4 vectors of that dtype and may be the same vector. A later call for the same op and dtype replaces the kernel.

Install kernel, a procedure (size a b out) computing out[i] = op(a[i], b[i]), as the implementation of binary op on dtype arrays. out may be a or b.

Install kernel, a procedure (rows cols src out), as the reduction rop (sum, mean, max or min) over axis (0 or 1) of a row-major rows x cols array. out receives cols results for axis 0 and rows results for axis 1.

Install kernel, a procedure (src offset shape strides dst), as the copy of strided dtype arrays of rank at most 4 into row-major order.

List of the (op . dtype) pairs that backend has unary kernels for.

Make backend the active activation backend. #f removes it.

The active activation backend, or #f.

Examples

Basic Operations

(import array-morphisms-core
        array-morphisms-basic-ops
        array-morphisms-structural-ops
        array-morphisms-realization)

;; Create arrays
(define a (morph-from-list '(1.0 2.0 3.0 4.0 5.0) '(5) 'f64))
(define b (morph-from-list '(2.0 4.0 6.0 8.0 10.0) '(5) 'f64))

;; Build lazy computation
(define squared (morph-map (lambda (x) (* x x)) a))
(define complex (morph+ squared (morph-map (lambda (x) (* x 2)) b)))

;; Force computation
(morph->list (realize complex))

Layer Normalization

(define (layer-norm x eps)
  (let* ((mean (morph-reduce 'mean x '(0)))
         (centered (morph- x mean))
         (variance (morph-reduce 'mean 
                               (morph* centered centered) 
                               '(0)))
         (std (morph-sqrt (morph+ variance eps))))
    (morph/ centered std)))

;; Create 4x6 activation matrix
(define activations
  (morph-from-list (map (lambda (i) (+ i 1.0)) (iota 24))
                   #(4 6) 'f64))

;; Apply layer normalization
(define normalized (layer-norm activations 1e-5))
(morph->list (realize normalized))

Memory Reuse with Context

(import array-morphisms-context)

;; Create context for repeated computation
(define ctx (make-morphism-context))

;; Trace phase
(define result1 (realize/ctx ctx expensive-morphism))
(finalize-context! ctx)

;; Replay phase - reuses buffers
(reset-context! ctx)
(define result2 (realize/ctx ctx expensive-morphism))

Batched Operations

(import array-morphisms-batch-ops)

;; Create batch of matrices
(define m1 (morph-from-list '(1.0 2.0 3.0 4.0) '(2 2) 'f64))
(define m2 (morph-from-list '(2.0 4.0 6.0 8.0) '(2 2) 'f64))
(define batch (morph-stack (list m1 m2) 0))  ; shape (2 2 2)

;; Batch map: apply function to each matrix
(define doubled (batch-map 
                  (lambda (m) (morph-map (lambda (x) (* x 2)) m))
                  batch))

;; Batch reduce: sum across batch dimension
(define summed (batch-reduce 'sum doubled))

Zero-Copy Structural Chain

(define x (morph-from-list '(0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 
                             8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0)
                           '(16) 'f64))

;; Stride-2 downsample
(define x-even (morph-slice x '(0) '(16) 2))      ; shape (8)

;; Stride-2 again - composed becomes stride-4
(define x-quarter (morph-slice x-even '(0) '(8) 2)) ; shape (4)

;; Reshape and transpose - all zero-copy views
(define as-2x2 (morph-reshape x-quarter #(2 2)))
(define transposed (morph-transpose as-2x2 '(1 0)))

;; Only when we realize is data accessed
(morph->list (realize transposed))

BLAS Integration

(import array-morphisms-blas-exec)

;; Matrix multiply automatically uses BLAS when available
(define A (morph-from-list '(1.0 2.0 3.0 4.0 5.0 6.0) '(2 3) 'f64))
(define B (morph-from-list '(7.0 8.0 9.0 10.0 11.0 12.0) '(3 2) 'f64))

;; Dispatches to the GEMM of the registered backend; the built-in
;; microBLAS backend is registered by default
(define C (morph-matmul A B))

;; Optionally replace it with the system-BLAS backend
;; (egg array-morphisms-blas) ...
(import array-morphisms-blas-egg-backend)
(register-blas-backend! (make-blas-egg-backend))

;; ... or with the crunch-compiled backend (egg array-morphisms-crunch)
;; (import array-morphisms-crunch-conv-backend)
;; (register-blas-backend! (make-crunch-blas-backend))

(morph->list (realize C))

Element-wise Kernels

(import array-morphisms-activation-exec
        array-morphisms-crunch-activations   ; egg array-morphisms-crunch
        array-morphisms-crunch-threads)

;; Activations, their derivatives, add/sub/mul/div, 2-D reductions and
;; strided copies now run as whole-array kernels split across threads
(register-activation-backend! (make-crunch-threaded-activation-backend))

;; Threading settings, read at every call
(crunch-thread-count 8)          ; threads requested
(crunch-thread-min-chunk 16384)  ; fewest elements per thread

;; Back to the Scheme code, with identical results
(register-activation-backend! #f)

Performance Notes

• Lazy Evaluation: Operations build expression trees, materialize once
• Zero-Copy Views: Structural operations (reshape, transpose, slice) are free
• Memory Reuse: Use contexts for repeated computations to enable buffer planning
• Batch Operations: Process multiple arrays together for efficiency
• BLAS Dispatch: Automatically uses optimized kernels for large operations
• Element-wise Kernels: Register an activation backend to run activations, binary arithmetic and reductions as compiled, optionally threaded, whole-array kernels

Author

Ivan Raikov

Repository

https://github.com/iraikov/array-morphisms

Requirements

• datatype
• matchable
• random-mtzig
• srfi-1
• srfi-4
• srfi-69

Optional companion eggs:

• array-morphisms-blas: system-BLAS backend (needs a BLAS library)
• array-morphisms-crunch: crunch-compiled BLAS, convolution and element-wise kernels (needs the crunch egg)

License

LGPL-3

References

Mullin, Lenore M.R. A Mathematics of Arrays. Ph.D. Thesis, Syracuse University, 1988.

The MoA formalism provides the theoretical foundation for index transformations as compositional abstractions.