pca-web

Client-side PCA and IncrementalPCA for the browser and Node.js — a zero-runtime-dependency TypeScript reimplementation of scikit-learn 1.9's sklearn.decomposition.PCA and IncrementalPCA, with optional WebGPU acceleration.

• scikit-learn is the specification. Every solver (full, covariance_eigh, arpack, randomized), every nComponents mode (int, fraction, 'mle', null), whitening, copy semantics, the auto-solver heuristic, and the deterministic sign convention (svd_flip(u_based_decision=False)) are ported from the sklearn 1.9.0 source and verified against fixtures generated by the real scikit-learn (see Numerical parity).
• Ships as ESM with TypeScript declarations. Node ≥ 20 and evergreen browsers.
• Zero runtime dependencies. The WebGPU path lives behind the pca-web/webgpu subpath so the core stays GPU-free for bundlers.
• Even the random number generator is a bit-exact numpy RandomState (MT19937) replica: with the same seed, svdSolver: 'randomized' reproduces sklearn's output to floating-point accuracy — not just statistically.
• Fits can run off the main thread (pca-web/client + pca-web/worker), report live progress with intermediate models, be aborted mid-fit, and be serialized for IndexedDB/postMessage/JSON — see the live demo (source in examples/demo/).

import { PCA } from 'pca-web';

const pca = new PCA({ nComponents: 2, svdSolver: 'randomized', randomState: 42 });
const embedding = pca.fitTransform(X); // X: number[][] or a Matrix

pca.components;               // principal axes (Matrix, k×p)
pca.explainedVarianceRatio;   // Float64Array
pca.singularValues;
pca.transform(Xnew);
pca.inverseTransform(embedding);

Installation

npm install pca-web

import { PCA, IncrementalPCA, Matrix } from 'pca-web';          // CPU core
import { WebGPUPCA, isWebGPUSupported } from 'pca-web/webgpu';  // GPU frontend

Inputs, outputs and dtypes

fit/transform accept either number[][] or a Matrix — a thin row-major wrapper around a typed array:

import { Matrix } from 'pca-web';
const X = new Matrix(new Float64Array(n * p), n, p); // zero-copy
const X32 = new Matrix(new Float32Array(n * p), n, p);

The result dtype follows the input, like sklearn: Float32Array in → float32 attributes and transforms out; Float64Array or number[][] in → float64 out. One deliberate divergence: float32 inputs are computed in float64 internally and rounded to float32 on output, whereas sklearn runs a native-f32 LAPACK pipeline — this port is therefore slightly more accurate for f32 data, and agrees with sklearn within the documented f32 tolerance class (~2e-3) rather than bit-for-bit. Transform results are Matrix (use .toArray() for number[][], or .data/.rows/.cols directly).

Like sklearn, copy: false lets fit overwrite the training data: the full solver centers it in place, and the truncated solvers (arpack, randomized) center and square it in place — the input is destroyed. The default copy: true never mutates inputs.

Fitting requires at least 2 samples (n − 1 is the variance denominator). This is a deliberate, documented deviation from sklearn, which accepts a single sample and returns an all-NaN model alongside a Python RuntimeWarning — a channel that does not exist in JS. The same applies to IncrementalPCA's first batch; later partialFit batches may be any size, including single rows.

API ↔ scikit-learn mapping

Names are idiomatic TypeScript camelCase; each maps 1:1 to its sklearn counterpart.

