Skip to content
Computational Lab.

Least Squares Curve Fitting

One exact solve, no descent required.

Beginner / linear-algebra / approximation / statistics

The least-squares line has a closed form: solve the 2x2 normal equations exactly by Cramer's rule and you are done - no iterations, no learning rate. The residual stems show what the fit absorbs and what it cannot, and the outlier slider demonstrates exactly what unweighted exactness costs.

Tagged: normal equations, OLS, outliers, R-squared

Live simulation

Raise the outlier rate and watch unweighted exactness pay for it.

Speed

Raise the outlier rate and watch unweighted exactness pay for it.

Mathematics

Fitting y = ax + b by least squares has a closed form: minimize Σ(yᵢ − a xᵢ − b)², set both partial derivatives to zero, and a 2×2 linear system (the normal equations) falls out - solved exactly by Cramer's rule.

[ Sxx Sx ] [a] [Sxy]

[ Sx N ] [b] = [Sy ]

a = (Sxy·N − Sx·Sy) / (Sxx·N − Sx²), b = (Sxx·Sy − Sx·Sxy) / (Sxx·N − Sx²)

This is the algebraic counterpart to iterative fitting: one exact solve instead of a descent. That exactness is also the lesson the outlier slider teaches - a single far point drags the whole line, because least squares squares its residuals and unweighted exactness has no defenses. The condition number of the normal matrix shows when the data itself (nearly constant x) makes the problem fragile.

How it works

Data is drawn from a seeded stream: x uniform, y from the true line plus Gaussian noise, plus an outlier jump for a controllable fraction. The fit, residuals, RMS, R², and the normal-matrix condition number all recompute on any parameter change; residual stems visualize each point's distance to the line.

With noise at zero the fit recovers the true coefficients to float precision - which the test suite asserts.

Implementation

