hand-rolling ols regression in rust, and the collinearity trap that ate my coefficients

lights out, again.

last time out i argued that a chart is just rectangles and lines, drew stintlab’s entire ui in rust on canvas2d to prove it, and promised — right at the bottom of that post — to come back for the degradation model itself: the math the pictures are pictures of. this is that post.

the pictures draw a number. how much slower a tire gets, lap after lap, as it wears out. pin that number down and you can start to guess when a car should box. so i tried to fit it with a multiple linear regression i wrote by hand — no nalgebra, no ndarray, just the normal equations and a determinant — and the thing handed me back confident, well-formed garbage. it took a comment in one of my own tests to explain why.

stintlab ↗ is a personal side project: i pull historical formula 1 data off a free public api ↗ , fit tire-degradation models in rust, and render the results in the browser. no users, no deadline, just me learning what regression actually feels like when you build it from the arithmetic up instead of calling .fit() on something someone else wrote. it’s the same paper-to-shipping-tool grind i wrote about in the drain post , except this time the “paper” is a first-year stats identity and the shipping tool is a pit-wall calculator that, for one glorious afternoon, lied to me with a straight face.

the model is two slopes fighting

here’s the whole physical claim, compressed into one line of doc comment:

        rust
        
    
/// Statistical degradation model fitted per (compound, circuit) pair.
///
/// Model: `lap_time = base + slope * tire_age - fuel_correction * (laps_total - lap_number)`
///
/// Computed in Rust at ingest time via OLS linear regression.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DegradationModel {
    pub compound: Compound,
    pub circuit_key: String,
    pub base_lap_time_ms: u32,
    pub degradation_slope_ms: f64,
    pub fuel_correction_ms: f64,
    pub r_squared: f64,
    pub sample_count: u32,
}

read it as two slopes pulling in opposite directions. a tire gets slower as it ages, so degradation_slope_ms is the milliseconds-per-lap penalty for tire_age. meanwhile the car burns fuel and gets lighter, so it gets faster as the race goes on — that’s the fuel_correction_ms term riding on laps_total - lap_number, the laps of fuel still onboard. a lap time is a base pace plus tire punishment minus fuel relief. two predictors, one intercept. that’s a textbook multiple linear regression, and i wanted to fit it without pulling a linear-algebra crate into stintlab-core.

ordinary least squares, no libraries

ordinary least squares has a closed form, and it’s short enough to type. you stack every lap into a design matrix X — a column of ones for the intercept, a column of tire ages, a column of fuel-remaining — stack the lap times into y, and the best-fit coefficients beta are the solution to the normal equations ↗ , (Xᵀ X) beta = Xᵀ y. that’s it. build two small things, solve one linear system.

        rust
        
    
// Build OLS matrices: y = X * beta
// X columns: [1 (intercept), tire_age, fuel_remaining]
// where fuel_remaining = laps_total - lap_number
let n = filtered.len();
let mut y = vec![0.0_f64; n];
let mut x = vec![[0.0_f64; 3]; n];

for (i, lap) in filtered.iter().enumerate() {
    y[i] = f64::from(lap.lap_time_ms.unwrap_or(0));
    x[i][0] = 1.0; // intercept
    x[i][1] = f64::from(lap.tire_age); // tire age
    x[i][2] = f64::from(laps_total.saturating_sub(lap.lap_number)); // fuel remaining
}

// Solve normal equations: (X^T X) beta = X^T y
let beta = solve_normal_equations_3x3(&x, &y)?;

three columns means Xᵀ X is a 3×3 matrix and Xᵀ y is a length-3 vector, no matter how many thousand laps you feed in. a fixed, tiny system. which is the entire reason i felt licensed to solve it the way i’m about to.

solving 3×3 by hand: cramer’s rule

for a system this small you don’t need gaussian elimination and you certainly don’t need a crate. cramer’s rule ↗ says each coefficient is a ratio of determinants: replace column k of the matrix with the right-hand side, take its determinant, divide by the determinant of the original. so i build Xᵀ X and Xᵀ y by hand and turn the crank:

        rust
        
    
/// Solve the 3x3 normal equations system `(X^T X) beta = X^T y` using Cramer's rule.
///
/// Returns `None` if the system is singular (determinant ~ 0).
fn solve_normal_equations_3x3(x: &[[f64; 3]], y: &[f64]) -> Option<[f64; 3]> {
    // Compute X^T X (3x3 symmetric matrix)
    let mut xtx = [[0.0_f64; 3]; 3];
    for i in 0..3 {
        for j in 0..3 {
            xtx[i][j] = x.iter().map(|row| row[i] * row[j]).sum();
        }
    }

    // Compute X^T y (3x1 vector)
    let mut xty = [0.0_f64; 3];
    for i in 0..3 {
        xty[i] = x.iter().zip(y.iter()).map(|(row, &yi)| row[i] * yi).sum();
    }

    // Solve via Cramer's rule for 3x3
    let det = det3x3(&xtx);
    if det.abs() < 1e-12 {
        return None;
    }

    let mut result = [0.0_f64; 3];
    for col in 0..3 {
        let mut m = xtx;
        for (row_idx, row) in m.iter_mut().enumerate() {
            row[col] = xty[row_idx];
        }
        result[col] = det3x3(&m) / det;
    }

    Some(result)
}

