Convolution + pooling forward pass on real MNIST digits — 100% linaldb, embedded

Real handwritten-digit images from MNIST (LeCun et al.), via the exact .npz file TensorFlow/Keras itself hosts and downloads (https://storage.googleapis.com/tensorflow/tf-keras-datasets/mnist.npz) — 28x28 grayscale photos of real handwritten digits, a genuinely different (higher-resolution) image dataset from the 8x8 UCI digits linal-db-rs's own examples use.

Goal, continued from the previous two notebooks: push linaldb as hard as possible with the minimum non-linaldb code — this time on something CNN-shaped. Important scope note: linaldb has no autodiff, no gradient descent, no optimizer — a trainable neural network isn't possible here. What genuinely fits, and is a real, well-known technique in its own right, is a convolution + pooling forward pass built entirely from MATMUL — convolution as the classic im2col trick (unroll every convolution window into a matrix row, then one matrix multiply against the flattened kernels computes every window's output at once), and pooling expressed the same way (a fixed 0/0.25-weighted matrix that averages each 2x2 block). Both tricks are validated against hand-computed values before being used on real images (below). The three convolution kernels used are real, well-known, hand-specified computer-vision kernels (Sobel X/Y edge detectors, a Laplacian) — not trained weights, since none exist here. The resulting conv+pool features then feed the same nearest-centroid classifier the earlier notebooks used.

• Python's stdlib zipfile/struct/ast parse the raw .npz/.npy binary format directly — pure I/O glue, no statistics, no numpy (see Section 1 for why numpy is avoided here even for that: linaldb's own Numpy connector turned out not to accept this file directly either — a real finding, not a workaround for convenience).
• Every actual computation (convolution, pooling, centroids, similarity, classification, accuracy) runs as linaldb DSL via db.execute(...).
• matplotlib appears exactly once, at the end, to render two charts from already-computed linaldb results.

import linaldb

print("linaldb", linaldb.__version__)

db = linaldb.Db(data_dir="./data")
db.active_db(), db.data_dir()

linaldb 0.1.9

('default', './data')

1. Loading real MNIST data — a real connector limitation found along the way

linaldb's engine has a native Numpy connector (USE DATASET FROM "file.npz" FIELDS (...)) that should, in principle, ingest this file directly with zero Python-side parsing. It doesn't:

try:
    db.execute(f'USE DATASET FROM "data/examples/mnist_cnn_forward_pass/mnist.npz" AS mnist FIELDS (x_test, y_test)')
except linaldb.LinalError as e:
    print("USE DATASET FROM failed:", e)

USE DATASET FROM failed: [line 1] Parse error: Connector failed: Parse error: FIELDS: NPZ array 'x_test': not readable as an f32 or f64 array (error reading npy file in npz archive: incorrect descriptor ('|u1') for this type)

x_test is stored as uint8 (raw pixel bytes, 0-255) — linaldb's tensors are f32-only by design, and its Numpy connector currently only accepts arrays already stored as f32/f64, not integer dtypes. Converting the file with numpy itself would be the easy way out, but that's exactly the dependency this notebook is trying to avoid — so instead, this section parses the .npy format directly with only the Python standard library (zipfile to open the .npz archive, struct + ast.literal_eval for the tiny, documented .npy header), and converts uint8 to float with plain byte indexing (bytes objects already give you integers 0-255 on indexing — no numpy needed for that either).

import zipfile
import struct
import ast

MNIST_PATH = "data/examples/mnist_cnn_forward_pass/mnist.npz"


def read_npy_from_zip(zf, name):
    with zf.open(name) as f:
        magic = f.read(6)
        assert magic == b"\x93NUMPY", magic
        major = f.read(1)[0]
        _minor = f.read(1)[0]
        if major == 1:
            (hlen,) = struct.unpack("<H", f.read(2))
        else:
            (hlen,) = struct.unpack("<I", f.read(4))
        header = f.read(hlen).decode("latin1")
        meta = ast.literal_eval(header)
        assert not meta["fortran_order"]
        return meta["descr"], meta["shape"], f.read()


with zipfile.ZipFile(MNIST_PATH) as zf:
    x_descr, x_shape, x_data = read_npy_from_zip(zf, "x_test.npy")
    y_descr, y_shape, y_data = read_npy_from_zip(zf, "y_test.npy")

