# MOSAIC Online

Andrew Kiruluta, UC Berkeley, CA. Aug 2026

Reference implementation of the refactored **MOSAIC** algorithm from the manuscript
*“MOSAIC: Closed-Form Riemannian Streaming Updates on a Low-Rank Covariance Bundle.”*

MOSAIC maintains a fixed-memory streaming state consisting of

- an orthonormal frame `U ∈ St(d,r)` representing the retained subspace;
- an SPD covariance fiber `S ∈ SPD(r)` and its inverse;
- a running mean `mu ∈ R^d`;
- a residual-noise estimate `sigma²`.

The geometric state is interpreted through the gauge relation

```text
(U, S) ~ (U Q, Q^T S Q),     Q ∈ O(r),
```

rather than as a globally independent Cartesian product. The implementation follows the corrected novelty framing of the manuscript: the Grassmann geometry and AIRM geometry are established tools; the algorithmic contribution is their samplewise specialization into a fixed-memory streaming recursion with an exact rank-one AIRM likelihood update and a stable precision-preconditioned subspace rule.

## 1. Algorithm implemented

For a centered sample `z = x - mu`, MOSAIC computes

```text
c   = U^T z
rho = z - U c
e   = ||rho||
```

and updates the residual noise floor by an EMA of `e²/(d-r)`.

### 1.1 Precision-preconditioned Grassmann step

The coordinate direction is

```text
xi = (S + sigma² I)^(-1) c.
```

This is intentionally described as a **noise-regularized precision preconditioner**, not as the exact Fisher natural gradient. After normalizing

```text
v = xi / ||xi||,        w = rho / ||rho||,
```

the rank-one Grassmann exponential rotates only the plane spanned by `U v` and `w`:

```text
theta = eta_u * atan2(||rho||, c^T v)

U+ = U + [ (cos(theta)-1) Uv + sin(theta) w ] v^T.
```

The update costs `O(d r)` and preserves orthonormality analytically.

### 1.2 Exact rank-one AIRM covariance step

For the one-sample retained-coordinate likelihood

```text
f_c(S) = 1/2 [ log det S + c^T S^(-1) c ],
```

the affine-invariant Riemannian gradient is

```text
grad f_c(S) = 1/2 (S - c c^T).
```

The exact exponential map simplifies to

```text
q = c^T S^(-1) c

g(q) = [exp(eta_s q / 2) - 1] / q

S+ = exp(-eta_s/2) [ S + g(q) c c^T ].
```

No matrix square root or generic matrix exponential is needed. `S^-1` is updated by Sherman--Morrison, so this part costs `O(r²)`.

`q` is simultaneously a retained-space Mahalanobis novelty score. The coefficient of `c c^T` is surprise-adaptive: low-`q` observations can be weighted more conservatively than the Euclidean first-order update, whereas sufficiently high-`q` observations receive more weight.

## 2. Repository layout

```text
mosaic-online-refactored/
├── README.md
├── LICENSE
├── pyproject.toml
├── requirements.txt
├── src/mosaic_online/
│   ├── __init__.py
│   ├── geometry.py        # exact Grassmann and AIRM primitives
│   ├── model.py           # streaming MOSAIC estimator
│   └── metrics.py
├── tests/
│   ├── test_geometry.py   # equation/invariant checks
│   └── test_model.py
├── scripts/
│   ├── smoke_digits.py    # standard benchmark smoke test
│   └── smoke_synthetic.py # controlled subspace recovery smoke test
├── examples/
│   └── minimal_stream.py
└── results/
    └── digits_smoke.json  # generated reference run
```

## 3. Installation

Python 3.10+ is recommended.

```bash
python -m venv .venv
source .venv/bin/activate
python -m pip install --upgrade pip
pip install -e '.[test]'
```

For only the core algorithm, NumPy is the sole runtime dependency:

```bash
pip install -e .
```

The benchmark additionally needs scikit-learn:

```bash
pip install -e '.[benchmark]'
```

## 4. Minimal use

```python
import numpy as np
from mosaic_online import MOSAIC

model = MOSAIC(
    n_features=64,
    rank=8,
    eta_u=0.08,
    eta_s=0.03,
    eta_mu=0.01,
    seed=0,
)

for x in stream:
    model.partial_fit(x)

coords = model.transform_one(x)
score = model.novelty_score(x)
reconstruction = model.reconstruct(x)
cov = model.covariance(include_residual=True)
```

The retained state does **not** grow with the number of observations. Its main floating-point storage is `d*r + r*r + d + O(1)`.

## 5. Standard benchmark smoke test: handwritten Digits

The default end-to-end smoke test uses `sklearn.datasets.load_digits`, a standard 1,797-example handwritten-digit benchmark with 64 pixel features. It is bundled with scikit-learn, so **no network download is required**.

The protocol is deliberately simple and reproducible:

1. scale pixels from `[0,16]` to `[0,1]`;
2. stratified 75/25 train/test split with seed 42;
3. shuffle the training split once;
4. stream the 1,347 training observations through MOSAIC **once**;
5. evaluate rank-16 reconstruction on the 450 held-out examples;
6. verify that the Grassmann frame remains orthonormal, the covariance fiber remains SPD, the stored inverse remains consistent, and all state values stay finite;
7. report batch PCA at rank 16 only as a familiar non-streaming reference ceiling—not as a like-for-like competitor.

Run:

```bash
python scripts/smoke_digits.py
```

or save machine-readable output:

```bash
python scripts/smoke_digits.py --json results/digits_smoke.json
```

The smoke test exits non-zero if any invariant fails. Its pass criteria are intentionally platform-tolerant:

