vitali

Multiscale geometric analysis of point clouds, in Rust, with Python bindings and a MATLAB bridge.

Two constructions live here, both from Mauro Maggioni and coauthors, and neither with a maintained implementation in any package ecosystem before this crate.

Geometric multi-resolution analysis builds a multiscale, piecewise-affine approximation of a point cloud: a tree partition of the data, a local PCA scaling basis and centre at every node, and a geometric wavelet basis relating each node to its parent. Allard, Chen and Maggioni, Multi-scale geometric methods for data sets II: geometric multi-resolution analysis, Applied and Computational Harmonic Analysis, 2012; Chen and Maggioni, Multiscale geometric wavelets for the analysis of point clouds, CISS, 2010.

Diffusion wavelets build an orthonormal multiresolution directly from a diffusion operator, rather than from translates and dilates of a fixed mother wavelet: rank-revealing Gram-Schmidt on the operator's columns gives a scaling basis, the operator is compressed onto it, the compressed operator is squared to represent the next dyadic power, and the process repeats. Coifman and Maggioni, Diffusion wavelets, Applied and Computational Harmonic Analysis, 2006.

The crate is named for the Vitali covering lemma. The greedy net that builds every scale of the GMRA partition is a Vitali covering argument.

Prior art in August 2026

Geometric multi-resolution analysis has three GitHub repositories, none above one star, one a 2015 Python toolkit and one an R package last pushed in February 2021. Nothing on PyPI, in the Julia General registry, or on crates.io.

Diffusion wavelets have nothing in any ecosystem. Wavelets.jl, ContinuousWavelets.jl and PyWavelets are classical wavelets on regular grids. The construction that builds an orthonormal multiresolution from a diffusion operator exists in the 2007 MATLAB sources and nowhere else.

Install

cargo add vitali     # Rust
pip install vitali   # Python

Rust

use vitali::gmra;

let points = gmra::datasets::d_sphere(4000, 4, 100, 0);
let tree = gmra::build(&points, &gmra::Options::default());

println!("{} nodes, depth {}", tree.nodes.len(), tree.depth());

let scales = gmra::reconstruct_all_scales(&tree, &points);
for (j, approx) in scales.iter().enumerate() {
    let err = (approx - &points).norm() / points.norm();
    println!("scale {}: relative error {:.3e}", j + 1, err);
}

Diffusion wavelets take a sparse operator:

use vitali::diffusion_wavelets::{Options, build_tree, circle_diffusion};

let t = circle_diffusion(256);
let tree = build_tree(&t, &Options::default());

for level in &tree.scaling {
    println!("level {}: dimension {}", level.level, level.ext_basis.ncols());
}

Point clouds are D x N in Rust, one column per point, following the reference MATLAB distribution.

Python

Point clouds are (n_points, n_features) in Python, one row per point, following scikit-learn. The boundary transposes.

import numpy as np
import vitali

X = vitali.datasets.d_sphere(4000, d=4, n_features=100, seed=0)
tree = vitali.gmra(X, min_leaf_size=50)

tree.median_scaling_dimension()   # 4.0, the manifold's true dimension
tree.depth
tree.num_nodes

scales = tree.reconstruct_all_scales(X)
for j, approx in enumerate(scales, start=1):
    err = np.linalg.norm(approx - X) / np.linalg.norm(X)
    print(f"scale {j}: relative error {err:.3e}")

A node exposes its own geometry:

node = tree.node(tree.leaves()[0])
node.centre          # (n_features,)
node.scaling_basis   # (local_dim, n_features), one direction per row
node.sigmas          # every singular value, before truncation
node.wavelet_basis   # the basis relating this node to its parent
node.radius

Diffusion wavelets take a dense array or any scipy sparse matrix, which is converted to CSR without ever being densified:

import scipy.sparse as sp
import vitali

T = sp.csr_matrix(vitali.circle_diffusion(256))
tree = vitali.diffusion_wavelets(T, precision=1e-10)

tree.scaling_dimensions()          # dimension at every level
tree.scaling()[0].ext_basis        # Phi_1 in the original delta basis
tree.wavelet()[1].ext_basis        # Psi_2 in the original delta basis

MATLAB

MATLAB reaches vitali through its Python interface, so the Rust runs underneath without a MEX build. Point pyenv at an interpreter your MATLAB release supports, once per machine:

pyenv(Version = "/path/to/python", ExecutionMode = "OutOfProcess");

Out-of-process execution keeps the extension's libraries away from the ones MATLAB loads itself, and is the sensible default.

Two conversions do all the work, and both are in matlab/ in this repository. tonumpy respects MATLAB's column-major order, which is where a hand-written conversion usually goes wrong:

addpath("matlab");

X = rand(4000, 100);              % n_points by n_features
tree = py.vitali.gmra(tonumpy(X), pyargs("min_leaf_size", int64(50)));

k = double(tree.median_scaling_dimension())
d = double(tree.depth)

scales = cell(tree.reconstruct_all_scales(tonumpy(X)));
approx = fromnumpy(scales{end});   % back to a MATLAB double matrix
err = norm(approx - X, "fro") / norm(X, "fro")

Three points decide whether a call works:

• Keyword arguments go through pyargs. Every option is keyword-only.
• Integers need an explicit cast. MATLAB numeric literals are doubles, and min_leaf_size, max_dim, levels, stop_fcns, manifold_dimension and every seed are integers. Pass int64(50), never 50.
• Orientation is transposed against the MATLAB toolboxes. GMRA.m and DWPTree.m take D x N, one column per point. vitali's Python surface takes (n_points, n_features). Data prepared by GenerateData_and_SetParameters.m passes as tonumpy(X.').

A diffusion operator goes across as a dense matrix:

T = fromnumpy(py.vitali.circle_diffusion(int64(256)));   % or your own operator
tree = py.vitali.diffusion_wavelets(tonumpy(T), pyargs("precision", 1e-10));

levels = cell(tree.scaling());
Phi1 = fromnumpy(levels{1}.ext_basis);

matlab/vitali_demo.m runs the whole sequence end to end and prints the recovered dimension and the per-scale reconstruction error.

Differences from the reference distribution

Partitioning. The reference distribution partitions with METIS-based nested dissection or a cover tree, both reached only through precompiled MEX binaries with no source in the distribution: nesdis.mexa64 and covertree.mexa64 do not load under GNU Octave ('covertree' undefined near line 4, checked directly). This crate recursively splits each node's points at the median of their projection onto the node's top local PCA direction, which is a principal-direction divisive partition tree. Node for node, reconstruction error and dimension counts here are not directly comparable against a run using either reference partitioner. The comparisons below are against independently generated data of known dimension, and against the reference distribution's own reported figures as a target range.

Reconstruction. The reference distribution's inverse transform accumulates wavelet detail along the root-to-leaf path, corrected at every step by a tangential-correction mechanism (addTangentialCorrections, its own default) that keeps the accumulated sum aligned with each node's own direct local PCA reconstruction. Without that correction the accumulated sum drifts, because two nodes' independently computed local PCA planes are not nested in general, unlike diffusion wavelets' scaling spaces: an early build of the wavelet-sum reconstruction here measured 26x worse error at the finest scale than the same node's own direct reconstruction (0.66 against 0.025 relative error) on curved test data, worsening with depth. Reproducing that bookkeeping (ScalBasisChange, ScalWavBasisChange, ScalCenterConsts) is out of scope. reconstruct_all_scales instead reconstructs each scale directly from that scale's own node, centre + Phi' Phi (x - centre), which is the value the tangential correction exists to make the accumulated sum equal. Wavelet bases and translation constants are still constructed and available on every node, for compression-oriented use.

Precision. This crate is f64 throughout. See the D-Cube result below.

op_threshold. DWPTree.m's docstring states that OpThreshold defaults to 0; its shipped code defaults it to Precision when left unset. Options::op_threshold follows the shipped code as an explicit field with no implicit default, so a caller wanting the docstring's behaviour passes 0.0 and a caller wanting the shipped behaviour passes precision itself.

Validation

35 tests. Diffusion wavelets are checked against the four structural identities the construction is stated by, and against MATLAB-exported fixtures; GMRA is checked against manifolds of known intrinsic dimension.

Structural identities of the diffusion wavelets

On the circle at N = 256:

Dimension recovery under GMRA

Median scaling dimension across the whole tree, against each manifold's true dimension, fully adaptive (no fixed manifold_dimension, internal_threshold = 0.5, leaf_threshold = 1e-2, both reference defaults), N = 4000:

Reconstruction error per scale under GMRA

The reference distribution's own double-precision run on D-Cube reconstructs to machine precision at every scale with no rise anywhere. Its single-precision run shows a non-monotone error rising to about 0.04 near the deepest scale. With manifold_dimension fixed at the true dimension (matching the reference distribution's own driver for these synthetic cases) and a leaf size comfortably above the dimension being estimated locally (min_leaf_size = 50 for 6-dimensional data), D-Cube at N = 4000, D = 50 gives:

[3.55e-16, 3.22e-16, 2.93e-16, 2.99e-16, 2.73e-16, 2.45e-16, 2.26e-16, 2.01e-16]

Machine precision at every one of the eight scales, and no rise between consecutive scales above 1e-12, which is the threshold the test asserts. Variation between neighbouring scales sits at the 1e-18 level, which is rounding rather than a trend. That reproduces the reference distribution's double-precision finding independently and through a different partitioner, and supports reading the reported single-precision spike as an artefact of that floating-point path rather than a property of the construction.

A smaller leaf size (min_leaf_size = 20, twenty points for a 6-dimensional local PCA) reproduces a related, smaller artefact: one leaf under-samples and keeps rank 5 rather than 6, leaving every dimension-fixed internal scale at machine precision and the finest scale at 0.37% relative error. Both runs are kept as tests. The effect is small, explained, and specific to under-sampled local PCA.

Curved data reconstructs approximately, with error decreasing by scale and settling above machine precision at the finest one. N = 2000, D = 50, min_leaf_size = 20, manifold_dimension fixed at the true dimension, seed 1:

D-Sphere (d = 4): [0.431, 0.246, 0.213, 0.172, 0.126, 0.074, 0.060, 0.025]
SwissRoll (d = 2): [0.390, 0.240, 0.179, 0.151, 0.102, 0.067, 0.035, 0.004]

cargo test -p vitali --release --test reconstruction -- --nocapture prints all of these, and crates/vitali/tests/ states the assertion behind each one.

Layout

crates/vitali/      the library, published to crates.io
crates/vitali-py/   PyO3 bindings, published to PyPI as vitali
matlab/             tonumpy, fromnumpy and a demo script

Licence

Apache-2.0.