From 5641ac39443925fb66e0bf7609f16a349a890a02 Mon Sep 17 00:00:00 2001
From: Hideki <izumi.hideki.36w@st.kyoto-u.ac.jp>
Date: Tue, 9 Jun 2026 17:07:02 +0900
Subject: [PATCH 1/6] Add optional numba solver

---
 .gitignore                     |   2 +
 pyproject.toml                 |   7 +
 speedtests/.gitignore          |   2 +
 speedtests/README.md           |  86 +++++
 speedtests/local.sh.example    |  12 +
 speedtests/speedrun.py         |  35 ++
 speedtests/speedrun.sh         |  73 ++++
 speedtests/speedrun_legacy.py  |  34 ++
 speedtests/speedrun_legacy.sh  |  56 +++
 speedtests/speedrun_mkl.sh     |  79 +++++
 speedtests/speedtest_common.py | 608 +++++++++++++++++++++++++++++++++
 src/fastl2lir/_numba.py        | 157 +++++++++
 src/fastl2lir/fastl2lir.py     |  75 ++--
 tests/test__numba.py           | 122 +++++++
 tests/test_fastl2lir.py        |   7 +-
 15 files changed, 1332 insertions(+), 23 deletions(-)
 create mode 100644 speedtests/.gitignore
 create mode 100644 speedtests/README.md
 create mode 100644 speedtests/local.sh.example
 create mode 100644 speedtests/speedrun.py
 create mode 100755 speedtests/speedrun.sh
 create mode 100644 speedtests/speedrun_legacy.py
 create mode 100755 speedtests/speedrun_legacy.sh
 create mode 100755 speedtests/speedrun_mkl.sh
 create mode 100644 speedtests/speedtest_common.py
 create mode 100644 src/fastl2lir/_numba.py
 create mode 100644 tests/test__numba.py

diff --git a/.gitignore b/.gitignore
index bdce9d0..a1c2a5a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,6 +2,8 @@
 ._*
 
 *.pyc
+*.nbc
+*.nbi
 *.egg-info
 build
 dist
diff --git a/pyproject.toml b/pyproject.toml
index 461cc25..9e48369 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,6 +17,13 @@ dependencies = [
     "threadpoolctl>=2.1.0",
     "tqdm>=4.64.1",
 ]
+
+[project.optional-dependencies]
+numba = [
+    "numba>=0.59",
+    "scipy>=1.8.0",
+]
+
 [project.urls]
 Homepage = "https://github.com/KamitaniLab/PyFastL2LiR"
 Repository = "https://github.com/KamitaniLab/PyFastL2LiR"