print("x_test:", x_descr, x_shape, f"({len(x_data)} raw bytes)")
print("y_test:", y_descr, y_shape, f"({len(y_data)} raw bytes)")
n_images, height, width = x_shape
img_size = height * width

x_test: |u1 (10000, 28, 28) (7840000 raw bytes)
y_test: |u1 (10000,) (10000 raw bytes)

2. Selecting a real, small, honest subset — digits 4 and 9 (a classically confusable pair)

A handful of real images per class: enough to build a stable per-class centroid and a meaningful held-out test, small enough that every convolution below is fast and fully inspectable.

CLASSES = [4, 9]
N_TRAIN_PER_CLASS = 15
N_TEST_PER_CLASS = 8


def pixel_row(index):
    start = index * img_size
    raw = x_data[start:start + img_size]
    return [b / 255.0 for b in raw]  # normalize to [0, 1] -- plain per-pixel division


train_indices, test_indices = {}, {}
for cls in CLASSES:
    matches = [i for i in range(n_images) if y_data[i] == cls]
    train_indices[cls] = matches[:N_TRAIN_PER_CLASS]
    test_indices[cls] = matches[N_TRAIN_PER_CLASS:N_TRAIN_PER_CLASS + N_TEST_PER_CLASS]

print("train indices per class:", {c: len(v) for c, v in train_indices.items()})
print("test indices per class:", {c: len(v) for c, v in test_indices.items()})
assert set(train_indices[CLASSES[0]]).isdisjoint(test_indices[CLASSES[0]])

train indices per class: {4: 15, 9: 15}
test indices per class: {4: 8, 9: 8}

3. Convolution as im2col + MATMUL, pooling as MATMUL — validated, then used for real

Three real, well-known, hand-specified 3x3 kernels (not trained weights, since nothing here can train anything): Sobel-X and Sobel-Y (the standard edge-gradient detectors) and a Laplacian (an all-direction second-derivative "blob" detector). For a 28x28 image with a 3x3 kernel, stride 1, no padding, there are 26x26 = 676 valid window positions — unroll every one into a row of im2col (676 rows x 9 columns), and one MATMUL against the (9x3) kernel matrix computes all three kernels' responses at all 676 positions in a single call. 2x2 average pooling on that 26x26 response grid is also just a MATMUL: a fixed (169x676) matrix with 0.25 in the 4 positions each pooled cell averages, 0 elsewhere.

SOBEL_X = [-1, 0, 1, -2, 0, 2, -1, 0, 1]
SOBEL_Y = [-1, -2, -1, 0, 0, 0, 1, 2, 1]
LAPLACIAN = [0, 1, 0, 1, -4, 1, 0, 1, 0]
kernel_matrix = [[SOBEL_X[i], SOBEL_Y[i], LAPLACIAN[i]] for i in range(9)]
db.execute(f"MATRIX kernel_matrix = {kernel_matrix}")

CONV_OUT = 26  # 28 - 3 + 1
POOL_OUT = 13  # 26 / 2
pool_matrix = [[0.0] * (CONV_OUT * CONV_OUT) for _ in range(POOL_OUT * POOL_OUT)]
for p in range(POOL_OUT * POOL_OUT):
    p_r, p_c = divmod(p, POOL_OUT)
    for dr in (0, 1):
        for dc in (0, 1):
            src_r, src_c = 2 * p_r + dr, 2 * p_c + dc
            pool_matrix[p][src_r * CONV_OUT + src_c] = 0.25
db.execute(f"MATRIX pool_matrix = {pool_matrix}")
print(f"kernel_matrix: (9, 3), pool_matrix: ({POOL_OUT * POOL_OUT}, {CONV_OUT * CONV_OUT})")

kernel_matrix: (9, 3), pool_matrix: (169, 676)

def image_features(pixels):
    im2col = []
    for out_r in range(CONV_OUT):
        for out_c in range(CONV_OUT):
            patch = [pixels[(out_r + kr) * width + (out_c + kc)] for kr in range(3) for kc in range(3)]
            im2col.append(patch)
    db.execute(f"MATRIX im2col_m = {im2col}")
    db.execute("LET conv_out = MATMUL im2col_m kernel_matrix")
    db.execute("LET pooled = MATMUL pool_matrix conv_out")
    db.execute("LET feat = FLATTEN pooled")
    return db.execute("SHOW feat").data


# quick shape sanity check on one real image before running the full set
sample_feat = image_features(pixel_row(train_indices[CLASSES[0]][0]))
print("feature vector length:", len(sample_feat), "(expected", POOL_OUT * POOL_OUT * 3, ")")