the determinant itself is the cofactor expansion you memorised and then forgot:

        rust
        
    
/// Compute the determinant of a 3x3 matrix.
fn det3x3(m: &[[f64; 3]; 3]) -> f64 {
    m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1])
        - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0])
        + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])
}

stare at the guard in the solver, because it’s the hero of this whole post: if det.abs() < 1e-12 { return None; }. cramer’s rule divides by the determinant. when the determinant is zero the system has no unique solution, and dividing by a number a hair off zero doesn’t give you the answer — it gives you a wildly amplified one, a coefficient blown up by numerical noise. the guard says: if the matrix is singular, don’t pretend. return None. i wrote that line as defensive boilerplate, the kind of thing you add because dividing by zero is embarrassing. i did not yet understand that it was about to catch a real, structural bug in my data — and that catching it was the point.

the trap: two slopes the data can’t tell apart

i wrote the fit. i wrote a test that generated fake laps from a known model and checked that the regression recovered the known slopes. it didn’t. sometimes it returned None; when it did fit, the two slopes came out as nonsense that summed to something plausible but individually meant nothing. the fit was “good” — high r-squared — and completely wrong.

the comment i left for future-me when i finally understood is the most important thing in the repo:

        rust
        
    
/// Generate synthetic lap data with known degradation characteristics.
///
/// Simulates multiple drivers with different stint start laps to break
/// the collinearity between `tire_age` and `laps_total - lap_number`.
/// Without this, OLS cannot separate the two effects (perfect linear dependence).
fn make_test_laps(/* ... */) -> Vec<Lap> {
    // Generate data from 3 "drivers" with different stint starts
    // to break the collinearity between tire_age and fuel_remaining.
    let stint_configs: &[(u16, &str)] = &[
        (1, "DR1"),  // stint starting lap 1
        (15, "DR2"), // stint starting lap 15
        (25, "DR3"), // stint starting lap 25
    ];
    // ...
}

here is the trap, and it’s pure algebra. take a single stint that starts on lap s. on any lap in that stint, tire_age = lap_number - s + 1 and fuel_remaining = laps_total - lap_number. add them:

text

tire_age + fuel_remaining = laps_total - s + 1

a constant. for one stint, fuel_remaining is just (a constant) - tire_age. my two predictor columns aren’t independent — one is an exact affine function of the other, and both live in the span of the intercept column. the design matrix loses a rank, Xᵀ X goes singular, its determinant collapses to zero, and the guard fires. this is textbook perfect multicollinearity ↗ : when predictors have an exact linear relationship the moment matrix can’t be inverted, and the coefficients aren’t just hard to estimate, they’re not identifiable — infinitely many (slope, fuel_correction) pairs fit the data equally well. the regression wasn’t broken. the question i was asking it was unanswerable.

and now the fix reads itself. three drivers, three different stint starts, three different constants laps_total - s + 1. across the pooled data fuel_remaining is no longer one global line in tire_age — a car that bolted on fresh tires at lap 25 has old rubber and a light fuel load at the same time a lap-1 starter has old rubber and no fuel left. that decoupling is the only thing that lets ols pull the two slopes apart. the stint_configs array isn’t test scaffolding; it’s the experimental design that makes the coefficients exist at all. real f1 data gives you this for free, because drivers pit on different laps. my first synthetic generator, one tidy stint, did not — and so it asked an impossible question and the matrix, correctly, refused to answer.

with varied stint starts, the same fit recovers the planted parameters cleanly:

        rust
        
    
#[test]
fn fit_model_recovers_known_parameters() {
    // 60 laps across 3 drivers (20 each) with varying stint starts
    let laps = make_test_laps(60, 90_000, 60.0, 30.0, 57, Compound::Medium);
    let model = fit_model("test_circuit", Compound::Medium, &laps, 57)
        .expect("should fit with sufficient multi-driver data");
    // ...
    assert!(
        (model.degradation_slope_ms - 60.0).abs() < 5.0,
        "slope {}, expected ~60",
        model.degradation_slope_ms
    );

    assert!(
        (model.fuel_correction_ms - 30.0).abs() < 5.0,
        "fuel_correction {}, expected ~30",
        model.fuel_correction_ms
    );

    assert!(
        model.r_squared > 0.95,
        "r_squared {} should be > 0.95",
        model.r_squared
    );
}

before you even get to fit: filtering

the fit is only as good as the laps you hand it, so there’s a gate in front. a hard minimum sample count, and a three-sigma outlier pass to throw out safety-car laps, traffic, and the lap someone spun:

        rust
        
    
