# skoro
> Skoro is a Rust library for writing a probabilistic model once, generatively, and running different inference algorithms on it, built for models that grow, branch and are estimated over time.
---
---
title: What is skoro?
summary: A library for writing a probabilistic model once and running any inference algorithm on it, built for models that live over time.
order: 1
---
# What is skoro?
Skoro is a Rust library for writing probabilistic models once, generatively, and running different
inference algorithms on the same model. This page says what that means, who it is for, and the
ideas the rest of the library follows from. The argument for each of them is in
[the charter](/archive/charter/writeup.md).
## One model, written generatively
A model is written the way it is thought of: each uncertain quantity, and the conditional
distribution it has given the quantities it depends on. Consider a range measurement $z$ between a
vehicle at $x$ and a target at $y$:
$$
z \sim \mu(\mathrm{d}z \mid x, y)
$$
That one statement serves three situations. The vehicle and the target can be estimated jointly
from $z$; the vehicle alone can be estimated from $z$ when the target is known; and $z$ can be
simulated from given $x$ and $y$. In skoro these are not three models. The model holds the joint
distribution, and which quantities are observed is part of each question asked of it, so the same
model answers all three, from different threads at once.
Factor-graph libraries, by contrast, describe a model by the residuals of one particular posterior.
Each of the three situations needs a different set of factors, and the third, simulation, cannot be
expressed at all.
## Inference reads the model
The library records a model's structure as the model is built, much as automatic differentiation
records a computation: which quantities depend on which, what each density term scores, which of
its parts are a residual or a whitening factor, which operations can be differentiated to second
order. An inference algorithm reads what it needs from that record. If the model meets the
algorithm's assumptions, the algorithm runs; if not, it refuses before it starts and names what is
missing and where. The user never describes the model a second time for the algorithm's benefit.
> [!DECISION] No capability flags
> What a model can do is read off its structure, never declared by its user. See
> [position 6 of the charter](/archive/charter/writeup.md#6-algorithms-refuse-by-structure).
## Models that live over time
Skoro is built for sequential estimation, where a model is never finished:
- **It grows.** Each tick appends new states and measurements, and each kind of measurement can be
different.
- **It forgets.** States that stop mattering are marginalized onto what remains and retired.
- **It branches.** When it is uncertain which target a measurement came from, the model is held
under several hypotheses at once, as cheap branches of one another, and an algorithm weighs them.
- **It is shared.** One writer builds while estimators run on immutable snapshots of the model,
without waiting for each other.
Sub-models are ordinary functions of quantities, placed under namespaces by the code that calls
them, so the model of one tick is written once and reused for every tick.
## Who it is for
1. **People developing inference algorithms.** Skoro's first purpose is to make new algorithms easy
to write against real models: what an algorithm can discover about a model matters as much as
how a model is built.
2. **People building real-time estimators** for robots and autonomous systems, on CPU or GPU, from
Rust, often under several hypotheses at once.
3. **People doing Bayesian inference** with quantities that live on manifolds, such as rotations
and poses, or that need second-order derivatives.
## What it is not
- **Not a general probabilistic programming language.** Each model has a fixed structure, built
forward. Structure that is itself random is handled one level up, over a family of models, rather
than by programs whose shape changes as they run, as in Gen.jl.
- **Not a replacement for mature solvers.** It does not compete with GTSAM or Ceres on the maturity
of its least-squares solvers, nor with Stan or NumPyro on the breadth of its distributions. What
it offers is the combination: one generative model, a structure algorithms can read,
manifold-aware derivatives of second order, and a model that lives across threads and over time.
## The ideas everything follows from
The [charter](/archive/charter/writeup.md) argues each of these in full:
- **Conditioning belongs to the query**, never to the model.
- **Derivatives come from composition.** A few primitive operations carry hand-written derivatives,
and everything else, distributions and manifold retractions included, is built from them; each
operation reports the order of derivative it supports.
- **Static at the edges, dynamic inside.** Models are written with typed handles the compiler
checks, while everything stored inside is nameable, so models can be serialized and driven from
other languages.
- **Randomness is explicit.** Every random draw comes from a key, so one seed reproduces a run.
- **Precision belongs to the query.** A model commits to no float precision; each estimator chooses
its own, and the solves where conditioning matters run on the host in double precision.
---
---
title: Reading these docs
summary: Two tracks that say what skoro is, notes of the mathematics, and an archive of how it came to be.
order: 90
---
# Reading these docs
These docs are kept in four parts, and each answers a different question.
## What skoro is
Two tracks say what the library is today, and are rewritten whenever it changes.
- **The [guide](/guide.md) is top-down.** It starts from what a user writes (a model, a
distribution, a density to evaluate) and reaches into the machinery only as far as an example
needs. Read it to use skoro.
- **The [concepts](/concepts.md) are bottom-up.** They build the library from its primitive
objects upward, each idea resting only on the ones before it. Read them to understand why skoro
works the way it does, or before changing it.
The [reference](/reference.md) holds what is looked up rather than read: the crates, the
glossary, the notation.
## What it rests on
The [mathematics](/math.md) is a library of notes, one idea to a note: a retraction, a tangent
space, the Fisher information. A note says what it builds on, and the map on that page draws the
whole network, together with the archive's writeups that use it.
## How it came to be
The [archive](/archive.md) keeps the record of every piece of development. A feature's entry
holds the mathematics worked out for it, the implementation companion that turned the mathematics
into a design, and the plan it was built by, each kept as it was approved. A decision's entry holds
the argument for it.
The archive is history, and history is not rewritten. Where the concepts track _tells_ you what is,
an entry _discusses_ what might be, as it stood when it was written. When the library later
changes, a new entry revises the old one, and the old one says so at its top. Every page lists the
entries that made it what it is, under "How this came to be", and every entry lists the pages it
changed.
> [!NOTE]
> Every document is also served as Markdown (follow "View as Markdown" at its foot), and
> [llms.txt](../../llms.txt) indexes all of them for language models.
---
---
title: How skoro is developed
summary: Mathematics first, then the design, a plan, the code, and last the pages that say what is.
order: 91
---
# How skoro is developed
Every change to skoro is an entry in the [archive](/archive.md), worked on in a branch of its
own and moved through a fixed series of stages. Each stage ends when a person approves what it
produced; nothing moves on by itself.
## Features
A feature (an algorithm, a distribution, an example) goes through five stages.
1. **Mathematics.** A writeup of the mathematics at play, revised until it is right. Ideas that
stand on their own become [math notes](/math.md), which later writeups build on.
2. **Implementation.** A companion that turns the mathematics into a design: the types and traits,
sketched in code, where the neighbouring interfaces constrain them, and the tests that will show
the design works. It is meant to be complete enough that someone who knows Rust, the mathematics
and those interfaces could build it without designing anything of their own.
3. **Plan.** The implementation's prose distilled into a checklist of tasks, each with its tests.
4. **Build.** The code, committed as the plan's boxes are ticked.
5. **Publish.** The pages of the guide and the concepts are brought up to date, quoting the code
that now exists, and the entry is linked to the pages it changed.
A **decision** (a position the library takes) is written up and published. A **chore** (a port, a
refactor) is planned, built and published.
## What the stages guarantee
The stages are kept by a small tool, `bin/wf`, and checked on every pull request:
- A branch works on one entry at a time.
- An approved document is frozen. Changing it means reopening its stage, which is recorded, with
its reason, in the entry's history. A design is allowed to change its story; it is not allowed to
forget it.
- Code changes only once a plan is approved, and the pages that say what skoro is change only at
the end, when there is code to quote.
- An entry is finished, or explicitly set aside, before its branch merges.
The procedure itself, written for the agents that help with it, is in the repository's
[`AGENTS.md`](https://github.com/mvarble/skoro/blob/main/AGENTS.md).
---
---
title: Automatic differentiation over frames
summary: How a straight-line program is differentiated by composing each step's derivative maps over its computation graph, forward with frames of tangents, backwards over its backpropagation graph with frames of cotangents, and forward-over-reverse for Hessian products.
---
# Automatic differentiation over frames
Automatic differentiation computes the derivatives of a program exactly, up to rounding, from the
derivatives of the steps it is made of. This note states the three ways of composing them that a
library needs (pushing tangents forward, pulling cotangents back, and differentiating the pull-back
once more for second order), each over a *frame* of several directions at once.
## Programs and their computation graphs
A *straight-line program* has slots $s_0, s_1, \dots$, each holding a vector, and runs a sequence of
steps
$$
s_{o_k} \leftarrow \map_k\bigl(s_{i_{k,1}}, \dots, s_{i_{k,m_k}}\bigr), @tag(step)
$$
each writing one slot from slots written earlier. Each step's map $\map_k$ is smooth, with partial
Jacobians $\partial_j \map_k$, one per input, of shape (output dimension) by (dimension of input
$j$). When an input is not differentiated, its partial Jacobian is simply never used.
A program is pictured by its *computation graph*. Its nodes are the slots, and it has an edge from
each input of a step to that step's output, so a node with incoming edges is computed by the step
that writes it, from the nodes those edges come from. Because each step reads only slots written
before it, the computation graph has no cycles, and the order of the steps is one of its
topological orders. Its *roots*, the nodes with no incoming edge and so no ancestor, are the
slots no step writes: the program's inputs, the variables $\point$ we differentiate with respect
to. Some nodes are
*targets*, the values to differentiate, $f(\point)$. A node with several outgoing edges is a slot
read by several steps.
**Figure: A computation graph.** The program $s_3 \leftarrow \point_1 \odot \point_2$, $s_4 \leftarrow \exp s_3$, $f \leftarrow s_3 + s_4$. Each node is a slot; each edge runs from an input of a step to its output. The roots $\point_1$ and $\point_2$ are the program's inputs, and $f$ is its target.
Nodes:
- $\point_1$: Root (an input)
- $\point_2$: Root (an input)
- $s_3$: Slot written by a step, $\map_1 = {\odot}$
- $s_4$: Slot written by a step, $\map_2 = \exp$
- $f$: Slot written by a step, $\map_3 = {+}$
Edges:
- $\point_1$ → $s_3$: input of a step
- $\point_2$ → $s_3$: input of a step
- $s_3$ → $s_4$: input of a step
- $s_3$ → $f$: input of a step
- $s_4$ → $f$: input of a step
A *frame* is a matrix whose columns are directions: a tangent frame at a slot of dimension
$\dimension$ is a $\dimension \times n$ matrix $\tangent$ of $n$ tangent vectors, and a cotangent
frame a $\dimension \times n$ matrix $\alpha$ whose columns are linear functionals, written as
columns. Below, $\tangent_s$ and $\alpha_s$ are the frames held at slot $s$.
## Forward: pushing tangent frames
Seed each root with a tangent frame $\tangent_{\point}$ (missing means zero), and run the steps in
order, each pushing its inputs' frames forward along the computation graph's edges by its
*Jacobian–vector product*:
$$
\tangent_{s_{o_k}} = \sum_j \partial_j \map_k \; \tangent_{s_{i_{k,j}}} . @tag(jvp)
$$
By the chain rule, the frame reaching a target is $\sum_l \partial f / \partial \point_l \cdot \tangent_{\point_l}$:
one pass with a frame of $n$ columns gives $n$ directional derivatives, and seeding the columns of
an identity gives the whole Jacobian.
## Reverse: pulling cotangent frames over the backpropagation graph
The reverse pass runs on the *backpropagation graph*: the same nodes as the computation graph, with
every edge reversed and labelled by the pull-back $\partial_j \map_k\transpose$ of the step it came
from. Since reversing the edges swaps incoming for outgoing, a node is a root of the computation
graph if and only if it is a *leaf* of the backpropagation graph, a node with no outgoing edge.
Cotangents enter at the targets, flow along its edges, and end at those leaves, the program's
inputs. A node with several incoming edges, a slot read by several steps, receives a cotangent
along each, and holds their sum.
**Figure: Its backpropagation graph.** The same slots with every edge reversed and labelled by the pull-back it applies. Cotangents enter at the target $f$ and end at the leaves $\point_1$ and $\point_2$, the computation graph's roots; $s_3$ has two incoming edges, and its cotangent is their sum.
Nodes:
- $\point_1$: Root (an input)
- $\point_2$: Root (an input)
- $s_3$: Slot written by a step, $\map_1 = {\odot}$
- $s_4$: Slot written by a step, $\map_2 = \exp$
- $f$: Slot written by a step, $\map_3 = {+}$
Edges:
- $s_3$ → $\point_1$: pull-back (“$\partial_1 \map_1\transpose$”)
- $s_3$ → $\point_2$: pull-back (“$\partial_2 \map_1\transpose$”)
- $s_4$ → $s_3$: pull-back (“$\partial \map_2\transpose$”)
- $f$ → $s_3$: pull-back (“$\partial_1 \map_3\transpose$”)
- $f$ → $s_4$: pull-back (“$\partial_2 \map_3\transpose$”)
Seed each target with a cotangent frame $\alpha_f$, every other slot with zero, and run the steps
backwards, each adding its *vector–Jacobian product*, the pull-back of its output's cotangent, into
each input:
$$
\alpha_{s_{i_{k,j}}} \mathrel{+}= \partial_j \map_k\transpose \, \alpha_{s_{o_k}} . @tag(vjp)
$$
The accumulation, which sums the cotangents arriving along every incoming edge of the
backpropagation graph, is the chain rule for a slot with several consumers. When the steps are done,
each root holds $\alpha_{\point_l} = \sum_f (\partial f / \partial \point_l)\transpose \alpha_f$. For a scalar
target seeded with $\alpha_f = 1$, that is the gradient $\nabla_{\point_l} f$, whatever the number
of roots.
Each column of a frame costs a small multiple of one evaluation of the program, in either
direction. A Jacobian with few rows is cheaper in reverse, one with few columns forward.
## Forward-over-reverse: Hessian products
For second derivatives, differentiate the reverse pass itself along a path $\point + t\tangent$. Run
the forward pass of [](eq:jvp) with velocities $\tangent$ seeded at the roots, so that every slot
holds its value $s$ and its velocity $\tangent_s$, then run the reverse pass over the
backpropagation graph carrying two cotangents per slot: $\alpha_s$, as in [](eq:vjp), and
$\beta_s$, its derivative along the path. One column of $\tangent$ is one path; a frame of $n$
columns runs $n$ paths at once, and then $\beta_s$ is a frame of $n$ columns while $\alpha_s$ is the
same for all of them. (A cotangent frame of several columns together with a velocity frame of
several columns would make $\beta$ a three-way array; the passes here take one cotangent column.)
Differentiating [](eq:vjp) along the path by the product rule gives
$$
\beta_{s_{i_{k,j}}} \mathrel{+}= \hvp_{k,j}(\alpha_{s_{o_k}}) + \partial_j \map_k\transpose \, \beta_{s_{o_k}},
\qquad
\hvp_{k,j}(\alpha) \defeq \frac{\dd}{\dd t}\Big|_{t=0} \partial_j \map_k(s + t\tangent_s)\transpose \alpha
= \sum_l \bigl(\partial_l \partial_j \map_k\,[\tangent_{s_{i_{k,l}}}]\bigr)\transpose \alpha . @tag(leibniz)
$$
The first term is the step's own second-order map, which depends on the velocities of its inputs;
the second is the step's ordinary pull-back, applied to $\beta$. A library writes the first for
each step and the program applies the rule. Two consequences are useful when writing $\hvp$ for a
step. If a step is linear, its partial Jacobians are constant and its $\hvp$ is zero. If the
pull-back into input $j$ does not involve input $l$, then the velocity of input $l$ contributes
nothing to $\hvp_{k,j}$.
> **Proposition (Forward-over-reverse computes Hessian products).** Let $f$ be a scalar target, seed its cotangent with $\alpha_f = 1$ and $\beta_f = 0$, and seed the
> roots with velocities $\tangent$. Then the reverse recurrences deliver, at every root $\point_i$,
>
> $$
> \alpha_{\point_i} = \nabla_{\point_i} f(\point), \qquad \beta_{\point_i} = \sum_l \nabla_{\point_i}\nabla_{\point_l} f(\point)\,\tangent_{\point_l} .
> $$
The cotangents $\alpha_s$ are functions of the roots: the reverse recurrence holds as an identity
in $\point$, with every $\partial_j \map_k$ evaluated at the slot values the program computes from
$\point$. Differentiate that identity along $\point + t\tangent$ for one column $\tangent$. By the
product rule, $\beta_s \defeq \frac{\dd}{\dd t}\alpha_s(\point + t\tangent)\big|_{t=0}$ satisfies
[](eq:leibniz), because the slot values move with the velocities the forward pass computed, and the
seeds $\alpha_f = 1$ and $\beta_f = 0$ are the value and the derivative of a constant. At a root,
$\alpha_{\point_i}(\point) = \nabla_{\point_i} f(\point)$ identically in $\point$, so its derivative
along the path is $\sum_l \nabla_{\point_i} \nabla_{\point_l} f(\point)\,\tangent_{\point_l}$.
With all roots stacked into one vector, $\beta_\point = \nabla^2 f(\point)\,\tangent$: a frame of
$n$ velocities gives $n$ columns of the Hessian for about $n$ times the cost of a gradient, and the
identity frame gives all of it.
---
---
title: Derivatives in tangent charts
summary: The derivative of a map between manifolds, its gradient and its Hessian product, all read in the tangent charts a retraction provides, and why automatic differentiation computes them unchanged.
depends_on: [retraction, automatic-differentiation]
---
# Derivatives in tangent charts
A [retraction](/math/retraction.md) gives every point of a manifold a chart, its tangent space. This note
reads derivatives in those charts: the derivative of a map between manifolds as a matrix, the
gradient of a scalar function as a vector, and its second derivative along a direction. Each is
what an optimizer or a sampler needs at a point, and none needs a metric.
## The derivative
Let $\map : \manifold \to \otherManifold$ be smooth, between manifolds of dimensions $\dimension$
and $e$, each with a retraction, and let $\otherPoint = \map(\point)$. The derivative of $\map$ at
$\point$ is the $e \times \dimension$ matrix of $\map$ read in the chart at $\point$ and the chart at
$\otherPoint$:
$$
\D \map(\point)\, \tangent \defeq \frac{\dd}{\dd t}\Big|_{t=0} \bigl(\map(\point \oplus t \tangent) \ominus \otherPoint\bigr),
\qquad \tangent \in \R^\dimension . @tag(derivative)
$$
Equivalently, $\D \map(\point)$ is the ordinary Jacobian at $0$ of the *local representative*
$\hat \map_\point = R_\otherPoint^{-1} \circ \map \circ R_\point$, a map from a neighbourhood of $0$
in $\R^\dimension$ to $\R^e$. Because each $\D R_\point(0)$ is the identity, it is also the
intrinsic differential of $\map$, written in the bases that identify the tangent spaces with
$\R^\dimension$ and $\R^e$.
> **Proposition (The chain rule in tangent charts).** For smooth $\map : \manifold \to \otherManifold$ and $G : \otherManifold \to \bbZ$,
> $\D(G \circ \map)(\point) = \D G(\map(\point))\,\D \map(\point)$, every derivative taken in tangent
> charts.
Near $0$, $\widehat{G \circ \map}_\point = R_z^{-1} \circ G \circ \map \circ R_\point
= (R_z^{-1} \circ G \circ R_\otherPoint) \circ (R_\otherPoint^{-1} \circ \map \circ R_\point)
= \hat G_\otherPoint \circ \hat \map_\point$, with $z = G(\otherPoint)$, because
$R_\otherPoint \circ R_\otherPoint^{-1}$ is the identity near $\otherPoint$. The ordinary chain rule
at $0$, with $\hat \map_\point(0) = 0$, gives the claim.
So derivatives in charts compose exactly as Jacobians do, and the passes of
[automatic differentiation](/math/automatic-differentiation.md) apply unchanged to a program whose steps
are maps between manifolds, as long as each step's derivative is taken in charts. To derive one by
hand:
1. perturb the input, $\point_t = \point \oplus t \tangent$;
2. push it through, $\otherPoint_t = \map(\point_t)$;
3. solve $\otherPoint_t = \otherPoint \oplus \otherTangent(t)$ for $\otherTangent(t)$;
4. the Jacobian–vector product is $\otherTangent'(0)$, and the vector–Jacobian product is its
transpose.
The transpose is taken in coordinates: a cotangent at $\otherPoint$ is a row of $e$ numbers acting
on tangent coordinates, and $\D \map(\point)\transpose$ maps it to a cotangent at $\point$. No inner
product on the tangent spaces is involved.
## The gradient
For a scalar function $f : \manifold \to \R$, its gradient in the chart at $\point$ is
$$
g(\point) \defeq \nabla (f \circ R_\point)(0) = \D f(\point)\transpose 1 \in \R^\dimension, @tag(gradient)
$$
the cotangent that a reverse pass seeded with $1$ delivers. It says how $f$ changes along a step:
$f(\point \oplus \tangent) = f(\point) + g(\point)\transpose \tangent + O(\lVert\tangent\rVert^2)$.
With no metric in play, it is the differential of $f$ written in chart coordinates, not a
Riemannian gradient.
## The Hessian product
The second derivative needs a choice that the first does not, because moving along a direction
moves the chart as well. The *chart Hessian* differentiates the gradient along a direction
$\tangent$, reading the gradient at each point in that point's own chart:
$$
H(\point)\, \tangent \defeq \frac{\dd}{\dd t}\Big|_{t=0} g(\point \oplus t \tangent) . @tag(chart-hessian)
$$
The other candidate is the Hessian of the pull-back $f \circ R_\point$, which keeps the chart at
$\point$ fixed. On a flat manifold the two are the same matrix, the ordinary Hessian. On a curved
one they differ by a term proportional to the gradient.
> **Proposition (The chart Hessian and the Hessian of the pull-back).** Let $\point_t = \point \oplus t \tangent$ and $\varphi_t = R_\point^{-1} \circ R_{\point_t}$, defined
> near $0$ for small $t$. Then
>
> $$
> H(\point)\,\tangent = \nabla^2 (f \circ R_\point)(0)\, \tangent + \Bigl(\frac{\dd}{\dd t}\Big|_{t=0} \D\varphi_t(0)\Bigr)\transpose g(\point),
> $$
>
> so $H(\point)\,\tangent = \nabla^2 (f \circ R_\point)(0)\, \tangent$ at every stationary point, where
> $g(\point) = 0$.
For small $t$ and $u$ near $0$, $R_{\point_t}(u) = R_\point(\varphi_t(u))$, with $\varphi_0$ the
identity and $\varphi_t(0) = R_\point^{-1}(\point_t) = t \tangent$ by
[](eq:math/retraction/roundtrip). So $f \circ R_{\point_t} = (f \circ R_\point) \circ \varphi_t$ near
$0$, and
$$
g(\point_t) = \D \varphi_t(0)\transpose\, \nabla (f \circ R_\point)(t \tangent) .
$$
Differentiating at $t = 0$, with $\D \varphi_0(0) = I$ and $\nabla (f \circ R_\point)(0) = g(\point)$,
gives the claim.
On the rotations, with the retraction $R_\point(u) = \point \operatorname{Exp}(u)$, the
Baker–Campbell–Hausdorff formula gives
$\varphi_t(u) = t \tangent + u + \tfrac{t}{2}\, \tangent \times u$, up to terms of order
$t^2 \lVert u\rVert$ and $t \lVert u\rVert^2$, so the correction is
$\tfrac12 [g(\point)]_\times \tangent$: the chart Hessian is the Hessian of the pull-back plus the
skew-symmetric matrix $\tfrac12 [g(\point)]_\times$, and it is not symmetric away from a stationary
point. On the positive reals, with $R_\point(u) = \point\, e^{u}$, the charts commute instead:
$\varphi_t(u) = t \tangent + u$, whose derivative does not move with $t$, so the correction vanishes
and the two Hessians agree everywhere. Which symmetric matrix to use where they differ, if one is
needed, is the algorithm's decision. At a
stationary point, which is where curvature is most often read (a Laplace approximation, a Newton
step's convergence), the two agree.
## Programs over manifolds
The chart Hessian is the one a forward-over-reverse pass computes, provided each step's
second-order map is read in charts the same way.
> **Proposition (The passes compute chart derivatives).** Let a program's steps be maps between manifolds, with partial derivatives $\partial_j \map$ taken in
> tangent charts, and define each step's second-order map at its inputs $s$, with velocities
> $\tangent_s$, by
>
> $$
> \hvp_j(\alpha) = \frac{\dd}{\dd t}\Big|_{t=0} \partial_j \map(s \oplus t\tangent_s)\transpose \alpha .
> $$
>
> Then for a scalar target $f$, the reverse pass delivers the chart gradient $g(\point)$ at the
> roots, and forward-over-reverse with velocities $\tangent$ delivers the chart Hessian product
> $H(\point)\,\tangent$.
First order is [%full](statement:chain-rule) applied step by step. For the second, the reverse
recurrence holds identically in $\point$, with every $\partial_j \map$ evaluated at the slot values
computed from $\point$ and read in the charts at those values. Differentiate it along
$\point_t = \point \oplus t \tangent$. Each slot traces a curve $s(\point_t)$ through $s(\point)$
whose velocity in the chart at $s(\point)$ is $\tangent_s = \D s(\point)\, \tangent$, the velocity
the forward pass computes. The derivative of the smooth function
$s \mapsto \partial_j \map(s)\transpose \alpha$ along a curve depends only on the curve's velocity,
so it equals the derivative along $s \oplus t \tangent_s$, which is $\hvp_j(\alpha)$. The product
rule then gives the forward-over-reverse recurrence, and at a root the differentiated quantity is
$g(\point_t)$, whose derivative is [](eq:chart-hessian).
---
---
title: Models built forward, for the measure theorist
summary: Why an abstract probability space serves modelling as an abstract vector space serves linear algebra, why the way back to coordinates is harder in probability, and how building a model forward from densities restores it, so that every random variable factors through the quantities and every expectation is an explicit integral.
---
# Models built forward, for the measure theorist
This note is for the reader who knows measure theory far better than statistics, and who wants to
know what a probabilistic model *is* before being told how to compute with one. Its argument is an
analogy with linear algebra, and its conclusion is that the models skoro builds sit in a corner of
probability theory where the analogy is exact: one can reason abstractly and compute concretely,
and move between the two without meeting the edge cases that usually make the move painful.
## The abstraction, and the way back
Linear algebra is done in an abstract vector space $V$. A mathematician working there never says
which coordinate of which $\R^\dimension$ holds which quantity of interest: vectors are added,
subspaces $W \subseteq V$ are taken, linear maps are composed, and the dimension can be anything.
When a computation is needed, a basis is chosen, $V$ becomes $\R^\dimension$, and the abstract
statements become statements about arrays of numbers. The way back is a single theorem: every
finite-dimensional vector space is isomorphic to a coordinate space, and every choice of basis
gives such an isomorphism.
Probability is done in an abstract probability space $(\Omega, \mathscr F, \mathrm P)$. A modeller
working there never says what $\Omega$ is, nor which sets are in $\mathscr F$, nor what
$\mathrm P$ assigns to them. Instead they introduce random variables $\randomPoint : \Omega \to \manifold$,
measurable maps into the spaces where the quantities of the model live, and speak only of those:
their laws, their joint laws, the independence of one from another, the conditional law of one
given another. Any number of them may be introduced, and new ones may be introduced later, provided
the space is rich enough to carry them. The abstraction does for the modeller what $V$ does for the
linear algebraist: it relieves them of a coordinate system they do not need.
The way back to coordinates is where the analogy usually breaks. Several things can go wrong at
once:
- **The space is underdetermined.** Many triples $(\Omega, \mathscr F, \mathrm P)$ carry random
variables with the same joint law, and $\mathscr F$ may hold far more than the model ever
mentions. Nothing in the model says which.
- **Conditioning needs regularity.** The conditional law of $\randomPoint$ given $\otherRandomPoint$, for $\randomPoint$ with values in
$\manifold$ and $\otherRandomPoint$ in $\otherManifold$, defined as a kernel $\kernel$ from $\otherManifold$ to
$\manifold$ characterized by
$$
\mathrm P(\randomPoint \in A,\ \otherRandomPoint \in B) = \int_B \kernel(\otherPoint, A)\, \mathrm P(\otherRandomPoint \in \dd\otherPoint)
\qquad \text{for all measurable } A \subseteq \manifold,\ B \subseteq \otherManifold, @tag(regular-conditional)
$$
a *regular conditional distribution*, exists only under assumptions on the spaces (standard Borel is the
usual one), and even then only up to null sets of $\otherRandomPoint$'s law. Conditioning on an event of
probability zero, which is what observing a continuous quantity does, depends on which version is
chosen; the Borel–Kolmogorov paradox is the classic example of two natural versions disagreeing.
- **Densities need a reference.** A density is a Radon–Nikodym derivative, and exists only when
the law is absolutely continuous with respect to a chosen measure. The law of a function of
quantities, such as $(\randomPoint, \randomPoint)$ or a projection, often has none.
- **Existence needs a theorem.** Writing down a family of random variables does not make a space
that carries them. Kolmogorov's extension theorem builds one from consistent finite-dimensional
laws, but asks for topological regularity of the spaces to do so.
So the modeller who wants to compute, to write a probability as an integral and evaluate it, has
to justify every step back from $(\Omega, \mathscr F, \mathrm P)$, and usually does not.
## A corner of probability where the way back is clean
Applied models are almost never as general as the theory allows. They are built from two
ingredients, and both remove an edge case.
**Seeds with densities.** Every quantity lives in a space with a *reference measure*: Lebesgue
measure on $\R^\dimension$ (or an open subset of it, such as the positive orthant), or counting
measure on a countable set, which may itself sit inside a vector space, such as the integers. Every
distribution a quantity is drawn from is absolutely continuous with respect to its reference
measure, so it is given by a function.
> **Theorem (Radon–Nikodym).** Let $\nu$ be a $\sigma$-finite measure on a measurable space $(\manifold, \mathscr X)$, and $\mu$ a
> finite measure on it with $\mu \ll \nu$, that is, $\nu(A) = 0$ implies $\mu(A) = 0$. Then there is
> a measurable $\density : \manifold \to [0, \infty)$ with
>
> $$
> \mu(A) = \int_A \density \,\dd\nu \qquad \text{for every } A \in \mathscr X,
> $$
>
> and any two such functions agree $\nu$-almost everywhere. One writes $\density = \dd\mu / \dd\nu$
> and $\mu(\dd\point) = \density(\point)\, \nu(\dd\point)$.
With respect to Lebesgue measure, $\density$ is a *probability density function*:
$\mu(\dd\point) = \density(\point)\, \dd\point$, its values are not probabilities, every single
point has measure zero, and $\density$ is determined only up to a Lebesgue-null set. With respect to
counting measure, nothing is null but the empty set, and the density is pinned down completely by
the mass it gives to points.
> **Proposition (A mass function is its measure on points).** Let $\manifold = \{\point_0, \point_1, \dots\}$ be countable, with every subset measurable, and let
> $\nu = \sum_k \delta_{\point_k}$ be counting measure. Every probability measure $\mu$ on
> $\manifold$ is absolutely continuous with respect to $\nu$, its density is unique, and it is
>
> $$
> \density(\point_k) = \mu(\{\point_k\}), \qquad \text{so that} \qquad
> \mu(\dd\point) = \sum_k \density(\point_k)\, \delta_{\point_k}(\dd\point) .
> $$
If $\nu(A) = 0$ then $A = \emptyset$, so $\mu \ll \nu$, and $\nu$ is $\sigma$-finite because
$\manifold$ is countable. Countable additivity gives
$\mu(A) = \sum_{\point_k \in A} \mu(\{\point_k\}) = \int_A \density \,\dd\nu$ with
$\density(\point_k) = \mu(\{\point_k\})$, and uniqueness almost everywhere is uniqueness
everywhere, since the only null set is empty.
This is why $\density$ is called a *probability mass function* in that case, and why the two cases
can be treated as one: a density is always with respect to a reference measure, and a mass function
is a density with respect to counting measure.
**Quantities introduced one at a time.** Each quantity is drawn from a distribution whose
parameters may depend on any quantity introduced before it, through any measurable map, however
nonlinear. A vehicle's position is drawn given its previous position, pushed through a motion
model; a range measurement is drawn given the vehicle and a landmark, through the distance between
them.
> **Definition (A model built forward).** A _model built forward_ is a finite sequence of _quantities_ $i = 1, \dots, n$, each with
>
> 1. a space $(\manifold_i, \mathscr X_i)$ and a $\sigma$-finite _reference measure_ $\nu_i$ on it:
> Lebesgue measure on $\R^{\dimension_i}$ or on an open subset of it, or counting measure on a
> countable set; and
> 2. a _conditional density_: a measurable function
> $\density_i : \manifold_i \times \manifold_1 \times \cdots \times \manifold_{i-1} \to [0, \infty)$
> such that $\density_i(\cdot \mid \point_{ for every $\point_{
> Its _joint density_ is
> $\density(\point) = \prod_{i=1}^n \density_i(\point_i \mid \point_{ product measure $\nu = \nu_1 \otimes \cdots \otimes \nu_n$ on
> $\manifold = \manifold_1 \times \cdots \times \manifold_n$.
A conditional density written as $\density_i(\point_i \mid \point_{ **Proposition (The joint law of a model built forward).** The measure $\mu(\dd\point) = \density(\point)\, \nu(\dd\point)$ is a probability measure on
> $\manifold$. For each $i$, its marginal on the first $i$ quantities has density
> $\prod_{j \le i} \density_j(\point_j \mid \point_{ $\nu_1 \otimes \cdots \otimes \nu_i$, and
>
> $$
> \kernel_i(\point_{ $$
>
> is a regular conditional distribution of the $i$-th quantity given the ones before it.
The function $\density$ is non-negative and measurable on the product. By Tonelli's theorem,
integrating over $\point_n$ first, then $\point_{n-1}$, and so on, each factor
$\density_i(\cdot \mid \point_{
The probability space can then be taken to be $\Omega = \manifold$, with its product
$\sigma$-algebra, $\mathrm P = \mu$, and the quantities the coordinate projections. This is the
analogue of choosing a basis: the quantities are the coordinates. Nothing is lost by the choice,
because every statement the model makes is about the joint law of its quantities, and that law is
$\mu$ on any space that carries them.
Nor is the restriction to finitely many quantities essential. The construction never used anything
about the spaces beyond measurability, and it extends to infinitely many quantities built forward,
with no topological assumption, by a theorem that is the right one for models built this way.
> **Theorem (Ionescu-Tulcea).** Let $(\manifold_i, \mathscr X_i)_{i \ge 1}$ be any measurable spaces, $\mu_1$ a probability measure
> on $\manifold_1$, and for each $i \ge 2$ let $\kernel_i$ be a probability kernel from
> $\manifold_1 \times \cdots \times \manifold_{i-1}$ to $\manifold_i$. There is a unique probability
> measure $\mu$ on $\prod_{i \ge 1} \manifold_i$, with its product $\sigma$-algebra, whose marginal on
> every finite initial product is $\mu_1 \otimes \kernel_2 \otimes \cdots \otimes \kernel_i$.
Kolmogorov's theorem builds a law from consistent marginals, and needs the spaces to be nice for
the marginals to determine one. Ionescu-Tulcea builds it from successive conditional laws, which is
exactly what a model built forward provides, and needs nothing. Appending a quantity is extending a
basis: the old coordinates keep their joint law, and a new one is added with its conditional law
given them.
The second thing the restriction buys is a guarantee about what the model can say. A model speaks
only of its quantities, so the random variables it can express are the $\sigma(\randomPoint)$-measurable ones,
$\randomPoint = (\randomPoint_1, \dots, \randomPoint_n)$, and those are exactly the functions of the quantities.
> **Lemma (Doob–Dynkin).** Let $\randomPoint : \Omega \to \manifold$ be a measurable map from a measurable space $(\Omega, \mathscr F)$ to
> any measurable space $(\manifold, \mathscr X)$, and $\otherRandomPoint : \Omega \to \R$. Then $\otherRandomPoint$ is
> $\sigma(\randomPoint)$-measurable if and only if $\otherRandomPoint = g \circ \randomPoint$ for some measurable $g : \manifold \to \R$.
If $\otherRandomPoint = g \circ \randomPoint$, then $\otherRandomPoint^{-1}(B) = \randomPoint^{-1}(g^{-1}(B)) \in \sigma(\randomPoint)$. Conversely, every set in
$\sigma(\randomPoint)$ is $\randomPoint^{-1}(A)$ for some $A \in \mathscr X$, so an indicator $\mathbf 1_{\randomPoint^{-1}(A)}$ is
$\mathbf 1_A \circ \randomPoint$, and by linearity every $\sigma(\randomPoint)$-measurable simple function is $g \circ \randomPoint$
for a simple $g$. A $\sigma(\randomPoint)$-measurable $\otherRandomPoint$ is a pointwise limit of such functions
$g_k \circ \randomPoint$; the set where $(g_k)$ converges is measurable, and $g = \lim_k g_k$ there and $0$
elsewhere is measurable with $\otherRandomPoint = g \circ \randomPoint$, because the sequence converges at every point of the
image of $\randomPoint$.
> **Corollary (Every expectation is an integral against the joint density).** Let a model built forward be carried by quantities $\randomPoint = (\randomPoint_1, \dots, \randomPoint_n)$ on a probability space
> $(\Omega, \mathscr F, \mathrm P)$, so that $\randomPoint$ has law $\mu$. Every random variable $\otherRandomPoint$ the model
> can express, that is every $\sigma(\randomPoint)$-measurable one, is $g \circ \randomPoint$ for a measurable $g$, and when
> it is integrable,
>
> $$
> \E[\otherRandomPoint] = \int_{\manifold} g(\point)\, \density(\point)\, \nu(\dd\point)
> = \int_{\manifold_1} \!\!\cdots \int_{\manifold_n} g(\point) \prod_{i=1}^n \density_i(\point_i \mid \point_{ $$
The first claim is [%full](statement:doob-dynkin). For the second, the law of $\randomPoint$ is $\mu$, so
$\E[g \circ \randomPoint] = \int g \,\dd\mu$ by the change of variables for image measures, and
$\int g \,\dd\mu = \int g \density \,\dd\nu$ because $\density = \dd\mu / \dd\nu$. The iterated form
is Fubini's theorem, which applies because $g \density$ is $\nu$-integrable.
This is the probabilistic half of the analogy made precise. In linear algebra, a vector of $V$ is
a combination of basis vectors and a linear functional is a combination of coordinates. Here, a
random variable is a function of the quantities, an event of the model is a measurable set of
their values, and its probability, an expectation, a marginal or a moment is an integral against
the product of the conditional densities, in the order the model was built. The integral is
explicit, though it is rarely available in closed form; estimating it is what inference
algorithms are for.
## Conditioning and likelihoods
Observing quantities is conditioning on them, and in this corner it needs no existence theorem:
the conditional law is written down from the densities.
> **Proposition (Conditioning is division of densities).** Split the quantities of a model built forward into observed ones $O$ and free ones $F$, and let
>
> $$
> \density_O(\point_O) = \int \density(\point_F, \point_O)\, \nu_F(\dd\point_F) .
> $$
>
> Then $\density_O$ is a density of the law of $\randomPoint_O$, finite $\nu_O$-almost everywhere, and the
> kernel
>
> $$
> \kernel(\point_O, \dd\point_F) = \frac{\density(\point_F, \point_O)}{\density_O(\point_O)}\, \nu_F(\dd\point_F)
> \qquad \text{where } 0 < \density_O(\point_O) < \infty,
> $$
>
> defined arbitrarily elsewhere, is a regular conditional distribution of $\randomPoint_F$ given $\randomPoint_O$.
By Tonelli, $\density_O$ is measurable and $\int \density_O \,\dd\nu_O = 1$, so it is a density of
the law of $\randomPoint_O$ and is finite $\nu_O$-almost everywhere. The set where it is $0$ has probability
$0$ under that law, and so does the set where it is infinite. For measurable $A$ and $B$,
$$
\mathrm P(\randomPoint_F \in A,\ \randomPoint_O \in B)
= \int_B \int_A \density(\point_F, \point_O)\, \nu_F(\dd\point_F)\, \nu_O(\dd\point_O)
= \int_B \kernel(\point_O, A)\, \density_O(\point_O)\, \nu_O(\dd\point_O),
$$
the two exceptional sets contributing nothing to either side, which is the defining property of a
regular conditional distribution.
The version is chosen, not discovered: changing the densities on a null set, or the reference
measures by a change of coordinates, changes $\kernel$ on a null set of observations, which is all the
freedom the Borel–Kolmogorov paradox exploits. In a model built forward the choice is made once,
by the densities and reference measures the model was written with, and every observation is
conditioned on by that same rule. That is also why two densities may only be compared, or divided,
when they are with respect to the same reference measure.
The same computation explains a factor that scores no quantity. Multiplying the joint density by a
non-negative function and normalizing is a change of measure, and it is always conditioning on data
that the model does not carry as a quantity.
> **Proposition (A likelihood factor is a change of measure).** Let $L : \manifold \to [0, \infty)$ be measurable with $0 < Z = \int L \,\dd\mu < \infty$. Then
> $\mu_L(\dd\point) = L(\point)\, \mu(\dd\point) / Z$ is a probability measure with density
> $L \density / Z$ with respect to $\nu$. If $L(\point) = h(\otherPoint_0 \mid \point)$ for a
> conditional density $h$ of a further quantity $\otherRandomPoint$ appended to the model, then $\mu_L$ is the
> conditional law of $\randomPoint$ given $\otherRandomPoint = \otherPoint_0$, as conditioning by densities gives it.
$\mu_L$ is a probability measure because $\int L \,\dd\mu = Z$, and its density with respect to
$\nu$ is $L \density / Z$. With $\otherRandomPoint$ appended, the joint density is $\density(\point)\, h(\otherPoint \mid \point)$,
so [%full](statement:conditioning), with $\randomPoint$ free and $\otherRandomPoint$ observed at $\otherPoint_0$, gives the
conditional density $\density(\point)\, h(\otherPoint_0 \mid \point) / \density_{\otherRandomPoint}(\otherPoint_0)$,
where $\density_{\otherRandomPoint}(\otherPoint_0) = \int h(\otherPoint_0 \mid \point)\, \mu(\dd\point) = Z$.
So a *likelihood factor* is exactly what its name says: the likelihood of an observation that is
not a quantity of the model, already conditioned on.
## What is left outside
Two things a measure theorist might expect are deliberately absent. A function of quantities, a
*mapped value*, has a law, the image of $\mu$, but in general no density with respect to any
reference measure, so a model built forward never scores one: densities are given for quantities,
which were introduced with them, and for nothing else. And an *improper* seed, a "density" with
infinite integral such as a flat density on a line, is not a probability measure at all, and the
joint law above is then a $\sigma$-finite measure rather than a probability. Conditioning by
division of densities still makes sense when the normalizing integral that results, the $Z$ or
$\density_O$ above, is finite and positive, and then yields a probability measure. Whether it is
depends on what is observed: it is a condition on the query, not on the model.
## How skoro is built on this
Skoro is designed by this philosophy. A skoro model is not a function $\density$ to be evaluated,
but a probability space built forward, quantity by quantity, and the library sets out to let its
users work with that space in every way the measure theory above allows, of which evaluating
$\density$ is only one.
- **Quantities are specified by their augmenting kernels.** Each quantity $\randomPoint_i$ enters the model
with the kernel that augments the space built so far,
$$
\kernel_i(\point_{
---
title: Product manifolds
summary: A product of manifolds is a manifold whose tangent space stacks its factors', and whose retraction acts factor by factor; derivatives in its charts are block matrices of the factors'.
depends_on: [retraction, derivatives-in-charts]
---
# Product manifolds
Several quantities taken together, such as a pose and a velocity, or every state in a window of an
estimator, are one point of a product manifold. This note shows that the product inherits a
retraction from its factors, and how derivatives in its charts break into blocks.
## The product and its tangent space
Let $\manifold_1, \dots, \manifold_k$ be manifolds of dimensions $\dimension_1, \dots, \dimension_k$, each with a
[retraction](/math/retraction.md). Their product $\manifold = \manifold_1 \times \cdots \times \manifold_k$
is a manifold of dimension $\dimension = \sum_i \dimension_i$, whose points are tuples $\point = (\point_1, \dots, \point_k)$. Its
tangent space at $\point$ is the direct sum of the factors', identified with $\R^\dimension$ by stacking:
$$
\T_\point \manifold = \T_{\point_1}\manifold_1 \oplus \cdots \oplus \T_{\point_k}\manifold_k \;\cong\; \R^\dimension,
\qquad \tangent = (\tangent_1; \dots; \tangent_k) . @tag(stack)
$$
(The $\oplus$ of vector spaces here is the direct sum, not the retraction.)
> **Proposition (The product retraction).** If each $R^i$ is a retraction on $\manifold_i$, then
>
> $$
> R_\point(\tangent) = \bigl(R^1_{\point_1}(\tangent_1), \dots, R^k_{\point_k}(\tangent_k)\bigr)
> $$
>
> is a retraction on $\manifold = \manifold_1 \times \cdots \times \manifold_k$, its inverse is
> $\otherPoint \ominus \point = (\otherPoint_1 \ominus \point_1; \dots; \otherPoint_k \ominus \point_k)$, and it is flat when every factor is.
$R_\point(0) = (\point_1, \dots, \point_k) = \point$. The derivative of $R_\point$ at $0$ is block diagonal, its $i$-th
block the derivative of $R^i_{\point_i}$ at $0$, which is $I_{\dimension_i}$; so it is $I_\dimension$. The inverse and
flatness follow factor by factor.
## Derivatives in product charts
Let $\map : \manifold \to \otherManifold$ be smooth, with $\otherPoint = \map(\point)$. Because a step in the product chart
moves each factor in its own chart, the [derivative in charts](/math/derivatives-in-charts.md) of $\map$ is
a row of blocks,
$$
\D \map(\point) = \bigl[\, \D_1 \map(\point) \;\; \cdots \;\; \D_k \map(\point) \,\bigr],
\qquad
\D_i \map(\point)\, \tangent_i = \frac{\dd}{\dd t}\Big|_{t=0} \bigl(\map(\point_1, \dots, \point_i \oplus t \tangent_i, \dots, \point_k) \ominus \otherPoint\bigr),
@tag(blocks)
$$
each the derivative in the $i$-th factor's chart with the other factors held fixed. A map into a
product is a column of blocks in the same way, one per factor of the output. For a scalar $f : \manifold \to \R$, the
chart gradient stacks the factors' partial gradients, $g = (g_1; \dots; g_k)$, and the chart
Hessian is the block matrix whose $(i, j)$ block differentiates $g_i$ along the $j$-th factor:
$H_{ij}(\point)\, \tangent_j = \frac{\dd}{\dd t} g_i(\point_1, \dots, \point_j \oplus t \tangent_j, \dots, \point_k)$ at
$t = 0$.
In particular, the retraction of the product is, as a map of the pair $(\point, \tangent)$, the factors'
retractions side by side, and its derivatives of every order are block diagonal in the factors.
Nothing about a product's derivatives requires more than its factors' derivatives and the
bookkeeping of [](eq:stack).
---
---
title: Retractions and tangent charts
summary: A retraction moves a point of a manifold along a tangent vector; it makes each tangent space a chart around its point, and gives the operations ⊕ and ⊖.
---
# Retractions and tangent charts
Estimation moves points. On $\R^\dimension$ a step is an addition, $\point + \tangent$; on a curved manifold, such
as the rotations or the unit sphere, there is no addition, and a *retraction* takes its place. This
note defines it, shows how it turns every tangent space into a chart, and introduces the notation
$\oplus$ and $\ominus$ used throughout these notes.
## The definition
Let $\manifold$ be a smooth manifold of dimension $\dimension$. In computation, a tangent space
$\T_\point\manifold$ is never handled abstractly: each is identified with $\R^\dimension$ by a basis chosen
smoothly in $\point$, and a tangent vector is a column of $\dimension$ numbers. We assume such an identification
throughout. Some manifolds have none that is smooth everywhere: by the hairy ball theorem, no
smooth choice of tangent bases covers the whole sphere $S^2$. There the identification is smooth on
an open part of the manifold, and a computation stays inside it.
> **Definition (Retraction).** A _retraction_ on a manifold $\manifold$ of dimension $\dimension$ is a smooth map
> $R : \T\manifold \to \manifold$, written $(\point, \tangent) \mapsto R_\point(\tangent)$, such that for every
> $\point \in \manifold$, with $\T_\point\manifold$ identified with $\R^\dimension$,
>
> 1. $R_\point(0) = \point$, and
> 2. the derivative of $R_\point$ at $0$ is the identity: $\D R_\point(0) = I_\dimension$.
The definition is the one used in optimization on manifolds, where retractions replace the
exponential map of a Riemannian manifold. The exponential map is a retraction, but
so is any map that agrees with it to first order, which is usually cheaper. Nothing here needs a
metric.
Two examples. On $\R^\dimension$, $R_\point(\tangent) = \point + \tangent$. On the positive reals, with
$\T_\point(0, \infty)$ identified with $\R$ by measuring relative change,
$R_\point(\tangent) = \point\,e^{\tangent}$: it is defined for every $\tangent$, never leaves the positive reals,
and satisfies $\frac{\dd}{\dd\tangent} \point e^\tangent \big|_{0} = \point$, which is the identity
in the chosen identification (a tangent at $\point$ of coordinate $1$ is the velocity $\point$).
## Tangent charts
Because $\D R_\point(0)$ is invertible, the inverse function theorem makes $R_\point$ a diffeomorphism
from a neighbourhood of $0 \in \R^\dimension$ onto a neighbourhood of $\point$. So $R_\point$ is a *chart* of
$\manifold$ around $\point$, the tangent chart at $\point$, and every point near $\point$ has coordinates in it.
As is common in robotics, we write
$$
\point \oplus \tangent \defeq R_\point(\tangent), \qquad \otherPoint \ominus \point \defeq R_\point^{-1}(\otherPoint), @tag(oplus)
$$
the second for $\otherPoint$ in the neighbourhood where $R_\point$ is invertible. In words, $\otherPoint \ominus \point$ is where
$\otherPoint$ sits in the chart at $\point$. By construction,
$$
(\point \oplus \tangent) \ominus \point = \tangent @tag(roundtrip)
$$
for every $\tangent$ in the chart's domain, and $\point \oplus (\otherPoint \ominus \point) = \otherPoint$ for every $\otherPoint$ in its
image. On $\R^\dimension$, $\oplus$ and $\ominus$ are addition and subtraction; on the positive reals above,
$\otherPoint \ominus \point = \log(\otherPoint / \point)$.
## The mirror identity
On $\R^\dimension$, $\point \ominus (\point \oplus \tangent) = -\tangent$ as well: the step back from $\point \oplus \tangent$ to
$\point$ is the reverse of the step forward. On a curved manifold it generally is not, because the left
side is measured in the chart at the *moved* point $\point \oplus \tangent$, not at $\point$.
> **Lemma (The mirror identity holds to first order).** For $\tangent$ small enough, $\point \ominus (\point \oplus \tangent) = -\tangent + O(\lVert\tangent\rVert^2)$.
Write $\otherPoint = \point \oplus \tangent$ and $\varphi(\tangent, c) = R_\point^{-1}(R_{\point \oplus \tangent}(c))$, the chart
at $\otherPoint$ read in the chart at $\point$. It is smooth near $(0, 0)$, with $\varphi(\tangent, 0) = \tangent$
and, since $R_{\point \oplus 0} = R_\point$, $\varphi(0, c) = c$. Hence
$\varphi(\tangent, c) = \tangent + c + O(\lVert\tangent\rVert\,\lVert c\rVert)$. The coordinate
$c = \point \ominus \otherPoint$ is the one with $\varphi(\tangent, c) = R_\point^{-1}(\point) = 0$, so
$c = -\tangent + O(\lVert\tangent\rVert\,\lVert c\rVert)$, and $\lVert c \rVert = O(\lVert\tangent\rVert)$
gives the claim.
The identity holds exactly in some cases: on $\R^\dimension$, on the positive reals above, and on any Lie
group whose retraction is $R_\point(\tangent) = \point \operatorname{Exp}(\tangent)$, since then
$\point \ominus (\point \oplus \tangent) = \operatorname{Log}(\operatorname{Exp}(\tangent)^{-1}) = -\tangent$. It
fails at second order for other retractions, and a computation that needs "the point whose
coordinate is $-\tangent$" must not assume it.
A manifold whose retraction is $\point + \tangent$, so that both identities hold exactly and every chart is
the same affine chart, is called *flat* in these notes.
---
---
title: "The core"
pages: []
---
# The core
The core of skoro, built to [the charter](/archive/charter.md): the API for building models, with
quantity types (the Euclidean manifolds, products, and the discrete types), the primitive
operations and their derivatives of first and second order, composition, the compilation graph
and its compiler, the probabilistic graph, terms and bindings, the writer with its snapshots,
branches and building scopes, addresses, explicit random keys, and the layer over burn. It ends where models can be built and
differentiated but no distribution exists yet, except private ones in tests. Much of it is ported
from the `archive` branch, reshaped to the charter. When it is done, `REVISIT.md` lists every
feature of the archive not yet brought back, each to return as an entry of its own.
---
---
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.
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).
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] .
> $$
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$.
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.
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)$.
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?
---
---
title: "What skoro is"
pages: [start/what-is-skoro]
---
# What skoro is
The library's charter: what skoro is for, who it serves, and the positions on which everything
else in it rests. It was written when the repository started over, from an interrogation of the
design that came before (kept on the `archive` branch), so that each later entry has a stated
purpose to build towards and to be measured against. It settled skoro's purpose and its users, and
twelve positions on the shape of its core; the home page and
[What is skoro?](/start/what-is-skoro.md) were written from it.
---
---
summary: A model is written once, generatively, and any inference algorithm reads what it needs from the model's structure, refusing by structure when the model does not fit.
---
# Context
## Context
### The problem skoro exists for
Estimation software for robots and autonomous systems is, for the most part, factor-graph
software. In GTSAM and its relatives a model is a set of *factors*: residual functions whose
weighted squares are summed and minimized. That formulation has two costs.
The first is that it is not generative. A factor is one particular reading of a conditional
density (its negative logarithm, seen as a least-squares residual), and the conditional itself is
never written down. Consider a range measurement $z$ between a vehicle at $x$ and a target at $y$:
$$
z \sim \mu(\mathrm{d}z \mid x, y)
$$
That one statement covers three situations: estimating $x$ and $y$ jointly from $z$, estimating
$x$ from $z$ when $y$ is known, and simulating $z$ from given $x$ and $y$. A factor-graph library
needs a different factor class for the first two, different key management for each, and has no
account of the third at all, since a residual cannot be sampled from. Even a Kalman filter does
not fit naturally: the filter is a statement about a hierarchy of conditional distributions, and
the factor graph holds only the residuals of the posterior that hierarchy implies.
The second cost is that a factor graph serves one algorithm. It is built for maximum a posteriori
estimation of a nonlinear least-squares problem, and an algorithm that needs anything else (a
sampler, a filter, a variational approximation) has to be told the model's structure again, in its
own terms.
Probabilistic programming languages make the opposite choice. Gen.jl in particular separates
*modeling* from *inference*: a model is a generative program, and inference algorithms are written
against an interface to it. That separation is the one skoro wants. What those languages are not
built for is the model that lives over time. In sequential estimation the model grows as
measurements arrive, its pieces are not alike (different sensors, different measurement models at
different times), old states are retired as they stop mattering, and estimation runs concurrently
with all of it. A language that builds a program once and then runs inference on it does not
manage that well, and is rarely embedded in a real-time Rust system.
### What came before
The repository started over on 2026-10-09. The library as it stood is kept on the
[`archive`](https://github.com/mvarble/skoro/tree/archive) branch: fourteen crates, a model client
over a shared graph, a hand-written engine for first- and second-order derivatives on manifolds,
the Normal distribution, SO(3) and SE(3), and a block-sparse Gauss–Newton machinery with
marginalization for sliding windows. Most of it worked, and much of it is coming back. What it
taught, and what this document responds to, is a short list:
- **Every operation hand-wrote four derivative maps**, which made each new distribution expensive
and was where wrong derivatives hid.
- **Conditioning was a mutation.** Observing a quantity changed the model, for every reader,
which made "simulate, then condition, then infer" and posterior predictive checks awkward.
- **One lock guarded the model**, and every handle carried access to it.
- **Type erasure lost information.** A manifold assembled at run time could retract a point but
could not differentiate the retraction, so product manifolds were out of reach.
- **Implementations were duplicated** across compile-time and run-time dimensioned types.
- **Randomness was global**: no seed, and no way to sample one quantity on its own.
- **Fourteen crates** kept useful internals out of reach of their siblings.
## Purpose
> [!DECISION] What skoro is
> Skoro is a library for writing probabilistic models once, generatively, and running different
> inference algorithms on the same model. The model is built forward, the way a user thinks of it:
> each uncertain quantity, and the conditional distribution it has given the quantities it depends
> on. The library records the structure as it is built, much as automatic differentiation records
> a computation, together with whatever later computations need. An inference algorithm reads that
> structure. If the model meets the algorithm's assumptions, the algorithm runs; if not, it says so
> with an error. The user never describes the model a second time for the algorithm's benefit.
Three kinds of user follow from this, in order of priority:
1. **People developing inference algorithms.** The reading surface (what an algorithm can discover
about a model, and how cheaply) matters as much as the building surface. Exploring new
algorithms is skoro's first purpose.
2. **People building real-time estimators** for robots and autonomous systems: heterogeneous
sequential models that grow, retire old states, and are estimated while measurements stream in,
on CPU or GPU, from Rust, often under several hypotheses at once about what the measurements
came from.
3. **People doing Bayesian inference** with quantities that live on manifolds, such as rotations
and poses, or that need second-order derivatives.
Skoro is not a general probabilistic programming language. Each of its models has a fixed
structure, and structure that is itself random is handled one level up, over a family of models
(see below). It does not compete with Stan or NumPyro on breadth of distributions or with GTSAM on
the maturity of its least-squares solvers. What it offers that they do not is the combination: one
generative model, a structure that algorithms can read, manifold-aware second-order derivatives,
and a model that lives across threads and over time.
## Positions
### 1. A model is the joint; conditioning belongs to the query
A model holds the generative structure and nothing else: the joint distribution of its quantities,
never a posterior. Which quantities are observed, and their values, is part of each *query*: an
inference call, a density evaluation, a generation. `observe` builds such an observation set; it
does not change the model.
The same model can then serve several questions at once, from several threads: one estimator
treats the target's position as known, another estimates it, a third simulates measurements. Two
uses are worth naming, because they motivated the position:
- **Checking an inference algorithm against simulation.** Generate $(x, y)$ from the joint
$\mu(\mathrm{d}x, \mathrm{d}y)$, condition on the generated $y$, infer
$\hat\kappa(\cdot \mid y)$, and compare it with the generated $x$ (the error of a point estimate,
or the concentration of the approximate posterior around $x$).
- **Checking an estimator's consistency.** Generate predicted measurements
$z \sim \kappa(\mathrm{d}z \mid \hat{x})$ from an estimate and compare them with the
measurements actually received. With conditioning in the model, this needed a way to pretend the
model was joint again; with conditioning in the query, it is one more query.
### 2. Models are built forward, change by appending and retiring, and random structure lives one level up
A model grows by adding quantities, each together with its conditional distribution. The one
shrinking change is *retirement*: in sequential estimation, quantities that were once estimated
eventually stop mattering, so they are marginalized onto what remains and then removed. Appending
and retiring are the only mutations a model sees.
Within one model, structure is fixed: it is whatever the builder created, and a random draw does
not add or remove nodes. Random structure is nevertheless central to skoro's purpose (which
measurement came from which target, whether a target was born or has died, which motion model is
in force), and it is handled in two ways.
**Bounded structure, inside one model.** When the structure has a known bound, it is expressed
with discrete quantities and the selection primitive (position 7). Data association is a
categorical choice per measurement among target slots and clutter; a random number of objects up
to a maximum is a set of slots, each with an existence variable gating its terms; a switching model
is a discrete mode selecting among parameter sets; a nonparametric prior is truncated. The graph
never changes, and samplers move the discrete quantities. This costs computation on empty slots
and needs a bound, and asks nothing new of the core.
**Unbounded structure, over a family of models.** In general, the random thing is *which model*.
A structural hypothesis $S$ (this measurement belongs to that target; a target was born here) lives
outside any graph, with a prior $p(S)$, and each hypothesis is a model built from it by ordinary
appending and retiring. Inference over structure (multiple-hypothesis tracking, sequential Monte
Carlo over hypotheses, reversible-jump Markov chain Monte Carlo with birth and death moves) keeps
one model per live hypothesis, runs ordinary inference inside each, and weighs the hypotheses by
$p(S)$ and an estimate of each model's evidence. A birth is an append on one hypothesis's model,
a death a retirement, and pruning a hypothesis is dropping its model.
Multiple-hypothesis estimation is a first-class use of skoro, and the core is shaped for it in two
ways. Hypotheses that differ by one association share almost everything, so models *branch*
cheaply from one another (position 3). And the same quantity must be recognizable across
hypotheses, so that they can be compared, merged and pruned, which is what *addresses* are for
(position 4). The algorithms that infer structure are later entries; the core only has to make
them possible.
This is close to what Gen.jl does, without what made it impractical here. A Gen trace records which
random choices a program made, which fixes its structure, together with their values. In skoro the
structure is an explicit hypothesis and the values live in a typed graph, so nothing has to reflect
over typed programs.
### 3. One writer per branch, many snapshots
A model has one writer, which appends and retires, and publishes immutable *snapshots* that share
structure with earlier ones. Inference runs on a snapshot, without a lock: a long sampling run
never blocks the thread adding measurements, and the thread adding measurements never changes a
model out from under an estimator. Memory belonging to retired quantities is reclaimed when the
last snapshot holding it is dropped.
A snapshot can also be *forked*: turned into a new writer whose appends and retirements form a
branch of their own, copying nothing that the branch does not change. Branches are how a model
under several structural hypotheses is held (position 2): each hypothesis is a branch, a new
measurement associated three ways is three cheap forks, and a pruned hypothesis is a branch nobody
holds any longer. There is one writer per branch, never two writers on one.
This replaces the archive's design, in which a shared graph sat behind one lock, every handle
carried a client that could change it, and drop guards fed a mark-and-sweep collector. It also
removes the need for a separate read-only access level: a snapshot is read-only by construction.
**Building reads like ordinary code.** A model is written as functions of quantities, with no
builder in their signatures. An operator and a function a user writes are the same kind of thing:
```rust
// The range from a vehicle at `x` to a target at `y`: an ordinary function of quantities.
fn range(x: Q>, y: Q>) -> Q {
(y - x).norm()
}
```
(The names in these sketches are illustrative; the core's implementation settles them.) `y - x`
and `range(x, y)` both add nodes to the model, and neither names the writer. They can do
so because building happens inside a *building scope* that a writer opens, and every operation on
a quantity adds its node to the writer whose scope is open on the current thread:
```rust
let x = writer.build(|| {
let x = vehicle_prior.sample_as("vehicle/0");
let z = Normal::new(range(x, target), sigma).sample_as("range/0");
x
})?;
```
The writer is the only thing that changes a model, and a scope is how it lends itself to the code
that builds; outside a scope, an operation on a quantity is refused. A `Q` is then a plain typed
identifier, cheap to copy: it holds no reference to a writer, so it can be sent to other threads and used in
queries against any snapshot that contains its quantity. The same `Q` is valid on every branch
that shares the quantity, which multiple-hypothesis building depends on: each hypothesis's writer
opens its own scope, and the same function of the same quantities builds into each.
```rust
for (hypothesis, writer) in branches.iter_mut() {
writer.build(|| associate(measurement, hypothesis.target(), hypothesis.clutter()))?;
}
```
Errors cannot be returned from `y - x`, and an operator that panicked on a shape mismatch would
break the rule that errors are values. So an operation that fails inside a scope still returns a
`Q`, marked failed: anything built on it is failed too, and the scope reports the first failure,
with where it happened, as the error of `build`. Nothing a failed scope built is published.
**How a scope works.** The exclusive access a writer needs is taken once, when the scope opens,
not by each operation: `build` borrows the writer mutably for the length of the closure, so the
compiler guarantees nothing else touches it meanwhile. Opening the scope stages a copy-on-write
draft of the latest snapshot and registers that draft with the current thread. An operation finds
the draft through the thread and borrows it only for as long as it takes to insert its own node,
never while user code or a distribution runs, so a distribution's `sample_as` can build its
log-density with ordinary operations inside the same scope. The draft belongs to one thread and is
invisible to readers, who hold immutable snapshots; when the scope ends without a failure, the
draft becomes the latest snapshot in a single atomic step, so a reader sees the model as it was
before the scope or as it is after it, never half of it. If the closure panics, the draft is
unregistered and discarded. Because a thread-local cannot be generic over the backend, the draft is
held type-erased and recovered by each operation; an operation that finds no draft, or a draft for
another backend, returns a failed `Q`. Scopes nest, and an operation builds into the innermost.
### 4. Static at the edges, dynamic inside, and everything nameable
The building API is typed. A quantity is a `Q`, where `T` names its type (a manifold or a
discrete set), and distributions and operations are generic, so the compiler checks as much of a
model as it can while the model is written. Inside, the graph is type-erased: it holds quantities
of different types, added at run time from streams of data, and that needs dynamic dispatch. The
boundary between the two is deliberate, not accidental.
Everything the graph stores type-erased is *nameable*: each quantity type, operation and
distribution has a registered identifier and plain parameters. That is what lets a model be
serialized (an autonomous system can write its model out when it fails, to be inspected later),
and what will let Python bindings drive a fixed registry of distributions and operations while new
ones are written in Rust. Neither the format nor the bindings are part of the core; the principle
is, because retrofitting it is far costlier than respecting it from the start.
Quantities are nameable too. Each has an *address* chosen by its builder, such as "the position of
target 3 at tick 40", which stays the same across snapshots, branches and serialization, beside
the internal identifier the graph uses. Addresses are what let two branches agree that they hold
the same quantity, what random keys are derived from (position 10), and what a serialized model is
read back by.
**Addresses are hierarchical, and sub-models are functions.** An address is a path of segments,
each a name or an integer (tick `t` and measurement `j` are numbers, and stay numbers), displayed
as `40/target/3`. A *namespace* is a frame on the same stack as a building scope: inside
`scope::namespace(t, || ...)`, a quantity added as `x` is addressed `t/x`. A sub-model is then an
ordinary function of quantities that names its own quantities locally and knows nothing about where
it is placed; the code that calls it places it:
```rust
fn step(prev: Q, u: &Odometry, zs: &[Range]) -> Q {
let x = Normal::new(prev.compose(u.delta()), q).sample_as("x");
for (j, z) in zs.iter().enumerate() {
scope::namespace(j, || Normal::new(range(x, landmark(z.id)), sigma).sample_as("z"));
}
x
}
let mut x = writer.build(|| prior.sample_as("x"))?;
for (t, data) in stream.enumerate() {
x = writer.build(|| scope::namespace(t, || step(x, &data.u, &data.zs)))?;
}
// "x", then "0/x", "0/0/z", "0/1/z", "1/x", "1/0/z", ...
```
Sub-models connect through their arguments and results, as Gen.jl's generative functions do, never
by looking each other's quantities up by name; looking a quantity up by address, relative to the
open namespaces, remains for code that resumes building from a snapshot. Namespaces nest freely.
They are also the natural unit of retirement: when tick 0 leaves a sliding window,
`writer.retire("0")` marginalizes everything under `0/` onto its Markov blanket and removes it,
which makes declarative the bookkeeping the archive did by hand on every tick. A namespace is part
of an address within one model, and is unrelated to branching: `1/x` names the same quantity in
every branch that contains it.
Three things are deliberately absent: code generated and compiled ahead of time (SymForce's
approach does not fit models built at run time), macro languages for models (hard to reason about
and unfriendly to editors), and panics across the public API (errors are values, which bindings
will need).
### 5. Derivatives come from composition, and their order is reported
Skoro differentiates for itself. burn's `Autodiff` backend was the obvious candidate and was ruled
out for two reasons: it has no second order (Hessian–vector products, which some samplers need,
require differentiating the backward pass), and it pulls back one cotangent vector per pass, while
a Jacobian-based optimizer wants a whole *frame* of them in one. Tangent-space derivatives on
non-Euclidean manifolds looked possible inside burn, but nothing guaranteed them deep in its
implementation.
A small set of *primitive* operations carries hand-written derivative maps (forward, JVP, VJP and
HVP, over frames), written once and tested hard. Everything else, distributions' log-densities and
manifolds' retractions included, is a *composite*: a subgraph of primitives, differentiated by the
chain rule. A composite is stored as a named, reusable subgraph and expanded when the model is
compiled, so a serialized model reflects the calls its user made. A hand-written fused version of a
composite is an optimization, tested against the composite, and some operations will need one for
numerical care (the logarithm of SO(3) near angle $\pi$, say).
Each operation states the highest order of derivative it supports, and a compiled program's order
is the least over its steps. A pass that needs more is refused before it runs, naming the operation
responsible. The case that motivates this is the bridge to neural networks: a model may contain an
operation that runs burn's `Autodiff` internally, translating between skoro's tensors and burn's at
its boundary. Such an operation supports first order only, and a second-order algorithm on that
model must fail with that explanation rather than with zeros.
### 6. Algorithms refuse by structure
The rule that runs through every position here: *what a model can do is read off its structure,
never declared by its user.* There are no capability flags and no capability traits. Instead,
each fact an algorithm depends on is recorded where it naturally lives, and the algorithm checks it
before it starts:
- the order of derivative each operation supports (above);
- whether a quantity is discrete, and whether it is observed (below);
- whether each density term carries the binding an algorithm reads (below);
- whether each distribution has a sampler, and whether that sampler is exact (below).
A refusal names what is missing and where. An algorithm that silently skips what it cannot read
(the archive's least-squares readers dropped any term without a residual) must be paired with a
check that nothing was skipped, since a well-conditioned solve of the wrong system is the worst
outcome an estimator can have.
### 7. Quantities: runtime dimensions, typed wrappers, two kinds, and products
**One implementation per idea.** Quantity types are dimensioned at run time (`Euclidean(n)`, a
matrix of `m` rows and `n` columns), and operations are written once against them. Fixed-size
forms such as `Euclidean<3>` are thin typed wrappers that convert to the runtime ones, so the
compiler catches shape errors where sizes are known and nothing is implemented twice. Shape
mismatches between runtime quantities are errors when the node is added.
**Two kinds of quantity.** A continuous quantity lives on a manifold: it has a tangent space, a
retraction, and an inverse retraction, and estimation moves it through its tangent space. A
discrete quantity takes values in a set and has none of these. Values of both kinds are stored as
tensors (floating-point and integer), and an operation's inputs may be of either kind. The rules
that follow:
- derivatives stop at discrete quantities: an operation is differentiated in its continuous inputs
only, and a discrete quantity is never a coordinate of a gradient, Jacobian or Hessian;
- a gradient-based algorithm on a model with an unobserved discrete quantity refuses, saying so;
- *selection*, choosing among continuous values by a discrete index, is a primitive, differentiable
in the values and not in the index, because mixtures, switching models and data association are
built from it;
- `Categorical(n)` and an unbounded `Integer` are the discrete types, and `Boolean` is a typed name
for `Categorical(2)`;
- a quantity with a continuous part and a discrete part, such as a pose with a mode, is two
quantities.
Distributions over discrete quantities, and samplers that move them, arrive later through the
workflow. The kinds themselves are in the core because adding them later would touch everything.
**Products, provided the manifold interface builds subgraphs.** A product of quantity types is a
quantity type: its tangent space is the direct sum of its factors', and it retracts factor by
factor. Products make a joint distribution over several quantities ordinary (it scores one quantity
of a product type, and the quantities a user works with are projections of it), and they give a
marginal over a Markov blanket of run-time size a type to live on.
The archive could not build products, because a manifold *returned* its retraction as a typed
operation, which type erasure made opaque. Here a manifold *builds* its retraction into the graph:
```rust
trait DynManifold {
fn dimension(&self) -> usize; // the tangent space's
fn storage_size(&self) -> usize; // the stored point's
// Adds the nodes computing `point ⊕ tangent`, and returns the result.
fn build_retraction(&self, g: &mut dyn GraphBuilder, point: Node, tangent: Node) -> Node;
fn build_inverse_retraction(&self, g: &mut dyn GraphBuilder, a: Node, b: Node) -> Node;
}
```
A product's retraction is then a subgraph like any other: slice the point and the tangent into
blocks, let each factor build its retraction on its block, and concatenate. Its derivatives, of
both orders, follow from composition, and nothing about the factors' types needs to be recovered.
The position is accepted on that condition: the core's implementation must show that the
retraction of a runtime product of two Euclidean factors and a real line, and its first and second
derivatives, agree with the factor-wise computation. If it cannot, this position is reopened.
### 8. Terms have a fixed shape; roles are optional and arrive with their readers
A density term is the log-density node, the quantity it scores (a joint factor over existing
quantities scores none), and the quantities it is conditioned on. The model records all three when
the term is added, and readers use them for structural questions: which terms touch a quantity,
what its Markov blanket is, what retiring it disturbs.
Anything further a term offers is a *binding*: a labelled reference from the term to a node that
plays a *role* for it. A Gaussian term binds its whitened residual, its mean and its covariance
factor; a term with a closed-form expected information binds a root of it. A least-squares reader
finds the residuals, and a Fisher-scoring reader the roots, without the user saying anything: the
distribution declared them as it was added. Bindings count as references, so what a term binds
lives as long as the term.
The core defines no role at all. Each role arrives with the distribution or reader that gives it
meaning. A role is a typed key: it carries its name (its stable identifier) and the multiplicity it
expects (exactly one, at most one, or any number), and a lookup returns typed handles accordingly,
so the check that a term bound what it should happens once, in the lookup, rather than in every
reader.
### 9. A distribution has a density, and may have a sampler with a density of its own
A distribution provides the log-density of its term, and may provide a sampler. A sampler that
draws exactly from the distribution says so. One that draws from something else reports its own
log-density, which is how an improper prior (a flat density on a line) can still be generated
from, by a broad Gaussian say, with the mismatch accounted for. The importance weight of a
generated point is then a sum over terms:
$$
\log w(x) = \sum_t \bigl(\log p_t(x) - \log q_t(x)\bigr), @tag(weight)
$$
where $q_t = p_t$, and the summand vanishes, for a term whose sampler is exact. A distribution with
no sampler at all is legal; generating from a model that contains one is refused, by structure.
This replaces the archive's two flags on each term ("in the target", "in the proposal"), and the
rule that a quantity generated but not scored had to be scored by a separate factor.
### 10. Randomness is explicit
Every sampling call takes a random key, which splits deterministically into independent keys, as in
JAX. One seed reproduces a whole run; each chain, particle set or thread has its own stream; no
hidden state couples them. The key a quantity's sampler receives is derived from the run's key and
the quantity's address, never from the order in which nodes were added, so that sampling one
quantity disturbs no other, retiring old states does not change the draws of new ones, and two
branches holding the same quantity draw it alike.
burn exposes no control over its random-number generation, so a sampler draws its uniform variates
from the key on the host and moves them to the device when it must. That costs a transfer, never
reproducibility.
### 11. burn behind a layer of skoro's own; precision chosen per query; solves on the host
burn was chosen for performant tensor operations that run on CPU and GPU from one Rust codebase,
with cubecl making custom kernels simple to write. The backend trait is sealed to burn's CPU backend
and its cubecl GPU backends: `Autodiff` never runs a model (it appears only inside the bridging
operation), and skoro supports only the backends it tests.
Operations reach burn through a small module of skoro's own, holding the tensor shapes skoro uses
(batched points, tangent frames, integer values) and the tensor operations its primitives need.
burn is not yet at version 1 and breaks between releases; an upgrade should touch that module and
the primitives' kernels, never a distribution or an estimator. Composition keeps that module small,
since only the primitives touch tensors.
**Precision belongs to the query.** burn's backend type does not fix float precision: the element
type it names is largely ignored, and precision is decided at run time, by a process-wide default
fixed at first use. Precision is also a choice that differs by estimator, not by model:
Levenberg–Marquardt wants double precision, because its conditioning demands it, while a sampler
on a GPU wants single precision, for speed. So:
- **A model commits to no precision.** Constants, observed values and initial points are plain
data in the graph, stored on the host in double precision. They must be nameable and serializable
anyway (position 4), and storing them as device tensors would fix a precision the moment the
model is built.
- **The estimator chooses precision when it compiles or evaluates.** Compiling a program for a
query takes a device and a precision, and that is when stored values become device tensors of
the requested type. Compiled programs are cached by snapshot, targets, device and precision, so
one snapshot can serve a double-precision solver on the CPU and a single-precision sampler on a
GPU at once, as it serves different conditionings at once (position 1).
- **Skoro's tensor layer manages precision, and never defers to burn's default.** Every tensor
skoro creates is created at an explicit precision, and every conversion is an explicit cast.
Primitives' kernels are written for both precisions and tested at both. A device that cannot
compute in double precision, as many GPUs through `wgpu` cannot, is reported when a query asks it
to, before anything runs, like any other unmet requirement (position 6).
- **Random draws do not depend on precision.** A sampler draws its uniform variates from its key in
double precision on the host and casts them, so the same key gives the same draw at either
precision, up to rounding.
**The solves where conditioning lives run on the host, in double precision.** Dense
factorizations (the QR of a least-squares step, Cholesky, triangular solves, the elimination behind
marginalization) run on the host in double precision, through a dense linear-algebra library
(faer, as in the archive), whatever precision the model was evaluated at. burn has no QR
factorization, and the alternative, forming normal equations on the device, squares the condition
number. The arrangement also serves the mixed case: residuals and Jacobians evaluated in single
precision on a GPU, for speed, and the step solved in double precision on the host, where the
conditioning matters. Such a Jacobian arrives already rounded, which is the estimator's trade to
make knowingly; the solve itself is never the weak point. This is a principle for the estimators'
entries to build on; how each one uses it is theirs to argue.
### 12. Two crates and an umbrella
The core is two crates: `skoro-tensor`, the sealed backend and the layer over burn, and
`skoro-core`, everything that makes a model a model (identities and errors, quantity types and
manifolds, the primitives and the differentiation engine, the graph and its compiler, terms and
bindings, the writer and its snapshots, and the typed building API), separated by modules rather
than crates. The `skoro` crate re-exports both. Each later feature (a distribution, a manifold
such as SO(3), an estimator) is a crate of its own, built only on the core's public API, which
keeps that API honest: whatever a feature needs, the core must expose.
## Consequences
What becomes easy:
- running several algorithms, and several conditionings, on one model, concurrently;
- holding a model under many structural hypotheses at once, as cheap branches of one another;
- composing a model from sub-models written as functions, placed under namespaces, and retiring a
whole namespace (a tick leaving a window) in one call;
- simulation-based checks of inference algorithms and consistency checks of estimators;
- adding a distribution or a manifold operation, which is mostly composition of existing pieces,
with derivatives of both orders for free;
- writing an algorithm that discovers what it needs and refuses clearly when it is absent;
- reproducing any run that involved randomness, from one seed;
- writing a model out for later inspection, and, later, driving one from Python.
What becomes hard, or costs something:
- **Random structure is not automatic.** A model's structure is fixed; structure that is itself
random needs either a bound (slots and discrete switches) or an algorithm over a family of
branches, and those algorithms are later work. Programs whose shape changes as they run, in the
manner of Gen.jl, are not supported.
- **Builders choose addresses.** Every quantity needs a meaningful, stable address, which is one
more thing to name while building.
- **Compiled programs are longer.** Composition yields more, smaller steps than hand-written fused
operations, until profiling earns a fusion.
- **Every reader takes an observation set**, and "is this quantity free?" becomes a question about
the model and the query together.
- **Snapshots and branches are a new design**, untested in the archive.
- **Building scopes are implicit state.** Which writer an operation builds into is decided by the
scope open on the current thread, not by anything in its signature; an operation outside a scope
is refused when it runs, not when it compiles, and a failure inside one surfaces only when the
scope ends.
- **Registration has a cost.** Every type-erased thing must be named and registered, built-in or
not.
- **Products are conditional**, and the core must prove them before they are relied upon.
- **Explicit keys** add an argument to every sampling call.
- **Precision costs a conversion.** A stored value becomes a device tensor, at the precision a
query asks for, once per compiled program; and every kernel must be correct, and tested, at two
precisions.
## Alternatives
**Factors as the model** (GTSAM). Rejected for the reasons in the context: a factor is one reading
of a conditional density, serves one algorithm, and cannot be sampled from. Factors remain
available as what a least-squares reader *derives* from a generative model.
**Programs as the model** (Gen.jl). Gen's separation of modeling and inference is adopted; its
representation is not. A model that is a program re-run whenever a random choice changes its shape
needs efficient re-execution, which is hard, and reflection over typed programs, which Rust makes
impractical; and a program held as a closure cannot be serialized. Random structure is instead
handled over families of fixed-structure models held as branches (position 2), which keeps each
model cheap for algorithms to read.
**Random structure as out of scope.** Simpler, and tempting, since a graph built forward cannot
express it directly. Rejected because it would abandon multiple-hypothesis estimation, which
real-time estimation in cluttered environments depends on, when forkable snapshots and addresses
make it affordable.
**burn's `Autodiff`.** Rejected for the lack of second order and of frame pullbacks. Kept, through
the bridging operation, for neural networks inside a model.
**Hand-written derivatives for every operation** (the archive). Replaced by composition, which
removes most derivative code and the bugs that lived in it.
**Code generated ahead of time** (SymForce). Rejected: skoro's models are built at run time, from
streams of data, and must not need a compilation step.
**Conditioning by mutation, and one lock around a shared graph** (the archive). Replaced by
query-time conditioning and snapshots, for the reasons in positions 1 and 3.
**Compile-time dimensions only, or both kinds implemented separately** (the archive). Replaced by
runtime dimensions with typed wrappers, so nothing is implemented twice.
**Handles that carry their writer** (as the archive's carried its client). Operators work
anywhere a handle is, but a handle becomes bound to one thread and one branch: it cannot be sent to
an estimator, and building the same function into several hypotheses means rebinding every handle
to each branch. Building scopes keep the sugar and leave handles free.
**Lifetimes on handles instead of a thread-local** (`Q<'s, B, T>`, holding a reference to the
draft). Using a handle outside its scope would then fail to compile rather than fail when run, and
no thread-local would be needed. Declined, for now, because the lifetime appears in every function
a user writes, which is the burden building scopes exist to remove; it is the fallback if the
run-time checks prove error-prone.
**Importing a built model under a namespace.** Possible, since built models are nameable, but an
imported model arrives with its own priors fixed, while a sequential model needs each step's prior
to be the previous step's result. Sub-models as functions with arguments express that; imports do
not, and are left out until a use for them appears.
**An explicit builder in every signature** (`builder.add_mapped_quantity(...)`, or
`fn f(b: &mut Builder, x: Q) -> Q`). Rejected: it makes the code that describes a model read like
code that manipulates a graph, and an operator cannot take a builder at all.
**Precision fixed by the backend or the model** (the archive's precision policy). A backend-wide
policy cannot serve a double-precision solver and a single-precision sampler on one model, and a
model that stores device tensors fixes precision before any estimator has been chosen. Precision
was moved to the query.
**Factorizations on the device.** Rejected while burn has no QR and normal equations square the
condition number; dense solves stay on the host in double precision.
**Target and proposal flags on terms** (the archive). Replaced by samplers that carry their own
densities.
**Many small crates** (the archive). Replaced by two, since crate boundaries kept useful internals
apart without making the design any clearer.