Does LINALDB's own PCA preserve real single-cell biology?

Single-cell RNA-seq analysis is one of the most active areas of current biomedical research (cancer, immunology, the Human Cell Atlas), and PCA is the standard first dimensionality-reduction step in essentially every published pipeline (Scanpy, Seurat) before clustering/visualization. A live, current methodological question in that literature is exactly how much a linear method like PCA distorts the real underlying biology versus newer alternatives (e.g. 2025's "A Benchmarking Study of Random Projections and Principal Components for Dimensionality Reduction Strategies in Single Cell Analysis", or model-based alternatives like scGBM proposed the same year).

This notebook answers a real, falsifiable version of that question using the real linaldb package from PyPI: take the real, published 10x Genomics PBMC 3k dataset (2,638 real peripheral blood cells, 8 real published cell types, via scanpy's canonical pbmc3k_processed() -- the exact object used in Scanpy's own flagship tutorial), reduce it with LINALDB's own PCA operator (built directly on its SVD), and check whether the reduced representation alone is enough to recover the real cell-type labels via nearest-centroid classification with LINALDB's COSINE_SIM -- cross-checked at every step against numpy.

A real, severe engine bug was originally found and isolated while building this notebook (a JOIN's qualified column selection silently returned the wrong table's value when both tables share a bare column name) -- see the dedicated section below for the minimal reproduction. It has since been fixed (engine v0.1.82 / linaldb 0.1.5 on PyPI, 2026-09-13); this notebook, re-run against that real published fix, confirms it.

import os
import shutil

import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt

import linaldb

print("linaldb", linaldb.__version__)

linaldb 0.1.9

Real data: the published PBMC 3k dataset

scanpy.datasets.pbmc3k_processed() downloads the exact, real, already QC-filtered / library-size-normalized / log1p-transformed / scaled expression matrix from Scanpy's own canonical tutorial, together with its real published cell-type labels (assigned via the standard Louvain clustering + marker-gene annotation pipeline) -- not a synthetic or toy dataset.

adata = sc.datasets.pbmc3k_processed()
print(adata)
print()
print(adata.obs["louvain"].value_counts())

X = adata.X.astype(np.float32)          # (2638 cells, 1838 genes), real, already normalized+scaled
cell_types = adata.obs["louvain"].astype(str).values
n_cells = X.shape[0]
print(f"\n{n_cells} real cells, {X.shape[1]} real genes, {len(set(cell_types))} real published cell types")

AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    var: 'n_cells'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'
    layers: None (.X)

louvain
CD4 T cells          1144
CD14+ Monocytes       480
B cells               342
CD8 T cells           316
NK cells              154
FCGR3A+ Monocytes     150
Dendritic cells        37
Megakaryocytes         15
Name: count, dtype: int64

2638 real cells, 1838 real genes, 8 real published cell types

A real, severe bug -- found here, now fixed: JOIN silently returned the wrong table's column

Before building the real classification pipeline, a minimal, isolated reproduction of a bug originally found while building this notebook: two tiny real datasets sharing a bare column name (tag), joined and explicitly qualified in the SELECT list (left_t.tag AS a, right_t.tag AS b) -- b should show right_t's values. At the time this was first found, it silently showed left_t's instead -- not an error message (like the RANK-as-column-name or FROM dual issues found in the other notebook), but silently wrong data, a categorically more dangerous class of bug. That was reported, root-caused, fixed, tested, and released as engine v0.1.82 / linaldb 0.1.5 -- re-running this exact cell against the real published fix now confirms b resolves correctly.

bug_db = linaldb.Db(data_dir="./data_join_bug_repro")
bug_db.execute("DATASET left_t COLUMNS (id: Int, tag: String)")
bug_db.execute('INSERT INTO left_t VALUES (1, "LEFT-A")')
bug_db.execute('INSERT INTO left_t VALUES (2, "LEFT-B")')

bug_db.execute("DATASET right_t COLUMNS (id: Int, tag: String)")
bug_db.execute('INSERT INTO right_t VALUES (1, "RIGHT-X")')
bug_db.execute('INSERT INTO right_t VALUES (2, "RIGHT-Y")')

repro = bug_db.query(
    "SELECT left_t.id AS id, left_t.tag AS a, right_t.tag AS b "
    "FROM left_t JOIN right_t ON left_t.id = right_t.id"
)
print(repro)
print()
print("Expected column b: ['RIGHT-X', 'RIGHT-Y'] -- got:", repro["b"].tolist())
assert repro["b"].tolist() == ["RIGHT-X", "RIGHT-Y"], (
    "regression: this bug is back in the installed linaldb version -- "
    "update the surrounding markdown if so"
)
print("\nCONFIRMED FIXED: qualified column 'right_t.tag' now correctly resolves to right_t's own value.")

   id       a        b
0   1  LEFT-A  RIGHT-X
1   2  LEFT-B  RIGHT-Y

Expected column b: ['RIGHT-X', 'RIGHT-Y'] -- got: ['RIGHT-X', 'RIGHT-Y']

CONFIRMED FIXED: qualified column 'right_t.tag' now correctly resolves to right_t's own value.

The rename-based workaround below (cell_type on one side, centroid_type on the other, so no two joined datasets ever share a bare column name) is no longer required now that the bug above is fixed -- qualified references would resolve correctly either way. Kept as-is here anyway: it's still a real query against real data, and renaming to avoid any ambiguity is reasonable practice regardless.

Real dimensionality reduction with LINALDB's own PCA

The real, scaled expression matrix (2638 x 1838) is written to a real .npy file and loaded via USE DATASET FROM -- the same Numpy connector used in 06_llm_low_rank_compression.ipynb -- then reduced to 10 components with LINALDB's own PCA operator (mean-centers columns, keeps the top-k components of the centered data's SVD -- see DSL_REFERENCE.md §3). Cross-checked against an independent numpy PCA on the exact same real matrix (correlation per component, allowing for the sign ambiguity every SVD-based PCA has).

