Skip to content

Structure and laziness

Two ideas do most of the work in linox: arithmetic builds a description rather than a result, and that description is rewritten into whatever structure is most useful.

Laziness

A + B does not add anything:

import jax.numpy as jnp
import linox

a = linox.Diagonal(jnp.arange(1.0, 5.0))
b = linox.Diagonal(jnp.ones(4))

total = a + b
assert isinstance(total, linox.Diagonal)   # rewritten, not deferred

Two diagonals can be added cheaply, so they are. Where no such rule exists, the sum is held as a composite and evaluated only on application:

import jax
import jax.numpy as jnp
import linox

key = jax.random.PRNGKey(0)
m = linox.Matrix(jax.random.normal(key, (4, 4)))
k = linox.Kronecker(linox.Matrix(jnp.eye(2)), linox.Matrix(jnp.eye(2)))

composite = m @ k
assert type(composite).__name__ == "ProductLinearOperator"

# Applied right to left; no 4x4 product is formed.
assert (composite @ jnp.ones(4)).shape == (4,)

Rewriting

Sums are inspected for patterns worth specialising. The most important is an isotropic shift — a matrix plus a multiple of the identity, which is what regularisation and GP jitter look like:

import jax
import jax.numpy as jnp
import linox

key = jax.random.PRNGKey(0)
dense = jax.random.normal(key, (4, 4))
spd = linox.Matrix(dense @ dense.T + 4 * jnp.eye(4))

shifted = spd + 0.1 * linox.Identity(4)
assert type(shifted).__name__ == "IsotropicAdditiveLinearOperator"

That is not cosmetic. The rewritten operator caches an eigendecomposition of A and serves inverse, square root, log, determinant and powers from it, since the eigenvectors of A and s·I + A coincide and only the eigenvalues shift.

A diagonal plus a symmetric low-rank term is recognised similarly and solved by the Woodbury identity rather than by forming the sum.

The rewrites are meaning-preserving:

import jax
import jax.numpy as jnp
import linox

key = jax.random.PRNGKey(0)
dense = jax.random.normal(key, (4, 4))
raw = dense @ dense.T + 4 * jnp.eye(4)

rewritten = linox.Matrix(raw) + 0.5 * linox.Identity(4)
plain = linox.Matrix(raw + 0.5 * jnp.eye(4))

assert jnp.allclose(linox.todense(rewritten), linox.todense(plain))

Structure is a promise about the operator

Some structure is intrinsic — a Diagonal is diagonal by construction. Some is asserted by you:

import jax
import jax.numpy as jnp
import linox

key = jax.random.PRNGKey(0)
dense = jax.random.normal(key, (4, 4))
spd = dense @ dense.T + 4 * jnp.eye(4)

asserted = linox.PSD(linox.Matrix(spd))
assert asserted.is_psd
assert asserted.is_symmetric

Sym, PSD and SPD attach a property without checking it. They are promises, and wrapping does not copy or materialise anything — the wrapper delegates every matvec to the operator underneath.

When laziness ends

A composite collapses to numbers when you apply it, solve with it, or ask for its dense form. Everything before that is bookkeeping — which is why building a large expression costs nothing until it is used.