| pca-web | scikit-learn | | --- | --- | | new PCA({ nComponents, copy, whiten, svdSolver, tol, iteratedPower, nOversamples, powerIterationNormalizer, randomState }) | PCA(n_components, copy, whiten, svd_solver, tol, iterated_power, n_oversamples, power_iteration_normalizer, random_state) | | fit(X) / fitTransform(X) / transform(X) / inverseTransform(Y) | fit / fit_transform / transform / inverse_transform | | getCovariance() / getPrecision() | get_covariance() / get_precision() | | scoreSamples(X) / score(X) | score_samples / score | | getFeatureNamesOut() | get_feature_names_out() | | fitAsync(X, opts?) / fitTransformAsync(X, opts?) | — (non-blocking fit, see below) | | toModel() / PCA.fromModel(m) / modelToJSON / modelFromJSON | — (serialization, see below) | | components | components_ | | explainedVariance / explainedVarianceRatio | explained_variance_ / explained_variance_ratio_ | | singularValues / mean / noiseVariance | singular_values_ / mean_ / noise_variance_ | | nComponents / nSamples / nFeaturesIn | n_components_ / n_samples_ / n_features_in_ | | resolvedSvdSolver | _fit_svd_solver (the solver 'auto' dispatched to) | | new IncrementalPCA({ nComponents, whiten, copy, batchSize }) | IncrementalPCA(n_components, whiten, copy, batch_size) | | partialFit(X) / batchSize / nSamplesSeen / variance | partial_fit / batch_size_ / n_samples_seen_ / var_ |

nComponents semantics match sklearn exactly: an integer ≥ 0; a fraction in (0, 1) selecting the smallest k whose cumulative explained-variance ratio exceeds it (full/covariance_eigh); the string 'mle' for Minka's MLE (full/covariance_eigh, requires n ≥ p); or null (default) for min(nSamples, nFeatures) (arpack: min − 1). Validation errors mirror sklearn's messages and check order.