diff --git a/speedtests/.gitignore b/speedtests/.gitignore
new file mode 100644
index 0000000..e5d4295
--- /dev/null
+++ b/speedtests/.gitignore
@@ -0,0 +1,2 @@
+local.sh
+*.csv
\ No newline at end of file
diff --git a/speedtests/README.md b/speedtests/README.md
new file mode 100644
index 0000000..95d7c2d
--- /dev/null
+++ b/speedtests/README.md
@@ -0,0 +1,86 @@
+# Speed Tests
+
+This directory contains speed-run scripts for measuring PyFastL2LiR fitting
+runtime. These scripts are not the package test suite; correctness tests live
+under `tests/`.
+
+The speed runs generate synthetic fMRI-like regression data, fit FastL2LiR on
+several target shapes, compare predictions between runners when more than one
+runner is used, print timing summaries, and save CSV summaries in this
+directory.
+
+## Files
+
+- `speedrun.sh`: Run the current working-tree package with the project Python
+  environment. It compares `solver="numpy"` and `solver="numba"`.
+- `speedrun.py`: Python entry point used by `speedrun.sh`.
+- `speedrun_legacy.sh`: Run the legacy installed `fastl2lir` package from a
+  configured legacy Python environment.
+- `speedrun_legacy.py`: Python entry point used by `speedrun_legacy.sh`.
+- `speedrun_mkl.sh`: Run the current working-tree package from a configured
+  MKL-oriented Python environment.
+- `speedtest_common.py`: Shared benchmark setup, timing, prediction checks,
+  summaries, and CSV writing.
+- `local.sh.example`: Template for machine-local environment paths.
+
+## Local Setup
+
+Copy `local.sh.example` to `local.sh` and edit paths for the machine:
+
+```bash
+cp speedtests/local.sh.example speedtests/local.sh
+```
+
+`local.sh` is ignored by git. Relative paths are resolved from the repository
+root.
+
+## Running
+
+Current environment:
+
+```bash
+./speedtests/speedrun.sh
+```
+
+Legacy environment:
+
+```bash
+./speedtests/speedrun_legacy.sh
+```
+
+MKL-oriented environment:
+
+```bash
+./speedtests/speedrun_mkl.sh
+```
+
+For a quick smoke run, restrict the synthetic data size and case list:
+
+```bash
+BENCH_N=24 BENCH_P=12 BENCH_K=4 BENCH_REPEATS=1 BENCH_CASES=fc8 ./speedtests/speedrun.sh
+```
+
+## Controls
+
+- `BENCH_N`: Number of samples. Default: `1000`.
+- `BENCH_P`: Number of input features or voxels. Default: `15000`.
+- `BENCH_K`: Number of selected features. Default: `500`.
+- `BENCH_REPEATS`: Repeats per benchmark case. Default: `3`.
+- `BENCH_CASES`: Comma-separated subset of `fc8`, `fc6`, `conv5`, and
+  `conv5_chunk10`.
+- `BENCH_RUNNERS`: Comma-separated subset of available runner labels. The
+  current speed run supports `numpy` and `numba`; the legacy speed run supports
+  `legacy`.
+- `BENCH_SERVER`: Label written to the terminal output and CSV. Defaults to the
+  hostname.
+
+## Outputs
+
+CSV files are written under `speedtests/`:
+
+- `speedtest_<server>.csv` for the current and MKL-oriented speed runs.
+- `speedtest_legacy_<server>.csv` for the legacy speed run.
+
+The CSV columns include the server label, benchmark case, matrix dimensions,
+repeat count, median runtime columns, speedup columns when multiple runners are
+present, and the maximum prediction difference on a small validation slice.
diff --git a/speedtests/local.sh.example b/speedtests/local.sh.example
new file mode 100644
index 0000000..184728d
--- /dev/null
+++ b/speedtests/local.sh.example
@@ -0,0 +1,12 @@
+# Copy this file to speedtests/local.sh and edit paths for your machine.
+# speedtests/local.sh is intentionally ignored by git.
+
+# Python environment used by speedtests/speedrun.sh. Relative paths are resolved from the
+# repository root.
+VENV_PATH=.venv
+
+# Legacy environment used by speedtests/speedrun_legacy.sh.
+LEGACY_CONDA=.venv-legacy
+
+# MKL-oriented environment used by speedtests/speedrun_mkl.sh.
+MKL_CONDA=.venv-mkl
diff --git a/speedtests/speedrun.py b/speedtests/speedrun.py
new file mode 100644
index 0000000..e87fe5b
--- /dev/null
+++ b/speedtests/speedrun.py
@@ -0,0 +1,35 @@
+"""Current-environment speed test for PyFastL2LiR implementations.
+
+This script compares L2-regularized linear regression fits on synthetic
+fMRI-like matrices in the active project Python environment. The mathematical
+objects are a Gaussian input matrix X, Gaussian target arrays Y with dense and
+convolutional feature shapes, ridge weights W, bias b, and predictions on a
+small validation slice. The execution stages are environment reporting,
+synthetic data generation, repeated fitting for the package numpy solver, the
+package numba solver, prediction-difference checks, median runtime summaries,
+and CSV writing. The saved output is speedtests/speedtest_<server>.csv.
+
+Run directly with:
+    speedtests/speedrun.sh
+"""
+
+import os
+import socket
+from pathlib import Path
+
+import fastl2lir
+
+from speedtest_common import run_speedtest
+
+server_name = os.environ.get("BENCH_SERVER") or socket.gethostname()
+
+run_speedtest(
+    fastl2lir_module=fastl2lir,
+    runners=[
+        {"label": "numpy", "fit_kwargs": {"solver": "numpy"}},
+        {"label": "numba", "fit_kwargs": {"solver": "numba"}},
+    ],
+    output_csv=Path("speedtests") / f"speedtest_{server_name}.csv",
+    comparison="current FastL2LiR.fit numpy and numba solver comparison",
+    module_names=("numba",),
+)
diff --git a/speedtests/speedrun.sh b/speedtests/speedrun.sh
new file mode 100755
index 0000000..f8a42d5
--- /dev/null
+++ b/speedtests/speedrun.sh
@@ -0,0 +1,73 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+script_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+repo_root="$(cd "${script_dir}/.." && pwd)"
+local_config="${script_dir}/local.sh"
+
+if [[ -f "$local_config" ]]; then
+    # shellcheck source=/dev/null
+    source "$local_config"
+fi
+
+resolve_path() {
+    local path_value="$1"
+
+    if [[ "$path_value" = /* ]]; then
+        printf "%s\n" "$path_value"
+    else
+        printf "%s/%s\n" "$repo_root" "$path_value"
+    fi
+}
+
+cd "$repo_root"
+
+# Use the project Python environment only. Thread handling is intentionally
+# left to fastl2lir and the runtime environment under test.
+export PYTHONUNBUFFERED="${PYTHONUNBUFFERED:-1}"
+export BENCH_SERVER="${BENCH_SERVER:-$(hostname)}"
+
+VENV_PATH="${VENV_PATH:-.venv}"
+python_bin="$(resolve_path "$VENV_PATH")/bin/python"
+
+if [[ ! -x "$python_bin" ]]; then
+    if [[ "$VENV_PATH" != ".venv" ]]; then
+        echo "ERROR: Python is not executable: $python_bin" >&2
+        echo "Create that environment or update VENV_PATH in speedtests/local.sh." >&2
+        exit 1
+    fi
+    if ! command -v uv >/dev/null 2>&1; then
+        echo "ERROR: .venv is missing and uv is not available." >&2
+        echo "Install uv or set VENV_PATH in speedtests/local.sh." >&2
+        exit 1
+    fi
+    uv sync --extra numba
+    python_bin="$(resolve_path "$VENV_PATH")/bin/python"
+fi
+
+if ! "$python_bin" - <<'PY'
+import sys
+
+if sys.version_info < (3, 10):
+    raise SystemExit("Python >= 3.10 is required")
+
+import fastl2lir  # noqa: F401
+import numba  # noqa: F401
+import numpy  # noqa: F401
+import scipy  # noqa: F401
+PY
+then
+    if [[ "$VENV_PATH" != ".venv" ]]; then
+        echo "ERROR: Python environment is incomplete: $python_bin" >&2
+        echo "Install fastl2lir, numpy, scipy, and numba there, or update VENV_PATH." >&2
+        exit 1
+    fi
+    if ! command -v uv >/dev/null 2>&1; then
+        echo "ERROR: Python environment is incomplete and uv is not available." >&2
+        echo "Set VENV_PATH in speedtests/local.sh to a complete environment." >&2
+        exit 1
+    fi
+    uv sync --extra numba
+fi
+
+exec "$python_bin" -u speedtests/speedrun.py "$@"
diff --git a/speedtests/speedrun_legacy.py b/speedtests/speedrun_legacy.py
new file mode 100644
index 0000000..ba08d2f
--- /dev/null
+++ b/speedtests/speedrun_legacy.py
@@ -0,0 +1,34 @@
+"""Legacy-environment speed test for PyFastL2LiR.
+
+This script measures L2-regularized linear regression runtime for the legacy
+FastL2LiR implementation in the configured legacy Python environment. The
+mathematical objects are synthetic Gaussian input matrices X, synthetic target
+arrays Y, ridge weights W, bias b, and validation predictions used to check
+repeat consistency. The execution stages are environment reporting, benchmark
+case construction, repeated legacy fitting, median runtime summarization, and
+CSV writing. The saved output is speedtests/speedtest_legacy_<server>.csv.
+
+Run directly with:
+    speedtests/speedrun_legacy.sh
+"""
+
+import os
+import socket
+from pathlib import Path
+
+import fastl2lir
+
+from speedtest_common import run_speedtest
+
+
+server_name = os.environ.get("BENCH_SERVER") or socket.gethostname()
+
+run_speedtest(
+    fastl2lir_module=fastl2lir,
+    runners=[
+        {"label": "legacy", "fit_kwargs": {}},
+    ],
+    output_csv=Path("speedtests") / f"speedtest_legacy_{server_name}.csv",
+    comparison="legacy FastL2LiR.fit runtime",
+    module_names=("scipy",),
+)
diff --git a/speedtests/speedrun_legacy.sh b/speedtests/speedrun_legacy.sh
new file mode 100755
index 0000000..0b3d237
--- /dev/null
+++ b/speedtests/speedrun_legacy.sh
@@ -0,0 +1,56 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+script_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+repo_root="$(cd "${script_dir}/.." && pwd)"
+local_config="${script_dir}/local.sh"
+
+if [[ -f "$local_config" ]]; then
+    # shellcheck source=/dev/null
+    source "$local_config"
+fi
+
+resolve_path() {
+    local path_value="$1"
+
+    if [[ "$path_value" = /* ]]; then
+        printf "%s\n" "$path_value"
+    else
+        printf "%s/%s\n" "$repo_root" "$path_value"
+    fi
+}
+
+cd "$repo_root"
+
+export PYTHONUNBUFFERED="${PYTHONUNBUFFERED:-1}"
+export BENCH_SERVER="${BENCH_SERVER:-$(hostname)}"
+
+LEGACY_CONDA="${LEGACY_CONDA:-.venv-legacy}"
+conda_root="$(resolve_path "$LEGACY_CONDA")"
+python_bin="${conda_root}/bin/python"
+
+if [[ ! -x "$python_bin" ]]; then
+    echo "ERROR: legacy Python is not executable: $python_bin" >&2
+    echo "Set LEGACY_CONDA in speedtests/local.sh." >&2
+    exit 1
+fi
+
+if [[ -f "${conda_root}/etc/profile.d/conda.sh" ]]; then
+    source "${conda_root}/etc/profile.d/conda.sh"
+    conda activate base
+else
+    export PATH="${conda_root}/bin:${PATH}"
+fi
+
+"$python_bin" - <<'PY'
+import sys
+
+import fastl2lir  # noqa: F401
+import numpy  # noqa: F401
+import scipy  # noqa: F401
+
+if sys.version_info[:2] != (3, 8):
+    print(f"WARNING: expected Python 3.8, got {sys.version.split()[0]}")
+PY
+
+exec "$python_bin" -u speedtests/speedrun_legacy.py "$@"
diff --git a/speedtests/speedrun_mkl.sh b/speedtests/speedrun_mkl.sh
new file mode 100755
index 0000000..199c6f5
--- /dev/null
+++ b/speedtests/speedrun_mkl.sh
@@ -0,0 +1,79 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+script_dir="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+repo_root="$(cd "${script_dir}/.." && pwd)"
+local_config="${script_dir}/local.sh"
+
+if [[ -f "$local_config" ]]; then
+    # shellcheck source=/dev/null
+    source "$local_config"
+fi
+
+resolve_path() {
+    local path_value="$1"
+
+    if [[ "$path_value" = /* ]]; then
+        printf "%s\n" "$path_value"
+    else
+        printf "%s/%s\n" "$repo_root" "$path_value"
+    fi
+}
+
+cd "$repo_root"
+
+# Use the requested MKL-oriented environment while importing this working
+# tree's src package. Override BENCH_RUNNERS or MKL_BENCH_RUNNERS to change the
+# solver subset for a specific local environment.
+export PYTHONUNBUFFERED="${PYTHONUNBUFFERED:-1}"
+export BENCH_SERVER="${BENCH_SERVER:-$(hostname)_mkl}"
+MKL_BENCH_RUNNERS="${MKL_BENCH_RUNNERS:-numpy,numba}"
+export BENCH_RUNNERS="${BENCH_RUNNERS:-$MKL_BENCH_RUNNERS}"
+export PYTHONPATH="${repo_root}/src${PYTHONPATH:+:${PYTHONPATH}}"
+
+MKL_CONDA="${MKL_CONDA:-.venv-mkl}"
+mkl_root="$(resolve_path "$MKL_CONDA")"
+python_bin="${mkl_root}/bin/python"
+
+if [[ ! -x "$python_bin" ]]; then
+    echo "ERROR: MKL Python is not executable: $python_bin" >&2
+    echo "Set MKL_CONDA in speedtests/local.sh." >&2
+    exit 1
+fi
+
+"$python_bin" - <<'PY'
+import sys
+
+if sys.version_info < (3, 10):
+    raise SystemExit("Python >= 3.10 is required")
+
+missing = []
+for module_name in ("numpy", "scipy", "threadpoolctl"):
+    try:
+        __import__(module_name)
+    except ImportError:
+        missing.append(module_name)
+
+if missing:
+    raise SystemExit(
+        "Missing MKL environment dependencies: "
+        + ", ".join(missing)
+        + "\nInstall them in the MKL environment, then rerun speedtests/speedrun_mkl.sh."
+    )
+
+import numpy
+import scipy
+import threadpoolctl
+
+import fastl2lir
+
+print("MKL speed-test environment")
+print(f"  python: {sys.version.split()[0]}")
+print(f"  numpy: {numpy.__version__}")
+print(f"  scipy: {scipy.__version__}")
+print(f"  fastl2lir file: {fastl2lir.__file__}")
+print(f"  threadpoolctl: {threadpoolctl.__version__}")
+print(f"  threadpool info: {threadpoolctl.threadpool_info()}")
+PY
+
+exec "$python_bin" -u speedtests/speedrun.py "$@"
diff --git a/speedtests/speedtest_common.py b/speedtests/speedtest_common.py
new file mode 100644
index 0000000..82adce2
--- /dev/null
+++ b/speedtests/speedtest_common.py
@@ -0,0 +1,608 @@
+"""Shared speed-test routines for PyFastL2LiR solver experiments.
+
+This module benchmarks L2-regularized linear regression fits on synthetic
+fMRI-like design matrices. The mathematical objects are an input matrix X
+(samples by voxels), a target tensor Y (samples by model units), selected
+feature counts for ridge regression, fitted weights W, fitted bias b, and
+short-slice predictions used to check numerical agreement between runners.
+
+Execution proceeds by reading benchmark dimensions from environment variables,
+building named target-shape cases, generating fresh Gaussian X and Y arrays for
+each repeat, fitting each configured runner in alternating order, comparing
+predictions against the first runner, summarizing median runtimes, and writing a
+CSV table. Saved outputs are speed-test CSV files under speedtests/ whose columns
+include the server name, case dimensions, per-runner median fit times,
+speedups relative to the baseline runner, and maximum prediction differences.
+"""
+
+import importlib
+import os
+import socket
+import statistics
+import sys
+import time
+from pathlib import Path
+
+import numpy as np
+
+
+RNG = np.random.default_rng()
+
+
+def read_positive_int(name, default):
+    """Read a positive integer setting from the environment.
+
+    Inputs:
+        name: Environment variable name.
+        default: Integer value used when the variable is unset.
+    Output:
+        The parsed positive integer.
+    This function validates that benchmark dimensions and repeat counts are
+    positive so the experiment fails before allocating arrays.
+    """
+    value = os.environ.get(name)
+    if value is None:
+        return default
+
+    try:
+        parsed = int(value)
+    except ValueError as exc:
+        raise SystemExit(f"{name} must be a positive integer") from exc
+
+    if parsed < 1:
+        raise SystemExit(f"{name} must be a positive integer")
+
+    return parsed
+
+
+def get_benchmark_settings():
+    """Collect benchmark dimensions and repeat counts.
+
+    Inputs:
+        Environment variables BENCH_N, BENCH_P, BENCH_K, and BENCH_REPEATS.
+    Output:
+        Tuple of sample count n, voxel count p, selected feature count k, and
+        repeat count.
+    This function centralizes the size controls used by all speed-test cases.
+    """
+    n = read_positive_int("BENCH_N", 1000)
+    p = read_positive_int("BENCH_P", 15000)
+    k = read_positive_int("BENCH_K", 500)
+    repeats = read_positive_int("BENCH_REPEATS", 3)
+    return n, p, k, repeats
+
+
+def make_cases(k, repeats):
+    """Build the target-shape cases for the speed test.
+
+    Inputs:
+        k: Number of selected input features used by ridge fitting.
+        repeats: Number of repeats for each case.
+        Optional BENCH_CASES environment variable with comma-separated case
+        names.
+    Output:
+        List of tuples containing case name, Y shape after the sample axis,
+        fit parameters, and repeat count.
+    This function defines dense and chunked target layouts that mimic neural
+    network feature layers.
+    """
+    cases = [
+        ("fc8", (1000,), dict(alpha=100.0, n_feat=k), repeats),
+        ("fc6", (4096,), dict(alpha=100.0, n_feat=k), repeats),
+        (
+            "conv5_chunk5",
+            (5, 14, 14),
+            dict(alpha=100.0, n_feat=k, chunk_size=196),
+            repeats,
+        ),
+        (
+            "conv4_chunk5",
+            (5, 28, 28),
+            dict(alpha=100.0, n_feat=k, chunk_size=784),
+            repeats,
+        ),
+    ]
+
+    case_filter = os.environ.get("BENCH_CASES")
+    if not case_filter:
+        return cases
+
+    selected_cases = {name.strip() for name in case_filter.split(",") if name.strip()}
+    known_cases = {name for name, *_ in cases}
+    unknown_cases = selected_cases - known_cases
+    if unknown_cases:
+        raise SystemExit(
+            "Unknown BENCH_CASES entries: " + ", ".join(sorted(unknown_cases))
+        )
+
+    return [case for case in cases if case[0] in selected_cases]
+
+
+def filter_runners(runners):
+    """Select benchmark runners from BENCH_RUNNERS when requested.
+
+    Inputs:
+        runners: Ordered list of available runner dictionaries.
+        Optional BENCH_RUNNERS environment variable with comma-separated runner
+        labels.
+    Output:
+        Ordered list of selected runner dictionaries.
+    This function lets a researcher run the same speed-test script in
+    environments that lack optional solvers such as numba while keeping the
+    experiment definition explicit in the terminal log and CSV columns.
+    """
+    runner_filter = os.environ.get("BENCH_RUNNERS")
+    if not runner_filter:
+        return runners
+
+    selected_labels = {
+        label.strip() for label in runner_filter.split(",") if label.strip()
+    }
+    known_labels = {runner["label"] for runner in runners}
+    unknown_labels = selected_labels - known_labels
+    if unknown_labels:
+        raise SystemExit(
+            "Unknown BENCH_RUNNERS entries: " + ", ".join(sorted(unknown_labels))
+        )
+
+    selected_runners = [
+        runner for runner in runners if runner["label"] in selected_labels
+    ]
+    if not selected_runners:
+        raise SystemExit("BENCH_RUNNERS did not select any runners")
+
+    return selected_runners
+
+
+def module_version(module_name):
+    """Return a module version string for the environment report.
+
+    Inputs:
+        module_name: Importable Python module name.
+    Output:
+        Version string, '(unknown)', or '(not installed)'.
+    This function records the numerical stack used for a speed-test run.
+    """
+    try:
+        module = importlib.import_module(module_name)
+    except ImportError:
+        return "(not installed)"
+    return getattr(module, "__version__", "(unknown)")
+
+
+def package_file(module):
+    """Return the resolved file path for an imported module.
+
+    Inputs:
+        module: Imported Python module object.
+    Output:
+        Absolute path to the module file.
+    This function helps distinguish installed package code from local source
+    modules during benchmark reporting.
+    """
+    return str(Path(module.__file__).resolve())
+
+
+def make_data(n, p, y_shape):
+    """Generate synthetic regression data for one benchmark repeat.
+
+    Inputs:
+        n: Number of samples.
+        p: Number of input features or voxels.
+        y_shape: Target-unit shape after the sample axis.
+    Output:
+        Tuple (X, Y) of float64 Gaussian arrays.
+    This function creates fresh comparable input and target arrays for all
+    runners in a single repeat.
+    """
+    X = RNG.normal(size=(n, p)).astype(np.float64)
+    Y = RNG.normal(size=(n,) + y_shape).astype(np.float64)
+    return X, Y
+
+
+def runner_module(fastl2lir_module, runner):
+    """Resolve the FastL2LiR module used by one runner.
+
+    Inputs:
+        fastl2lir_module: Default module used when a runner has no override.
+        runner: Runner dictionary.
+    Output:
+        Module object that provides FastL2LiR.
+    This function lets a speed test compare runners that may use different
+    modules while sharing the same benchmark data.
+    """
+    return runner.get("module", fastl2lir_module)
+
+
+def fit_once(fastl2lir_module, X, Y, runner, params):
+    """Fit one runner on one generated dataset.
+
+    Inputs:
+        fastl2lir_module: Default module that provides FastL2LiR.
+        X: Input matrix for ridge fitting.
+        Y: Target matrix or tensor for ridge fitting.
+        runner: Runner dictionary with label, optional module, optional
+        model_kwargs, and fit_kwargs.
+        params: Case-level fit parameters such as alpha, n_feat, and chunk_size.
+    Output:
+        Tuple of elapsed seconds and fitted model.
+    This function constructs the requested FastL2LiR implementation, merges
+    case and runner fit arguments, runs fit, and prints the elapsed time.
+    """
+    label = runner["label"]
+    module = runner_module(fastl2lir_module, runner)
+    model_kwargs = dict(runner.get("model_kwargs", {}))
+    fit_kwargs = dict(params)
+    fit_kwargs.update(runner.get("fit_kwargs", {}))
+
+    print(f"  start {label} fit", flush=True)
+    start = time.perf_counter()
+    model = module.FastL2LiR(**model_kwargs).fit(X, Y, **fit_kwargs)
+    elapsed = time.perf_counter() - start
+    print(f"  done {label} fit: {elapsed:.4f} s", flush=True)
+    return elapsed, model
+
+
+def ordered_runners(runners, repeat_index):
+    """Return runner order for one repeat.
+
+    Inputs:
+        runners: Ordered list of runner dictionaries.
+        repeat_index: Zero-based repeat index.
+    Output:
+        Runners in forward order for even repeats and reverse order for odd
+        repeats.
+    This function reduces systematic timing bias from always running one
+    implementation first.
+    """
+    if repeat_index % 2 == 0:
+        return runners
+    return list(reversed(runners))
+
+
+def measure_case(fastl2lir_module, n, p, y_shape, params, repeats, runners):
+    """Measure all runners for one benchmark case.
+
+    Inputs:
+        fastl2lir_module: Default module that provides FastL2LiR.
+        n: Number of samples.
+        p: Number of input features or voxels.
+        y_shape: Target-unit shape after the sample axis.
+        params: Case-level fit parameters.
+        repeats: Number of repeated fits.
+        runners: Runner dictionaries to compare.
+    Output:
+        Tuple of per-label elapsed-time lists and maximum prediction difference
+        versus the first runner.
+    This function generates data, fits all runners, and checks prediction
+    agreement on the first samples for every repeat.
+    """
+    times_by_label = {runner["label"]: [] for runner in runners}
+    max_pred_diffs = []
+
+    for repeat_index in range(repeats):
+        print(f"  repeat {repeat_index + 1}/{repeats}")
+
+        # Generate fresh data for every repeat. Within a repeat, all runners use
+        # the same data so correctness and timing remain comparable.
+        X, Y = make_data(n, p, y_shape)
+
+        models = {}
+        for runner in ordered_runners(runners, repeat_index):
+            elapsed, model = fit_once(fastl2lir_module, X, Y, runner, params)
+            times_by_label[runner["label"]].append(elapsed)
+            models[runner["label"]] = model
+
+        pred_slice = X[: min(32, X.shape[0])]
+        baseline_label = runners[0]["label"]
+        baseline_pred = models[baseline_label].predict(pred_slice)
+        for runner in runners[1:]:
+            pred = models[runner["label"]].predict(pred_slice)
+            max_pred_diffs.append(float(np.max(np.abs(baseline_pred - pred))))
+
+    max_pred_diff = max(max_pred_diffs) if max_pred_diffs else 0.0
+    return times_by_label, max_pred_diff
+
+
+def median_by_label(times_by_label):
+    """Compute median runtime for each runner label.
+
+    Inputs:
+        times_by_label: Mapping from runner label to elapsed-time list.
+    Output:
+        Mapping from runner label to median elapsed seconds.
+    This function provides the primary speed-test statistic for noisy repeated
+    runs.
+    """
+    return {label: statistics.median(times) for label, times in times_by_label.items()}
+
+
+def speedup_values(runners, medians):
+    """Compute speedups relative to the first runner.
+
+    Inputs:
+        runners: Ordered list of runner dictionaries.
+        medians: Mapping from runner label to median elapsed seconds.
+    Output:
+        Mapping from non-baseline runner label to baseline_time / runner_time.
+    This function gives each comparison an explicit name when more than two
+    implementations are benchmarked.
+    """
+    if len(runners) < 2:
+        return {}
+    baseline = runners[0]["label"]
+    return {
+        runner["label"]: medians[baseline] / medians[runner["label"]]
+        for runner in runners[1:]
+    }
+
+
+def summary_header(runners):
+    """Build the CSV header for the benchmark summary.
+
+    Inputs:
+        runners: Ordered list of runner dictionaries.
+    Output:
+        Comma-separated CSV header string.
+    This function includes per-runner median time columns and named speedup
+    columns relative to the first runner.
+    """
+    time_columns = ",".join(f"{runner['label']}_median_s" for runner in runners)
+    columns = f"server,case,n,p,k,q,repeats,{time_columns}"
+    if len(runners) >= 2:
+        baseline = runners[0]["label"]
+        speedup_columns = ",".join(
+            f"{runner['label']}_speedup_vs_{baseline}" for runner in runners[1:]
+        )
+        columns += f",{speedup_columns}"
+    columns += ",max_pred_diff"
+    return columns
+
+
+def summary_line(server_name, row, runners):
+    """Build one CSV data line for a benchmark case.
+
+    Inputs:
+        server_name: Host or user-specified benchmark server label.
+        row: Case summary dictionary.
+        runners: Ordered list of runner dictionaries.
+    Output:
+        Comma-separated CSV row string.
+    This function serializes dimensions, medians, speedups, and prediction
+    agreement into the saved summary table.
+    """
+    base = [
+        server_name,
+        row["case"],
+        str(row["n"]),
+        str(row["p"]),
+        str(row["k"]),
+        str(row["q"]),
+        str(row["repeats"]),
+    ]
+    times = [f"{row['medians'][runner['label']]:.6f}" for runner in runners]
+    values = base + times
+    for runner in runners[1:]:
+        values.append(f"{row['speedups'][runner['label']]:.2f}")
+    values.append(f"{row['max_pred_diff']:.3e}")
+    return ",".join(values)
+
+
+def print_environment(
+    fastl2lir_module, module_names, n, p, k, repeats, cases, server_name
+):
+    """Print Python, package, and benchmark-size metadata.
+
+    Inputs:
+        fastl2lir_module: Default FastL2LiR package module.
+        module_names: Extra import names whose versions should be printed.
+        n: Number of samples.
+        p: Number of input features or voxels.
+        k: Number of selected features.
+        repeats: Repeat count.
+        cases: Benchmark case definitions.
+        server_name: Host or user-specified benchmark server label.
+    Output:
+        None.
+    This function reports enough runtime context to interpret saved speed-test
+    numbers later.
+    """
+    print("Environment")
+    print(f"  python: {sys.version.split()[0]}")
+    print(f"  numpy: {np.__version__}")
+    for module_name in module_names:
+        print(f"  {module_name}: {module_version(module_name)}")
+    print(f"  fastl2lir: {getattr(fastl2lir_module, '__version__', '(unknown)')}")
+    print(f"  fastl2lir file: {package_file(fastl2lir_module)}")
+    print(f"  server: {server_name}")
+    for env_name in (
+        "OPENBLAS_NUM_THREADS",
+        "MKL_NUM_THREADS",
+        "OMP_NUM_THREADS",
+        "NUMBA_NUM_THREADS",
+    ):
+        print(f"  {env_name} env: {os.environ.get(env_name, '(unset)')}")
+    try:
+        numba = importlib.import_module("numba")
+    except ImportError:
+        pass
+    else:
+        print(f"  numba.get_num_threads(): {numba.get_num_threads()}")
+    print(f"  benchmark size: n={n}, p={p}, n_feat={k}, repeats={repeats}")
+    print(f"  benchmark cases: {', '.join(name for name, *_ in cases)}")
+
+
+def print_intro(comparison, runners):
+    """Print the experiment description before benchmark cases run.
+
+    Inputs:
+        comparison: Human-readable comparison label.
+        runners: Ordered list of runner dictionaries.
+    Output:
+        None.
+    This function describes the primary metric, runner labels, data freshness,
+    and ordering policy.
+    """
+    print("\nSpeed test results")
+    print(f"  comparison: {comparison}")
+    print("  primary metric: median runtime")
+    print(f"  runners: {', '.join(runner['label'] for runner in runners)}")
+    print("  execution: runners are fit sequentially, one implementation at a time")
+    print("  tqdm: not patched")
+    print("  data: regenerated for every repeat")
+    if len(runners) >= 2:
+        print("  order: alternates between forward and reverse runner order")
+    print("")
+
+
+def print_runner_sources(fastl2lir_module, runners):
+    """Print module and argument details for each runner.
+
+    Inputs:
+        fastl2lir_module: Default module that provides FastL2LiR.
+        runners: Ordered list of runner dictionaries.
+    Output:
+        None.
+    This function records the module, model arguments, and fit arguments that
+    define each computation.
+    """
+    print("Runner sources")
+    for runner in runners:
+        module = runner_module(fastl2lir_module, runner)
+        model_kwargs = runner.get("model_kwargs", {})
+        fit_kwargs = runner.get("fit_kwargs", {})
+        print(f"  {runner['label']} module: {module.__name__}")
+        print(f"  {runner['label']} file: {package_file(module)}")
+        print(f"  {runner['label']} model_kwargs: {model_kwargs}")
+        print(f"  {runner['label']} fit_kwargs: {fit_kwargs}")
+    print("")
+
+
+def run_speedtest(
+    fastl2lir_module,
+    runners,
+    output_csv,
+    comparison,
+    module_names=(),
+):
+    """Run the full speed-test workflow and save the CSV summary.
+
+    Inputs:
+        fastl2lir_module: Default module that provides FastL2LiR.
+        runners: Ordered runner dictionaries.
+        output_csv: Path where the summary CSV will be saved.
+        comparison: Human-readable comparison label.
+        module_names: Extra import names whose versions should be printed.
+    Output:
+        None.
+    This function orchestrates environment reporting, case measurement,
+    human-readable summaries, CSV writing, and saved-path reporting.
+    """
+    runners = filter_runners(runners)
+    n, p, k, repeats = get_benchmark_settings()
+    cases = make_cases(k, repeats)
+    server_name = os.environ.get("BENCH_SERVER") or socket.gethostname()
+
+    print_environment(
+        fastl2lir_module, module_names, n, p, k, repeats, cases, server_name
+    )
+    print_intro(comparison, runners)
+    print_runner_sources(fastl2lir_module, runners)
+
+    rows = []
+    for name, y_shape, params, case_repeats in cases:
+        print(name)
+        q = int(np.prod(y_shape))
+        print(
+            f"  hyperparameters: n={n}, p={p}, k={params['n_feat']}, "
+            f"q={q}, repeats={case_repeats}"
+        )
+        times_by_label, max_pred_diff = measure_case(
+            fastl2lir_module, n, p, y_shape, params, case_repeats, runners
+        )
+        medians = median_by_label(times_by_label)
+        rows.append(
+            {
+                "case": name,
+                "xshape": (n, p),
+                "yshape": (n,) + y_shape,
+                "n": n,
+                "p": p,
+                "k": params["n_feat"],
+                "q": q,
+                "repeats": case_repeats,
+                "times_by_label": times_by_label,
+                "medians": medians,
+                "speedups": speedup_values(runners, medians),
+                "max_pred_diff": max_pred_diff,
+            }
+        )
+        print("")
+
+    print_case_summaries(rows, runners)
+    write_summary(output_csv, rows, runners, server_name)
+    print(f"\nSaved CSV: {output_csv}")
+
+
+def print_case_summaries(rows, runners):
+    """Print per-case timing and correctness summaries.
+
+    Inputs:
+        rows: Case summary dictionaries produced by run_speedtest.
+        runners: Ordered list of runner dictionaries.
+    Output:
+        None.
+    This function presents raw repeated times, median times, named speedups,
+    and maximum prediction differences for researcher inspection.
+    """
+    for row in rows:
+        print(row["case"])
+        print(f"  X={row['xshape']}, Y={row['yshape']}")
+        print(
+            f"  hyperparameters: n={row['n']}, p={row['p']}, "
+            f"k={row['k']}, q={row['q']}, repeats={row['repeats']}"
+        )
+        for runner in runners:
+            label = runner["label"]
+            rounded_times = [round(t, 4) for t in row["times_by_label"][label]]
+            print(f"  {label} times: {rounded_times} s")
+        median_text = ", ".join(
+            f"{runner['label']}={row['medians'][runner['label']]:.4f}s"
+            for runner in runners
+        )
+        for runner in runners[1:]:
+            baseline = runners[0]["label"]
+            label = runner["label"]
+            median_text += (
+                f", {label}_speedup_vs_{baseline}={row['speedups'][label]:.2f}x"
+            )
+        print(f"  median: {median_text}")
+        print(f"  max prediction diff on first samples: {row['max_pred_diff']:.3e}")
+        print("")
+
+
+def write_summary(output_csv, rows, runners, server_name):
+    """Write and print the CSV benchmark summary.
+
+    Inputs:
+        output_csv: Path where the summary CSV will be saved.
+        rows: Case summary dictionaries produced by run_speedtest.
+        runners: Ordered list of runner dictionaries.
+        server_name: Host or user-specified benchmark server label.
+    Output:
+        None.
+    This function creates the temp output directory when needed and writes the
+    same summary rows that are printed to the terminal.
+    """
+    output_csv = Path(output_csv)
+
+    print("Summary table")
+    print(summary_header(runners))
+    for row in rows:
+        print(summary_line(server_name, row, runners))
+
+    output_csv.parent.mkdir(parents=True, exist_ok=True)
+    with output_csv.open("w", encoding="utf-8") as f:
+        f.write(summary_header(runners) + "\n")
+        for row in rows:
+            f.write(summary_line(server_name, row, runners) + "\n")
diff --git a/src/fastl2lir/_numba.py b/src/fastl2lir/_numba.py
new file mode 100644
index 0000000..c3a2535
--- /dev/null
+++ b/src/fastl2lir/_numba.py
@@ -0,0 +1,157 @@
+"""Numba implementations for PyFastL2LiR."""
+
+import numpy as np
+from numba import get_num_threads, njit, prange, set_num_threads
+from threadpoolctl import threadpool_info
+
+
+_WARMED_UP_DTYPES = set()
+
+
+def _blas_backend_names():
+    """Return detected BLAS backend names from threadpoolctl."""
+    names = []
+    for info in threadpool_info():
+        if info.get("user_api") != "blas":
+            continue
+        internal_api = info.get("internal_api")
+        prefix = info.get("prefix")
+        filepath = info.get("filepath")
+        parts = [part for part in (internal_api, prefix, filepath) if part]
+        names.append(" ".join(parts))
+    return names
+
+
+def _uses_mkl():
+    """Return True when any detected BLAS backend is Intel MKL."""
+    names = _blas_backend_names()
+    return any("mkl" in name.lower() for name in names)
+
+
+def _load_blas_backend_for_detection(dtype=np.float64):
+    """Run a tiny BLAS call so threadpoolctl can see the active backend."""
+    x = np.ones((1, 1), dtype=dtype)
+    np.matmul(x, x)
+
+
+def check_numba_solver_environment(dtype=np.float64):
+    """Raise RuntimeError when solver="numba" is running with Intel MKL."""
+    _load_blas_backend_for_detection(dtype)
+    if not _uses_mkl():
+        return
+
+    detected = ", ".join(_blas_backend_names()) or "unknown"
+    raise RuntimeError(
+        "solver='numba' is not supported when the detected "
+        f"BLAS backend is Intel MKL. Detected backend: {detected}"
+    )
+
+
+def validate_numba_solver(numba_num_threads, dtype=np.float64):
+    """Validate optional numba solver dependencies, threads, and BLAS backend."""
+    if numba_num_threads is not None and numba_num_threads < 1:
+        raise ValueError("numba_num_threads must be a positive integer or None")
+
+    check_numba_solver_environment(dtype)
+
+
+def fit_selected_ridge_numba(X, C, W0, W1, n_feat, dtype, numba_num_threads=4):
+    """Fit selected-feature ridge regression with the optional numba kernel."""
+    validate_numba_solver(numba_num_threads, dtype)
+
+    previous_num_threads = None
+    if numba_num_threads is not None:
+        previous_num_threads = get_num_threads()
+        set_num_threads(min(numba_num_threads, previous_num_threads))
+
+    try:
+        warmup_numba(dtype)
+        W = np.zeros((C.shape[0], X.shape[1] - 1), dtype=dtype)
+        b = np.zeros((1, C.shape[0]), dtype=dtype)
+        return _fit_selected_ridge_numba(X, C, W, b, W0, W1, n_feat)
+    finally:
+        if previous_num_threads is not None:
+            set_num_threads(previous_num_threads)
+
+
+@njit(cache=True, parallel=True)
+def _fit_selected_ridge_numba(X, C, W, b, W0, W1, n_feat):
+    n_outputs = C.shape[0]
+    bias_index = X.shape[1] - 1
+
+    for index_outputDim in prange(n_outputs):
+        C0 = np.abs(C[index_outputDim, :])
+        feat_idx = np.argsort(C0)[::-1]
+        feat_idx = feat_idx[:n_feat]
+
+        I_with_bias = np.empty(feat_idx.size + 1, dtype=feat_idx.dtype)
+        for i in range(feat_idx.size):
+            I_with_bias[i] = feat_idx[i]
+        I_with_bias[feat_idx.size] = bias_index
+
+        W0_sub = np.zeros((I_with_bias.size, I_with_bias.size), dtype=W0.dtype)
+        for i in range(I_with_bias.size):
+            for j in range(I_with_bias.size):
+                W0_sub[i, j] = W0[I_with_bias[i], I_with_bias[j]]
+
+        rhs = np.zeros((I_with_bias.size, 1), dtype=W1.dtype)
+        for i in range(I_with_bias.size):
+            rhs[i, 0] = W1[index_outputDim, I_with_bias[i]]
+
+        Wb = np.linalg.solve(W0_sub.astype(np.float64), rhs.astype(np.float64)).astype(
+            W.dtype
+        )
+
+        for i in range(feat_idx.size):
+            W[index_outputDim, feat_idx[i]] = Wb[i, 0]
+        b[0, index_outputDim] = Wb[-1, 0]
+
+    return W.T, b
+
+
+def warmup_numba(dtype=np.float64):
+    """Compile or load the numba ridge kernel before fitting real data.
+
+    The feature-selection fit path calls ``_fit_selected_ridge_numba`` once
+    per fit. With ``cache=True``, numba can reuse compiled code across Python
+    processes, but the first call in a process still pays dispatcher/cache-load
+    and sometimes compilation overhead. If that first call happens on the real
+    training matrix, benchmark and notebook users see a large first-fit spike.
+
+    Warm up the same kernel once per dtype with tiny arrays so the one-time
+    cost is paid before the real solve loop. This does not remove numba's
+    startup cost, but it keeps that cost independent of the user's data size
+    and makes subsequent feature-selection fits reflect steady-state speed.
+    """
+    dtype = np.dtype(dtype)
+    if dtype in _WARMED_UP_DTYPES:
+        return
+
+    X = np.array(
+        [
+            [1.0, 0.0, 1.0],
+            [0.0, 1.0, 1.0],
+            [1.0, 1.0, 1.0],
+        ],
+        dtype=dtype,
+    )
+    C = np.array(
+        [
+            [0.3, 0.2],
+            [0.1, 0.4],
+        ],
+        dtype=dtype,
+    )
+    W = np.zeros((2, 2), dtype=dtype)
+    b = np.zeros((1, 2), dtype=dtype)
+    W0 = np.matmul(X.T, X) + np.eye(X.shape[1], dtype=dtype)
+    W1 = np.array(
+        [
+            [1.0, 0.5, 1.5],
+            [0.5, 1.0, 1.5],
+        ],
+        dtype=dtype,
+    )
+
+    _fit_selected_ridge_numba(X, C, W, b, W0, W1, 1)
+    _WARMED_UP_DTYPES.add(dtype)
diff --git a/src/fastl2lir/fastl2lir.py b/src/fastl2lir/fastl2lir.py
index decb9f6..87f3a18 100644
--- a/src/fastl2lir/fastl2lir.py
+++ b/src/fastl2lir/fastl2lir.py
@@ -8,6 +8,11 @@
 from threadpoolctl import threadpool_limits
 from tqdm import tqdm
 
+try:
+    from ._numba import fit_selected_ridge_numba
+except ImportError:
+    fit_selected_ridge_numba = None
+
 
 class FastL2LiR(object):
     """Fast L2-regularized linear regression class."""
@@ -53,6 +58,8 @@ def fit(
         chunk_size=0,
         cache_dir="./cache",
         dtype=np.float64,
+        solver="numpy",
+        numba_num_threads=4,
     ):
         """Fit the L2-regularized linear model with the given data.
 
@@ -80,6 +87,15 @@ def fit(
             because this is an operation for each unit.
             (The sample selection operation itself does not essentially need
             to record the selected voxel.)
+        solver: str ('numpy' or 'numba')
+            Solver used for fitting with feature selection.
+            'numpy' preserves the original implementation and is the default.
+            'numba' uses the experimental numba implementation only when
+            feature selection is enabled and save_select_feat is False.
+        numba_num_threads: int or None
+            Maximum number of numba threads used while solver='numba'. The
+            previous numba thread count is restored after fitting. Set to None
+            to leave the current numba thread count unchanged.
 
         Returns
         -------
@@ -148,6 +164,8 @@ def fit(
                         n_feat=n_feat,
                         use_all_features=no_feature_selection,
                         dtype=dtype,
+                        solver=solver,
+                        numba_num_threads=numba_num_threads,
                     )
                 w_list.append(W)
                 b_list.append(b)
@@ -179,6 +197,8 @@ def fit(
                     n_feat=n_feat,
                     use_all_features=no_feature_selection,
                     dtype=dtype,
+                    solver=solver,
+                    numba_num_threads=numba_num_threads,
                 )
 
         self.__W = W
@@ -262,7 +282,15 @@ def predict(self, X, dtype=np.float64, save_select_feat=False, spatial_norm=None
         return Y
 
     def __sub_fit(
-        self, X, Y, alpha=0, n_feat=0, use_all_features=True, dtype=np.float64
+        self,
+        X,
+        Y,
+        alpha=0,
+        n_feat=0,
+        use_all_features=True,
+        dtype=np.float64,
+        solver="numpy",
+        numba_num_threads=4,
     ):
         if use_all_features:
             # Without feature selection
@@ -298,24 +326,33 @@ def __sub_fit(
             W1 = np.matmul(Y.T, X)
             C = C.T
 
-            with threadpool_limits(limits=1, user_api="blas"):
-                for index_outputDim in tqdm(range(Y.shape[1])):
-                    C0 = abs(C[index_outputDim, :])
-                    feat_idx = np.argsort(C0)
-                    feat_idx = feat_idx[::-1]
-                    feat_idx = feat_idx[0:n_feat]
-                    feat_idx = np.hstack((feat_idx, X.shape[1] - 1))
-                    W0_sub = (
-                        W0.ravel()[
-                            (
-                                feat_idx + (feat_idx * W0.shape[1]).reshape((-1, 1))
-                            ).ravel()
-                        ]
-                    ).reshape(feat_idx.size, feat_idx.size)
-                    Wb = np.linalg.solve(W0_sub, W1[index_outputDim, feat_idx])
-                    W[index_outputDim, feat_idx[:-1]] = Wb[:-1]
-                    b[0, index_outputDim] = Wb[-1]
-                W = W.T
+            if solver == "numpy":
+                with threadpool_limits(limits=1, user_api="blas"):
+                    for index_outputDim in tqdm(range(Y.shape[1])):
+                        C0 = abs(C[index_outputDim, :])
+                        feat_idx = np.argsort(C0)
+                        feat_idx = feat_idx[::-1]
+                        feat_idx = feat_idx[0:n_feat]
+                        feat_idx = np.hstack((feat_idx, X.shape[1] - 1))
+                        W0_sub = (
+                            W0.ravel()[
+                                (
+                                    feat_idx + (feat_idx * W0.shape[1]).reshape((-1, 1))
+                                ).ravel()
+                            ]
+                        ).reshape(feat_idx.size, feat_idx.size)
+                        Wb = np.linalg.solve(W0_sub, W1[index_outputDim, feat_idx])
+                        W[index_outputDim, feat_idx[:-1]] = Wb[:-1]
+                        b[0, index_outputDim] = Wb[-1]
+                    W = W.T
+            elif solver == "numba":
+                if fit_selected_ridge_numba is None:
+                    raise ImportError("solver='numba' requires numba to be installed")
+                W, b = fit_selected_ridge_numba(
+                    X, C, W0, W1, n_feat, dtype, numba_num_threads
+                )
+            else:
+                raise ValueError("Unknown solver specified:", solver)
 
         return W, b
 
diff --git a/tests/test__numba.py b/tests/test__numba.py
new file mode 100644
index 0000000..a6aba9f
--- /dev/null
+++ b/tests/test__numba.py
@@ -0,0 +1,122 @@
+"""Tests for the optional numba solver."""
+
+from unittest import TestCase
+from unittest.mock import patch
+
+import numpy as np
+
+import fastl2lir
+import fastl2lir._numba as fastl2lir_numba
+
+
+class TestNumbaSolver(TestCase):
+    """Tests for optional numba feature-selection fitting."""
+
+    def skip_if_mkl_backend(self):
+        """Skip solver execution tests when the current BLAS backend is MKL."""
+        fastl2lir_numba._load_blas_backend_for_detection()
+        if fastl2lir_numba._uses_mkl():
+            self.skipTest("numba solver is not supported with Intel MKL")
+
+    def test_numba_solver_matches_numpy_solver(self):
+        """Test that the numba feature-selection solver matches numpy."""
+
+        try:
+            import numba  # noqa: F401
+        except ImportError:
+            self.skipTest("numba is not installed")
+        self.skip_if_mkl_backend()
+
+        rng = np.random.default_rng(0)
+        X = rng.normal(size=(24, 12))
+        Y = rng.normal(size=(24, 5))
+        X_test = rng.normal(size=(7, 12))
+
+        model_numpy = fastl2lir.FastL2LiR()
+        model_numba = fastl2lir.FastL2LiR()
+
+        model_numpy.fit(X, Y, alpha=0.5, n_feat=4, solver="numpy")
+        model_numba.fit(X, Y, alpha=0.5, n_feat=4, solver="numba")
+
+        np.testing.assert_allclose(model_numba.W, model_numpy.W)
+        np.testing.assert_allclose(model_numba.b, model_numpy.b)
+        np.testing.assert_allclose(
+            model_numba.predict(X_test), model_numpy.predict(X_test)
+        )
+
+    def test_numba_solver_restores_thread_count(self):
+        """Test that the numba solver does not leak thread-count changes."""
+
+        try:
+            from numba import get_num_threads, set_num_threads
+        except ImportError:
+            self.skipTest("numba is not installed")
+        self.skip_if_mkl_backend()
+
+        rng = np.random.default_rng(1)
+        X = rng.normal(size=(24, 12))
+        Y = rng.normal(size=(24, 5))
+        original_num_threads = get_num_threads()
+
+        try:
+            fastl2lir.FastL2LiR().fit(X, Y, n_feat=4, solver="numba")
+            self.assertEqual(get_num_threads(), original_num_threads)
+
+            set_num_threads(1)
+            fastl2lir.FastL2LiR().fit(
+                X, Y, n_feat=4, solver="numba", numba_num_threads=None
+            )
+            self.assertEqual(get_num_threads(), 1)
+        finally:
+            set_num_threads(original_num_threads)
+
+    def test_numba_solver_rejects_mkl_backend(self):
+        """Test that numba fitting is rejected with an Intel MKL backend."""
+
+        try:
+            import numba  # noqa: F401
+        except ImportError:
+            self.skipTest("numba is not installed")
+
+        rng = np.random.default_rng(2)
+        X = rng.normal(size=(24, 12))
+        Y = rng.normal(size=(24, 5))
+
+        with patch.object(
+            fastl2lir_numba,
+            "threadpool_info",
+            return_value=[
+                {
+                    "user_api": "blas",
+                    "internal_api": "mkl",
+                    "prefix": "libmkl_rt",
+                    "filepath": "/example/libmkl_rt.so",
+                }
+            ],
+        ):
+            with self.assertRaisesRegex(RuntimeError, "Intel MKL"):
+                fastl2lir.FastL2LiR().fit(X, Y, n_feat=4, solver="numba")
+
+    def test_numba_solver_allows_non_mkl_backend(self):
+        """Test that numba environment validation allows non-MKL BLAS."""
+
+        with (
+            patch.object(
+                fastl2lir_numba,
+                "_load_blas_backend_for_detection",
+                return_value=None,
+            ),
+            patch.object(
+                fastl2lir_numba,
+                "threadpool_info",
+                return_value=[
+                    {
+                        "user_api": "blas",
+                        "internal_api": "blis",
+                        "prefix": "libblis",
+                        "filepath": "/example/libblis.so",
+                    }
+                ],
+            ),
+        ):
+            fastl2lir_numba.check_numba_solver_environment()
diff --git a/tests/test_fastl2lir.py b/tests/test_fastl2lir.py
index cd144e3..cfb6c83 100644
--- a/tests/test_fastl2lir.py
+++ b/tests/test_fastl2lir.py
@@ -154,12 +154,11 @@ def test_reshape(self):
         np.testing.assert_array_equal(model_test.b.shape, (1,) + Y_shape[1:])
         np.testing.assert_array_equal(pred_test.shape, Y_shape)
 
-
     def test_save_select_feat_default_select_sample(self):
-        '''save_select_feat=True with default select_sample=None should not raise.'''
-        data = np.load('./tests/testdata_nfeat.npz')
+        """save_select_feat=True with default select_sample=None should not raise."""
+        data = np.load("./tests/testdata_nfeat.npz")
         model = fastl2lir.FastL2LiR()
-        model.fit(data['x_tr'], data['y_1d'], n_feat=20, save_select_feat=True)
+        model.fit(data["x_tr"], data["y_1d"], n_feat=20, save_select_feat=True)
 
 
 if __name__ == "__main__":

From 2b69f81319365ef1722144a145799df61c3a7dc9 Mon Sep 17 00:00:00 2001
From: Hideki <izumi.hideki.36w@st.kyoto-u.ac.jp>
Date: Tue, 9 Jun 2026 17:20:52 +0900
Subject: [PATCH 2/6] docs: update numba solver documentation to clarify dtype
 behavior

---
 src/fastl2lir/fastl2lir.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/fastl2lir/fastl2lir.py b/src/fastl2lir/fastl2lir.py
index 87f3a18..9b8ae6b 100644
--- a/src/fastl2lir/fastl2lir.py
+++ b/src/fastl2lir/fastl2lir.py
@@ -92,6 +92,11 @@ def fit(
             'numpy' preserves the original implementation and is the default.
             'numba' uses the experimental numba implementation only when
             feature selection is enabled and save_select_feat is False.
+            Note: the numba solver always performs the per-unit linear solve
+            in float64 and casts the result back to ``dtype``. With
+            ``dtype=np.float32`` its weights/bias therefore differ slightly
+            from the numpy solver (which solves in float32); the two agree
+            only up to float32 precision, not bit-for-bit.
         numba_num_threads: int or None
             Maximum number of numba threads used while solver='numba'. The
             previous numba thread count is restored after fitting. Set to None

From 197cfa3ea701108ac937189067f8995d8fffecd9 Mon Sep 17 00:00:00 2001
From: Hideki <izumi.hideki.36w@st.kyoto-u.ac.jp>
Date: Tue, 9 Jun 2026 17:26:16 +0900
Subject: [PATCH 3/6] fix: use mergesort for stable feature index sorting in
 numba solver

---
 src/fastl2lir/_numba.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/fastl2lir/_numba.py b/src/fastl2lir/_numba.py
index c3a2535..115d1d9 100644
--- a/src/fastl2lir/_numba.py
+++ b/src/fastl2lir/_numba.py
@@ -81,7 +81,7 @@ def _fit_selected_ridge_numba(X, C, W, b, W0, W1, n_feat):
 
     for index_outputDim in prange(n_outputs):
         C0 = np.abs(C[index_outputDim, :])
-        feat_idx = np.argsort(C0)[::-1]
+        feat_idx = np.argsort(C0, kind="mergesort")[::-1]
         feat_idx = feat_idx[:n_feat]
 
         I_with_bias = np.empty(feat_idx.size + 1, dtype=feat_idx.dtype)

From 490d86f6ca7f692cde7ebcc2bb2ddda9e709ea05 Mon Sep 17 00:00:00 2001
From: Hideki <izumi.hideki.36w@st.kyoto-u.ac.jp>
Date: Tue, 9 Jun 2026 17:26:23 +0900
Subject: [PATCH 4/6] fix: validate solver input and improve numba feature
 selection handling

---
 src/fastl2lir/fastl2lir.py | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/fastl2lir/fastl2lir.py b/src/fastl2lir/fastl2lir.py
index 9b8ae6b..df27a65 100644
--- a/src/fastl2lir/fastl2lir.py
+++ b/src/fastl2lir/fastl2lir.py
@@ -107,6 +107,8 @@ def fit(
         self
             Returns an instance of self.
         """
+        if solver not in ("numpy", "numba"):
+            raise ValueError("solver must be 'numpy' or 'numba'")
 
         if X.dtype != dtype:
             X = X.astype(dtype)
@@ -137,6 +139,13 @@ def fit(
             if (spatial_norm is not None) or (select_sample is not None):
                 save_select_feat = True
 
+        if solver == "numba" and (save_select_feat or no_feature_selection):
+            warnings.warn(
+                "solver='numba' is only used for feature-selection fitting "
+                "when save_select_feat is False. Falling back to the numpy "
+                "fit path for this call."
+            )
+
         # Chunking
         if chunk_size > 0:
             chunks = self.__get_chunks(range(Y.shape[1]), chunk_size)
@@ -335,7 +344,7 @@ def __sub_fit(
                 with threadpool_limits(limits=1, user_api="blas"):
                     for index_outputDim in tqdm(range(Y.shape[1])):
                         C0 = abs(C[index_outputDim, :])
-                        feat_idx = np.argsort(C0)
+                        feat_idx = np.argsort(C0, kind="mergesort")
                         feat_idx = feat_idx[::-1]
                         feat_idx = feat_idx[0:n_feat]
                         feat_idx = np.hstack((feat_idx, X.shape[1] - 1))

From d9912bc6081213fd47fa46a72299c19090f1a5f6 Mon Sep 17 00:00:00 2001
From: Hideki <izumi.hideki.36w@st.kyoto-u.ac.jp>
Date: Tue, 9 Jun 2026 17:26:33 +0900
Subject: [PATCH 5/6] test: add comprehensive tests for numba solver
 functionality and validation

---
 tests/test__numba.py | 85 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 85 insertions(+)

diff --git a/tests/test__numba.py b/tests/test__numba.py
index a6aba9f..6a65e18 100644
--- a/tests/test__numba.py
+++ b/tests/test__numba.py
@@ -44,6 +44,91 @@ def test_numba_solver_matches_numpy_solver(self):
             model_numba.predict(X_test), model_numpy.predict(X_test)
         )
 
+    def test_numba_solver_matches_numpy_solver_float32(self):
+        """Test that float32 numba fitting agrees with numpy within float32 tolerance."""
+
+        try:
+            import numba  # noqa: F401
+        except ImportError:
+            self.skipTest("numba is not installed")
+        self.skip_if_mkl_backend()
+
+        rng = np.random.default_rng(3)
+        X = rng.normal(size=(32, 14)).astype(np.float32)
+        Y = rng.normal(size=(32, 6)).astype(np.float32)
+        X_test = rng.normal(size=(9, 14)).astype(np.float32)
+
+        model_numpy = fastl2lir.FastL2LiR()
+        model_numba = fastl2lir.FastL2LiR()
+
+        model_numpy.fit(X, Y, alpha=0.5, n_feat=5, dtype=np.float32, solver="numpy")
+        model_numba.fit(X, Y, alpha=0.5, n_feat=5, dtype=np.float32, solver="numba")
+
+        self.assertEqual(model_numpy.W.dtype, np.float32)
+        self.assertEqual(model_numba.W.dtype, np.float32)
+        np.testing.assert_allclose(model_numba.W, model_numpy.W, rtol=1e-4, atol=1e-5)
+        np.testing.assert_allclose(model_numba.b, model_numpy.b, rtol=1e-4, atol=1e-5)
+        np.testing.assert_allclose(
+            model_numba.predict(X_test, dtype=np.float32),
+            model_numpy.predict(X_test, dtype=np.float32),
+            rtol=1e-4,
+            atol=1e-5,
+        )
+
+    def test_numba_solver_matches_numpy_solver_with_tied_correlations(self):
+        """Test stable feature selection when correlations tie at the cutoff."""
+
+        try:
+            import numba  # noqa: F401
+        except ImportError:
+            self.skipTest("numba is not installed")
+        self.skip_if_mkl_backend()
+
+        h1 = np.array([1, -1, 1, -1, 1, -1, 1, -1.0])
+        h2 = np.array([1, 1, -1, -1, 1, 1, -1, -1.0])
+        h3 = np.array([1, 1, 1, 1, -1, -1, -1, -1.0])
+        h4 = np.array([1, -1, -1, 1, 1, -1, -1, 1.0])
+        X = np.column_stack([h1, h2, h3, h4])
+        Y = np.column_stack([2 * h1 + h2 + h3, -h1 + h2 + h3])
+        X_test = X[[0, 2, 4], :]
+
+        model_numpy = fastl2lir.FastL2LiR()
+        model_numba = fastl2lir.FastL2LiR()
+
+        model_numpy.fit(X, Y, alpha=0.5, n_feat=2, solver="numpy")
+        model_numba.fit(X, Y, alpha=0.5, n_feat=2, solver="numba")
+
+        np.testing.assert_allclose(model_numba.W, model_numpy.W)
+        np.testing.assert_allclose(model_numba.b, model_numpy.b)
+        np.testing.assert_allclose(
+            model_numba.predict(X_test), model_numpy.predict(X_test)
+        )
+
+    def test_solver_is_validated_before_fit_path_selection(self):
+        """Test that invalid solver names are rejected even without feature selection."""
+
+        rng = np.random.default_rng(4)
+        X = rng.normal(size=(16, 5))
+        Y = rng.normal(size=(16, 3))
+
+        with self.assertRaisesRegex(ValueError, "solver must be"):
+            fastl2lir.FastL2LiR().fit(X, Y, n_feat=0, solver="typo")
+
+    def test_numba_solver_warns_when_falling_back_to_numpy_fit_path(self):
+        """Test that numba requests warn when the numba path is not used."""
+
+        rng = np.random.default_rng(5)
+        X = rng.normal(size=(16, 5))
+        Y = rng.normal(size=(16, 3))
+
+        with self.assertWarnsRegex(UserWarning, "solver='numba' is only used"):
+            fastl2lir.FastL2LiR().fit(X, Y, n_feat=0, solver="numba")
+
+        with self.assertWarnsRegex(UserWarning, "solver='numba' is only used"):
+            fastl2lir.FastL2LiR().fit(
+                X, Y, n_feat=3, solver="numba", save_select_feat=True
+            )
+
     def test_numba_solver_restores_thread_count(self):
         """Test that the numba solver does not leak thread-count changes."""
 

From b776b981c4ad720f52a959851f6497b278946f0f Mon Sep 17 00:00:00 2001
From: Hideki Izumi <hidekizumi@gmail.com>
Date: Tue, 9 Jun 2026 17:58:10 +0900
Subject: [PATCH 6/6] Potential fix for pull request finding

fix my mistakes...

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
---
 speedtests/README.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/speedtests/README.md b/speedtests/README.md
index 95d7c2d..48b3ec0 100644
--- a/speedtests/README.md
+++ b/speedtests/README.md
@@ -66,8 +66,8 @@ BENCH_N=24 BENCH_P=12 BENCH_K=4 BENCH_REPEATS=1 BENCH_CASES=fc8 ./speedtests/spe
 - `BENCH_P`: Number of input features or voxels. Default: `15000`.
 - `BENCH_K`: Number of selected features. Default: `500`.
 - `BENCH_REPEATS`: Repeats per benchmark case. Default: `3`.
-- `BENCH_CASES`: Comma-separated subset of `fc8`, `fc6`, `conv5`, and
-  `conv5_chunk10`.
+- `BENCH_CASES`: Comma-separated subset of `fc8`, `fc6`, `conv5_chunk5`, and
+  `conv4_chunk5`.
 - `BENCH_RUNNERS`: Comma-separated subset of available runner labels. The
   current speed run supports `numpy` and `numba`; the legacy speed run supports
   `legacy`.