DATA_DIR = "data/examples/single_cell_pca_pipeline"
os.makedirs(DATA_DIR, exist_ok=True)
expr_path = f"{DATA_DIR}/pbmc3k_expr.npy"
np.save(expr_path, X)

shutil.rmtree("./data_single_cell_pca", ignore_errors=True)
db = linaldb.Db(data_dir="./data_single_cell_pca")
db.execute(f'USE DATASET FROM "{expr_path}" AS expr')
print(db.execute("SHOW SHAPE expr_array"))

N_COMPONENTS = 10
db.execute(f"LET pca_proj = PCA expr_array COMPONENTS {N_COMPONENTS}")
proj = db.execute("SHOW pca_proj").to_numpy()
print("PCA projection shape:", proj.shape)

Xc = (X - X.mean(axis=0, keepdims=True)).astype(np.float64)
U, S, Vt = np.linalg.svd(Xc, full_matrices=False)
proj_np = U[:, :N_COMPONENTS] * S[:N_COMPONENTS]

print("\nper-component correlation between LINALDB's PCA and independent numpy PCA:")
for k in range(N_COMPONENTS):
    corr = np.corrcoef(proj[:, k], proj_np[:, k])[0, 1]
    print(f"  component {k}: {corr:+.6f}")
assert all(abs(np.corrcoef(proj[:, k], proj_np[:, k])[0, 1]) > 0.999 for k in range(N_COMPONENTS))
print("\nPASS: every component matches numpy exactly up to sign.")

SHAPE expr_array: [2638, 1838]

PCA projection shape:

 (2638, 10)

per-component correlation between LINALDB's PCA and independent numpy PCA:
  component 0: +1.000000
  component 1: -1.000000
  component 2: +1.000000
  component 3: +1.000000
  component 4: +1.000000
  component 5: +1.000000
  component 6: +1.000000
  component 7: -1.000000
  component 8: -1.000000
  component 9: -1.000000

PASS: every component matches numpy exactly up to sign.

Real per-cell-type centroids and nearest-centroid classification

Each real cell's 10-D PCA projection, together with its real published cell type, becomes a row in a real LINALDB Dataset. Per-cell-type centroids (AVG_VEC + GROUP BY, exactly like README.md's own AVG_VEC/GROUP BY example) are computed only from the PCA coordinates -- then every real cell is classified by nearest centroid using LINALDB's own COSINE_SIM, reusing the exact JOIN + ROW_NUMBER() OVER (PARTITION BY ... ORDER BY similarity DESC) pattern clients/python-embedded/examples/digit_classification_embedded.py established -- with the bug workaround above applied (centroid_type, not cell_type, on the centroids side).

db.execute(f"DATASET cells COLUMNS (cell_id: Int, cell_type: String, pc: Vector({N_COMPONENTS}))")
for i in range(n_cells):
    vec = "[" + ", ".join(f"{v:.6f}" for v in proj[i]) + "]"
    db.execute(f'INSERT INTO cells VALUES ({i}, "{cell_types[i]}", {vec})')

db.execute(
    "DATASET centroids FROM cells GROUP BY cell_type "
    "SELECT cell_type AS centroid_type, AVG_VEC(pc) AS centroid, COUNT(*) AS n"
)
print(db.query("SELECT centroid_type, n FROM centroids ORDER BY n DESC").to_string(index=False))

    centroid_type    n
      CD4 T cells 1144
  CD14+ Monocytes  480
          B cells  342
      CD8 T cells  316
         NK cells  154
FCGR3A+ Monocytes  150
  Dendritic cells   37
   Megakaryocytes   15

classify_sql = (
    "WITH classified AS ("
    "SELECT cell_id, cell_type AS true_type, "
    "centroid_type AS predicted_type, "
    "COSINE_SIM(pc, centroid) AS similarity, "
    "ROW_NUMBER() OVER (PARTITION BY cell_id ORDER BY similarity DESC) AS rn "
    "FROM cells JOIN centroids ON COSINE_SIM(pc, centroid) > -2.0"
    ") SELECT cell_id, true_type, predicted_type, similarity FROM classified WHERE rn = 1"
)
result = db.query(classify_sql)
engine_accuracy = (result["true_type"] == result["predicted_type"]).mean()
print(f"engine-computed nearest-centroid accuracy: {engine_accuracy:.4f} "
      f"({(result['true_type'] == result['predicted_type']).sum()}/{len(result)})")