tol is accepted for API parity; the built-in Lanczos solver always converges the requested triplets to machine precision (equivalent to sklearn's default tol=0).

Solver selection

| solver | use when | notes | | --- | --- | --- | | 'auto' (default) | in doubt | sklearn 1.9's exact heuristic: covariance_eigh when p ≤ 1000 and n ≥ 10·p; else full when max(n, p) ≤ 500 or nComponents is 'mle'; else randomized when 1 ≤ k < 0.8·min(n, p); else full | | 'full' | small/medium data, all components | Golub–Reinsch SVD of the centered data | | 'covariance_eigh' | many samples, few features (n ≫ p) | eigendecomposition of the p×p covariance; fastest for tall data, GPU-accelerated | | 'randomized' | large data, few components | Halko et al. randomized SVD, seed-compatible with sklearn; GPU-accelerated | | 'arpack' | k strictly < min(n, p), sparse-ish spectra | Golub–Kahan–Lanczos with full reorthogonalization (scipy svds equivalent) |

RandomState (exported) is the numpy-compatible RNG; pass an integer seed (< 2³²) or an instance as randomState.

IncrementalPCA

import { IncrementalPCA } from 'pca-web';

const ipca = new IncrementalPCA({ nComponents: 8, batchSize: 200 });
ipca.fit(X);                    // internally batched over gen_batches(n, batchSize)
// or stream without holding the dataset in memory:
for (const chunk of chunks) {
  ipca.partialFit(chunk);
}
ipca.transform(Xnew);

Semantics ported from sklearn: batchSize defaults to 5 * nFeatures in fit; small tail batches are merged like gen_batches; the incremental mean/variance uses float64 accumulators (so mean/variance are always float64, while float32 components stay float32 after the first partialFit and become float64 once batches are stacked — numpy's vstack promotion, faithfully replicated); noiseVariance follows sklearn's k ∉ {batch n, p} rule. The parity suite verifies the complete fitted state after every partial_fit step against sklearn.

WebGPU acceleration

import { WebGPUPCA, isWebGPUSupported } from 'pca-web/webgpu';

const pca = new WebGPUPCA({ nComponents: 16, svdSolver: 'randomized', randomState: 42 });
await pca.fit(X);          // async: GPU when available and worthwhile, CPU otherwise
pca.backend;               // 'webgpu' | 'cpu' — what the last fit actually used
pca.gpuAdapterInfo;        // e.g. "apple / metal-3"
await pca.transform(Xnew); // GPU projection for large inputs
pca.dispose();             // release GPU resources

• What runs on the GPU: the large-GEMM hotspots — the Gram product of covariance_eigh, the power-iteration panel products of randomized (the data matrix stays resident on the device across iterations), and large transform projections. Everything else (eigh, panel QR/LU, small SVDs, sign flips, validation, bookkeeping) runs through the same code as the CPU class, so semantics cannot drift.
• Fallback: without WebGPU (or for full/arpack, or inputs below minGpuElements, default 2¹⁸ elements), WebGPUPCA delegates to the CPU PCA — results are bit-identical to new PCA(...) by construction (verified in the Node test suite).
• Device injection: pass your own device: GPUDevice (it is never destroyed by dispose()); powerPreference is forwarded to requestAdapter.
• TypeScript: the pca-web/webgpu declarations use the standard WebGPU globals (GPUDevice, …). With lib: ["dom"] (TypeScript ≥ 5.9 ships WebGPU in the DOM lib) they resolve with no extra packages. Node-only lib targets instead add the types-only package @webgpu/types (declared here as an optional peer dependency) with "types": ["@webgpu/types"] and skipLibCheck in tsconfig. The runtime has zero dependencies either way, and the core pca-web entry needs nothing.
• Precision: WGSL has no f64, and (measured on Chrome/Metal) the shader compiler's fast-math destroys classic compensated-arithmetic tricks. The GEMM kernels therefore use fma-exact products of double-single operands with exact integer-binned accumulation (5×13-bit i32 bins per element; integer adds cannot be fast-mathed), recombined in float64. The engine self-checks precision on creation and refuses the GPU (falling back to CPU) if the driver cannot reach near-f64 accuracy — measured 5.3e-13 worst-case GEMM relative error on Apple Metal, vs ~1e-5 for a plain f32 kernel.
• Measured CPU↔GPU equivalence (browser harness, Apple M-series, apple / metal-3): components/singular values/variances agree to ~1e-10 absolute / ~1e-8 relative worst-case on float64 fits (dominated by eigenvector gap sensitivity, not GEMM error); float32 fits round to identical f32 outputs.

Benchmarks from this machine (npm run bench:browser, median of 3 fits, GPU time includes upload + readback):

| case | size | solver | CPU (ms) | GPU (ms) | speedup | | --- | --- | --- | ---: | ---: | ---: | | cov_50000x200 | 50000×200 | covariance_eigh | 730 | 237 | 3.1× | | cov_100000x100 | 100000×100 | covariance_eigh | 396 | 149 | 2.7× | | rand_20000x500_nc16 | 20000×500 | randomized | 3332 | 637 | 5.2× | | rand_5000x2000_nc32 | 5000×2000 | randomized | 4843 | 631 | 7.7× | | cov_f32_50000x200 (f32) | 50000×200 | covariance_eigh | 1250 | 235 | 5.3× |

Progress, non-blocking fits and cancellation

Every estimator has four fit entry points with identical numerics — fit and fitTransform run the solver to completion synchronously; fitAsync and fitTransformAsync run the same solver steps time-sliced through the event loop (default budget 12 ms per slice), so the UI keeps painting. Observer-on vs observer-off, sync vs async: the fitted models are bit-identical (enforced by tests, including that snapshots draw nothing from the RNG stream).

const pca = new PCA({ nComponents: 16, svdSolver: 'randomized', randomState: 42 });
const controller = new AbortController();

await pca.fitAsync(X, {
  budgetMs: 12,                      // solver time per event-loop slice
  signal: controller.signal,         // abort mid-fit → rejects with name === 'AbortError'
  snapshot: { scores: true, every: 1 },
  onProgress(e) {
    render(e.fraction);              // [0,1] monotone, or null while indeterminate
    if (e.snapshot) plot(e.snapshot.scores);   // evolving embedding of the training rows
  },
});

onProgress fires synchronously from inside the fit with {estimator, solver, phase, step, totalSteps, fraction, snapshot?, detail?}. Cadence and snapshot availability per solver:

| solver | phases | one event per… | fraction | intermediate snapshots | | --- | --- | --- | --- | --- | | full | decompose → finalize | decomposition (indeterminate) | null → 1 | finalize only | | covariance_eigh | gram → decompose → finalize | Gram chunk (~2²² elements) | 0–0.85 → pinned → 1 | finalize only | | randomized | power-iteration → finalize | power iteration | 0–0.9 → 1 | every iteration (opt-in) | | arpack | lanczos-step → finalize | convergence check; detail: {basisSize, jmax, maxResidual} | null → 1 | at checkpoints | | IncrementalPCA | batch → finalize | partial_fit batch | linear → 1 | every batch |

Snapshots are opt-in (snapshot: { components: true } or { scores: true }): fresh float64 copies in sklearn's sign convention, safe to keep or transfer. On the randomized solver a snapshot costs roughly one extra solver step per iteration — every: 2 halves that.

Semantics: an onProgress exception propagates and fails the fit. An aborted or failed PCA fit leaves the estimator unfitted (an abort before the first step leaves a previous model intact); IncrementalPCA keeps the completed-batch model, matching partial_fit semantics. Abort accepts any {aborted, reason?} object — a real AbortSignal works, none is required. Concurrent fits on one instance throw. WebGPUPCA.fit(X, options) takes the same options; abort during a GPU fit releases device buffers and does not fall back to a CPU refit.

Web Workers

pca-web/client exports WorkerPCA and WorkerIncrementalPCA — async proxies that run the estimator in a worker and stream progress back. The client bundle contains no solver code (enforced by an import-graph test): apps that fit only in the worker ship the numerics once, inside the worker bundle.

import { WorkerPCA } from 'pca-web/client';

const pca = new WorkerPCA({
  nComponents: 16,
  svdSolver: 'randomized',
  randomState: 42,            // numeric seeds only — a live RandomState can't cross threads
  backend: 'webgpu',          // optional: WebGPU inside the worker, CPU fallback
});
await pca.fit(X, {
  signal: controller.signal,  // aborts land mid-fit (the worker fit is time-sliced)
  transfer: true,             // optional: move X's buffer instead of copying it
  onProgress: (e) => render(e.fraction, e.snapshot),
  progress: { minIntervalMs: 33, snapshot: { scores: true } },  // worker-side throttle
});

pca.components;               // sync — every fit piggybacks the model back to a client mirror
pca.explainedVarianceRatio;   // sync
await pca.transform(Xnew);    // compute methods stay in the worker
const model = pca.exportModel();     // sync snapshot of the mirror (see Model serialization)
await pca.dispose();          // frees the estimator; terminates the worker it owns

Progress events are throttled worker-side (minIntervalMs, default 33 ms, latest-wins within a phase; phase boundaries and finalize always delivered, seq monotone). transfer: true detaches the input's buffer from the caller; views are tight-sliced first so a subarray's parent buffer is never touched. With backend: 'webgpu', await pca.info() reports which backend actually executed plus the adapter ({backend, gpuAdapterInfo, webgpuAvailable}).

Where the worker comes from. By default WorkerPCA spawns the packaged entry via new Worker(new URL('./worker.js', import.meta.url), { type: 'module' }) — webpack 5, Vite build, and native ESM all resolve that statically. If your toolchain can't, pass { worker }:

| toolchain | recipe | | --- | --- | | webpack 5, native ESM | default works — or { worker: new Worker(new URL('pca-web/worker', import.meta.url), { type: 'module' }) } | | Vite | make a one-line entry worker-entry.ts containing import 'pca-web/worker', pass { worker: () => new Worker(new URL('./worker-entry.ts', import.meta.url), { type: 'module' }) }, and set worker: { format: 'es' } in vite.config.ts (what examples/demo does); or import PcaWorker from 'pca-web/worker?worker' → { worker: () => new PcaWorker() }. In dev also add 'pca-web' to optimizeDeps.exclude. | | esbuild | bundle pca-web/worker as its own entry point and pass its URL: { worker: () => new Worker('/assets/pca-worker.js', { type: 'module' }) } | | Node ≥ 20 | worker_threads — bridge a MessagePort (below) |

Node (worker_threads Worker is an EventEmitter, not an EventTarget, so hand the client a MessagePort, which is one):

// pca.worker.mjs — runs in the thread
import { workerData } from 'node:worker_threads';
import { attachPCAWorker } from 'pca-web/worker';
attachPCAWorker(workerData.port);

// main thread
import { MessageChannel, Worker } from 'node:worker_threads';
import { WorkerPCA } from 'pca-web/client';
const { port1, port2 } = new MessageChannel();
new Worker(new URL('./pca.worker.mjs', import.meta.url), {
  workerData: { port: port2 },
  transferList: [port2],
});
const pca = new WorkerPCA({ nComponents: 3, worker: port1 });

attachPCAWorker(port) (from pca-web/worker) attaches the request handler to any port-shaped endpoint; the packaged worker entry just calls it on the worker's global scope. Requests on one worker execute strictly in order; aborts are handled out-of-band. Errors are marshalled back instanceof-correct where it matters (NotFittedError, AbortError), and terminate() rejects in-flight calls with WorkerTerminatedError.

Model serialization

toModel() captures the complete fitted state as a plain object of typed arrays — cheap, structured-clone-friendly (IndexedDB, postMessage, structuredClone), and validated on rehydration. PCA.fromModel(model) / IncrementalPCA.fromModel(model) restore an estimator whose every method output is bit-identical to the original's; a rehydrated IncrementalPCA continues partialFit streams exactly as if never interrupted.

// IndexedDB (structured clone — no JSON, typed arrays stored directly)
store.put(pca.toModel(), 'my-model');
const pca2 = PCA.fromModel(await get(store, 'my-model'));

// JSON when you need text (bit-exact round-trip for f64 and f32)
localStorage.setItem('m', modelToJSON(pca.toModel()));
const pca3 = PCA.fromModel(modelFromJSON(localStorage.getItem('m')!) as PCAModel);

The worker proxies speak the same format: exportModel() (synchronous, from the client-side mirror) and WorkerPCA.fromModel(model, options), which hydrates the worker without a refit. Models carry a formatVersion; assertValidModel rejects malformed or future-versioned input.

Numerical parity and tolerances

Parity is defined by tests, not intent. python/generate_fixtures.py runs real scikit-learn 1.9.0 (numpy 2.5.1) over 96 PCA/IncrementalPCA cases — every solver, every nComponents mode, whiten on/off, float32/float64, rank-deficient and badly scaled data, per-step partial_fit state — plus RNG-stream and LAPACK-kernel fixtures, and commits the outputs (fixtures/, ~11 MB). The TypeScript suite fits the same inputs and compares every fitted attribute and method output elementwise — including signs; comparisons are never up-to-sign.

Assertion bounds by class (all cases pass; see tests/pca-parity.test.ts for the full table):

| class | components | explained variance (rel) | transforms | notes | | --- | --- | --- | --- | --- | | float64 full/arpack | 1e-9 abs / 1e-7 rel | 1e-9 | 1e-8 abs / 1e-7 rel | observed component diffs ≤ ~1e-12 on unit-scale data | | float64 covariance_eigh | 1e-9 abs / 1e-7 rel | 1e-9 | 1e-7 abs / 1e-6 rel | Gram route squares the condition number; tail singular values carry a ~√eps noise floor | | float64 randomized (LU/QR/auto) | 1e-7 abs / 1e-5 rel | 1e-7 | 1e-6 abs / 1e-5 rel | same-seed comparison; observed ≤ ~1e-10 | | float32 (all solvers) | 2e-3 | 2e-3 | 5e-3 | f32 rounding dominates | | IncrementalPCA float64 | 1e-8 abs / 1e-6 rel | 1e-8 | 1e-7 | observed ≤ ~6e-13; mean/var bit-exact |

Documented caveats, discovered and verified during porting:

• noiseVariance is the cancellation (totalVar − Σev)/(min − k); its relative error is the explained-variance error amplified by ~totalVar/noiseVar. Bounds: 1e-11 abs (roundoff- scale zeros) / 1e-9 rel deterministic, 1e-4 rel randomized.
• powerIterationNormalizer: 'none' (explicitly requested, ≥3 iterations) is documented by sklearn itself as numerically unstable: trailing directions carry (σᵢ/σ₁)^(2·nIter+1) ≈ 1e-13 relative weight, and sklearn's own components move ~2e-4 under a 1-ulp input perturbation (measured). The suite compares that configuration in a dedicated loose class (components 5e-3 abs); no implementation can do better than the algorithm's conditioning.
• Rank deficiency: components beyond the numerical rank span an arbitrary orthonormal basis of the null space — LAPACK's choice and this library's legitimately differ. The suite masks comparisons beyond the fixture rank (their variances are still compared, near zero), and skips inverseTransform comparisons of out-of-span test data when rank < k.
• Score/log-likelihood values are large negatives; they agree to ~1e-10 relative (deterministic f64).
• Deliberate deviations: fits with fewer than 2 samples throw a clear error (sklearn returns an all-NaN model with a RuntimeWarning). All-constant (zero-variance) data keeps sklearn's exact behavior: a well-formed model whose explainedVarianceRatio is NaN (0/0).

Demo app

examples/demo/ is a Vite app exercising the whole feature surface — worker vs main-thread execution, live progress with snapshot embeddings, abort, WebGPU-in-worker with adapter readout, IncrementalPCA partialFit streaming, eigen-digit/reconstruction/outlier panels, and model persistence (JSON + IndexedDB restore-without-refit). It consumes the library strictly through its public package exports, doubling as a packaging test.

npm run demo:dev       # build the library, then vite dev
npm run demo:build     # static build (examples/demo/dist)

The deploy workflow (.github/workflows/deploy-pages.yml) publishes it to GitHub Pages on every push to main; one-time repo setup: Settings → Pages → Source: "GitHub Actions". The bundled digits dataset is regenerated by npm run demo:data (same venv as fixtures).

Testing and development

npm test                     # Node suite: parity vs sklearn fixtures + API behavior,
                             # async/observer bit-equivalence, abort, model round-trips,
                             # worker protocol/parity over MessageChannel (275 tests)
npm run test:coverage        # the same suite with a V8 coverage report (report-only)
npm run test:browser         # real-GPU harness (Playwright Chromium); reports exactly which
                             # cases executed on which adapter, and fails if the GPU didn't run
npm run test:browser:worker  # real Worker in Chromium: bit-equal results, progress on the
                             # main thread, mid-fit abort, transfer semantics, WebGPU-in-worker
npm run typecheck && npm run lint && npm run build
npm run bench                # CPU timings (Node)
npm run bench:browser        # CPU vs GPU timings (browser)

Regenerating fixtures requires Python with the pinned reference stack (scikit-learn 1.9.0, numpy 2.5.1, scipy 1.18.0 — python/requirements.txt):

python -m venv .venv && .venv/bin/pip install -r python/requirements.txt
npm run fixtures       # rewrites fixtures/ deterministically

The numerics core (src/numeric/) — Golub–Reinsch SVD, symmetric tridiagonal eigensolver, Householder QR with LAPACK dlarfg semantics, partial-pivot LU (scipy permute_l compatible), Golub–Kahan–Lanczos, Minka MLE, and the MT19937 RandomState replica — is dependency-free TypeScript, unit-tested against numpy/scipy fixtures.

License

MIT