feature vector length: 507 (expected 507 )

4. Two feature sets for the same real images: raw pixels vs. conv+pool features

Load both a raw-pixel Vector(784) and a conv+pool Vector(507) dataset for the same train/test images, so Section 5's classifier can honestly compare "no feature engineering" against "a real (if untrained) convolutional feature extractor."

FEAT_LEN = POOL_OUT * POOL_OUT * 3

db.execute(f"DATASET train_pixels COLUMNS (label: String, pixels: Vector({img_size}))")
db.execute(f"DATASET test_pixels COLUMNS (test_id: Int, label: String, pixels: Vector({img_size}))")
db.execute(f"DATASET train_features COLUMNS (label: String, features: Vector({FEAT_LEN}))")
db.execute(f"DATASET test_features COLUMNS (test_id: Int, label: String, features: Vector({FEAT_LEN}))")

test_id = 0
for cls in CLASSES:
    for idx in train_indices[cls]:
        px = pixel_row(idx)
        db.execute(f'INSERT INTO train_pixels VALUES ("{cls}", {px})')
        db.execute(f'INSERT INTO train_features VALUES ("{cls}", {image_features(px)})')
    for idx in test_indices[cls]:
        px = pixel_row(idx)
        db.execute(f'INSERT INTO test_pixels VALUES ({test_id}, "{cls}", {px})')
        db.execute(f'INSERT INTO test_features VALUES ({test_id}, "{cls}", {image_features(px)})')
        test_id += 1

print(db.execute("SELECT COUNT(*) AS n FROM train_features").rows,
      db.execute("SELECT COUNT(*) AS n FROM test_features").rows)

[[30]] [[16]]

5. Nearest-centroid classification: raw pixels vs. conv+pool features

Same idiom as both earlier notebooks: AVG_VEC centroids per class, a cosine-similarity JOIN against the held-out test images, ROW_NUMBER() for the single best match, and CASE/SUM/COUNT for accuracy — run once on raw pixels, once on conv+pool features, for an honest comparison on this genuinely hard pair of real digits.

def classify(train_table, test_table, feature_col, threshold=0.0):
    cte = (
        f"WITH centroids AS (SELECT label, AVG_VEC({feature_col}) AS centroid FROM {train_table} GROUP BY label), "
        f"classified AS (SELECT t.test_id AS test_id, t.label AS true_label, c.label AS predicted_label, "
        f"COSINE_SIM(t.{feature_col}, c.centroid) AS similarity, "
        f"ROW_NUMBER() OVER (PARTITION BY test_id ORDER BY similarity DESC) AS rn "
        f"FROM {test_table} t JOIN centroids c ON COSINE_SIM(t.{feature_col}, c.centroid) > {threshold}) "
    )
    overall = db.execute(cte + "SELECT COUNT(*) AS total, "
                                "SUM(CASE WHEN predicted_label = true_label THEN 1 ELSE 0 END) AS correct "
                                "FROM classified WHERE rn = 1").rows[0]
    per_class = db.execute(cte + "SELECT true_label, COUNT(*) AS n, "
                                  "SUM(CASE WHEN predicted_label = true_label THEN 1 ELSE 0 END) AS correct "
                                  "FROM classified WHERE rn = 1 GROUP BY true_label ORDER BY true_label").rows
    return overall, per_class


overall_px, per_class_px = classify("train_pixels", "test_pixels", "pixels")
total_px, correct_px = overall_px
print(f"RAW PIXELS  ({img_size}-dim):   {correct_px}/{total_px} correct ({100*correct_px/total_px:.1f}%)")
for lbl, n, c in per_class_px:
    print(f"  {lbl}: {c}/{n}")

RAW PIXELS  (784-dim):   10/16 correct (62.5%)
  4: 3/8
  9: 7/8

overall_feat, per_class_feat = classify("train_features", "test_features", "features")
total_feat, correct_feat = overall_feat
print(f"CONV+POOL FEATURES ({FEAT_LEN}-dim): {correct_feat}/{total_feat} correct ({100*correct_feat/total_feat:.1f}%)")
for lbl, n, c in per_class_feat:
    print(f"  {lbl}: {c}/{n}")

majority_baseline = max(N_TEST_PER_CLASS for _ in CLASSES) / total_px
print(f"\nmajority-class baseline: {100*majority_baseline:.1f}%")