/// Minimum number of valid laps required to fit a degradation model.
const MIN_SAMPLES: usize = 30;

        rust
        
    
let n_f = times.len() as f64;
let mean = times.iter().sum::<f64>() / n_f;
let variance = times.iter().map(|t| (t - mean).powi(2)).sum::<f64>() / n_f;
let std_dev = variance.sqrt();

// Exclude outliers (>3 standard deviations from mean)
let filtered: Vec<&Lap> = valid
    .into_iter()
    .filter(|l| {
        let t = f64::from(l.lap_time_ms.unwrap_or(0));
        (t - mean).abs() <= 3.0 * std_dev
    })
    .collect();

and after the fit, r-squared — the fraction of lap-time variance the model explains — as a one-number sanity read on whether the line is worth anything:

        rust
        
    
let y_mean = y.iter().sum::<f64>() / n as f64;
let ss_tot: f64 = y.iter().map(|yi| (yi - y_mean).powi(2)).sum();
let ss_res: f64 = y
    .iter()
    .enumerate()
    .map(|(i, yi)| {
        let y_hat = beta[0] + beta[1] * x[i][1] + beta[2] * x[i][2];
        (yi - y_hat).powi(2)
    })
    .sum();

let r_squared = if ss_tot > 0.0 {
    1.0 - ss_res / ss_tot
} else {
    0.0
};

(all those as f64 casts on sample counts are exactly what clippy’s cast_precision_loss ↗ pedantic lint yells about — a u64 past 2⁵³ loses bits in an f64. i suppress it locally with a note, because a lap count is never getting anywhere near that, and a blanket #[allow] would silence the one place it might someday matter.)

the honest part: what i cut, and what i’m side-eyeing

the blog rule around here is that i don’t get to pretend. so:

“confidence” is just r-squared wearing a lanyard. the pit-window predictor reports a confidence number, and here’s where it comes from:

        rust
        
    
// Confidence derived from R-squared of the models used
let confidence = best_r_squared_min;

that’s the smaller r-squared of the two models involved in a stop. which means my “confidence” is goodness-of-fit on the data i already saw, not a calibrated interval that says anything about a lap that hasn’t happened. a model can fit its own history at r-squared 0.97 and still be a terrible prophet for lap 40. i warned in the last post that r-squared was doing more emotional labor in this predictor than it should, and this is the line where it clocks in. it’s a plausible-looking number standing in for real predictive uncertainty, and it is not the same thing.

the fuel sign convention is a hypothesis, not a law. look at how i pull the fuel term out:

        rust
        
    
let degradation_slope_ms = beta[1];
// fuel_correction is negated: model is base + slope*age - fuel_corr*remaining
// OLS fits: y = b0 + b1*age + b2*remaining, so fuel_correction = -b2
let fuel_correction_ms = -beta[2];

the name fuel_correction_ms, positive, subtracted, encodes a physical story: more fuel onboard, faster the car sheds it, laps get quicker over a stint. but ols is free to fit whatever sign the data supports, and on real telemetry — where a heavy car is a slow car — a correctly-signed fuel effect could come out with the opposite sign to what the variable’s name promises. the regression will happily hand me a coefficient either way; nothing in the math checks it against physics. i haven’t validated the sign against a full season yet, and until i do, that minus sign is an assertion i should treat as a bug candidate, not a fact.

cramer’s rule is a party trick that doesn’t scale. it’s genuinely fine for 3×3 — fast, exact, no dependencies. but cofactor-expansion determinants blow up factorially and get numerically shaky as the matrix grows, so the moment i want a fourth predictor (track temp? a tire-warmup term?) this whole approach is the wrong tool and i rewrite the solver around a proper decomposition. i’m not shipping cramer’s rule past three columns, and i knew that when i wrote it.

none of these are secret. the guard turned the identifiability bug into an honest None instead of a confident wrong answer, and that’s the pattern i’d defend: make the code refuse rather than lie.

lessons learned

a singular matrix is the honest error message. i wrote det.abs() < 1e-12 -> None as defensive noise and it turned out to be the load-bearing line — it converted “these two slopes are mathematically inseparable” from silent garbage into a None i had to deal with. when the math can’t answer, make the code say so.
collinearity is a property of your design, not your solver. no amount of clever numerics separates tire_age from fuel_remaining inside a single stint, because inside a single stint they are the same variable up to a constant. the fix lived in the data — three drivers, three stint starts — not the regression.
name a thing and you’ve made a claim. calling a coefficient fuel_correction asserts a sign and a story the regression never agreed to. write the units, write the sign convention, and go check it against reality before you trust it.

next up: the pit-window strategy engine — a brute-force search over every candidate pit lap and every compound, why “just try them all” is the right call at this scale, and my ongoing argument with that r-squared-as-confidence number, which deserves a real prediction interval and doesn’t have one yet.

the code’s on github ↗ , singular matrices and unverified sign conventions and all. it’s a work in progress, like everything here.