```text
all state finite                    True
||U^T U - I||_F                     < 1e-8
lambda_min(S)                       > 1e-12
||S S^-1 - I||_F                   < 1e-5
held-out centered distortion        < 0.60
```

The reconstruction threshold is a **sanity criterion**, not a state-of-the-art claim. Batch PCA should normally reconstruct better because it sees the full training matrix jointly and is used here only to contextualize the rank-16 compression level.

## 6. Geometry/unit smoke tests

Run the complete test suite:

```bash
pytest
```

The tests cover:

- exact orthonormality of the rank-one Grassmann exponential;
- numerical agreement between the closed-form AIRM update and the generic `sqrtm/expm` affine-invariant exponential;
- positive-definiteness under large positive AIRM step sizes;
- Sherman--Morrison inverse consistency;
- finite novelty scores and reconstruction shapes;
- fixed-memory state accounting.

A controlled synthetic stream also checks whether the learned subspace improves relative to its random initialization:

```bash
python scripts/smoke_synthetic.py
```

## 7. Hyperparameters

### `eta_u`
Controls the Grassmann rotation. Start small on stationary data and raise it when the subspace must track drift.

Typical exploratory range:

```text
0.02, 0.05, 0.1, 0.2, 0.4
```

### `eta_s`
Controls the exact AIRM covariance step. Because the exponential update remains SPD for every positive mathematical step size, it does not have the same positivity boundary as the Euclidean truncation. Very large values can still create poor numerical/statistical behavior, so practical tuning is still required.

Typical exploratory range:

```text
0.005, 0.02, 0.05, 0.1
```

### `eta_mu`
EMA rate for the running location. Set `0` if inputs are already centered. For stationary uncentered streams, values around `0.001--0.02` are a reasonable starting range.

### `noise_beta`
EMA rate for the residual variance estimate. The manuscript uses `0.005`.

### `rank`
Choose according to the compression budget and reconstruction objective. Because the current package intentionally implements the **fixed-rank core algorithm**, automatic rank growth/shrinkage from exploratory manuscript code is not enabled in this reference release. That keeps the central geometric mechanism easy to test and avoids conflating the novelty claim with a heuristic rank controller.

## 8. Numerical safeguards

The exact AIRM expression contains `exp(eta_s*q/2)`. Extremely surprising samples can overflow ordinary floating point even though the underlying symbolic formula is valid. The implementation therefore clips the exponent by default with `exp_cap=40`.

This is a numerical guard, not an alternative geometry. For equation-verification tests it is disabled with `exp_cap=None`.

The rank-one Grassmann exponential preserves orthogonality analytically. Optional periodic QR cleanup is available through `reorthogonalize_every`, but it is off by default. If QR is enabled, the corresponding gauge transformation is applied to `S` so the ambient covariance representation remains coherent.

## 9. What this repository does and does not claim

This package implements the **refactored** manuscript claims:

- Grassmann rank-one geodesics are established prior art.
- Affine-invariant SPD geometry is established prior art.
- Fixed-rank covariance quotient/bundle geometry is established prior art.
- The implemented subspace direction is a stable precision preconditioner, **not** asserted to be the exact Fisher natural gradient.
- The algorithmic centerpiece is the exact square-root-free rank-one AIRM likelihood recursion combined with a frame-coherent streaming subspace update.

The Digits run is a smoke benchmark for reproducibility and numerical health. It is not presented as evidence of state-of-the-art recognition accuracy or as a comprehensive empirical comparison. A publication-grade evaluation should add multiple real streaming datasets, tuned online baselines, concept drift, anisotropy sweeps, repeated seeds, confidence intervals, and wall-clock/memory profiling.

## 10. Reproducing the packaged reference result

From the repository root:

```bash
python -m pip install -e '.[test]'
pytest
python scripts/smoke_synthetic.py
python scripts/smoke_digits.py --json results/digits_smoke.json
```

All three commands should exit with status 0.

## 11. Citation / manuscript traceability

The source code is organized so the principal equations map directly to the manuscript:

- `grassmann_exp_rank1`: rank-one Grassmann exponential;
- `MOSAIC.partial_fit`: precision-preconditioned subspace update and frame-coherent coordinate refresh;
- `spd_geodesic_step`: exact rank-one AIRM fiber exponential;
- `novelty_score`: `q = c^T S^-1 c`.

For scientific use, cite the accompanying MOSAIC manuscript and clearly distinguish the corrected bundle/precision-preconditioning formulation from earlier draft terminology.

Manuscript can be found here: [10.13140/RG.2.2.17951.52646](https://doi.org/10.13140/RG.2.2.17951.52646/1)





## 12. Packaged smoke-test result

The repository was preflighted before packaging with Python/NumPy/SciPy/scikit-learn in the build environment. The packaged seed-42 Digits result in `results/digits_smoke.json` is:

```text
train/test samples                  1347 / 450
ambient dimension / retained rank  64 / 16
MOSAIC centered distortion         0.16156
batch PCA reference distortion     0.15405
||U^T U - I||_F                    1.65e-15
lambda_min(S)                      5.29e-02
||S S^-1 - I||_F                   7.96e-15
MOSAIC state                       1,345 floats
all smoke checks                   PASS
```

The synthetic subspace test reduced normalized projector error from approximately `0.940` at random initialization to `0.0147` after one pass through 2,500 samples.

A five-seed Digits preflight (`0, 1, 7, 42, 99`) produced held-out MOSAIC distortions between approximately `0.159` and `0.166`; every run passed the numerical invariants. Exact throughput depends strongly on CPU/BLAS and should not be compared across machines without a controlled setup.