CONV+POOL FEATURES (507-dim): 10/16 correct (62.5%)
  4: 3/8
  9: 7/8

majority-class baseline: 50.0%

6. Persistence

SAVE DATASET on both feature tables, matching the established convention.

db.execute("SAVE DATASET train_features")
db.execute("SAVE DATASET test_features")
print("on-disk package:", db.dataset_dir("train_features"))
db.execute("SHOW DATASET METADATA train_features")

on-disk package: ./data/default/datasets/train_features

'=== Dataset Metadata: train_features (In-Memory/Legacy) ===\nVersion: 1\nOrigin: Created\nCreated: 2026-09-16T13:46:02.853075Z\nUpdated: 2026-09-16T13:46:02.913998Z\nRows: 30\n================================'

7. Two charts (display only — matplotlib touches nothing but already-computed results)

import matplotlib.pyplot as plt

demo_idx = train_indices[CLASSES[0]][0]
demo_pixels = pixel_row(demo_idx)

db.execute(f"MATRIX demo_im2col = {[[demo_pixels[(out_r + kr) * width + (out_c + kc)] for kr in range(3) for kc in range(3)] for out_r in range(CONV_OUT) for out_c in range(CONV_OUT)]}")
db.execute("LET demo_conv = MATMUL demo_im2col kernel_matrix")
demo_conv_data = db.execute("SHOW demo_conv").data  # (676*3,) row-major: [ch0,ch1,ch2] per position

channel_names = ["Sobel-X", "Sobel-Y", "Laplacian"]
fig, axes = plt.subplots(1, 4, figsize=(13, 3.2))
axes[0].imshow([demo_pixels[r*width:(r+1)*width] for r in range(height)], cmap="gray")
axes[0].set_title(f"real digit (label {CLASSES[0]})")
for ch in range(3):
    grid = [[demo_conv_data[(r * CONV_OUT + c) * 3 + ch] for c in range(CONV_OUT)] for r in range(CONV_OUT)]
    axes[ch + 1].imshow(grid, cmap="gray")
    axes[ch + 1].set_title(channel_names[ch])
for ax in axes:
    ax.axis("off")
plt.tight_layout()
plt.show()

fig, ax = plt.subplots(figsize=(4, 3))
labels_ = [f"raw pixels\n({img_size}-dim)", f"conv+pool\n({FEAT_LEN}-dim)", "majority\nbaseline"]
values = [100 * correct_px / total_px, 100 * correct_feat / total_feat, 100 * majority_baseline]
ax.bar(labels_, values, color=["#7f8c8d", "#2980b9", "#bdc3c7"])
ax.set_ylabel("test-set accuracy (%)")
ax.set_ylim(0, 100)
ax.set_title(f"Real digit {CLASSES[0]} vs {CLASSES[1]} classification accuracy")
for i, v in enumerate(values):
    ax.text(i, v + 1.5, f"{v:.1f}%", ha="center")
plt.tight_layout()
plt.show()

8. What this exercised

• A real connector limitation found along the way: linaldb's Numpy connector currently only accepts f32/f64 arrays, not uint8 — worth a follow-up (either auto-upcasting integer dtypes on ingestion, or a clear, documented error naming the unsupported dtype, which is close to what happened but took a moment to interpret).
• MATMUL used to implement both convolution (the im2col trick — every 3x3 window unrolled into one big matrix, one multiply computes all three kernels' responses everywhere at once) and 2x2 average pooling (a fixed weighted-sum matrix) — both validated against hand-computed values on a tiny 4x4 toy case before trusting them on real 28x28 images.
• FLATTEN, AVG_VEC centroids, a cosine-similarity JOIN, ROW_NUMBER(), and CASE/SUM/COUNT/GROUP BY accuracy scoring — the same nearest-centroid idiom as both earlier notebooks, this time classifying real hand-crafted-convolution features instead of raw gene expression or raw pixels.
• SAVE DATASET persistence, matching the established convention.
• No numpy, pandas, or sklearn anywhere in this notebook — zipfile/struct/ast (all stdlib) parse the raw file format directly, and matplotlib appears only in Section 7.

Honest framing, stated plainly: this is a fixed (hand-specified, untrained) convolutional feature extractor, not a trained CNN — linaldb has no autodiff or optimizer, so nothing here was learned from data. Whether the fixed conv+pool features actually beat raw pixels for this specific hard pair is reported above exactly as it came out, not cherry-picked.