Lenia Gradients in Depth 
¶
The companion to 65 - Lenia Gradients. That notebook uses the gradient — to make a creature faster, to grow a picture, to kill a soliton. This one asks how the gradient behaves:
- Is it real? Autodiff differentiates the program it is given, and a program can be discontinuous while still returning a plausible number.
- How far back can it see before chaos destroys it?
- Where does it go blind, and why?
- What can it never distinguish, no matter how much data it is given?
- Where does it point in space?
- What does colour cost?
None of this is needed to optimize a Lenia rule. All of it explains why those optimizations behave the way they do.
Installation¶
You will need Python 3.12 or later, and a working JAX installation. For example, you can install JAX with:
%pip install -U "jax[cuda]"
/home/faldor_google_com/dev/cax/.venv/bin/python3: No module named pip
Note: you may need to restart the kernel to use updated packages.
Then, install CAX from PyPi:
%pip install -U "cax[examples]"
/home/faldor_google_com/dev/cax/.venv/bin/python3: No module named pip
Note: you may need to restart the kernel to use updated packages.
Import¶
import contextlib
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import optax
from cax.cs.lenia import (
FreeKernelParams,
Lenia,
LeniaGrowthParams,
LeniaKernelParams,
LeniaRuleParams,
exponential_growth_fn,
exponential_kernel_fn,
free_kernel_fn,
load_pattern,
)
Configuration¶
seed = 0
spatial_dims = (96, 96)
channel_size = 1
R = 13
T = 10
state_scale = 1
key = jax.random.key(seed)
Setup¶
The same helpers as notebook 54: build a system and place a pattern. Simulations run
through the system's driver, cs(state, num_steps=...).
def make_lenia(rule_params, **kwargs):
"""Build a Lenia system, with keyword overrides for the shared configuration."""
options = {
"spatial_dims": spatial_dims,
"channel_size": channel_size,
"R": R,
"T": T,
"state_scale": state_scale,
}
options.update(kwargs)
return Lenia(rule_params=rule_params, **options)
def sample_state(pattern, num_channels=channel_size, dims=spatial_dims):
"""Place a pattern at the center of an empty grid."""
mid = tuple(dim // 2 for dim in dims)
slices = tuple(
slice(m - c // 2, m + c - c // 2) for m, c in zip(mid, pattern.shape[:-1])
)
return jnp.zeros((*dims, num_channels)).at[slices].set(pattern)
def is_inexact(leaf):
"""Check whether a leaf holds floating point values."""
return jnp.issubdtype(jnp.asarray(leaf).dtype, jnp.inexact)
def partition(tree):
"""Split a pytree into its float leaves and its remaining leaves."""
return (
jax.tree.map(lambda x: x if is_inexact(x) else None, tree),
jax.tree.map(lambda x: None if is_inexact(x) else x, tree),
)
def combine(params, static):
"""Reassemble a pytree split by `partition`."""
return jax.tree.map(
lambda x, y: y if x is None else x, params, static, is_leaf=lambda x: x is None
)
orbium_pattern, orbium_rule_params = load_pattern("Orbium")
state_init = sample_state(orbium_pattern)
params, static = partition(orbium_rule_params)
1. Is the gradient real?¶
jax.grad always returns a number. It differentiates the program, and a program can be
discontinuous while still producing a plausible derivative at every point it is asked
about. The way to check is the definition: for a smooth function the central difference
$$\frac{L(\theta + \varepsilon) - L(\theta - \varepsilon)}{2\varepsilon}$$
converges to the derivative as $\varepsilon \to 0$. If it converges to something else, or does not converge, the program is not what autodiff thinks it is.
This section is the one place in either notebook that needs float64: in float32 the rounding error swamps the difference quotient long before $\varepsilon$ gets small enough to matter. It is scoped to a context manager so nothing else inherits it.
@contextlib.contextmanager
def float64():
"""Run a block in float64 — finite differences need more precision than float32."""
jax.config.update("jax_enable_x64", True)
try:
yield
finally:
jax.config.update("jax_enable_x64", False)
def as_float64(tree):
"""Cast every floating point leaf of a pytree to float64."""
return jax.tree.map(lambda x: x.astype(jnp.float64) if is_inexact(x) else x, tree)
def mean_state(params, static, state, *, num_steps=32, **kwargs):
"""Compute the mean of the state after a rollout, as a scalar loss."""
cs = make_lenia(combine(params, static), **kwargs)
return jnp.mean(cs(state, num_steps=num_steps))
def finite_difference(loss_fn, params, index, eps):
"""Take a central difference of `loss_fn` along a single leaf of `params`."""
leaves, treedef = jax.tree.flatten(params)
def bump(sign):
bumped = list(leaves)
bumped[index] = bumped[index] + sign * eps
return jax.tree.unflatten(treedef, bumped)
return float((loss_fn(bump(1)) - loss_fn(bump(-1))) / (2 * eps))
def audit(rule_params, state, *, eps_values=(1e-4, 1e-6, 1e-8), **kwargs):
"""Compare autodiff against finite differences for every float parameter."""
params, static = partition(as_float64(rule_params))
state = state.astype(jnp.float64)
def loss_fn(params):
return mean_state(params, static, state, **kwargs)
grads = jax.grad(loss_fn)(params)
names = [
jax.tree_util.keystr(p)
for p, _ in jax.tree_util.tree_flatten_with_path(params)[0]
]
columns = "".join(f"{f'FD {eps:g}':>14}" for eps in eps_values)
print(f"{'parameter':<24}{'autodiff':>14}{columns}")
for index, (name, leaf) in enumerate(zip(names, jax.tree.leaves(grads))):
row = f"{name:<24}{float(np.asarray(leaf).ravel()[0]):>14.6e}"
for eps in eps_values:
row += f"{finite_difference(loss_fn, params, index, eps):>14.6e}"
print(row)
with float64():
audit(orbium_rule_params, state_init)
parameter autodiff FD 0.0001 FD 1e-06 FD 1e-08
.weight 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 .kernel_params.r 1.354891e-02 2.210893e-02 8.612956e-01 8.477837e+01 .kernel_params.beta 1.682684e-19 0.000000e+00 0.000000e+00 -8.673617e-11 .growth_params.mean -8.474663e-04 -1.786997e-03 -8.474811e-04 -8.474664e-04
.growth_params.std 2.861261e-01 3.102742e-01 2.863225e-01 2.861261e-01
Three different behaviours in one table.
growth_params.mean and growth_params.std are exact. The finite difference walks onto
the autodiff value and stays there. These are honestly differentiable.
weight and kernel_params.beta are zero — not small, zero. That is the normalization
gauge: Orbium has one kernel with one ring, so the only direction those parameters have is
the scale direction, and the loss is scale-free in it.
kernel_params.r diverges. Autodiff reports about 1.4e-02; the finite difference grows
without bound as $\varepsilon$ shrinks. A difference quotient that grows like
$1/\varepsilon$ means the numerator is not shrinking at all — the loss steps. Multiplying
it back out shows the step size directly.
with float64():
params_64, static_64 = partition(as_float64(orbium_rule_params))
state_64 = state_init.astype(jnp.float64)
names = [
jax.tree_util.keystr(p)
for p, _ in jax.tree_util.tree_flatten_with_path(params_64)[0]
]
index_r = names.index(".kernel_params.r")
def loss_64(params):
"""Compute the float64 rollout loss as a function of the float parameters."""
return mean_state(params, static_64, state_64)
print(f"{'eps':>10}{'FD':>16}{'FD x 2 eps':>16}")
for eps in [1e-5, 1e-6, 1e-7, 1e-8]:
fd = finite_difference(loss_64, params_64, index_r, eps)
print(f"{eps:>10.0e}{fd:>16.4e}{fd * 2 * eps:>16.4e}")
eps FD FD x 2 eps
1e-05 9.8413e-02 1.9683e-06
1e-06 8.6130e-01 1.7226e-06
1e-07 8.4901e+00 1.6980e-06
1e-08 8.4778e+01 1.6956e-06
FD x 2 eps settles at about 1.7e-06 and stays there as $\varepsilon$ shrinks. That is the
signature of a jump: the loss really does step by a fixed amount however small a
neighbourhood you look in.
Notebook 54 shows the cause — the Gaussian core does not vanish where the support mask cuts
it, so grid cells enter the kernel discontinuously as r grows. Under the exponential core
they enter with weight zero, and the same audit converges.
with float64():
audit(orbium_rule_params, state_init, kernel_fn=exponential_kernel_fn)
parameter autodiff FD 0.0001 FD 1e-06 FD 1e-08
.weight 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 .kernel_params.r 5.178198e-03 5.694372e-03 5.234086e-03 5.178199e-03 .kernel_params.beta -8.570026e-18 1.301043e-14 -1.301043e-12 0.000000e+00 .growth_params.mean 2.277975e-02 2.015727e-02 2.200741e-02 2.277975e-02
.growth_params.std 5.704332e-02 1.309604e-01 6.748625e-02 5.704332e-02
Every parameter now agrees to the digits printed, r included.
Notice how small $\varepsilon$ had to get. At 1e-4 even the growth parameters are visibly
off, because the loss is only piecewise smooth: Lenia clips its state to [0, 1], and in
this run 98% of cells sit exactly on that boundary, so a finite difference taken over too
wide a step straddles cells crossing it. Section 3 returns to what clipping does to the
gradient.
The lesson generalizes past Lenia: a gradient that exists is not the same as a gradient that is correct, and the difference costs one finite-difference sweep.
2. How far back can a gradient see?¶
Backpropagating through $n$ steps multiplies $n$ Jacobians together, so the gradient
inherits whatever the dynamics do to small perturbations. The cheapest way to measure that
is forward mode: jax.jvp propagates one tangent vector alongside the trajectory, and since
the scan returns every state, a single call gives the tangent at every step,
$$\delta_k = \frac{\partial\,\mathrm{state}_k}{\partial\,\mathrm{state}_0}\,\delta_0 .$$
If $\|\delta_k\| \sim e^{\lambda k}$ then $\lambda$ is a Lyapunov exponent — measured by autodiff rather than by running two simulations and subtracting.
def tangent_norms(state, rule_params, key, *, num_steps):
"""Propagate one tangent vector and return its norm at every step."""
tangent = jax.random.normal(key, state.shape)
tangent = tangent / jnp.linalg.norm(tangent)
def states_fn(state):
return make_lenia(rule_params)(state, num_steps=num_steps, return_states=True)[
1
]
_, tangents = jax.jvp(states_fn, (state,), (tangent,))
return jnp.linalg.norm(tangents.reshape(num_steps, -1), axis=-1)
def fit_exponent(norms, *, start):
"""Fit log-norm against step count, returning the slope and its R-squared."""
steps = np.arange(start, len(norms)) + 1
log_norms = np.log(np.asarray(norms[start:]))
slope, intercept = np.polyfit(steps, log_norms, 1)
residual = log_norms - (slope * steps + intercept)
return float(slope), float(1 - np.var(residual) / np.var(log_norms))
num_horizon_steps = 256
key, subkey_soup, subkey_a, subkey_b = jax.random.split(key, 4)
state_soup = jax.random.uniform(subkey_soup, (*spatial_dims, channel_size))
norms_orbium = tangent_norms(
state_init, orbium_rule_params, subkey_a, num_steps=num_horizon_steps
)
norms_soup = tangent_norms(
state_soup, orbium_rule_params, subkey_b, num_steps=num_horizon_steps
)
for name, norms in [("Orbium", norms_orbium), ("soup", norms_soup)]:
rate, r_squared = fit_exponent(norms, start=32)
print(
f"{name:<8} |delta_1| = {float(norms[0]):.2e}"
f" |delta_256| = {float(norms[-1]):.2e}"
f" lambda = {rate:+.4f}/step R2 = {r_squared:.3f}"
)
Orbium |delta_1| = 1.52e-01 |delta_256| = 9.15e-01 lambda = +0.0024/step R2 = 0.102 soup |delta_1| = 9.50e-01 |delta_256| = 3.93e+05 lambda = +0.0288/step R2 = 0.947
plt.figure(figsize=(6, 3.5))
plt.semilogy(
np.arange(1, num_horizon_steps + 1), norms_orbium, label="Orbium (a soliton)"
)
plt.semilogy(
np.arange(1, num_horizon_steps + 1), norms_soup, label="random soup (chaotic)"
)
plt.xlabel("step")
plt.ylabel(r"$\|\delta_k\|$")
plt.title("Growth of a perturbation, measured by autodiff")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
The same rule, two initial states, six orders of magnitude apart.
Orbium is a soliton. Its tangent wanders around $10^0$ for 256 steps without settling into a line, and the exponential fit returns an $R^2$ near zero — the fit reporting that there is no exponential trend to find. Backpropagating through the whole rollout is safe, which is why notebook 54 can optimize over a 128-step budget without gradient clipping heroics.
The soup is chaotic, with a clean exponential fit. Its e-folding time is the horizon: a few hundred steps of backpropagation through a soup produces a gradient whose magnitude says more about the dynamics than about the loss.
Reverse mode inherits the same rate.
print(f"{'steps':>8}{'|grad| Orbium':>18}{'|grad| soup':>16}")
for horizon in [1, 4, 16, 64, 256]:
row = f"{horizon:>8}"
for state in [state_init, state_soup]:
grads = jax.grad(mean_state)(params, static, state, num_steps=horizon)
grad_norm = jnp.sqrt(sum(jnp.sum(x**2) for x in jax.tree.leaves(grads)))
row += f"{float(grad_norm):>18.3e}"
print(row)
steps |grad| Orbium |grad| soup
1 9.700e-02 0.000e+00
4 2.030e-01 5.436e-04
16 2.859e-01 2.107e+02
64 3.699e-01 6.019e+03
256 6.286e-01 1.241e+06
Note the soup's first rows: the gradient is not merely small, it is gone. Section 3 explains why.
3. Where the gradient goes blind¶
Two mechanisms remove whole regions of state space from the gradient's view, and they have opposite symptoms.
Clipping blocks one path out of two¶
The Lenia update ends with
state = jnp.clip(state + self._growth(perception) / self.T, 0.0, 1.0)
and jnp.clip has derivative zero outside [0, 1]. T is the time resolution and each step
adds growth / T, so a small T means large steps. 5N7KKM is catalogued at T = 1/3 —
growth × 3 per step — which pins almost the whole grid to the boundary.
multi_pattern, multi_rule_params = load_pattern("5N7KKM")
multi_dims, multi_channels, multi_radius, multi_scale = (128, 128), 3, 12, 2
scaled = multi_pattern.repeat(multi_scale, axis=0).repeat(multi_scale, axis=1)
state_multi = sample_state(scaled, num_channels=multi_channels, dims=multi_dims)
def clip_diagnostics(t_value, *, num_steps=8):
"""Report the clipped fraction and how much of the state gradient it costs."""
def loss_fn(state):
cs = Lenia(
spatial_dims=multi_dims,
channel_size=multi_channels,
R=multi_radius,
T=t_value,
state_scale=multi_scale,
rule_params=multi_rule_params,
)
final_state, states = cs(state, num_steps=num_steps, return_states=True)
return jnp.mean(final_state), states
(_, states), grads = jax.value_and_grad(loss_fn, has_aux=True)(state_multi)
return float(jnp.mean((states <= 0.0) | (states >= 1.0))), float(
jnp.mean(grads == 0.0)
)
print(f"{'T':>6}{'clipped cells':>16}{'dL/dstate0 exactly zero':>26}")
for t_value in [1 / 3, 2.0, 10.0]:
clipped, blind = clip_diagnostics(t_value)
print(f"{t_value:>6.2f}{clipped:>16.1%}{blind:>26.1%}")
T clipped cells dL/dstate0 exactly zero
0.33 97.4% 0.0%
2.00 97.1% 0.0%
10.00 97.1% 0.0%
97% of cells sit on the clip boundary at every step — and ∂L/∂state₀ is exactly zero on
only about 3% of them.
Clipping does not disconnect a cell, because a cell has two paths out and the clip blocks
only one. The direct path, state → next state, goes through the clip. The
convolutional path does not: perception is a plain convolution, so a cell still moves
its neighbours' potentials whether or not its own update was clipped. That is why Lenia
stays differentiable in a regime that looks completely saturated.
An off-band potential removes the rule entirely¶
The other blind region has nothing to do with clipping. The growth mapping is a Gaussian
bell, and Orbium's is narrow: mean = 0.15, std = 0.015. Far enough from the mean, the
bell is numerically flat and its derivative vanishes.
growth_params = LeniaGrowthParams(mean=jnp.array(0.15), std=jnp.array(0.015))
growth_derivative = jax.grad(lambda u: exponential_growth_fn(u, growth_params))
print(f"{'u':>6}{'distance (std)':>16}{'growth':>12}{'d growth / d u':>18}")
for u in [0.15, 0.18, 0.20, 0.25, 0.30, 0.50]:
print(
f"{u:>6.2f}{abs(u - 0.15) / 0.015:>16.1f}"
f"{float(exponential_growth_fn(jnp.array(u), growth_params)):>12.6f}"
f"{float(growth_derivative(jnp.array(u))):>18.3e}"
)
u distance (std) growth d growth / d u
0.15 0.0 1.000000 -0.000e+00 0.18 2.0 -0.729329 -3.609e+01 0.20 3.3 -0.992268 -1.718e+00 0.25 6.7 -1.000000 -1.985e-07 0.30 10.0 -1.000000 -2.572e-19 0.50 23.3 -1.000000 -0.000e+00
potential = make_lenia(orbium_rule_params).perceive(state_soup)
distance = (jnp.mean(potential) - 0.15) / 0.015
print(
f"a soup's potential: mean = {float(jnp.mean(potential)):.4f}, "
f"{float(distance):.1f} standard deviations from the growth mean"
)
print()
grads_rule = jax.grad(mean_state)(params, static, state_soup, num_steps=1)
grads_state = jax.grad(lambda state: mean_state(params, static, state, num_steps=1))(
state_soup
)
print("after one step from a soup:")
print(
f" |dL/dtheta| = "
f"{float(jnp.sqrt(sum(jnp.sum(x**2) for x in jax.tree.leaves(grads_rule)))):.3e}"
)
print(f" |dL/dstate0| = {float(jnp.linalg.norm(grads_state)):.3e}")
a soup's potential: mean = 0.5025, 23.5 standard deviations from the growth mean after one step from a soup: |dL/dtheta| = 0.000e+00 |dL/dstate0| = 9.887e-03
A soup sits more than twenty standard deviations from the growth band, so every cell is in
the flat region and the rule parameters vanish from the loss — they reach it only
through the bell. The state gradient survives, because state + growth / T adds growth to
the state a cell already had, and that path never passes through the bell.
Two blind spots with opposite symptoms, then. Clipping costs a cell one of its two paths and leaves the parameters visible. An off-band potential costs the parameters everything and leaves the state gradient intact — which is exactly why a rule search started from a random soup does not move, and why every optimization in notebook 54 starts from a living creature.
4. What the gradient cannot distinguish¶
A gradient can only separate parameters that change the loss. Some do not — not because they are unimportant, but because another parameter can undo them exactly.
free_kernel_fn builds its kernel from Gaussian bumps under a soft mask:
mask * jnp.sum(b * bell(radius / r, a, w), axis=-1)
The bump enters only through radius / r compared against a with width w, so replacing
$(r, a, w)$ by $(sr,\; a/s,\; w/s)$ leaves the expression algebraically unchanged. The
kernel really has a bump at absolute radius $a\,r$ with absolute width $w\,r$: two degrees
of freedom carried by three parameters.
free_rule_params = LeniaRuleParams(
channel_source=jnp.array([0]),
channel_target=jnp.array([0]),
weight=jnp.array([1.0]),
kernel_params=FreeKernelParams(
r=jnp.array([1.0]),
b=jnp.array([[1.0]]),
a=jnp.array([[0.5]]),
w=jnp.array([[0.15]]),
),
growth_params=LeniaGrowthParams(mean=jnp.array([0.15]), std=jnp.array([0.015])),
)
def free_trajectory(log_params, *, num_steps=8):
"""Roll a single-bump free-kernel rule forward from the Orbium seed."""
rule_params = LeniaRuleParams(
channel_source=jnp.array([0]),
channel_target=jnp.array([0]),
weight=jnp.array([1.0]),
kernel_params=FreeKernelParams(
r=jnp.exp(log_params["log_r"]),
b=jnp.array([[1.0]]),
a=jnp.exp(log_params["log_a"]),
w=jnp.exp(log_params["log_w"]),
),
growth_params=LeniaGrowthParams(
mean=jnp.exp(log_params["log_growth_mean"]),
std=jnp.exp(log_params["log_growth_std"]),
),
)
cs = make_lenia(rule_params, kernel_fn=free_kernel_fn)
return cs(state_init, num_steps=num_steps, return_states=True)[1]
truth = {
"log_r": jnp.log(jnp.array([1.0])),
"log_a": jnp.log(jnp.array([[0.5]])),
"log_w": jnp.log(jnp.array([[0.15]])),
"log_growth_mean": jnp.log(jnp.array([0.15])),
"log_growth_std": jnp.log(jnp.array([0.015])),
}
for scale in [0.8, 1.3]:
rescaled = dict(truth)
rescaled["log_r"] = truth["log_r"] + jnp.log(scale)
rescaled["log_a"] = truth["log_a"] - jnp.log(scale)
rescaled["log_w"] = truth["log_w"] - jnp.log(scale)
deviation = jnp.max(jnp.abs(free_trajectory(rescaled) - free_trajectory(truth)))
print(f"s = {scale}: max trajectory deviation = {float(deviation):.3e}")
s = 0.8: max trajectory deviation = 9.388e-07 s = 1.3: max trajectory deviation = 5.960e-07
Bitwise-level agreement, so the two rules are the same dynamics. Now the consequence: fit a
recorded trajectory and see what comes back. The starting point is 20% off in every
parameter and deliberately 30% high in r.
observed = free_trajectory(truth)
key, subkey = jax.random.split(key)
subkeys = jax.random.split(subkey, len(truth))
log_params = {
name: value + 0.2 * jax.random.normal(subkey, value.shape)
for (name, value), subkey in zip(truth.items(), subkeys)
}
log_params["log_r"] = log_params["log_r"] + jnp.log(1.3)
def identification_loss(log_params):
"""Compute the mean squared error between the simulated and observed trajectory."""
return jnp.mean(jnp.square(free_trajectory(log_params) - observed))
optimizer = optax.chain(optax.clip_by_global_norm(1.0), optax.adam(3e-2))
opt_state = optimizer.init(log_params)
@jax.jit
def identification_step(log_params, opt_state):
"""Take one Adam step on the identification loss."""
loss, grads = jax.value_and_grad(identification_loss)(log_params)
updates, opt_state = optimizer.update(grads, opt_state, log_params)
return optax.apply_updates(log_params, updates), opt_state, loss
losses = []
for _ in range(400):
log_params, opt_state, loss = identification_step(log_params, opt_state)
losses.append(float(loss))
print(f"loss {losses[0]:.3e} -> {losses[-1]:.3e}")
loss 2.496e-03 -> 1.110e-14
print(f"{'parameter':<18}{'truth':>12}{'fitted':>12}{'rel. error':>14}")
for name in truth:
true_value = float(jnp.exp(truth[name]).ravel()[0])
fitted_value = float(jnp.exp(log_params[name]).ravel()[0])
print(
f"{name[4:]:<18}{true_value:>12.5f}{fitted_value:>12.5f}"
f"{(fitted_value - true_value) / true_value:>13.2%}"
)
print()
print(f"{'combination':<18}{'truth':>12}{'fitted':>12}{'rel. error':>14}")
for name, factors in [("a * r", ("log_a", "log_r")), ("w * r", ("log_w", "log_r"))]:
true_value = float(
jnp.exp(truth[factors[0]]).ravel()[0] * jnp.exp(truth[factors[1]]).ravel()[0]
)
fitted_value = float(
jnp.exp(log_params[factors[0]]).ravel()[0]
* jnp.exp(log_params[factors[1]]).ravel()[0]
)
print(
f"{name:<18}{true_value:>12.5f}{fitted_value:>12.5f}"
f"{(fitted_value - true_value) / true_value:>13.3%}"
)
parameter truth fitted rel. error r 1.00000 0.82827 -17.17% a 0.50000 0.60367 20.73% w 0.15000 0.18110 20.73% growth_mean 0.15000 0.15000 0.00% growth_std 0.01500 0.01500 0.00% combination truth fitted rel. error a * r 0.50000 0.50000 0.000% w * r 0.15000 0.15000 -0.000%
The trajectory is reproduced to the last bits of float32 — and r is still wrong by double
digits, with a and w wrong in the opposite direction to compensate. The identifiable
combinations a·r and w·r come back to five decimal places, and both growth parameters
are exact.
The fit is not bad; the model is non-identifiable, and there is no gradient along the
direction that would have brought r home. The practical consequence: fix r = 1 when
fitting a free kernel, and let a and w carry the shape. Notebook 54's soliton search
does exactly that.
5. Where the gradient points in space¶
∂L/∂state₀ has the shape of the grid, so it can be looked at as an image. It answers a
question that is otherwise hard to ask: where would a nudge have mattered?
def final_mass(state, rule_params, *, num_steps):
"""Compute the total mass left after a rollout."""
return jnp.sum(make_lenia(rule_params)(state, num_steps=num_steps))
sensitivity = jax.grad(final_mass)(state_init, orbium_rule_params, num_steps=64)
limit = float(jnp.max(jnp.abs(sensitivity)))
fig, axes = plt.subplots(1, 3, figsize=(11, 3.5))
axes[0].imshow(state_init[..., 0], cmap="viridis")
axes[0].set_title("initial state")
axes[1].imshow(sensitivity[..., 0], cmap="RdBu_r", vmin=-limit, vmax=limit)
axes[1].set_title(r"$\partial\,\mathrm{mass}_{64}\,/\,\partial\,\mathrm{state}_0$")
axes[2].imshow(jnp.abs(sensitivity[..., 0]) > 0.05 * limit, cmap="gray")
axes[2].set_title("above 5% of the peak")
for axis in axes:
axis.axis("off")
plt.tight_layout()
plt.show()
print(
f"nonzero on {float(jnp.mean(sensitivity != 0.0)):.1%} of the grid, "
"above 5% of its peak on "
f"{float(jnp.mean(jnp.abs(sensitivity) > 0.05 * limit)):.1%}"
)
nonzero on 97.8% of the grid, above 5% of its peak on 3.9%
The sensitivity is nonzero almost everywhere — 64 steps is long enough for influence to cross the grid — but very unevenly spread, and signed: red cells would have added mass 64 steps later, blue cells would have removed it. The strongest response sits on the creature's rim rather than its body, which is where notebook 54's adversarial attack ends up putting its perturbation.
6. What colour costs¶
LeniaUpdate normalizes the kernel weights globally:
self.normalized_weight = rule_params.weight / jnp.sum(rule_params.weight)
Growth is then a weighted average over all kernels. So the total growth budget is fixed at one, and channels share it: a creature spread over three channels gets roughly a third of the growth rate per channel that the same creature gets in one.
That is not a bug — it is Chan's multi-kernel convention, and it keeps growth bounded however many kernels a rule has. But it means a three-channel creature runs on a slower clock than its one-channel twin, and the two are not comparable. Here is the same creature, Orbium's rule and Orbium's pattern, in one channel and copied into three.
def speed_of(cs, state, *, num_steps=192):
"""Measure how far the center of mass travels, in R per unit of Lenia time."""
_, states = cs(state, num_steps=num_steps, return_states=True)
grid = jnp.array(states.shape[1:3], jnp.float32)
def center_of_mass(state):
mass_grid = jnp.sum(state, axis=-1)
centers = []
for axis in range(2):
mass_axis = jnp.sum(mass_grid, axis=1 - axis)
angles = 2 * jnp.pi * jnp.arange(state.shape[axis]) / state.shape[axis]
resultant = jnp.sum(mass_axis * jnp.exp(1j * angles))
centers.append(
(jnp.angle(resultant) % (2 * jnp.pi)) / (2 * jnp.pi) * state.shape[axis]
)
return jnp.stack(centers)
centers = jax.vmap(center_of_mass)(states)
steps = jnp.diff(centers, axis=0)
steps = steps - jnp.round(steps / grid) * grid
return float(jnp.linalg.norm(jnp.sum(steps, axis=0)) / R * T / num_steps)
mono_cs = make_lenia(orbium_rule_params)
colour_rule_params = LeniaRuleParams(
channel_source=jnp.arange(3, dtype=jnp.int32),
channel_target=jnp.arange(3, dtype=jnp.int32),
weight=jnp.ones((3,)),
kernel_params=LeniaKernelParams(
r=jnp.repeat(orbium_rule_params.kernel_params.r, 3),
beta=jnp.repeat(orbium_rule_params.kernel_params.beta, 3, axis=0),
),
growth_params=LeniaGrowthParams(
mean=jnp.repeat(orbium_rule_params.growth_params.mean, 3),
std=jnp.repeat(orbium_rule_params.growth_params.std, 3),
),
)
colour_cs = make_lenia(colour_rule_params, channel_size=3)
colour_state = jnp.repeat(state_init, 3, axis=-1)
mono_speed = speed_of(mono_cs, state_init)
colour_speed = speed_of(colour_cs, colour_state)
print(f"one channel : {mono_speed:.4f}")
print(f"three channels: {colour_speed:.4f}")
print(f"ratio : {mono_speed / colour_speed:.2f}x")
one channel : 0.4620 three channels: 0.1656 ratio : 2.79x
The same rule, the same pattern, a factor of about three. Nothing about the creature changed — only how many ways its growth budget is divided.
This is worth knowing before comparing creatures: a one-channel and a three-channel result belong on separate leaderboards, and a colour creature that looks slow may simply be running on a slower clock.
Takeaways¶
- A gradient that exists is not a gradient that is correct. One finite-difference sweep separates the two, and it is the only way to catch a program that is discontinuous where autodiff reports a finite slope.
- The usable rollout length is a property of the pattern, not of the optimizer. A soliton can be differentiated through indefinitely; a chaotic soup cannot.
- Two blind regions, opposite symptoms. Clipping takes one of a cell's two paths and leaves the parameters visible; an off-band potential takes the parameters entirely and leaves the state gradient intact.
- Some parameters are unrecoverable in principle. A perfect fit with wrong parameters is what non-identifiability looks like, and no amount of data fixes it — only a better parameterization does.
- Normalization has consequences beyond the gradient. Sharing one growth budget across channels puts a colour creature on a slower clock than its one-channel twin.