---
summary: The mathematics the core rests on, from quantities, their charts and their products, through operations and their derivatives of first and second order composed over programs, to queries, terms, importance weights, random keys and retirement.
depends_on: [retraction, automatic-differentiation, derivatives-in-charts, product-manifold, generative-models]
---

# The problem

## The problem

The core is the part of skoro that every later entry builds on: the quantities a model is made of,
the operations that relate them, the derivatives of those operations, the log-density a model
defines, and the changes a model goes through as it is built, branched and retired. It ships no
distribution, but it must be exactly right about what a distribution, an estimator or a sampler
will ask of it, because each of those will trust it without checking.

This document works out that mathematics, at the level of precision the implementation needs.
Most of it is standard, and the standard parts are in the math notes it depends on:
[retractions](/math/retraction.md), [automatic differentiation over frames](/math/automatic-differentiation.md),
[derivatives in tangent charts](/math/derivatives-in-charts.md) and
[product manifolds](/math/product-manifold.md). What is particular to skoro, and what this
document settles, is how those pieces meet [the charter](/archive/charter/writeup.md): which quantities
are differentiated and when, what order of derivative a computation needs from each step, what
selection and discrete quantities do to derivatives, what a query does to a model's density, how
random keys make draws reproducible, and what retiring a quantity means.

## Setting