src/engine/experiments/least-squares.ts · 177 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import { gaussian } from "../core/rng";
3import type { ParameterDef, ThemeColors } from "../core/types";
4
5/**
6 * Least squares by the normal equations: the 2x2 system [Sxx Sx; Sx N] is
7 * solved exactly by Cramer's rule - no iterations, no descent. Deliberately
8 * the algebraic counterpart to gradient-based fitting: the outlier slider
9 * shows exactly what unweighted exactness costs.
10 */
11export class LeastSquaresExperiment extends BaseExperiment {
12 readonly id = "least-squares";
13
14 private xs: number[] = [];
15 private ys: number[] = [];
16 private fitA = 0;
17 private fitB = 0;
18 private rms = 0;
19 private r2 = 1;
20
21 private readonly trueA = 1.5;
22 private readonly trueB = 0.3;
23
24 protected params(): ParameterDef[] {
25 return [
26 { key: "points", label: "Data points", min: 10, max: 200, step: 10, defaultValue: 80 },
27 { key: "noise", label: "Noise sigma", min: 0, max: 0.4, step: 0.01, defaultValue: 0.1 },
28 { key: "outliers", label: "Outlier rate", min: 0, max: 30, step: 1, defaultValue: 0, unit: "%" },
29 ];
30 }
31
32 protected onReset(): void {
33 const count = this.num("points");
34 const sigma = this.num("noise");
35 const outlierRate = this.num("outliers") / 100;
36 this.xs = [];
37 this.ys = [];
38 for (let i = 0; i < count; i++) {
39 const x = this.rng();
40 let y = this.trueA * x + this.trueB + gaussian(this.rng) * sigma;
41 if (this.rng() < outlierRate) {
42 y += (this.rng() < 0.5 ? -1 : 1) * (0.8 + this.rng() * 0.8);
43 }
44 this.xs.push(x);
45 this.ys.push(y);
46 }
47 this.fit();
48 }
49
50 protected onParameterChange(): void {
51 this.reset();
52 }
53
54 /** Static plot: state changes only through parameters. */
55 protected onUpdate(): void {}
56
57 private fit(): void {
58 const n = this.xs.length;
59 let sx = 0;
60 let sy = 0;
61 let sxx = 0;
62 let sxy = 0;
63 for (let i = 0; i < n; i++) {
64 sx += this.xs[i];
65 sy += this.ys[i];
66 sxx += this.xs[i] * this.xs[i];
67 sxy += this.xs[i] * this.ys[i];
68 }
69 // Normal equations [[Sxx, Sx], [Sx, N]] [a, b]^T = [Sxy, Sy]^T
70 const determinant = sxx * n - sx * sx;
71 if (Math.abs(determinant) < 1e-12) {
72 this.fitA = 0;
73 this.fitB = sy / n;
74 } else {
75 this.fitA = (sxy * n - sx * sy) / determinant;
76 this.fitB = (sxx * sy - sx * sxy) / determinant;
77 }
78 let rss = 0;
79 let tss = 0;
80 const meanY = sy / n;
81 for (let i = 0; i < n; i++) {
82 const residual = this.ys[i] - (this.fitA * this.xs[i] + this.fitB);
83 rss += residual * residual;
84 tss += (this.ys[i] - meanY) ** 2;
85 }
86 this.rms = Math.sqrt(rss / n);
87 this.r2 = tss > 0 ? 1 - rss / tss : 1;
88 }
89
90 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
91 ctx.fillStyle = theme.bg;
92 ctx.fillRect(0, 0, this.width, this.height);
93 const pad = 18;
94 const toX = (x: number) => pad + x * (this.width - 2 * pad);
95 const valueSpan = 2.2;
96 const toY = (y: number) => this.height - pad - ((y + 0.8) / valueSpan) * (this.height - 2 * pad);
97 // Residual stems
98 ctx.strokeStyle = theme.viz[1];
99 ctx.lineWidth = 1;
100 ctx.globalAlpha = 0.6;
101 ctx.beginPath();
102 for (let i = 0; i < this.xs.length; i++) {
103 const fitY = this.fitA * this.xs[i] + this.fitB;
104 ctx.moveTo(toX(this.xs[i]), toY(this.ys[i]));
105 ctx.lineTo(toX(this.xs[i]), toY(fitY));
106 }
107 ctx.stroke();
108 ctx.globalAlpha = 1;
109 // True line
110 ctx.strokeStyle = theme.fgTertiary;
111 ctx.setLineDash([4, 4]);
112 ctx.beginPath();
113 ctx.moveTo(toX(0), toY(this.trueA * 0 + this.trueB));
114 ctx.lineTo(toX(1), toY(this.trueA + this.trueB));
115 ctx.stroke();
116 ctx.setLineDash([]);
117 // Fitted line
118 ctx.strokeStyle = theme.accent;
119 ctx.lineWidth = 2;
120 ctx.beginPath();
121 ctx.moveTo(toX(0), toY(this.fitB));
122 ctx.lineTo(toX(1), toY(this.fitA + this.fitB));
123 ctx.stroke();
124 // Data
125 ctx.fillStyle = theme.fg;
126 for (let i = 0; i < this.xs.length; i++) {
127 ctx.beginPath();
128 ctx.arc(toX(this.xs[i]), toY(this.ys[i]), 2.2, 0, Math.PI * 2);
129 ctx.fill();
130 }
131 ctx.fillStyle = theme.fgTertiary;
132 ctx.font = "11px ui-monospace, monospace";
133 ctx.fillText("TRUE", pad + 4, pad + 12);
134 ctx.fillStyle = theme.accent;
135 ctx.fillText("FITTED", pad + 48, pad + 12);
136 }
137
138 getMetrics() {
139 const condition =
140 Math.abs(this.fitA) > 0 ? this.conditionNumber().toExponential(2) : "n/a";
141 return {
142 fittedA: this.fitA.toFixed(3),
143 fittedB: this.fitB.toFixed(3),
144 trueValues: `${this.trueA} / ${this.trueB}`,
145 rms: this.rms.toFixed(4),
146 r2: this.r2.toFixed(4),
147 condition: condition,
148 };
149 }
150
151 /** 2x2 condition number via the ratio of extreme eigenvalues of M^T M. */
152 private conditionNumber(): number {
153 const n = this.xs.length;
154 let sxx = 0;
155 let sx = 0;
156 for (let i = 0; i < n; i++) {
157 sxx += this.xs[i] * this.xs[i];
158 sx += this.xs[i];
159 }
160 // Normal matrix [[sxx, sx], [sx, n]]; eigenvalues closed-form.
161 const tr = sxx + n;
162 const det = sxx * n - sx * sx;
163 const disc = Math.sqrt(Math.max(0, tr * tr - 4 * det));
164 const lambdaMax = (tr + disc) / 2;
165 const lambdaMin = Math.max(1e-12, (tr - disc) / 2);
166 return lambdaMax / lambdaMin;
167 }
168
169 describe(): string {
170 return `Normal equations fitted y = ${this.fitA.toFixed(2)}x + ${this.fitB.toFixed(2)} against truth ${this.trueA}x + ${this.trueB}; RMS residual ${this.rms.toFixed(3)}, R2 ${this.r2.toFixed(3)}.`;
171 }
172
173 entities(): number {
174 return this.xs.length;
175 }
176}
177

Performance

Time per tick
O(N) for the fit - a single pass to accumulate five sums
Space
The data arrays only
Notes
Same problem, different philosophy than gradient descent: exact algebra versus iterative refinement, each with its own failure modes.

Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.