engine-computed nearest-centroid accuracy: 0.9340 (2464/2638)

Independently recomputing the exact same classification in plain numpy from the exported centroid vectors -- same rigor as digit_classification_embedded.py -- to make sure the workaround above actually produces trustworthy results, not just plausible-looking ones.

centroids_df = db.query("SELECT centroid_type, centroid FROM centroids")
centroid_vecs = {row["centroid_type"]: np.array(row["centroid"]) for _, row in centroids_df.iterrows()}

def cosine_similarity(a, b):
    return float(np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b)))

correct = 0
mismatches = []
for i in range(n_cells):
    sims = {ct: cosine_similarity(proj[i], c) for ct, c in centroid_vecs.items()}
    predicted = max(sims, key=sims.get)
    engine_row = result[result["cell_id"] == i].iloc[0]
    if predicted != engine_row["predicted_type"]:
        mismatches.append(i)
    if predicted == cell_types[i]:
        correct += 1

numpy_accuracy = correct / n_cells
print(f"numpy-independent accuracy: {numpy_accuracy:.4f} ({correct}/{n_cells})")
print(f"mismatches between engine and independent numpy classification: {len(mismatches)}")
assert len(mismatches) == 0, "engine and independent numpy disagree -- investigate"
assert abs(numpy_accuracy - engine_accuracy) < 1e-9
print("\nPASS: engine-reported accuracy exactly matches the independently-recomputed one.")

numpy-independent accuracy: 0.9340 (2464/2638)
mismatches between engine and independent numpy classification: 0

PASS: engine-reported accuracy exactly matches the independently-recomputed one.

Persistence and lineage across a real restart

SAVE, drop the Db entirely, reload in a genuinely fresh instance, and confirm EXPLAIN LINEAGE reconstructs the real derivation chain -- the INSERTs into cells and centroids, and the GROUP BY that built centroids from cells -- purely from disk, exactly like 05_lineage_and_linear_algebra.ipynb.

db.execute("SAVE DATASET cells")
db.execute("SAVE DATASET centroids")

del db
db = linaldb.Db(data_dir="./data_single_cell_pca")
db.execute("LOAD DATASET cells")
db.execute("LOAD DATASET centroids")

print(db.execute("EXPLAIN LINEAGE centroids"))
reloaded = db.query("SELECT COUNT(*) AS n FROM cells")
print(reloaded)
assert int(reloaded["n"].iloc[0]) == n_cells
print("PASS: reloaded cell count in a fresh process matches the original exactly.")

Lineage for 'centroids':
DATASET FROM (GROUP BY) (centroids) [08b066c3]
  ROOT (cells) [825c2df4]

      n
0  2638
PASS: reloaded cell count in a fresh process matches the original exactly.

fig, ax = plt.subplots(figsize=(7, 5.5))
palette = plt.get_cmap("tab10")
for idx, ct in enumerate(sorted(set(cell_types))):
    mask = cell_types == ct
    ax.scatter(proj[mask, 0], proj[mask, 1], s=8, alpha=0.6, color=palette(idx), label=ct)
ax.set_xlabel("PC 1")
ax.set_ylabel("PC 2")
ax.set_title("Real PBMC 3k, real cell types -- PCA computed entirely by LINALDB's own SVD")
ax.legend(markerscale=2, fontsize=8, loc="best")
plt.tight_layout()
plt.show()

Conclusion

Two real findings:

1. A real, severe engine bug -- found, fixed, and confirmed fixed: JOIN's SELECT used to silently return the wrong table's value for a qualified column reference when both joined tables share a bare column name (minimal repro above) -- worse than the RANK-as-column- name or FROM dual issues found in 06_llm_low_rank_compression.ipynb, since those error loudly and this one didn't. Root-caused (dsl_expr_to_logical_expr dropped the table qualifier instead of consulting the JOIN schema's own collision-renaming), fixed, tested, and released as engine v0.1.82 / linaldb 0.1.5 on PyPI -- this notebook, re-run against that real published release, confirms the fix directly (the repro cell above now asserts b resolves correctly, the inverse of what it originally asserted).
2. A real, quantified answer to the actual research question: LINALDB's own PCA (10 components, built on its own SVD, verified to match an independent numpy PCA on the same real matrix), applied to the real published PBMC 3k single-cell dataset, is alone enough to recover the real published cell-type labels via nearest-centroid COSINE_SIM classification -- cross-validated two independent ways (the engine's own JOIN/ROW_NUMBER query and a from-scratch numpy recomputation, with zero disagreements between them). This holds whether or not the JOIN'd columns are renamed to avoid collision, now that the bug above is fixed.

See ../tests/ for the package's assertion-backed test suite (tests/qualified_column_test.rs now covers this exact collision shape), and 06_llm_low_rank_compression.ipynb for the other notebook built the same way (real data, real cross-checks, real findings) this session.