A *model* is a set of *quantities* and a set of *terms*. Each quantity takes values in its *type*:
a manifold $\manifold$ of dimension $\dimension$ with a retraction, or a countable set. Each term
$t$ is a real-valued function $\logterm_t$ of some of the quantities, its *log-density*, and the
model's log-density is their sum (the measure theory underneath is in
[the appendix](#appendix-the-probability-space-of-a-model)),

$$
	\log \density(\point) = \sum_t \logterm_t(\point) . @tag(log-density)
$$

For a manifold-valued quantity, $\point \oplus \tangent$ and $\otherPoint \ominus \point$ are its
retraction and inverse retraction [](eq:math/retraction/oplus), and derivatives are read in tangent
charts [](eq:math/derivatives-in-charts/derivative). A quantity's *dimension* is that of its tangent
space; its *storage size* is the number of reals that represent a point, which may be larger (a
rotation has dimension $3$ and is stored as $9$ or $4$ numbers). Every flat vector of derivative
coordinates in this document is indexed by dimension, never by storage.

Computation is batched: a model is evaluated at $N$ points at once, one per row of a batch. Rows
never interact: an evaluation at $N$ points is $N$ independent evaluations, and every derivative
below is per row. Batching is therefore absent from the mathematics, and only the implementation
needs to keep it straight.

A model is two graphs over the same objects, and the rest of this document uses both. Each answers
different questions, and neither is ever called simply "the graph". The running example is a
vehicle at $\point_0$ and then $\point_1$, which measures a landmark at $m$ once from each position,
as $z_0$ and $z_1$.

### The compilation graph

The *compilation graph* records how values are computed. Its nodes are the model's values: its
quantities, its constants, and its *mapped values*, the outputs of operations, among which are the
terms' log-densities. It has an edge from each input of an operation to the operation's output, so
a node with incoming edges is a mapped value, computed from the nodes those edges come from.

- It is acyclic, because a model is built forward: an operation is added only after its inputs.
- Its *roots*, the nodes with no incoming edge and so no ancestor, are the quantities and the
  constants. They are the inputs of the map that the model computes, and every mapped value is a
  function of the roots it is reached from.
- A quantity is a root even when it was drawn from a distribution whose parameters are mapped
  values. That dependence belongs to the probabilistic graph below, not to this one, which is why
  no derivative ever flows through a sampler.

**Figure: A compilation graph.** The log-density of the term $t_4$, $\logterm_4 = -\tfrac12 \sum_i (z_0 - (m - \point_0))_i^2$, as it is computed. Each edge runs from an input of an operation to its output. The roots are the quantities $\point_0$, $m$, $z_0$ and the constant $-\tfrac12$; every other node is a function of them.

Nodes:

- $\point_0$: Quantity
- $m$: Quantity
- $z_0$: Quantity
- `h`: Mapped value, $m - \point_0$
- `r`: Mapped value, $z_0 - (\cdot)$
- `sq`: Mapped value, $r \odot r$
- `sum`: Mapped value, $\textstyle\sum_i$
- $-\tfrac12$: Constant
- $\logterm_4$: A term's log-density

Edges:

- $m$ → `h`: input of an operation
- $\point_0$ → `h`: input of an operation
- $z_0$ → `r`: input of an operation
- `h` → `r`: input of an operation
- `r` → `sq`: input of an operation
- `r` → `sq`: input of an operation
- `sq` → `sum`: input of an operation
- `sum` → $\logterm_4$: input of an operation
- $-\tfrac12$ → $\logterm_4$: input of an operation

Evaluating a model compiles a program from the compilation graph: the subgraph of the nodes from
which a requested target is reached, its operations taken in a topological order. That subgraph is
the program's [computation graph](/math/automatic-differentiation.md#programs-and-their-computation-graphs),
and its roots are the roots of the compilation graph that the targets reach.

### The probabilistic graph

The *probabilistic graph* records what depends on what in distribution. Its nodes are the model's
quantities and its terms. It has an edge from each quantity a term is *conditioned on* to the term,
and an edge from each term to the quantity it *scores*, if any.

- It is bipartite: every edge joins a quantity and a term.
- Each quantity is scored by at most one term: the term of the distribution it was drawn from. A
  term that scores nothing is a *likelihood factor* over the quantities it is conditioned on: in
  Bayesian terms, the likelihood of data that is not a quantity of the model, such as the
  measurements a retirement conditioned on. A term that scores an observed quantity is a likelihood
  too, but only in the queries that observe it; a likelihood factor is one in every query.
- It is acyclic, because a term is conditioned only on quantities that existed before the quantity
  it scores. A topological order of it is therefore an *ancestral* order, in which each quantity
  comes after everything its distribution is conditioned on.
- Its edges are read off the compilation graph rather than declared: a term is conditioned on the
  quantities among the roots its log-density is reached from, other than the one it scores.

**Figure: A probabilistic graph.** An arrow into a term comes from a quantity it is conditioned on; a dotted arrow out of a term points to the quantity it scores. Retiring $\point_0$ and $z_0$ replaces the terms of $\point_0$, $\point_1$ and $z_0$ by one term over the Markov blanket $\{\point_1, m\}$.

Nodes:

- $t_1$: Term, $\density(\point_0)$
- $t_2$: Term, $\density(m)$
- $\point_0$: Quantity
- $m$: Quantity
- $t_3$: Term, $\density(\point_1 \mid \point_0)$
- $t_4$: Term, $\density(z_0 \mid \point_0, m)$
- $\point_1$: Quantity
- $z_0$: Quantity
- $t_5$: Term, $\density(z_1 \mid \point_1, m)$
- $z_1$: Quantity

Edges:

- $t_1$ → $\point_0$: scores a quantity
- $t_2$ → $m$: scores a quantity
- $\point_0$ → $t_3$: conditions a term
- $t_3$ → $\point_1$: scores a quantity
- $\point_0$ → $t_4$: conditions a term
- $m$ → $t_4$: conditions a term
- $t_4$ → $z_0$: scores a quantity
- $\point_1$ → $t_5$: conditions a term
- $m$ → $t_5$: conditions a term
- $t_5$ → $z_1$: scores a quantity

## Quantities

**Euclidean quantities.** The core's flat types are the real line $\R$,
vectors $\R^n$ and matrices $\R^{m \times n}$, with $\point \oplus \tangent = \point + \tangent$ and
$\otherPoint \ominus \point = \otherPoint - \point$. A matrix is flattened column by column into
$\R^{mn}$ for its tangent coordinates, so its dimension equals its storage size. The dimensions
$n$, $m$ are fixed for each quantity when it is created, and every operation checks them then.

**Positive quantities.** Scales, rates and variances live on the positive orthant
$(0, \infty)^n$, with the retraction of [the example in retractions](/math/retraction.md#the-definition),
elementwise:

$$
	\point \oplus \tangent = \point \odot e^{\tangent}, \qquad
	\otherPoint \ominus \point = \log \otherPoint - \log \point . @tag(positive)
$$

A tangent coordinate is a relative change, and a point is stored as its $n$ values, so its
dimension equals its storage size. The positive orthant is not flat in the sense of the notes, since
$\oplus$ is not addition, but it is flat in disguise: $\log(\point \oplus \tangent) = \log \point +
\tangent$, so the logarithm is one chart of the whole space in which every retraction is a
translation. Three things follow. The retraction is defined for every $\tangent \in \R^n$ and never
leaves the space. The identities $(\point \oplus \tangent) \ominus \point = \tangent$ and
$\point \ominus (\point \oplus \tangent) = -\tangent$ both hold exactly. And the chart Hessian equals
the Hessian of the pull-back everywhere, not only at stationary points: in
[%full](statement:math/derivatives-in-charts/chart-hessian), $\varphi_t(u) = t\tangent + u$, whose
derivative does not move, so the correction vanishes. What is not trivial is reading a positive
quantity as a real number, through the inclusion $(0, \infty)^n \subset \R^n$: its derivative in
charts is $\operatorname{diag}(\point)$, which moves with the point.

**Discrete quantities.** The core's discrete types are the finite sets
$\mathrm{Categorical}(n) = \{0, \dots, n - 1\}$, with $\mathrm{Boolean} = \mathrm{Categorical}(2)$,
and the integers $\mathbb Z$. A discrete quantity has no tangent space, no retraction and no
dimension, and it is never a coordinate of a derivative. Its values are written
$\discretePoint$. A density over a discrete quantity is a probability mass, a density with respect
to the counting measure.

**Products.** A product of continuous types is a continuous type: a quantity of type
$\manifold_1 \times \cdots \times \manifold_n$ is one point $(\point_1, \dots, \point_n)$, stored as
its factors' storage concatenated, with the stacked tangent space
[](eq:math/product-manifold/stack) and the factor-wise retraction of
[%full](statement:math/product-manifold/product-retraction). The core has products of continuous
types only. A quantity with a continuous part and a discrete part is two quantities, as the charter
has it, and a joint density over several discrete quantities needs no product, since a term may be
conditioned on any number of quantities.

A product's factors may be of any continuous type, including products, and of dimensions known
only at run time. That is what makes a product the natural type of a marginal over a Markov
blanket, whose members are only known when a quantity is retired ([Retirement](#retirement)).

## Operations

An *operation* is a smooth map

$$
	\otherPoint = \map(\point_1, \dots, \point_m; \discretePoint_1, \dots, \discretePoint_r) @tag(operation)
$$

from continuous inputs $\point_j \in \manifold_j$ and discrete inputs $\discretePoint_i$ to an
output $\otherPoint \in \otherManifold$, smooth in the continuous inputs for each value of the
discrete ones. In the compilation graph it is a mapped node with an incoming edge from each input.
Its partial derivatives $\partial_j \map$ are taken in tangent charts, and it may provide four maps,
each over a *frame* of directions
([automatic differentiation over frames](/math/automatic-differentiation.md)):

| map | given | returns |
| --- | --- | --- |
| forward | the inputs | $\otherPoint$ |
| tangent | a frame $\tangent_j$ per continuous input | $\otherTangent = \sum_j \partial_j \map\, \tangent_j$ |
| cotangent | a frame $\alpha$ at the output | $\partial_j \map\transpose \alpha$ per continuous input |
| second order | velocities $\tangent_j$, and one cotangent $\alpha$ | $\hvp_j(\alpha) = \frac{\dd}{\dd t}\big\vert_{t=0} \partial_j \map(\point \oplus t\tangent)\transpose \alpha$ per input |

The second-order map takes a velocity frame of several columns, one Hessian direction each, and a
single cotangent column. Its definition in charts is the one under which a forward-over-reverse pass
computes the chart Hessian [](eq:math/derivatives-in-charts/chart-hessian), by
[%full](statement:math/derivatives-in-charts/passes).

Discrete inputs carry no derivatives, and a discrete output carries none either. A map that turns a
discrete input into a continuous output (an *embedding*, such as an integer read as a real) has no
partial derivative in it. A map that turns a continuous input into a discrete output (a
*discretization*, such as the indicator $\mathbb 1[\point \ge 0]$) passes no cotangent back. For a
step function that is the true derivative everywhere but at its steps, and the core makes no
attempt to do better: a model that needs a derivative across a discretization needs a relaxation,
which is a modelling choice, not the core's.

### Two gates

A pass runs over a program's computation graph, whose roots include the quantities it
differentiates with respect to. Within a pass, an operation's continuous inputs differ in two ways
that the implementation must keep apart:

- an input is *differentiable* if it is reached from a root the pass differentiates with respect
  to, which decides whether cotangents ($\alpha$, and in a Hessian pass $\beta$) are pulled back
  into it along the
  [backpropagation graph](/math/automatic-differentiation.md#reverse-pulling-cotangent-frames-over-the-backpropagation-graph);
- an input *carries a velocity* if it is reached from a root seeded with a velocity in this pass,
  which decides whether its velocity enters $\hvp$.

The second implies the first, and not conversely: a Hessian pass is often seeded at one quantity
while all of them are differentiated. An operation's $\hvp_j$ is absent whenever $\partial_j \map$
is not needed, and a missing velocity counts as zero.

### Order

Every operation states its *order* $\order(\map)$: $0$ if it provides only the forward map, $1$ if
it also provides the tangent and cotangent maps, $2$ if it also provides the second-order map. Most
primitives have order $2$. An operation that runs burn's automatic differentiation inside, such as a
neural network, has order $1$; one whose derivative does not exist or was not written has order
$0$. What a pass needs of a step depends on where the step sits in the program's computation graph.

> **Definition (The order a pass requires).** Fix a pass over a program: the roots $L$ it differentiates with respect to, the roots
> $S \subseteq L$ seeded with velocities (for a Hessian pass), and the targets $T$. A step $k$ with map
> $\map_k$ is *on the pass* when, in the program's computation graph, its output is continuous and
> reaches a target in $T$, and one of its continuous inputs is reached from a root in $L$. The order
> the pass requires of step $k$ is
>
> $$
> 	\req(k) = \begin{cases}
> 		0 & \text{if step } k \text{ is not on the pass,} \\
> 		2 & \text{if it is, the pass is a Hessian pass, and an input of step } k \text{ is reached from } S, \\
> 		1 & \text{otherwise.}
> 	\end{cases}
> $$
>
> The pass is admissible when $\order(\map_k) \ge \req(k)$ for every step $k$.

A pass that is not admissible is refused before it runs, naming a step $k$ with
$\order(\map_k) < \req(k)$. This refines the charter's "least order over the program's steps" in two
ways that matter in practice. A step that no differentiated root reaches, or that reaches no target,
needs nothing, so a model may contain a first-order neural network and still be differentiated
twice with respect to quantities that do not feed it. And a step that only pulls back cotangents in
a Hessian pass, because no seeded velocity reaches it, needs only first order: its $\hvp$ terms are
zero by [](eq:math/automatic-differentiation/leibniz).

## Selection

*Selection* chooses one of several continuous values by a discrete index:

$$
	\otherPoint = \operatorname{select}(\discretePoint;\, \point_0, \dots, \point_{n-1}) = \point_\discretePoint,
	\qquad \discretePoint \in \mathrm{Categorical}(n),\; \point_i \in \manifold . @tag(select)
$$

All the $\point_i$ and $\otherPoint$ are of one type. For each $\discretePoint$ the map is a
projection, so its derivatives in charts are as simple as derivatives get: since
$\otherPoint = \point_\discretePoint$, the chart at $\otherPoint$ is the chart at
$\point_\discretePoint$, and

$$
	\partial_{\point_i} \otherPoint = \begin{cases} I & i = \discretePoint, \\ 0 & i \ne \discretePoint, \end{cases}
	\qquad \hvp_{\point_i} = 0 .
$$

The cotangent at $\otherPoint$ goes to the selected value and nowhere else, and the index receives
nothing. Mixtures, switching models and data association are built on selection: a mixture whose
component is a discrete quantity reads that component's parameters through selection, and a
measurement associated by a discrete quantity selects the target it came from. The derivative with
respect to the values, at a fixed index, is exactly what a gradient-based algorithm conditioned on
the discrete quantities needs.

## The primitives

The core ships all of the primitive operations below, with the derivative maps they carry; those
on Euclidean quantities first. In
the table, $\odot$ is the elementwise product, $\tangent_a$ the velocity of an input $a$ in a
Hessian pass (or its tangent, in the tangent map), and $\alpha$ the cotangent at the output. Each
frame is handled column by column, so one column is shown.

| operation | tangent | cotangent per input | second order per input |
| --- | --- | --- | --- |
| $a + b$ | $\tangent_a + \tangent_b$ | $\alpha,\ \alpha$ | $0,\ 0$ |
| $a - b$ | $\tangent_a - \tangent_b$ | $\alpha,\ -\alpha$ | $0,\ 0$ |
| $-a$ | $-\tangent_a$ | $-\alpha$ | $0$ |
| $a \odot b$ | $\tangent_a \odot b + a \odot \tangent_b$ | $b \odot \alpha,\ a \odot \alpha$ | $\tangent_b \odot \alpha,\ \tangent_a \odot \alpha$ |
| $s\,a$, $s \in \R$ | $\tangent_s\, a + s\, \tangent_a$ | $a\transpose \alpha,\ s\, \alpha$ | $\tangent_a\transpose \alpha,\ \tangent_s\, \alpha$ |
| $A B$ | $\tangent_A B + A \tangent_B$ | $\alpha B\transpose,\ A\transpose \alpha$ | $\alpha \tangent_B\transpose,\ \tangent_A\transpose \alpha$ |
| $\sum_i a_i$ | $\sum_i (\tangent_a)_i$ | $\alpha \mathbf 1$ | $0$ |
| $e^{a}$ (elementwise) | $e^{a} \odot \tangent_a$ | $e^{a} \odot \alpha$ | $e^{a} \odot \tangent_a \odot \alpha$ |
| $\log a$ (elementwise) | $\tangent_a / a$ | $\alpha / a$ | $-\tangent_a \odot \alpha / a^2$ |
| slice, concatenate, transpose | the same map | its adjoint | $0$ |
| select | [](eq:select) | to the selected value | $0$ |

The operations on positive quantities read their tangents and cotangents in the charts
[](eq:positive), that is, in logarithmic coordinates, where all but the inclusion are linear. Here
$a$ and $b$ are positive, $c$ is Euclidean, and the products, quotients and powers are
elementwise.

| operation | from, to | tangent | cotangent per input | second order per input |
| --- | --- | --- | --- | --- |
| $a$ as a real | $(0, \infty)^n \to \R^n$ | $a \odot \tangent_a$ | $a \odot \alpha$ | $a \odot \tangent_a \odot \alpha$ |
| $e^{c}$ | $\R^n \to (0, \infty)^n$ | $\tangent_c$ | $\alpha$ | $0$ |
| $\log a$ | $(0, \infty)^n \to \R^n$ | $\tangent_a$ | $\alpha$ | $0$ |
| $a \odot b$ | $(0, \infty)^n \to (0, \infty)^n$ | $\tangent_a + \tangent_b$ | $\alpha,\ \alpha$ | $0,\ 0$ |
| $1 / a$ | $(0, \infty)^n \to (0, \infty)^n$ | $-\tangent_a$ | $-\alpha$ | $0$ |

The exponential and logarithm here are not the rows of the first table: those map $\R^n$ to
$\R^n$, read in the flat charts, where their derivatives are $e^{a}$ and $1/a$; these change type,
and in the charts of the positive side their derivatives are identities. Both are kept, because a
model may want either. The retraction [](eq:positive) is the composite $\point \odot e^{\tangent}$
of the second and fourth rows, and its inverse $\log \otherPoint - \log \point$ of the third and a
Euclidean subtraction, so every derivative of both, of either order, is an identity or zero.

Every row is checked the same way: the tangent map is the derivative of the forward map, the
cotangent map its transpose, and the second-order map the derivative of the cotangent map along the
velocities, holding $\alpha$ fixed. The bilinear operations ($a \odot b$, $s\,a$, $A B$) are the
only ones whose second-order maps mix inputs: the pull-back into one input involves the other, so
the other's velocity appears. The linear ones have none. The discrete operations the core needs
(integer arithmetic and comparison) and the embedding of an integer as a real carry no
derivatives, by the conventions of [Operations](#operations).

## Programs, composites and passes

A model is evaluated by compiling a program from the compilation graph, as in
[The compilation graph](#the-compilation-graph). Four passes run over a compiled program: the
forward pass, the tangent pass [](eq:math/automatic-differentiation/jvp) along the edges of its
computation graph, the cotangent pass [](eq:math/automatic-differentiation/vjp) along the edges of
its backpropagation graph, and the Hessian pass [](eq:math/automatic-differentiation/leibniz), which
runs both. Derivatives are collected at the program's roots only, which are the leaves of its
backpropagation graph, and [Queries](#queries) says which roots those are.

A *composite* operation is a named subgraph of the compilation graph, such as the retraction of a
product, or (in later entries) a distribution's log-density: a set of operations with designated
inputs and one output, stored under one name. Its derivative maps are *defined* as the composition
of its parts', so expanding a composite into its parts at compilation changes no value that any
pass computes, and its order is that of its parts under [%full](statement:order). A hand-written
fused form of a composite, introduced for speed or for numerical care, is a different operation
with the same mathematics, and it is correct only if every pass over it agrees with every pass over
the composite, to within rounding.

The product retraction is the composite the charter's acceptance test is about.

> **Proposition (The composite product retraction agrees with its factors).** Let a product's retraction be built as the composite of slicing the point and the tangent into
> blocks, applying each factor's retraction to its blocks, and concatenating the results. Then every
> pass over a program containing it (values, tangent frames, cotangent frames, and Hessian products)
> computes the same values as the factor-wise computation, block by block.

<Proof>

Slicing and concatenation are linear, so their tangent maps are themselves, their cotangent maps
are their adjoints (scattering into a zero block and gathering a block), and their second-order maps
are zero. In the tangent pass, slicing hands each factor's retraction exactly its blocks of the
frames, and concatenation stacks what they return. In the cotangent and Hessian passes, the
concatenation's adjoint hands each factor exactly its blocks of $\alpha$ and $\beta$, the factors
pull them back, and the slicing's adjoint places each result in its own block, adding only zeros. By
[](eq:math/automatic-differentiation/leibniz), with the linear steps' $\hvp$ zero, the product's
$\beta$ is the factors' $\beta$s, stacked. So every pass computes the block-diagonal assembly of the
factors' passes, which is the derivative structure of [](eq:math/product-manifold/blocks).

</Proof>

The proof introduces no arithmetic beyond additions of zero, so a test may hold the agreement to a
tolerance of a few units in the last place.

### The acceptance test

The charter accepts its position on products on one condition: that the retraction of a runtime
product, and its first and second derivatives, agree with the factor-wise computation. On the
core's own types that is a weak test at second order. In the charts of every one of them a
retraction is a translation, so every second-order map along it is zero, and the second
derivatives under test are zeros on both sides. The test is therefore run on a product that
includes a curved factor, and composed with a function whose second derivative does not vanish.

**A test-only sphere.** The core's tests carry a private type, never shipped: the unit sphere
$S^2 = \{\point \in \R^3 : \lVert \point \rVert = 1\}$, of dimension $2$ and storage size $3$. Its
tangent plane at $\point$ is $\point^\perp$, identified with $\R^2$ by a matrix $B(\point)$ whose
two orthonormal columns span it. No such $B$ is smooth on the whole sphere (a smooth field of
tangent bases on $S^2$ would give a nowhere-vanishing tangent vector field, which the hairy ball
theorem rules out), so the identification is smooth on an open part of the sphere, and the tests
stay inside it. Its retraction projects back onto the sphere:

$$
	\point \oplus \tangent = \frac{\point + B(\point)\,\tangent}{\lVert \point + B(\point)\,\tangent \rVert},
	\qquad
	\otherPoint \ominus \point = \frac{B(\point)\transpose \otherPoint}{\point\transpose \otherPoint}
	\quad (\point\transpose \otherPoint > 0) . @tag(sphere)
$$

The inverse follows from $\point\transpose B(\point) = 0$: the point $\otherPoint$ is
$(\point + B(\point)\tangent)/n$ with $n$ the norm, so $\point\transpose \otherPoint = 1/n$ and
$B(\point)\transpose \otherPoint = \tangent / n$. The derivative of the retraction at
$\tangent = 0$ is $B(\point)\tangent - \point\,(\point\transpose B(\point)\tangent) = B(\point)\tangent$,
whose coordinates are $\tangent$, so it is a retraction. Unlike the core's types, the sphere is
curved: its tangent charts do not commute, the correction term of
[%full](statement:math/derivatives-in-charts/chart-hessian) is not zero away from stationary
points, and the retraction's second-order maps are not zero. Its retraction is a composite of three
private primitives with hand-written derivatives in charts: the inclusion $S^2 \subset \R^3$, the
map $(\point, \tangent) \mapsto B(\point)\tangent$, and the normalization $w \mapsto w / \lVert w \rVert$
from $\R^3 \setminus \{0\}$ onto the sphere.

**What it compares.** The test builds a runtime product of two Euclidean factors, a real line, a
positive factor and the sphere, and a nonlinear scalar function $f$ on it that reads the positive
factor and the sphere as real numbers, through their inclusions. It checks, at points away from
the singular set of $B$:

1. that every pass over the product's retraction, values, tangent frames, cotangent frames and
   Hessian products, agrees with the factor-wise computation, block by block, as
   [%full](statement:product-agreement) says it must;
2. that the chart gradient and the chart Hessian products of $f$ through the product's retraction,
   computed by the passes, agree with the factor-wise assembly of
   [](eq:math/product-manifold/blocks), and with finite differences of their definitions
   [](eq:math/derivatives-in-charts/gradient) and [](eq:math/derivatives-in-charts/chart-hessian);
3. that on the sphere's block, away from a stationary point, the chart Hessian differs from the
   Hessian of the pull-back by the correction term, and on every other block it does not.

The private primitives' own derivatives are checked the same way first, against finite differences
of [](eq:math/derivatives-in-charts/derivative), so that a failure of the acceptance test points at
the composition rather than at a primitive.

## Queries

The charter puts conditioning in the query: the model holds the joint density
[](eq:log-density), and each question asked of it says which quantities are fixed. A *query*
partitions the quantities its targets reach in the compilation graph into

- the *observed* quantities $\point_O$, fixed at values the query gives;
- the *free continuous* quantities $\point_F$, the variables of the query;
- the *free discrete* quantities $\discretePoint_D$, which a generation draws, and whose values any
  other query gives and holds fixed.

The query's log-density is then a function of $\point_F$ alone,

$$
	\log \density(\point_F;\, \point_O, \discretePoint_D) = \sum_t \logterm_t(\point_F, \point_O, \discretePoint_D), @tag(query)
$$

which is, up to the constant $\log \density(\point_O, \discretePoint_D)$, the log-density of the
conditional distribution of $\point_F$ given the rest. One model answers many queries: with the
measurements observed it is a posterior; with nothing observed, a prior to simulate from; with the
states observed, a likelihood of measurements.

**Roots.** The roots of a query's program, the roots of its computation graph, are the
quantities and constants its targets reach. Of these, the free continuous quantities are the ones
every pass differentiates with respect to, in their tangent charts; the observed and discrete
quantities and the constants are fixed inputs. A quantity's value is an input of the query even
when the quantity was drawn from a distribution whose parameters depend on other quantities, since
it is a root of the compilation graph: $\partial \log \density / \partial \point_i$ counts
$\point_i$'s appearances in the terms, never its role in generating its children.

**Folding.** A mapped value whose inputs are all fixed is itself fixed, and need not be computed on
every evaluation: it can be computed once and put in place of the subgraph that computes it, as a
new root. Nothing a pass computes changes, since no derivative flows into a fixed value. A value is
fixed at one of two levels, and is folded at that level.

- **In the model.** A value reached only from constants is the same in every query, snapshot and
  branch, so it is folded once, when it is added to the compilation graph, and the model keeps the
  result, in double precision on the host like every value it stores. A query casts it to its own
  precision like any other constant. Conditioning in the query does not prevent this: a constant is
  not a quantity, so no query can observe or free it.
- **In a compiled query.** A value reached from observed quantities, and from constants, is fixed
  only for the queries that observe those quantities, so it is folded when such a query is
  compiled, and reused by every evaluation of the compiled program.

A value reached from a free quantity is not folded. The difference shows in a Gaussian term written
with a covariance factor $L$, whose log-density reads its inverse and the logarithms of its diagonal
(a later entry decides how the Gaussian takes its covariance). If $L$ is a constant, both are
folded once in the model; if it is an observed quantity, once per compiled query; if it is free,
they are computed on every evaluation, together with their derivatives with respect to $L$.

**Discrete quantities in a query.** The gradient of [](eq:query) at fixed $\discretePoint_D$ is the
gradient of the conditional density of $\point_F$ given $\discretePoint_D$. An algorithm that
targets the distribution of $\point_F$ alone, with $\discretePoint_D$ unobserved, needs
$\discretePoint_D$ summed or sampled out, which a gradient cannot do; that is why the charter has a
gradient-based algorithm refuse a query with a free discrete quantity, rather than silently
condition on the value it was given.

**Flattened coordinates.** An estimator works with vectors, not with quantities. Ordered as the
estimator chooses, the free quantities $\point_F = (\point_1, \dots, \point_n)$ are one point of the
product manifold $\manifold_1 \times \cdots \times \manifold_n$, a tangent of which is a vector
$\tangent \in \R^\dimension$, $\dimension = \sum_i \dimension_i$, stacked by
[](eq:math/product-manifold/stack). An estimator's step is that product's retraction,

$$
	\point_F \oplus \tangent = (\point_1 \oplus \tangent_1, \dots, \point_n \oplus \tangent_n), @tag(step)
$$

and a Jacobian, gradient or Hessian over $\point_F$ is the block matrix of
[](eq:math/product-manifold/blocks). So the layout an estimator uses and a product quantity type are
the same mathematical object, and the core should build them from one construction.

## Terms and the importance weight

A term $t$ is a node of the probabilistic graph. Its log-density $\logterm_t$ is a node of the
compilation graph, computed from the quantities it is reached from there. Following the charter, a
term records the quantity it *scores*, $s_t$, if any, and the quantities it is *conditioned on*,
$C_t$, so that $\logterm_t$ reads as $\log \density_t(s_t \mid C_t)$. The quantities a term
involves, $V_t = \{s_t\} \cup C_t$, its neighbours in the probabilistic graph, determine every
structural question the core answers: which terms touch a quantity, what a quantity's Markov
blanket is, and what retiring it disturbs. They are only right if they contain every quantity
$\logterm_t$ depends on, so the core reads $V_t$ off the compilation graph, as the quantities among
the roots $\logterm_t$ is reached from, rather than trusting a declaration; $C_t$ is then $V_t$
without $s_t$. In the example, $\logterm_4$ is reached from $\point_0$, $m$ and $z_0$, and the term
$t_4$ scores $z_0$, so it is conditioned on $\point_0$ and $m$.

**The generative reading.** A distribution may provide a sampler for the quantity it scores, which
draws $s_t$ given $C_t$, with a log-density $\log \proposal_t(s_t \mid C_t)$ of its own; an exact
sampler has $\proposal_t = \density_t$ (charter, position 9). A query generates its free
quantities, continuous and discrete, ancestrally: each from the sampler of the term that scores it,
in a topological order of the probabilistic graph. In this section, $\point_F$ stands for all of
them. Generation is possible exactly when every free quantity is scored by a term with a sampler;
it is refused otherwise, naming a quantity that is not. The importance weight of a generated point
is

$$
	\log \weight(\point_F) = \sum_t \logterm_t(\point_F, \point_O) - \sum_{t \,:\, s_t \in F} \log \proposal_t(s_t \mid C_t), @tag(weight)
$$

which is the charter's [](eq:archive/charter/writeup/weight), with the observed quantities' terms
and the likelihood factors, the likelihoods of this query, contributing $\logterm_t$ alone, and
each term whose sampler is exact cancelling.

> **Proposition (Importance weights).** Let a query leave the quantities $\point_F$ free and fix $\point_O$, and suppose that every free
> quantity is scored by exactly one term, whose sampler has density $\proposal_t$ with
> $\proposal_t > 0$ wherever $\density_t > 0$, each pair with respect to the same base measure. If
> $\point_F$ is drawn ancestrally from the samplers, then for every integrable $h$,
>
> $$
> 	\E_\proposal\bigl[\weight(\point_F)\, h(\point_F)\bigr] = \int h(\point_F)\, \density(\point_F, \point_O)\, \dd \point_F = \density(\point_O)\; \E\bigl[h(\point_F) \mid \point_O\bigr] .
> $$

<Proof>

Ancestral generation draws $\point_F$ with density
$\proposal(\point_F) = \prod_{t : s_t \in F} \proposal_t(s_t \mid C_t)$ with respect to the product
of the base measures, since each quantity is drawn given quantities drawn or fixed before it. Then
$\weight = \density(\point_F, \point_O) / \proposal(\point_F)$ wherever $\proposal > 0$, which covers
the support of $\density(\cdot, \point_O)$ by the positivity condition, and
$\E_\proposal[\weight h] = \int h\, \density(\point_F, \point_O)\, \dd \point_F$.

</Proof>

Two conditions in the statement are the core's to check or to state. The base measures must match.
In the core, a continuous quantity's density is with respect to Lebesgue measure on the space its
values lie in, $\R^n$ or the open subset $(0, \infty)^n$ of it, and a discrete quantity's with
respect to counting measure, so this holds as long as every distribution keeps to that convention.
The one place it is easy to break is a density written in the logarithmic chart of a positive
quantity, such as a log-normal written as a normal in $u = \log \point$: with respect to Lebesgue
measure on $(0, \infty)^n$, it carries the Jacobian of the chart,

$$
	\log \density_{\point}(\point) = \log \density_u(\log \point) - \sum_i \log \point_i . @tag(log-chart-density)
$$

Densities on manifolds without such a chart will need a base measure of their own when the first
of them arrives. Positivity is the distribution's responsibility, and the core cannot check it.

## Random keys

Every draw comes from a key (charter, position 10). Mathematically a key is an element of a finite
set $\mathcal K$, with two deterministic maps: $\operatorname{split} : \mathcal K \to \mathcal K^n$,
and $\operatorname{fold} : \mathcal K \times \Sigma \to \mathcal K$, which folds in one address
segment from the set $\Sigma$ of names and integers. The key of the quantity at address
$a = a_1 / \cdots / a_m$ in a run with key $\key$ is

$$
	\key_a = \operatorname{fold}(\cdots \operatorname{fold}(\operatorname{fold}(\key, a_1), a_2) \cdots, a_m), @tag(address-key)
$$

and its sampler turns $\key_a$ into the uniform variates it consumes. Because $\key_a$ depends
on $\key$ and $a$ alone, a quantity's draw does not depend on the order in which the model was
built, on what was retired before it, or on which branch holds it, and a namespace's keys share the
prefix of folds that the namespace's own address determines.

The independence of draws rests on an idealization: that $\operatorname{fold}$ and
$\operatorname{split}$ behave as a pseudorandom function, so that distinct addresses yield keys
indistinguishable from independent uniform ones. Two things the core must ensure for that to mean
anything. Distinct addresses must give distinct inputs: each segment's encoding must be injective
and tagged by its kind, so that the name `"3"` and the integer $3$ differ, and folding segment by
segment, rather than hashing a concatenated string, rules out collisions such as `ab/c` against
`a/bc`. And each batch row of a run must have its own run key, split from the run's, so that rows
are independent.

Because the uniform variates are drawn on the host from the key and then converted to the
precision a query computes in, one key gives the same draws at single and double precision, up to
the rounding of the conversion.

## Retirement

Retiring is the one change that shrinks a model: the quantities under a namespace are removed, and
what they said about the rest is kept. Retirement is read on the probabilistic graph. Let $E$ be
the quantities to retire, $\mathcal T_E$ the terms joined to any of them ($V_t \cap E \neq \emptyset$),
and

$$
	B = \Bigl(\bigcup_{t \in \mathcal T_E} V_t\Bigr) \setminus E @tag(blanket)
$$

their *Markov blanket*: the other quantities joined to those terms. In the example, retiring the
first tick, $E = \{\point_0, z_0\}$, touches the terms $t_1$, $t_3$ and $t_4$, and its blanket is
$\{\point_1, m\}$.

Retirement must condition as well as marginalize. The retired namespace of a tick holds that tick's
measurements as well as its states, and integrating the measurements out of the joint would discard
exactly the information a sliding window keeps them for. So the quantities of $E$ split into those
fixed at given values, $E_O$ (usually the measurements, with the values the estimator observed),
and those integrated out, $E_F$. Since conditioning belongs to the query, the values of $E_O$ are
given to the retirement itself, and once retired they become part of the model: what remains is the
joint of the remaining quantities given those values.

> **Proposition (Retirement leaves the rest of the model unchanged).** Let $E$ be the quantities to retire, of which $E_O$ are fixed at given values and $E_F$ are
> integrated out, let $\mathcal T_E$ be the terms joined to a quantity of $E$ in the probabilistic
> graph, and let $B$ be the other quantities joined to those terms. Replace the terms $\mathcal T_E$
> by the single term
>
> $$
> 	\logterm_B(\point_B) = \log \int \exp \Bigl(\sum_{t \in \mathcal T_E} \logterm_t\bigl(\point_{E_F}, \point_{E_O}, \point_B\bigr)\Bigr)\, \dd \point_{E_F} .
> $$
>
> Then, as a function of every quantity outside $E$, the new model's density equals the old model's
> density with $E_F$ integrated out and $E_O$ fixed at its values.

<Proof>

Split [](eq:log-density) into the terms of $\mathcal T_E$, which involve only $E$ and $B$, and the
rest, which involve no quantity of $E$. Fixing $E_O$ and integrating over $E_F$, the second group
comes out of the integral as a constant factor, by Tonelli's theorem for the non-negative
integrand, and what remains of the first group is $\exp \logterm_B(\point_B)$.

</Proof>

The replacement must be computed from the terms $\mathcal T_E$ alone. Computing the marginal of the
whole model on $B$ instead would fold in the terms outside $\mathcal T_E$, which the model still
holds, and count them twice. The integral defining $\logterm_B$ rarely has a closed form. It has one
for linear-Gaussian terms, where retiring the oldest state of a chain is the update-and-predict step
of a Kalman filter; otherwise it is approximated, by a Laplace approximation at the current
estimate in the archive's sliding windows. Either way, computing it needs the distributions and the
linear algebra of later entries; what the core needs is the structure: $E$, $\mathcal T_E$, $B$,
the removal, and a place for the replacement, a likelihood factor over $B$ whose type, when it is one
density over the blanket's tangent coordinates jointly, is the product of the blanket's types.

Two consequences for the rest of the model:

- **Mapped values.** In the compilation graph, a mapped value under the retired namespace that a
  remaining node is reached from, and that is reached only from remaining roots, is not a function
  of $E$ and must not be removed with it. One reached from $E$ can only reach terms in
  $\mathcal T_E$, which are replaced, and other mapped values reached from $E$, which go with it.
- **Generation after retirement.** A blanket quantity whose scoring term was in $\mathcal T_E$
  ($\point_1$ in the example, whose transition $t_3$ was conditioned on $\point_0$) is scored by no
  term afterwards, only joined to the replacement. Unless the replacement is written in a
  generative form, as a term scoring that quantity with a sampler, generation from the new model is
  refused for it, by structure.

## In skoro's terms

- **Quantities** are Euclidean (real, vector, matrix), positive ($(0, \infty)^n$, with the
  retraction [](eq:positive)), and discrete (categorical, boolean, integer), and products of
  continuous quantities. Dimensions are fixed when a quantity is created, and every
  derivative coordinate is indexed by tangent dimension.
- **The compilation graph** holds values and operations, and programs are compiled from it; **the
  probabilistic graph** holds quantities and terms, and is read off the compilation graph and each
  term's scored quantity.
- **Operations** are maps [](eq:operation) with up to four derivative maps over frames, each with
  an order. The core's primitives are those of [The primitives](#the-primitives), and every other
  operation, the product's retraction included, is a composite of them.
- **Folding** replaces a value reached only from fixed inputs with a root: in the model when they
  are constants, in a compiled query when they include observed quantities.
- **Passes** are the forward, tangent, cotangent and Hessian passes over a compiled program,
  differentiating with respect to a query's free continuous quantities in their charts, refused
  before running when not admissible by [%full](statement:order).
- **The density** is the sum of terms [](eq:log-density). A query fixes the observed quantities and
  the discrete ones, and its density is [](eq:query). Generation is ancestral from samplers, and its
  weight is [](eq:weight).
- **Random keys** are derived from addresses by [](eq:address-key).
- **Retirement** removes a namespace's quantities $E$, conditioning on the given values of $E_O$ and
  integrating out $E_F$, and replaces the terms $\mathcal T_E$ with one term over the blanket
  [](eq:blanket).

## Appendix: the probability space of a model

A skoro model is a *model built forward* in the sense of
[models built forward, for the measure theorist](/math/generative-models.md), and every object
in this document has a place in that note's measure theory.

- **Quantities are the coordinates.** Each quantity lives in a space with a reference measure:
  Lebesgue measure on $\R^n$ and $\R^{m \times n}$ and on the open subset $(0, \infty)^n$, counting
  measure on $\mathrm{Categorical}(n)$ and on $\mathbb Z$, and product measures on products. The
  probability space can be taken to be the product of the quantities' spaces, with the quantities as
  its coordinates, and the probabilistic graph's ancestral order as the order in which they were
  built.
- **Terms are conditional densities.** A term that scores a quantity has, as its $\exp \logterm_t$,
  that quantity's conditional density given those it is conditioned on, so [](eq:log-density) is the
  logarithm of the joint density of
  [%full](statement:math/generative-models/joint-law). A likelihood factor is the change of measure
  of [%full](statement:math/generative-models/likelihood-factor).
- **Mapped values are functions of the quantities,** and by
  [%full](statement:math/generative-models/doob-dynkin) they are all the random variables a model can
  express. Their laws are image measures, with no density in general, which is why no term scores
  one.
- **A query conditions by dividing densities.** Its log-density [](eq:query) is that of the
  regular conditional distribution of [%full](statement:math/generative-models/conditioning), up to
  the constant $\log \density_O$, which is why a query may drop it.
- **Every number an algorithm estimates is an integral** against the joint density, by
  [%full](statement:math/generative-models/integral): an expectation, a marginal, the evidence
  $\density_O$. Generation draws from the conditional densities in the order the model was built,
  which is Ionescu-Tulcea's construction, and importance weighting [](eq:weight) estimates such an
  integral from those draws.
- **Retirement integrates and conditions.** [%full](statement:retirement) integrates the retired
  free quantities out of the terms that touch them and conditions on the retired observed ones,
  which is the same marginalization and the same division of densities, carried out once and kept.

The test-only sphere has no reference measure in this sense, since it is not an open subset of a
coordinate space; it carries no density in the core's tests. Manifold-valued quantities beyond the
positive orthant will need reference measures of their own (surface measure on a sphere, Haar
measure on the rotations) when they arrive.

## Open questions

> [!TODO] What retirement does in the core
> Computing the replacement $\logterm_B$ needs distributions and linear algebra the core does not
> have. I propose that the core's retirement take the replacement from its caller (given $E$,
> $E_O$'s values, $\mathcal T_E$ and $B$, which the core computes), remove $E$, and attach it, all
> in one atomic change; and that the core's tests exercise it with a private linear-Gaussian dummy
> whose replacement is known in closed form. Computing the replacement (Gaussian, Laplace) is then a
> later entry's. Should retirement instead wait for that entry altogether?

> [!TODO] Re-scoring a blanket quantity
> For generation to survive a retirement, the replacement must be able to score blanket quantities
> that lost their scoring term, which makes retirement more than removal plus a likelihood factor. Is
> that in scope for the core, or is a model that has retired quantities simply not generated from
> until a later entry says how?