Concept · intermediate
Train and test sets
See why a model is judged on data it never saw, as polynomials of rising degree overfit.
A model that fits its training data perfectly can still be useless. It may have learned the noise in those particular points rather than the pattern behind them, and the only way to find out is to test it on data it never saw. That’s why several pages on this site hold back a test set. This page shows what goes wrong without one.
The dataset
data.csv holds 36 drives of one car, each at a steady speed in km/h,
with its fuel use in litres per 100 km. The data is synthetic, and we
know the truth: fuel use follows a gentle U, lowest near 65 km/h, plus
noise. The split column marks 12 drives as train and 24 as test.
We kept the training set small on purpose, because that’s when
overfitting is easiest to see.
The method
We fit a polynomial of each degree from 1 to 8 to the training drives alone. A degree- polynomial has coefficients:
Here is the speed mapped onto −1 to 1, using the training range of 36 to 116 km/h. Without that step, degree 8 would add up raw speeds to the 16th power, around , next to small numbers, and the arithmetic would lose accuracy.
Least squares picks the coefficients that make the squared errors on the training drives as small as possible. Setting the slope of that error to zero gives linear equations in unknowns, the normal equations, and each implementation solves them by Gaussian elimination: subtract multiples of each row from the rows below until the system is triangular, then solve from the bottom up.
Then we measure each fit by its root-mean-square error, in litres per 100 km, once on the training drives and once on the test drives:
Reading the results
| Degree | Train RMSE | Test RMSE |
|---|---|---|
| 1 | 0.529 | 0.779 |
| 2 | 0.284 | 0.391 |
| 3 | 0.278 | 0.405 |
| 4 | 0.272 | 0.390 |
| 5 | 0.220 | 0.546 |
| 6 | 0.186 | 0.630 |
| 7 | 0.149 | 0.549 |
| 8 | 0.074 | 1.153 |
The training error falls at every step. Each extra coefficient gives the curve one more way to bend toward the 12 points, so it can only fit them better. By degree 8, it nearly runs through every one.
The test error tells the real story. A straight line is too stiff to follow the U. Degree 2, the true shape, cuts the test error in half. Degrees 3 and 4 barely change it. From degree 5, the test error rises, reaching 1.153 at degree 8, the worst of all.
Look at where the training speeds fall: five between 36 and 54 km/h, five between 101 and 116, and only two, 71 and 81, in between. In that gap, the degree-8 curve is free to wiggle. At 70 km/h it predicts 3.90 litres, and at 90 it predicts 6.51, when the truth is 4.52 and 5.06. Four test drives also sit just outside the training range, where every polynomial has to extrapolate. But even the 20 test drives inside the range give degree 8 an error of 0.80, against 0.33 for degree 2.
Degrees 2 and 4 nearly tie on the test set. When they do, choose the simpler one: it has fewer ways to be wrong on data you haven’t seen yet.
What to watch for across languages
None of these implementations calls a fitter or a solver, though several
languages ship one: Octave’s polyfit and backslash, Julia’s \, APL’s
⌹, Ada’s Solve, and R’s lm(). NumPy’s polyfit would do it in
Python, but NumPy isn’t installed. Each implementation writes out the
elimination instead, and three things vary.
The first is scaling. In Julia, solving the normal equations on raw
speeds with the built-in \ gives a degree-8 training error of 0.073996
instead of 0.073821, a difference big enough to fail this site’s tests.
R’s lm() gets the right answer even from raw speeds, because it avoids
the normal equations and uses a more stable method, a QR decomposition.
The second is how a row swap works. Where a matrix is a list of rows, as
in Python, Go, Java, Lua, and OCaml, swapping trades two references and
copies no numbers. C can’t assign one array row to another, so it swaps
entry by entry, as do Common Lisp and Pascal with their 2-D arrays.
Octave, Julia, R, and APL swap two rows in one assignment by indexing
both at once. In Swift, Tcl, and PHP, arrays are values, so solve
works on its own copy and the caller’s equations stay untouched.
The third is the small traps along the way. APL reads K-1↓⍳N as
K-(1↓⍳N), since it evaluates right to left. In PostgreSQL, m[2] on a
2-D array gives NULL, not row 2; you need m[2:2]. SWI-Prolog’s /
returns an integer when it divides exactly, so the code converts each
scaled speed to a float. And COBOL works in fixed-point decimals, with 18
digits after the point, and still matches double precision to within
1e-11.
Data preview
| speed | litres | split |
|---|---|---|
| 42 | 4.87 | train |
| 106 | 5.78 | train |
| 39 | 4.9 | train |
| 108 | 5.98 | train |
| 45 | 4.47 | train |
| 71 | 3.89 | train |
| 101 | 5.71 | train |
| 54 | 4.7 | train |
First 8 of 36 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
train_count | 12 |
test_count | 24 |
train_rmse_d1 | 0.528506 |
train_rmse_d2 | 0.284284 |
train_rmse_d3 | 0.277800 |
train_rmse_d4 | 0.271889 |
train_rmse_d5 | 0.219795 |
train_rmse_d6 | 0.186111 |
train_rmse_d7 | 0.149015 |
train_rmse_d8 | 0.073821 |
test_rmse_d1 | 0.779162 |
test_rmse_d2 | 0.390931 |
test_rmse_d3 | 0.405384 |
test_rmse_d4 | 0.389787 |
test_rmse_d5 | 0.545971 |
test_rmse_d6 | 0.630430 |
test_rmse_d7 | 0.548756 |
test_rmse_d8 | 1.152968 |
Implementations
"""Train and test error for polynomials of degree 1 to 8, using only the
standard library.
Each polynomial is fitted by least squares to the training drives only:
we build the normal equations and solve them by Gaussian elimination. Then
we measure its root-mean-square error on the training drives and on the
test drives it never saw.
"""
import csv
import math
import sys
MAX_DEGREE = 8
def read(path):
train, test = [], []
with open(path, newline="") as f:
for row in csv.DictReader(f):
point = (float(row["speed"]), float(row["litres"]))
(train if row["split"] == "train" else test).append(point)
return train, test
def solve(a, b):
"""Solve a x = b by Gaussian elimination with partial pivoting."""
n = len(b)
m = [row[:] + [b[i]] for i, row in enumerate(a)] # augmented matrix
for col in range(n):
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small.
pivot = max(range(col, n), key=lambda r: abs(m[r][col]))
m[col], m[pivot] = m[pivot], m[col]
for r in range(col + 1, n):
factor = m[r][col] / m[col][col]
for c in range(col, n + 1):
m[r][c] -= factor * m[col][c]
# Back substitution, from the last unknown up.
x = [0.0] * n
for r in range(n - 1, -1, -1):
x[r] = (m[r][n] - sum(m[r][c] * x[c] for c in range(r + 1, n))) / m[r][r]
return x
def fit(points, degree, scale):
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
# of Z is 1, z_i, z_i^2, ..., z_i^degree.
size = degree + 1
a = [[0.0] * size for _ in range(size)]
b = [0.0] * size
for x, y in points:
z = scale(x)
for i in range(size):
b[i] += y * z ** i
for j in range(size):
a[i][j] += z ** (i + j)
return solve(a, b)
def rmse(coefs, points, scale):
total = 0.0
for x, y in points:
z = scale(x)
predicted = sum(c * z ** i for i, c in enumerate(coefs))
total += (predicted - y) ** 2
return math.sqrt(total / len(points))
def main():
train, test = read(sys.argv[1])
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power
# (the largest entry at degree 8) would swamp the small ones and make
# the equations hard to solve accurately.
low = min(x for x, _ in train)
high = max(x for x, _ in train)
mid, half = (low + high) / 2, (high - low) / 2
def scale(x):
return (x - mid) / half
print(f"train_count {len(train)}")
print(f"test_count {len(test)}")
errors = []
for degree in range(1, MAX_DEGREE + 1):
coefs = fit(train, degree, scale)
errors.append((degree, rmse(coefs, train, scale), rmse(coefs, test, scale)))
for degree, train_error, _ in errors:
print(f"train_rmse_d{degree} {train_error:.6f}")
for degree, _, test_error in errors:
print(f"test_rmse_d{degree} {test_error:.6f}")
if __name__ == "__main__":
main()Each fit builds the normal equations, a small square system whose answer
is the polynomial’s coefficients, and solve works through it by Gaussian
elimination. At each column it swaps up the row with the largest entry,
called partial pivoting, which keeps rounding error from growing.
Speeds are first mapped onto −1 to 1 using the training range only, so
powers up to the 16th stay a sensible size. Python’s standard library has
statistics.linear_regression for a straight line, but nothing we use
here fits a polynomial; NumPy’s polyfit would, but NumPy isn’t
installed.
-- Train and test error for polynomials of degree 1 to 8, using only Ada's
-- standard library.
--
-- Each polynomial is fitted by least squares to the training drives only:
-- we build the normal equations and solve them by Gaussian elimination.
-- Then we measure its root-mean-square error on the training drives and on
-- the test drives it never saw.
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Long_Float_Text_IO; use Ada.Long_Float_Text_IO;
with Ada.Command_Line;
with Ada.Strings.Fixed; use Ada.Strings.Fixed;
with Ada.Numerics.Long_Elementary_Functions;
use Ada.Numerics.Long_Elementary_Functions;
with Ada.Numerics.Long_Real_Arrays;
use Ada.Numerics.Long_Real_Arrays;
with Ada.Containers.Vectors;
procedure Main is
Max_Degree : constant := 8;
-- Real_Vector and Real_Matrix come from Annex G of the standard, with
-- "*" for matrix and vector products. Their index ranges are ours to
-- choose, so coefficient vectors count from 0, like the powers.
package Float_Vectors is new Ada.Containers.Vectors (Positive, Long_Float);
-- The text of field number Col (counting from 1) in a CSV line.
function Field (Line : String; Col : Positive) return String is
Start : Positive := Line'First;
begin
for I in 1 .. Col - 1 loop
Start := Index (Line (Start .. Line'Last), ",") + 1;
end loop;
declare
Stop : constant Natural := Index (Line (Start .. Line'Last), ",");
begin
return (if Stop = 0 then Line (Start .. Line'Last)
else Line (Start .. Stop - 1));
end;
end Field;
-- The number of the header field called Name.
function Column (Header, Name : String) return Positive is
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
if Field (Header, C) = Name then
return C;
end if;
end loop;
raise Constraint_Error with "column not found: " & Name;
end Column;
-- A growable list of numbers as a Real_Vector, counting from 1.
function To_Vector (List : Float_Vectors.Vector) return Real_Vector is
Result : Real_Vector (1 .. Natural (List.Length));
begin
for I in Result'Range loop
Result (I) := List (I);
end loop;
return Result;
end To_Vector;
-- The design matrix: row i is 1, z_i, z_i^2, ..., z_i^Degree.
function Powers (Z : Real_Vector; Degree : Natural) return Real_Matrix is
Result : Real_Matrix (Z'Range, 0 .. Degree);
begin
for I in Z'Range loop
for J in 0 .. Degree loop
Result (I, J) := Z (I) ** J;
end loop;
end loop;
return Result;
end Powers;
-- Solve A x = B by Gaussian elimination with partial pivoting. A and B
-- share the index range 0 .. N.
function Gauss_Solve (A : Real_Matrix; B : Real_Vector) return Real_Vector
is
Last : constant Integer := B'Last;
-- The augmented matrix: A with B as an extra last column.
M : Real_Matrix (B'Range, B'First .. Last + 1);
X : Real_Vector (B'Range) := [others => 0.0];
Pivot : Integer;
begin
for I in B'Range loop
for J in B'Range loop
M (I, J) := A (I, J);
end loop;
M (I, Last + 1) := B (I);
end loop;
for Col in B'Range loop
-- Swap up the row with the largest entry in this column, which
-- keeps the rounding error small. A 2-D array has no row objects
-- to swap, so the two rows swap one entry at a time.
Pivot := Col;
for R in Col + 1 .. Last loop
if abs M (R, Col) > abs M (Pivot, Col) then
Pivot := R;
end if;
end loop;
for C in M'Range (2) loop
declare
Held : constant Long_Float := M (Col, C);
begin
M (Col, C) := M (Pivot, C);
M (Pivot, C) := Held;
end;
end loop;
for R in Col + 1 .. Last loop
declare
Factor : constant Long_Float := M (R, Col) / M (Col, Col);
begin
for C in Col .. Last + 1 loop
M (R, C) := M (R, C) - Factor * M (Col, C);
end loop;
end;
end loop;
end loop;
-- Back substitution, from the last unknown up.
for R in reverse B'Range loop
declare
Known : Long_Float := 0.0;
begin
for C in R + 1 .. Last loop
Known := Known + M (R, C) * X (C);
end loop;
X (R) := (M (R, Last + 1) - Known) / M (R, R);
end;
end loop;
return X;
end Gauss_Solve;
-- Least squares: the coefficients c solve (Z^T Z) c = Z^T y.
function Fit (Z, Y : Real_Vector; Degree : Natural) return Real_Vector is
D : constant Real_Matrix := Powers (Z, Degree);
begin
return Gauss_Solve (Transpose (D) * D, Transpose (D) * Y);
end Fit;
function RMSE (Coefs, Z, Y : Real_Vector) return Long_Float is
Errors : constant Real_Vector := Powers (Z, Coefs'Last) * Coefs - Y;
begin
-- Errors * Errors is the inner product, the sum of squared errors.
return Sqrt ((Errors * Errors) / Long_Float (Y'Length));
end RMSE;
procedure Say (Name : String; Value : Long_Float) is
begin
Put (Name & " ");
Put (Value, Fore => 1, Aft => 6, Exp => 0); -- six decimals, no E
New_Line;
end Say;
-- N'Image is " 1", so the slice drops the leading space.
function Tag (N : Natural) return String is (N'Image (2 .. N'Image'Last));
File : File_Type;
Train_X, Train_Y, Test_X, Test_Y : Float_Vectors.Vector;
begin
Open (File, In_File, Ada.Command_Line.Argument (1));
declare
Header : constant String := Get_Line (File);
Speed_C : constant Positive := Column (Header, "speed");
Litres_C : constant Positive := Column (Header, "litres");
Split_C : constant Positive := Column (Header, "split");
begin
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
Speed : constant Long_Float :=
Long_Float'Value (Field (Line, Speed_C));
Litres : constant Long_Float :=
Long_Float'Value (Field (Line, Litres_C));
begin
if Field (Line, Split_C) = "train" then
Train_X.Append (Speed);
Train_Y.Append (Litres);
else
Test_X.Append (Speed);
Test_Y.Append (Litres);
end if;
end;
end loop;
end;
Close (File);
declare
Speeds : constant Real_Vector := To_Vector (Train_X);
Low : Long_Float := Speeds (Speeds'First);
High : Long_Float := Speeds (Speeds'First);
begin
for S of Speeds loop
Low := Long_Float'Min (Low, S);
High := Long_Float'Max (High, S);
end loop;
declare
-- Map the training speeds onto -1 to 1. Raw speeds to the 16th
-- power (the largest entry at degree 8) would swamp the small
-- ones and make the equations hard to solve accurately.
Mid : constant Long_Float := (Low + High) / 2.0;
Half : constant Long_Float := (High - Low) / 2.0;
-- An aggregate with a for loop inside builds each vector of z.
Train_Z : constant Real_Vector (Speeds'Range) :=
[for I in Speeds'Range => (Speeds (I) - Mid) / Half];
Test_Speeds : constant Real_Vector := To_Vector (Test_X);
Test_Z : constant Real_Vector (Test_Speeds'Range) :=
[for I in Test_Speeds'Range => (Test_Speeds (I) - Mid) / Half];
Y_Train : constant Real_Vector := To_Vector (Train_Y);
Y_Test : constant Real_Vector := To_Vector (Test_Y);
Train_Errors, Test_Errors : Real_Vector (1 .. Max_Degree);
begin
for D in 1 .. Max_Degree loop
declare
Coefs : constant Real_Vector := Fit (Train_Z, Y_Train, D);
begin
Train_Errors (D) := RMSE (Coefs, Train_Z, Y_Train);
Test_Errors (D) := RMSE (Coefs, Test_Z, Y_Test);
end;
end loop;
Put_Line ("train_count" & Natural'Image (Speeds'Length));
Put_Line ("test_count" & Natural'Image (Test_Speeds'Length));
for D in 1 .. Max_Degree loop
Say ("train_rmse_d" & Tag (D), Train_Errors (D));
end loop;
for D in 1 .. Max_Degree loop
Say ("test_rmse_d" & Tag (D), Test_Errors (D));
end loop;
end;
end;
end Main;The matrices are Annex G’s Real_Matrix and Real_Vector, whose index
ranges we choose, so coefficients count from 0 like the powers. We use
their Transpose and * to build the normal equations and the
predictions, since multiplying matrices is a general tool.
Annex G also has Solve, which would solve the system in one call; we
skip it and write out Gaussian elimination instead. A 2-D array has no row
objects, so partial pivoting swaps two rows one entry at a time. An array
aggregate with a for inside maps every speed onto −1 to 1.
⍝ Train and test error for polynomials of degree 1 to 8, from scratch in GNU
⍝ APL. Each polynomial is fitted by least squares to the training drives
⍝ only: we build the normal equations and solve them by Gaussian
⍝ elimination. Then we measure its root-mean-square error on the training
⍝ drives and on the test drives it never saw.
⍝ APL reads right to left, with no operator precedence: 2×3+4 is 2×(3+4).
⍝ ⍝ starts a comment; ← assigns; ⍵ is a function's right argument.
∇NAME SAY N
⍝ Print "name value" for a whole number. ⍕ formats a number as text.
⎕←NAME,' ',⍕N
∇
∇NAME SAY6 X;T
⍝ Print "name value" with six decimal places.
T←0 6⍕|X ⍝ |X is the size of X; 0 6⍕ gives 6 decimals
T←(∨\T≠' ')/T ⍝ drop the leading blank ⍕ adds
T←((T[1]='.')/'0'),T ⍝ ⍕ writes .94 for 0.94, so put the 0 back
⎕←NAME,' ',((X<0)/'-'),T ⍝ sign as -, since APL's own minus is ¯
∇
MAX_DEGREE←8
⍝ Read the file. ⎕ARG is the command line; ↑¯1↑ takes its last item, the
⍝ data path. ⎕FIO[26] reads the whole file as one string.
TEXT←⎕FIO[26]↑¯1↑⎕ARG
LINES←(TEXT≠⎕UCS 10)⊂TEXT ⍝ ⊂ splits at the 0s of the mask: one item per line
⍝ Split each (¨) line at commas; ⊃ stacks them into a table.
ROWS←⊃{(⍵≠',')⊂⍵}¨LINES
⍝ A column by header name (⍳), minus the header (1↓).
COLUMN←{1↓ROWS[;ROWS[1;]⍳⊂⍵]}
SPEED←⍎¨COLUMN 'speed' ⍝ ⍎ turns text to a number
LITRES←⍎¨COLUMN 'litres'
⍝ 1 for each training row (≡ matches whole strings); / keeps those rows.
TRAIN←(COLUMN 'split')≡¨⊂'train'
TX←TRAIN/SPEED ⋄ TY←TRAIN/LITRES
QX←(~TRAIN)/SPEED ⋄ QY←(~TRAIN)/LITRES ⍝ ~ is not: the test rows
⍝ Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
⍝ largest entry at degree 8) would swamp the small ones and make the
⍝ equations hard to solve accurately. ⌈/ is the largest, ⌊/ the smallest.
MID←((⌈/TX)+⌊/TX)÷2 ⋄ HALF←((⌈/TX)-⌊/TX)÷2
SCALE←{(⍵-MID)÷HALF}
∇X←SOLVE M;N;K;P;V;R;C
⍝ Solve the augmented system M (N rows, N + 1 columns, b last) by
⍝ Gaussian elimination with partial pivoting.
N←≢M ⋄ K←1
PIVOT:
⍝ Swap up the row with the largest entry in column K, at or below row K,
⍝ which keeps the rounding error small. Indexing with K,P on both sides
⍝ swaps the two rows in one assignment.
V←|M[(K-1)↓⍳N;K]
P←(K-1)+V⍳⌈/V ⍝ ⍳ finds the first largest
M[K,P;]←M[P,K;]
⍝ Every row below K subtracts a multiple of row K at once: the outer
⍝ product ∘.× makes the table of factor × row K entry.
R←K↓⍳N ⋄ C←(K-1)↓⍳N+1
M[R;C]←M[R;C]-(M[R;K]÷M[K;K])∘.×M[K;C]
K←K+1 ⋄ →(K≤N)/PIVOT
⍝ Back substitution, from the last unknown up.
X←N⍴0 ⋄ K←N
BACK:
R←K↓⍳N ⍝ the unknowns already found
X[K]←(M[K;N+1]-M[K;R]+.×X[R])÷M[K;K]
K←K-1 ⋄ →(K≥1)/BACK
∇
∇COEFS←FIT D;Z;E;A;B
⍝ Least squares: the coefficients solve (Z^T Z) c = Z^T y. E is the table
⍝ of exponents i + j; Z∘.*E raises every speed to every one of them, and
⍝ +⌿ sums over the speeds, so A[i;j] is the sum of z^(i+j).
Z←SCALE TX
E←(0,⍳D)∘.+0,⍳D
A←+⌿Z∘.*E
B←TY+.×Z∘.*0,⍳D ⍝ the sum of y × z^i, for each i
COEFS←SOLVE A,B ⍝ , joins B on as the last column
∇
∇E←COEFS RMSE XY;P
⍝ XY is speeds and litres; each prediction is the powers of z times the
⍝ coefficients. *0.5 is the square root.
P←((SCALE⊃XY[1])∘.*¯1+⍳≢COEFS)+.×COEFS
E←((+/(P-⊃XY[2])*2)÷≢P)*0.5
∇
'train_count' SAY ≢TX
'test_count' SAY ≢QX
DEGREES←⍳MAX_DEGREE
FITS←FIT¨DEGREES
⍝ ⊂ encloses the prefix so ,¨ joins it to each degree's text, and ⊂ around
⍝ each pair of vectors lets ¨ hand the same data to every fit.
((⊂'train_rmse_d'),¨⍕¨DEGREES) SAY6¨FITS RMSE¨⊂TX TY
((⊂'test_rmse_d'),¨⍕¨DEGREES) SAY6¨FITS RMSE¨⊂QX QY
)OFFWe keep the augmented system as an ordinary APL matrix, so whole rows are
one index away. The pivot swap is a single assignment,
M[K,P;]←M[P,K;], and each elimination step updates every lower row at
once with an outer product. The normal equations need no loop either:
Z∘.*E raises every speed to every exponent in the table E, and +⌿
sums over the speeds.
APL’s right-to-left rule bites here: K-1↓⍳N drops one item before
subtracting, so the code writes (K-1)↓⍳N. We skip ⌹, APL’s built-in
matrix divide, which would solve the system in one symbol.
# Train and test error for polynomials of degree 1 to 8, in POSIX awk.
#
# Each polynomial is fitted by least squares to the training drives only:
# we build the normal equations and solve them by Gaussian elimination. Then
# we measure its root-mean-square error on the training drives and on the
# test drives it never saw. Rows are stored as they are read; the work
# happens in END.
BEGIN { FS = ","; MAX_DEGREE = 8 }
# The header line: find each column by name.
NR == 1 {
for (i = 1; i <= NF; i++) {
if ($i == "speed") c_speed = i
if ($i == "litres") c_litres = i
if ($i == "split") c_split = i
}
if (!c_speed || !c_litres || !c_split) {
print "missing a column (speed, litres, split)" > "/dev/stderr"
failed = 1
exit 1
}
next
}
# Train and test drives go into separate arrays, numbered from 1.
$c_split == "train" { n_train++; train_x[n_train] = $c_speed; train_y[n_train] = $c_litres; next }
{ n_test++; test_x[n_test] = $c_speed; test_y[n_test] = $c_litres }
function scale(x) { return (x - mid) / half }
# Least squares: the coefficients solve (Z^T Z) c = Z^T y, where row i of Z
# is 1, z_i, z_i^2, ..., z_i^degree. awk has no 2D arrays, so m[i, j] is a
# single key "i SUBSEP j". Column size + 1 holds the right-hand side, making
# m the augmented matrix. The answer lands in the global array coef[1..size].
function fit(degree, size, i, j, p, z) {
size = degree + 1
for (i = 1; i <= size; i++)
for (j = 1; j <= size + 1; j++) m[i, j] = 0
for (p = 1; p <= n_train; p++) {
z = scale(train_x[p])
for (i = 1; i <= size; i++) {
m[i, size + 1] += train_y[p] * z ^ (i - 1)
for (j = 1; j <= size; j++) m[i, j] += z ^ (i + j - 2)
}
}
solve(size)
}
# Gaussian elimination with partial pivoting on the augmented matrix m.
function solve(n, col, r, c, pivot, t, factor, total) {
for (col = 1; col <= n; col++) {
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small. With no row objects, swap cell by cell.
pivot = col
for (r = col + 1; r <= n; r++)
if (abs(m[r, col]) > abs(m[pivot, col])) pivot = r
for (c = 1; c <= n + 1; c++) {
t = m[col, c]; m[col, c] = m[pivot, c]; m[pivot, c] = t
}
for (r = col + 1; r <= n; r++) {
factor = m[r, col] / m[col, col]
for (c = col; c <= n + 1; c++) m[r, c] -= factor * m[col, c]
}
}
# Back substitution, from the last unknown up.
for (r = n; r >= 1; r--) {
total = m[r, n + 1]
for (c = r + 1; c <= n; c++) total -= m[r, c] * coef[c]
coef[r] = total / m[r, r]
}
}
function abs(v) { return v < 0 ? -v : v }
# Arrays are passed by reference, so one function serves both sets.
function rmse(size, xs, ys, count, p, i, z, predicted, total) {
total = 0
for (p = 1; p <= count; p++) {
z = scale(xs[p])
predicted = 0
for (i = 1; i <= size; i++) predicted += coef[i] * z ^ (i - 1)
total += (predicted - ys[p]) ^ 2
}
return sqrt(total / count)
}
END {
if (failed) exit 1
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power
# (the largest entry at degree 8) would swamp the small ones and make
# the equations hard to solve accurately. The + 0 makes each a number.
low = high = train_x[1] + 0
for (p = 2; p <= n_train; p++) {
if (train_x[p] < low) low = train_x[p] + 0
if (train_x[p] > high) high = train_x[p] + 0
}
mid = (low + high) / 2
half = (high - low) / 2
printf "train_count %d\n", n_train
printf "test_count %d\n", n_test
for (d = 1; d <= MAX_DEGREE; d++) {
fit(d)
train_rmse[d] = rmse(d + 1, train_x, train_y, n_train)
test_rmse[d] = rmse(d + 1, test_x, test_y, n_test)
}
for (d = 1; d <= MAX_DEGREE; d++) printf "train_rmse_d%d %.6f\n", d, train_rmse[d]
for (d = 1; d <= MAX_DEGREE; d++) printf "test_rmse_d%d %.6f\n", d, test_rmse[d]
}mawk, the awk in this image, has no arrays of arrays, so the augmented matrix is m[i, j], a single
key that joins the two indices with SUBSEP. Column size + 1 holds the
right-hand side. With no row objects to exchange, the pivot swap moves the
two rows cell by cell through a temporary.
The coefficients land in a global array, coef, and rmse takes the
train or test arrays as arguments, since awk passes arrays by reference.
This image’s mawk raises a negative scaled speed to a whole power with ^
correctly, so no power function is needed.
10 REM Train and test error for polynomials of degree 1 to 8, in classic
20 REM BASIC.
30 REM
40 REM Each polynomial is fitted by least squares to the training drives
50 REM only: we build the normal equations and solve them by Gaussian
60 REM elimination. Then we measure its root-mean-square error on the
70 REM training drives and on the test drives it never saw.
80 REM
90 REM run.sh passes the data path in the DATA environment variable,
100 REM because bwbasic can't read command-line arguments.
110 MAXDEG = 8
120 REM TX, TY hold the training speeds (scaled) and litres; SX, SY the
130 REM test ones. M is the augmented matrix, the equations with their
140 REM right-hand side in column NS + 1, and CF the coefficients.
150 DIM TX(500), TY(500), SX(500), SY(500), FL$(20)
160 DIM M(9, 10), CF(9), PW(17), TR(8), TE(8)
170 GOSUB 1000 : REM read the drives into TX, TY and SX, SY
180 REM Map the training speeds onto -1 to 1. Raw speeds to the 16th
190 REM power (the largest entry at degree 8) would swamp the small ones
200 REM and make the equations hard to solve accurately.
210 LO = TX(1) : HI = TX(1)
220 FOR I = 1 TO NT
230 IF TX(I) < LO THEN LO = TX(I)
240 IF TX(I) > HI THEN HI = TX(I)
250 NEXT I
260 MD = (LO + HI) / 2 : HF = (HI - LO) / 2
270 FOR I = 1 TO NT
280 TX(I) = (TX(I) - MD) / HF
290 NEXT I
300 FOR I = 1 TO NE
310 SX(I) = (SX(I) - MD) / HF
320 NEXT I
330 PRINT "train_count"; NT
340 PRINT "test_count"; NE
350 FOR D = 1 TO MAXDEG
360 GOSUB 2000 : REM build the equations for degree D in M
370 GOSUB 5000 : REM solve them into CF
380 GOSUB 6000 : REM the two errors, into TR(D) and TE(D)
390 NEXT D
400 FOR D = 1 TO MAXDEG
410 L$ = "train_rmse_d" + MID$(STR$(D), 2) : Y = TR(D) : GOSUB 3000
420 NEXT D
430 FOR D = 1 TO MAXDEG
440 L$ = "test_rmse_d" + MID$(STR$(D), 2) : Y = TE(D) : GOSUB 3000
450 NEXT D
460 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read every drive, sorting it into the training or test arrays.
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find which field of the header holds each column
1050 FOR J = 1 TO NF
1060 IF FL$(J) = "speed" THEN CX = J
1070 IF FL$(J) = "litres" THEN CY = J
1080 IF FL$(J) = "split" THEN CS = J
1090 NEXT J
1100 NT = 0 : NE = 0
1110 IF EOF(1) THEN 1190
1120 LINE INPUT #1, L$
1130 GOSUB 4000
1140 X = VAL(FL$(CX)) : Y = VAL(FL$(CY))
1150 IF FL$(CS) = "train" THEN NT = NT + 1 : TX(NT) = X : TY(NT) = Y
1160 IF FL$(CS) <> "train" THEN NE = NE + 1 : SX(NE) = X : SY(NE) = Y
1170 GOTO 1110
1190 CLOSE #1
1200 RETURN
2000 REM Least squares for degree D: the coefficients solve NS equations,
2010 REM where M(I, J) is the sum of z^(I + J - 2) over the training drives
2020 REM and M(I, NS + 1) the sum of y z^(I - 1). (Rows count from 1, so
2030 REM row I is power I - 1.)
2040 NS = D + 1
2050 FOR I = 1 TO NS
2060 FOR J = 1 TO NS + 1
2070 M(I, J) = 0
2080 NEXT J
2090 NEXT I
2100 FOR P = 1 TO NT
2110 REM PW(K) is z^(K - 1), built by multiplying, from z^0 = 1 up.
2120 PW(1) = 1
2130 FOR K = 2 TO 2 * NS - 1
2140 PW(K) = PW(K - 1) * TX(P)
2150 NEXT K
2160 FOR I = 1 TO NS
2170 M(I, NS + 1) = M(I, NS + 1) + TY(P) * PW(I)
2180 FOR J = 1 TO NS
2190 M(I, J) = M(I, J) + PW(I + J - 1)
2200 NEXT J
2210 NEXT I
2220 NEXT P
2230 RETURN
3000 REM Print L$ and Y with six decimals.
3010 PRINT L$; " ";
3020 PRINT USING "#.######"; Y
3030 RETURN
4000 REM Split L$ at its commas into FL$(1) to FL$(NF).
4010 NF = 0 : P = 1
4020 Q = INSTR(P, L$, ",")
4030 NF = NF + 1
4040 IF Q = 0 THEN FL$(NF) = MID$(L$, P) : RETURN
4050 FL$(NF) = MID$(L$, P, Q - P) : P = Q + 1
4060 GOTO 4020
5000 REM Solve the NS equations in M by Gaussian elimination with partial
5010 REM pivoting, leaving the answer in CF(1..NS). A FOR loop with an
5020 REM empty range still runs once here, so the loops that can be empty
5030 REM are WHILE loops.
5040 FOR C = 1 TO NS
5050 REM Swap up the row with the largest entry in this column, which
5060 REM keeps the rounding error small.
5070 PV = C
5080 R = C + 1
5090 WHILE R <= NS
5100 IF ABS(M(R, C)) > ABS(M(PV, C)) THEN PV = R
5110 R = R + 1
5120 WEND
5130 FOR J = 1 TO NS + 1
5140 T = M(C, J) : M(C, J) = M(PV, J) : M(PV, J) = T
5150 NEXT J
5160 R = C + 1
5170 WHILE R <= NS
5180 FC = M(R, C) / M(C, C)
5190 FOR J = C TO NS + 1
5200 M(R, J) = M(R, J) - FC * M(C, J)
5210 NEXT J
5220 R = R + 1
5230 WEND
5240 NEXT C
5250 REM Back substitution, from the last unknown up.
5260 FOR R = NS TO 1 STEP -1
5270 T = M(R, NS + 1)
5280 J = R + 1
5290 WHILE J <= NS
5300 T = T - M(R, J) * CF(J)
5310 J = J + 1
5320 WEND
5330 CF(R) = T / M(R, R)
5340 NEXT R
5350 RETURN
6000 REM Root-mean-square error of the fit CF on the training drives, into
6010 REM TR(D), and on the test drives, into TE(D).
6020 T = 0
6030 FOR P = 1 TO NT
6040 Z = TX(P) : GOSUB 7000
6050 T = T + (YH - TY(P)) ^ 2
6060 NEXT P
6070 TR(D) = SQR(T / NT)
6080 T = 0
6090 FOR P = 1 TO NE
6100 Z = SX(P) : GOSUB 7000
6110 T = T + (YH - SY(P)) ^ 2
6120 NEXT P
6130 TE(D) = SQR(T / NE)
6140 RETURN
7000 REM YH = the polynomial CF at Z, by Horner's rule:
7010 REM c1 + z (c2 + z (c3 + ...)).
7020 YH = 0
7030 FOR K = NS TO 1 STEP -1
7040 YH = YH * Z + CF(K)
7050 NEXT K
7060 RETURNThere are no built-ins to skip here, so every step is a numbered block of
lines. The equations sit in a fixed M(9, 10) array, big enough for degree
8, with the right-hand side as the last column. Swapping two rows is a FOR
loop that trades one entry at a time through T. Powers come from repeated
multiplication, and predictions use Horner’s rule.
Watch the loops. A FOR with an empty range still runs once in bwbasic:
FOR R = 3 TO 2 runs its body with R = 3. So elimination and back
substitution, whose inner ranges are empty on the last row, use WHILE
loops instead.
/* Train and test error for polynomials of degree 1 to 8, using only the C
* standard library.
*
* Each polynomial is fitted by least squares to the training drives only:
* we build the normal equations and solve them by Gaussian elimination.
* Then we measure its root-mean-square error on the training drives and on
* the test drives it never saw. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAX_DEGREE 8
#define MAX_TERMS (MAX_DEGREE + 1)
#define MAX_DRIVES 1000
#define MAX_LINE 256
/* Training and test drives in separate parallel arrays: speed and litres. */
static double train_x[MAX_DRIVES], train_y[MAX_DRIVES];
static double test_x[MAX_DRIVES], test_y[MAX_DRIVES];
static int train_count, test_count;
/* The training range, used to map every speed onto -1 to 1. */
static double mid, half;
static double scale(double x) { return (x - mid) / half; }
/* Cut a line at its commas, in place, and return the number of fields.
* strtok replaces each comma with '\0', so fields[i] points into line. */
static int split(char *line, char *fields[], int max) {
int n = 0;
for (char *f = strtok(line, ",\r\n"); f && n < max; f = strtok(NULL, ",\r\n"))
fields[n++] = f;
return n;
}
/* Solve a x = b by Gaussian elimination with partial pivoting. The
* augmented matrix m is n rows by n + 1 columns; the last column is b. */
static void solve(int n, double m[MAX_TERMS][MAX_TERMS + 1], double x[]) {
for (int col = 0; col < n; col++) {
/* Swap up the row with the largest entry in this column, which keeps
* the rounding error small. C arrays can't be assigned whole, so the
* swap goes entry by entry. */
int pivot = col;
for (int r = col + 1; r < n; r++)
if (fabs(m[r][col]) > fabs(m[pivot][col])) pivot = r;
for (int c = 0; c <= n; c++) {
double t = m[col][c];
m[col][c] = m[pivot][c];
m[pivot][c] = t;
}
for (int r = col + 1; r < n; r++) {
double factor = m[r][col] / m[col][col];
for (int c = col; c <= n; c++) m[r][c] -= factor * m[col][c];
}
}
/* Back substitution, from the last unknown up. */
for (int r = n - 1; r >= 0; r--) {
double s = m[r][n];
for (int c = r + 1; c < n; c++) s -= m[r][c] * x[c];
x[r] = s / m[r][r];
}
}
/* Least squares: the coefficients solve (Z^T Z) c = Z^T y, where row i of
* Z is 1, z_i, z_i^2, ..., z_i^degree. */
static void fit(int degree, double coefs[]) {
int size = degree + 1;
double m[MAX_TERMS][MAX_TERMS + 1] = {{0}};
for (int k = 0; k < train_count; k++) {
double z = scale(train_x[k]);
for (int i = 0; i < size; i++) {
m[i][size] += train_y[k] * pow(z, i);
for (int j = 0; j < size; j++) m[i][j] += pow(z, i + j);
}
}
solve(size, m, coefs);
}
static double rmse(const double coefs[], int degree, const double xs[],
const double ys[], int n) {
double total = 0;
for (int k = 0; k < n; k++) {
double z = scale(xs[k]), predicted = 0;
for (int i = 0; i <= degree; i++) predicted += coefs[i] * pow(z, i);
total += (predicted - ys[k]) * (predicted - ys[k]);
}
return sqrt(total / n);
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
FILE *file = fopen(argv[1], "r");
if (!file) { perror(argv[1]); return 1; }
char line[MAX_LINE], *fields[3];
if (!fgets(line, sizeof line, file)) { fprintf(stderr, "empty file\n"); return 1; }
while (fgets(line, sizeof line, file)) {
if (split(line, fields, 3) != 3) continue;
double x = strtod(fields[0], NULL), y = strtod(fields[1], NULL);
if (strcmp(fields[2], "train") == 0 && train_count < MAX_DRIVES) {
train_x[train_count] = x;
train_y[train_count++] = y;
} else if (test_count < MAX_DRIVES) {
test_x[test_count] = x;
test_y[test_count++] = y;
}
}
fclose(file);
/* Map the training speeds onto -1 to 1. Raw speeds to the 16th power
* (the largest entry at degree 8) would swamp the small ones. */
double low = train_x[0], high = train_x[0];
for (int k = 1; k < train_count; k++) {
if (train_x[k] < low) low = train_x[k];
if (train_x[k] > high) high = train_x[k];
}
mid = (low + high) / 2;
half = (high - low) / 2;
double train_error[MAX_DEGREE + 1], test_error[MAX_DEGREE + 1];
for (int d = 1; d <= MAX_DEGREE; d++) {
double coefs[MAX_TERMS];
fit(d, coefs);
train_error[d] = rmse(coefs, d, train_x, train_y, train_count);
test_error[d] = rmse(coefs, d, test_x, test_y, test_count);
}
printf("train_count %d\n", train_count);
printf("test_count %d\n", test_count);
for (int d = 1; d <= MAX_DEGREE; d++) printf("train_rmse_d%d %.6f\n", d, train_error[d]);
for (int d = 1; d <= MAX_DEGREE; d++) printf("test_rmse_d%d %.6f\n", d, test_error[d]);
return 0;
}We keep the training and test drives in separate parallel arrays of
doubles, then fit each degree from the training arrays alone. fit adds
up the normal equations into one fixed-size augmented matrix, a 2D array
with b as its last column, and solve runs Gaussian elimination with
partial pivoting on it.
C can’t assign one array row to another, so the pivot swap trades the two
rows entry by entry. Everything beyond pow, fabs, and sqrt from math.h is a
plain loop.
// Train and test error for polynomials of degree 1 to 8, using only the C++
// standard library.
//
// Each polynomial is fitted by least squares to the training drives only:
// we build the normal equations and solve them by Gaussian elimination.
// Then we measure its root-mean-square error on the training drives and on
// the test drives it never saw.
#include <algorithm>
#include <cmath>
#include <fstream>
#include <functional>
#include <print>
#include <sstream>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
constexpr int MAX_DEGREE = 8;
struct Point {
double x; // speed
double y; // litres
};
// Each row of a Matrix is its own vector, so rows can differ in length and
// be swapped as whole objects.
using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;
// Split one line at every comma.
std::vector<std::string> split(const std::string& line) {
std::vector<std::string> fields;
std::stringstream stream(line);
std::string field;
while (std::getline(stream, field, ',')) {
fields.push_back(field);
}
return fields;
}
// The position of a named column in the header.
std::size_t column(const std::vector<std::string>& header, const std::string& name) {
auto found = std::find(header.begin(), header.end(), name);
if (found == header.end()) {
throw std::runtime_error("no column " + name);
}
return found - header.begin();
}
// Sort the rows into training and test points by their split column.
void read(const std::string& path, std::vector<Point>& train, std::vector<Point>& test) {
std::ifstream file(path);
std::string line;
std::getline(file, line);
if (!line.empty() && line.back() == '\r') line.pop_back();
std::vector<std::string> header = split(line);
std::size_t speed = column(header, "speed");
std::size_t litres = column(header, "litres");
std::size_t split_col = column(header, "split");
while (std::getline(file, line)) {
if (!line.empty() && line.back() == '\r') line.pop_back();
if (line.empty()) {
continue;
}
std::vector<std::string> fields = split(line);
Point p{std::stod(fields.at(speed)), std::stod(fields.at(litres))};
(fields.at(split_col) == "train" ? train : test).push_back(p);
}
}
// Solve a x = b by Gaussian elimination with partial pivoting.
Vector solve(Matrix a, const Vector& b) {
std::size_t n = b.size();
for (std::size_t i = 0; i < n; i++) a[i].push_back(b[i]); // augmented matrix
for (std::size_t col = 0; col < n; col++) {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small. Swapping two vectors only trades their
// internal pointers; no entries are copied.
std::size_t pivot = col;
for (std::size_t r = col + 1; r < n; r++) {
if (std::abs(a[r][col]) > std::abs(a[pivot][col])) pivot = r;
}
std::swap(a[col], a[pivot]);
for (std::size_t r = col + 1; r < n; r++) {
double factor = a[r][col] / a[col][col];
for (std::size_t c = col; c <= n; c++) a[r][c] -= factor * a[col][c];
}
}
// Back substitution, from the last unknown up.
Vector x(n, 0.0);
for (std::size_t r = n; r-- > 0;) {
double s = a[r][n];
for (std::size_t c = r + 1; c < n; c++) s -= a[r][c] * x[c];
x[r] = s / a[r][r];
}
return x;
}
// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
// Z is 1, z_i, z_i^2, ..., z_i^degree.
Vector fit(const std::vector<Point>& points, int degree,
const std::function<double(double)>& scale) {
std::size_t size = degree + 1;
Matrix a(size, Vector(size, 0.0));
Vector b(size, 0.0);
for (const Point& p : points) {
double z = scale(p.x);
for (std::size_t i = 0; i < size; i++) {
b[i] += p.y * std::pow(z, i);
for (std::size_t j = 0; j < size; j++) a[i][j] += std::pow(z, i + j);
}
}
return solve(a, b);
}
double rmse(const Vector& coefs, const std::vector<Point>& points,
const std::function<double(double)>& scale) {
double total = 0.0;
for (const Point& p : points) {
double z = scale(p.x), predicted = 0.0;
for (std::size_t i = 0; i < coefs.size(); i++) predicted += coefs[i] * std::pow(z, i);
total += (predicted - p.y) * (predicted - p.y);
}
return std::sqrt(total / points.size());
}
int main(int argc, char* argv[]) {
if (argc != 2) {
std::println(stderr, "usage: program data.csv");
return 1;
}
std::vector<Point> train, test;
read(argv[1], train, test);
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power
// (the largest entry at degree 8) would swamp the small ones and make
// the equations hard to solve accurately. minmax_element finds both
// ends in one pass.
auto [low, high] = std::minmax_element(
train.begin(), train.end(), [](const Point& p, const Point& q) { return p.x < q.x; });
double mid = (low->x + high->x) / 2, half = (high->x - low->x) / 2;
auto scale = [mid, half](double x) { return (x - mid) / half; };
std::println("train_count {}", train.size());
std::println("test_count {}", test.size());
std::vector<std::pair<double, double>> errors;
for (int degree = 1; degree <= MAX_DEGREE; degree++) {
Vector coefs = fit(train, degree, scale);
errors.emplace_back(rmse(coefs, train, scale), rmse(coefs, test, scale));
}
for (int d = 1; d <= MAX_DEGREE; d++) {
std::println("train_rmse_d{} {:.6f}", d, errors[d - 1].first);
}
for (int d = 1; d <= MAX_DEGREE; d++) {
std::println("test_rmse_d{} {:.6f}", d, errors[d - 1].second);
}
return 0;
}We store each drive as a Point struct and the normal equations as a
vector of vectors, one vector per row. solve appends b to each row to
make the augmented matrix, then runs Gaussian elimination with partial
pivoting.
The pivot swap is a single std::swap of two row vectors, which trades
their internal buffers rather than copying entries. The training range
comes from one std::minmax_element pass. The C++ standard library has
no linear solver or polynomial fit, so we skip nothing beyond the
arithmetic in <cmath>.
*> Train and test error for polynomials of degree 1 to 8, in
*> COBOL.
*>
*> Each polynomial is fitted by least squares to the training
*> drives only: we build the normal equations and solve them by
*> Gaussian elimination. Then we measure its root-mean-square
*> error on the training drives and on the test drives it never
*> saw.
*> Values are fixed-point decimals carried to 18 decimal places.
IDENTIFICATION DIVISION.
PROGRAM-ID. TRAIN-AND-TEST.
ENVIRONMENT DIVISION.
INPUT-OUTPUT SECTION.
FILE-CONTROL.
SELECT DATA-FILE ASSIGN TO DYNAMIC WS-FILE-NAME
ORGANIZATION IS LINE SEQUENTIAL.
DATA DIVISION.
FILE SECTION.
FD DATA-FILE.
01 DATA-LINE PIC X(200).
WORKING-STORAGE SECTION.
01 WS-FILE-NAME PIC X(256).
01 WS-EOF PIC X VALUE "N".
88 END-OF-FILE VALUE "Y".
*> One CSV line split into fields.
01 WS-FIELDS.
05 WS-FIELD PIC X(40) OCCURS 20 TIMES.
01 WS-FIELD-COUNT PIC 99.
01 WS-POINTER PIC 999.
01 WS-LINE-LENGTH PIC 999.
*> Which CSV column holds each value, found from the header.
01 WS-SPEED-COL PIC 99 VALUE 0.
01 WS-LITRES-COL PIC 99 VALUE 0.
01 WS-SPLIT-COL PIC 99 VALUE 0.
*> Every drive: its speed, scaled speed Z, litres, and split.
01 WS-N PIC 9(4) VALUE 0.
01 WS-DRIVES.
05 WS-DRIVE OCCURS 1000 TIMES.
10 WS-SPEED PIC S9(6)V9(18).
10 WS-Z PIC S9(3)V9(18).
10 WS-Y PIC S9(6)V9(18).
10 WS-IS-TRAIN PIC X.
01 WS-TRAIN-COUNT PIC 9(4) VALUE 0.
01 WS-TEST-COUNT PIC 9(4) VALUE 0.
*> The training range, which maps every speed onto -1 to 1.
01 WS-LOW PIC S9(6)V9(18).
01 WS-HIGH PIC S9(6)V9(18).
01 WS-MID PIC S9(6)V9(18).
01 WS-HALF PIC S9(6)V9(18).
01 MAX-DEGREE PIC 9 VALUE 8.
01 D PIC 9.
01 N PIC 99.
*> The augmented normal equations: up to 9 rows of 10 entries, the
*> last entry of each row being b. WS-A-ROW(I) is a group item,
*> so one MOVE copies a whole row.
01 WS-A-TABLE.
05 WS-A-ROW OCCURS 9 TIMES.
10 WS-A PIC S9(9)V9(18) OCCURS 10 TIMES.
01 WS-SPARE-ROW.
05 FILLER PIC S9(9)V9(18) OCCURS 10 TIMES.
*> The powers z^0 to z^16 of one drive's scaled speed: WS-POW(E)
*> holds z^(E - 1), because tables count from 1.
01 WS-POW-TABLE.
05 WS-POW PIC S9(3)V9(18) OCCURS 17 TIMES.
*> The fitted coefficients, and each degree's two errors.
01 WS-X-TABLE.
05 WS-X PIC S9(9)V9(18) OCCURS 9 TIMES.
01 WS-ERRORS.
05 WS-TRAIN-ERROR PIC S9(3)V9(18) OCCURS 8 TIMES.
05 WS-TEST-ERROR PIC S9(3)V9(18) OCCURS 8 TIMES.
01 WS-WHICH PIC X.
01 WS-RMSE PIC S9(3)V9(18).
01 WS-SUM PIC S9(9)V9(18).
01 WS-FACTOR PIC S9(9)V9(18).
01 WS-PREDICTED PIC S9(9)V9(18).
01 WS-POWER PIC S9(3)V9(18).
01 WS-COUNT PIC 9(4).
01 PIVOT PIC 99.
01 I PIC 99.
01 J PIC 99.
01 K PIC 99.
01 R PIC 99.
01 C PIC 99.
01 S PIC 9(4).
*> Output: a name, then the value edited to six decimal places.
01 WS-NAME PIC X(20).
01 WS-RESULT PIC S9(9)V9(18).
01 WS-RESULT-OUT PIC -(9)9.9(6).
01 WS-COUNT-OUT PIC Z(8)9.
PROCEDURE DIVISION.
MAIN.
ACCEPT WS-FILE-NAME FROM ARGUMENT-VALUE
OPEN INPUT DATA-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
END-READ
PERFORM SPLIT-LINE
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-FIELD-COUNT
EVALUATE WS-FIELD(I)
WHEN "speed" MOVE I TO WS-SPEED-COL
WHEN "litres" MOVE I TO WS-LITRES-COL
WHEN "split" MOVE I TO WS-SPLIT-COL
END-EVALUATE
END-PERFORM
IF WS-SPEED-COL = 0 OR WS-LITRES-COL = 0
OR WS-SPLIT-COL = 0
DISPLAY "missing column" UPON SYSERR
MOVE 1 TO RETURN-CODE
STOP RUN
END-IF
PERFORM UNTIL END-OF-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
NOT AT END PERFORM ADD-DRIVE
END-READ
END-PERFORM
CLOSE DATA-FILE
PERFORM SCALE-SPEEDS
PERFORM VARYING D FROM 1 BY 1 UNTIL D > MAX-DEGREE
COMPUTE N = D + 1
PERFORM BUILD-SYSTEM
PERFORM SOLVE
MOVE "Y" TO WS-WHICH
PERFORM MEASURE-RMSE
MOVE WS-RMSE TO WS-TRAIN-ERROR(D)
MOVE "N" TO WS-WHICH
PERFORM MEASURE-RMSE
MOVE WS-RMSE TO WS-TEST-ERROR(D)
END-PERFORM
MOVE WS-TRAIN-COUNT TO WS-COUNT-OUT
DISPLAY "train_count " FUNCTION TRIM(WS-COUNT-OUT)
MOVE WS-TEST-COUNT TO WS-COUNT-OUT
DISPLAY "test_count " FUNCTION TRIM(WS-COUNT-OUT)
PERFORM VARYING D FROM 1 BY 1 UNTIL D > MAX-DEGREE
MOVE SPACES TO WS-NAME
STRING "train_rmse_d" D DELIMITED BY SIZE INTO WS-NAME
MOVE WS-TRAIN-ERROR(D) TO WS-RESULT
PERFORM PRINT-RESULT
END-PERFORM
PERFORM VARYING D FROM 1 BY 1 UNTIL D > MAX-DEGREE
MOVE SPACES TO WS-NAME
STRING "test_rmse_d" D DELIMITED BY SIZE INTO WS-NAME
MOVE WS-TEST-ERROR(D) TO WS-RESULT
PERFORM PRINT-RESULT
END-PERFORM
STOP RUN.
*> Split DATA-LINE on commas into WS-FIELD(1), WS-FIELD(2), ...
SPLIT-LINE.
MOVE SPACES TO WS-FIELDS
MOVE 0 TO WS-FIELD-COUNT
MOVE 1 TO WS-POINTER
MOVE FUNCTION LENGTH(FUNCTION TRIM(DATA-LINE TRAILING))
TO WS-LINE-LENGTH
PERFORM UNTIL WS-POINTER > WS-LINE-LENGTH
ADD 1 TO WS-FIELD-COUNT
UNSTRING DATA-LINE DELIMITED BY ","
INTO WS-FIELD(WS-FIELD-COUNT)
WITH POINTER WS-POINTER
END-UNSTRING
END-PERFORM.
*> Store one drive and count it as training or test.
ADD-DRIVE.
PERFORM SPLIT-LINE
ADD 1 TO WS-N
COMPUTE WS-SPEED(WS-N) =
FUNCTION NUMVAL(WS-FIELD(WS-SPEED-COL))
COMPUTE WS-Y(WS-N) =
FUNCTION NUMVAL(WS-FIELD(WS-LITRES-COL))
IF WS-FIELD(WS-SPLIT-COL) = "train"
MOVE "Y" TO WS-IS-TRAIN(WS-N)
ADD 1 TO WS-TRAIN-COUNT
ELSE
MOVE "N" TO WS-IS-TRAIN(WS-N)
ADD 1 TO WS-TEST-COUNT
END-IF.
*> Map the training speeds onto -1 to 1. Raw speeds to the 16th
*> power (the largest entry at degree 8) would swamp the small
*> ones and need far more digits. The test speeds use the same
*> training range.
SCALE-SPEEDS.
MOVE 999999 TO WS-LOW
MOVE -999999 TO WS-HIGH
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-N
IF WS-IS-TRAIN(S) = "Y"
IF WS-SPEED(S) < WS-LOW
MOVE WS-SPEED(S) TO WS-LOW
END-IF
IF WS-SPEED(S) > WS-HIGH
MOVE WS-SPEED(S) TO WS-HIGH
END-IF
END-IF
END-PERFORM
COMPUTE WS-MID = (WS-LOW + WS-HIGH) / 2
COMPUTE WS-HALF = (WS-HIGH - WS-LOW) / 2
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-N
COMPUTE WS-Z(S) ROUNDED =
(WS-SPEED(S) - WS-MID) / WS-HALF
END-PERFORM.
*> Least squares: the coefficients solve (Z^T Z) c = Z^T y, so
*> entry (I, J) sums z^(I+J-2) over the training drives, and the
*> last column sums y z^(I-1). Each drive's powers come from
*> repeated multiplication.
BUILD-SYSTEM.
MOVE ZEROS TO WS-A-TABLE
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-N
IF WS-IS-TRAIN(S) = "Y"
MOVE 1 TO WS-POW(1)
PERFORM VARYING I FROM 2 BY 1 UNTIL I > 2 * N - 1
COMPUTE WS-POW(I) ROUNDED =
WS-POW(I - 1) * WS-Z(S)
END-PERFORM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > N
PERFORM VARYING J FROM 1 BY 1 UNTIL J > N
ADD WS-POW(I + J - 1) TO WS-A(I, J)
END-PERFORM
COMPUTE WS-A(I, N + 1) ROUNDED =
WS-A(I, N + 1) + WS-Y(S) * WS-POW(I)
END-PERFORM
END-IF
END-PERFORM.
*> Gaussian elimination with partial pivoting, then back
*> substitution into WS-X.
SOLVE.
PERFORM VARYING K FROM 1 BY 1 UNTIL K > N
*> Swap up the row with the largest entry in column K, which
*> keeps the rounding error small. Three group MOVEs swap
*> the rows through a spare one.
MOVE K TO PIVOT
PERFORM VARYING R FROM K BY 1 UNTIL R > N
IF FUNCTION ABS(WS-A(R, K))
> FUNCTION ABS(WS-A(PIVOT, K))
MOVE R TO PIVOT
END-IF
END-PERFORM
MOVE WS-A-ROW(K) TO WS-SPARE-ROW
MOVE WS-A-ROW(PIVOT) TO WS-A-ROW(K)
MOVE WS-SPARE-ROW TO WS-A-ROW(PIVOT)
*> Each lower row subtracts the multiple of row K that
*> clears its entry in column K.
PERFORM VARYING R FROM K BY 1 UNTIL R > N
IF R > K
COMPUTE WS-FACTOR ROUNDED =
WS-A(R, K) / WS-A(K, K)
PERFORM VARYING C FROM K BY 1 UNTIL C > N + 1
COMPUTE WS-A(R, C) ROUNDED =
WS-A(R, C) - WS-FACTOR * WS-A(K, C)
END-PERFORM
END-IF
END-PERFORM
END-PERFORM
*> Back substitution, from the last unknown up.
PERFORM VARYING R FROM N BY -1 UNTIL R < 1
MOVE WS-A(R, N + 1) TO WS-SUM
PERFORM VARYING C FROM R BY 1 UNTIL C > N
IF C > R
COMPUTE WS-SUM ROUNDED =
WS-SUM - WS-A(R, C) * WS-X(C)
END-IF
END-PERFORM
COMPUTE WS-X(R) ROUNDED = WS-SUM / WS-A(R, R)
END-PERFORM.
*> Root-mean-square error of the current fit over the drives
*> whose WS-IS-TRAIN matches WS-WHICH.
MEASURE-RMSE.
MOVE 0 TO WS-SUM
MOVE 0 TO WS-COUNT
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-N
IF WS-IS-TRAIN(S) = WS-WHICH
ADD 1 TO WS-COUNT
MOVE 0 TO WS-PREDICTED
MOVE 1 TO WS-POWER
PERFORM VARYING I FROM 1 BY 1 UNTIL I > N
COMPUTE WS-PREDICTED ROUNDED =
WS-PREDICTED + WS-X(I) * WS-POWER
COMPUTE WS-POWER ROUNDED = WS-POWER * WS-Z(S)
END-PERFORM
COMPUTE WS-SUM ROUNDED =
WS-SUM + (WS-PREDICTED - WS-Y(S)) ** 2
END-IF
END-PERFORM
COMPUTE WS-RMSE ROUNDED =
FUNCTION SQRT(WS-SUM / WS-COUNT).
*> Print "name value", rounded to six decimal places.
PRINT-RESULT.
COMPUTE WS-RESULT-OUT ROUNDED = WS-RESULT
DISPLAY FUNCTION TRIM(WS-NAME) " "
FUNCTION TRIM(WS-RESULT-OUT).We hold the augmented normal equations in a table of rows, each row a
group item of ten entries. That makes the pivot swap three MOVEs
through a spare row, since a group MOVE copies every entry at once.
Each drive’s powers of z come from repeated multiplication into a
small table, so we never use ** on the scaled speeds.
Values are fixed-point decimals with 18 places. The degree-8 system divides by pivots as small as 0.00006, yet all 16 errors still agree with Python’s doubles to within 1e-11. GnuCOBOL has no matrix or polynomial-fitting function to skip.
;;;; Train and test error for polynomials of degree 1 to 8, using only
;;;; standard Common Lisp.
;;;;
;;;; Each polynomial is fitted by least squares to the training drives only:
;;;; we build the normal equations and solve them by Gaussian elimination.
;;;; Then we measure its root-mean-square error on the training drives and
;;;; on the test drives it never saw. The equations sit in a 2-D array, read
;;;; and written with (aref m i j).
(defparameter *max-degree* 8)
;; Read every decimal as a double, in the data and in this file.
;; eval-when makes the setting apply while the compiler reads this file, too,
;; not only once the program runs.
(eval-when (:compile-toplevel :load-toplevel :execute)
(setf *read-default-float-format* 'double-float))
(defun split-on-comma (line)
"The fields of LINE, split at every comma."
(loop for start = 0 then (1+ end)
for end = (position #\, line :start start)
collect (subseq line start end)
while end))
(defun parse-number (text)
"The number written in TEXT, read by the Lisp reader. Binding *read-eval*
to nil stops the reader from running code written as #.(...)."
(let ((*read-eval* nil))
(coerce (read-from-string text) 'double-float)))
(defun read-drives (path)
"The rows of the CSV file at PATH as two lists of (speed . litres) pairs,
the training drives and the test drives, each in file order."
(with-open-file (in path)
(let* ((header (split-on-comma (read-line in)))
(column (lambda (name)
(or (position name header :test #'string=)
(error "no column ~a" name))))
(speed-at (funcall column "speed"))
(litres-at (funcall column "litres"))
(split-at (funcall column "split"))
(train '())
(test '()))
(loop for line = (read-line in nil)
while line
unless (string= line "")
do (let* ((fields (split-on-comma line))
(point (cons (parse-number (nth speed-at fields))
(parse-number (nth litres-at fields)))))
(if (string= (nth split-at fields) "train")
(push point train)
(push point test))))
;; push adds to the front, so reverse restores the file order.
(values (reverse train) (reverse test)))))
(defun solve (a b n)
"Solve A x = B, N equations in N unknowns, by Gaussian elimination with
partial pivoting. A is an N by N array and B a vector; x comes back as a
vector."
;; The augmented matrix: A with B as an extra last column.
(let ((m (make-array (list n (1+ n)))))
(dotimes (i n)
(dotimes (j n)
(setf (aref m i j) (aref a i j)))
(setf (aref m i n) (aref b i)))
(dotimes (col n)
;; Swap up the row with the largest entry in this column, which keeps
;; the rounding error small. A 2-D array has no row objects to swap,
;; so rotatef swaps the two rows one entry at a time.
(let ((pivot col))
(loop for r from (1+ col) below n
when (> (abs (aref m r col)) (abs (aref m pivot col)))
do (setf pivot r))
(dotimes (c (1+ n))
(rotatef (aref m col c) (aref m pivot c))))
(loop for r from (1+ col) below n
do (let ((factor (/ (aref m r col) (aref m col col))))
(loop for c from col to n
do (decf (aref m r c) (* factor (aref m col c)))))))
;; Back substitution, from the last unknown up.
(let ((x (make-array n :initial-element 0d0)))
(loop for r from (1- n) downto 0
do (setf (aref x r)
(/ (- (aref m r n)
(loop for c from (1+ r) below n
sum (* (aref m r c) (aref x c))))
(aref m r r))))
x)))
(defun fit (points degree)
"Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
Z is 1, z_i, z_i^2, ..., z_i^degree. POINTS holds (z . y) pairs."
(let* ((size (1+ degree))
(a (make-array (list size size) :initial-element 0d0))
(b (make-array size :initial-element 0d0)))
(loop for (z . y) in points
do (dotimes (i size)
(incf (aref b i) (* y (expt z i)))
(dotimes (j size)
(incf (aref a i j) (expt z (+ i j))))))
(solve a b size)))
(defun rmse (coefs points)
(let ((total (loop for (z . y) in points
sum (let ((predicted (loop for c across coefs
for i from 0
sum (* c (expt z i)))))
(expt (- predicted y) 2)))))
(sqrt (/ total (length points)))))
(multiple-value-bind (train test) (read-drives (second sb-ext:*posix-argv*))
;; Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
;; largest entry at degree 8) would swamp the small ones and make the
;; equations hard to solve accurately.
(let* ((low (reduce #'min (mapcar #'car train)))
(high (reduce #'max (mapcar #'car train)))
(mid (/ (+ low high) 2))
(half (/ (- high low) 2))
(scale (lambda (points)
(mapcar (lambda (p) (cons (/ (- (car p) mid) half) (cdr p)))
points)))
(train-z (funcall scale train))
(test-z (funcall scale test))
(fits (loop for degree from 1 to *max-degree*
collect (fit train-z degree))))
;; ~d prints an integer and ~,6f a number with six decimals.
(format t "train_count ~d~%" (length train))
(format t "test_count ~d~%" (length test))
(loop for coefs in fits
for degree from 1
do (format t "train_rmse_d~d ~,6f~%" degree (rmse coefs train-z)))
(loop for coefs in fits
for degree from 1
do (format t "test_rmse_d~d ~,6f~%" degree (rmse coefs test-z)))))The normal equations sit in a 2-D array built with make-array and read
with aref. A 2-D array has no row objects to swap, so partial pivoting
swaps the two rows one entry at a time with rotatef, which exchanges the
values in two places. Elimination updates entries in place with decf.
Powers come from expt, and the sums are loop ... sum clauses. We read
every decimal as a double, since Common Lisp would otherwise read them as
single floats. No solver or fitting function is called: elimination, back
substitution, and the error sums are all written out.
# Train and test error for polynomials of degree 1 to 8, using only
# Elixir's standard library.
#
# Each polynomial is fitted by least squares to the training drives only: we
# build the normal equations and solve them by Gaussian elimination. Then we
# measure its root-mean-square error on the training drives and on the test
# drives it never saw. A matrix is a list of row lists, and nothing is
# changed in place.
defmodule TrainAndTest do
@max_degree 8
def max_degree, do: @max_degree
# Every drive as {split, speed, litres}, with each column found by its
# header name. Float.parse returns {number, rest}; "42" parses as 42.0 too.
def read_drives(path) do
[header | rows] =
path
|> File.read!()
|> String.split("\n", trim: true)
|> Enum.map(&String.split(&1, ","))
number = fn text -> {value, ""} = Float.parse(text); value end
for fields <- rows do
row = header |> Enum.zip(fields) |> Map.new()
{row["split"], number.(row["speed"]), number.(row["litres"])}
end
end
# Solve a x = b by Gaussian elimination with partial pivoting. Each row of
# the augmented matrix is a row of a with its entry of b on the end.
def solve(a, b) do
a
|> Enum.zip_with(b, fn row, v -> row ++ [v] end)
|> eliminate()
|> back_substitute()
end
# Take out the row with the largest first entry, which keeps the rounding
# error small, and subtract a multiple of it from every other row so their
# first entries become 0. Those rows, without their first column, are a
# smaller system of the same kind, so we recurse. Each row of the triangle
# is one entry shorter than the last.
defp eliminate([]), do: []
defp eliminate(rows) do
[p | pivot_rest] = pivot = Enum.max_by(rows, fn [x | _] -> abs(x) end)
reduced =
for [r | row_rest] <- List.delete(rows, pivot) do
factor = r / p
Enum.zip_with(pivot_rest, row_rest, fn pc, rc -> rc - factor * pc end)
end
[pivot | eliminate(reduced)]
end
# Back substitution, from the last unknown up. Each triangle row is
# [coefficient, entries for the later unknowns ..., right-hand side].
defp back_substitute([]), do: []
defp back_substitute([[coefficient | rest] | rows]) do
later = back_substitute(rows)
{entries, [rhs]} = Enum.split(rest, -1)
[(rhs - dot(entries, later)) / coefficient | later]
end
# zip_reduce walks two lists side by side, carrying a running total.
def dot(u, v), do: Enum.zip_reduce(u, v, 0.0, fn a, b, acc -> acc + a * b end)
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
# of Z is 1, z_i, z_i^2, ..., z_i^degree. points holds {z, y} pairs.
def fit(points, degree) do
a =
for i <- 0..degree do
for j <- 0..degree, do: Enum.sum(for {z, _} <- points, do: z ** (i + j))
end
b = for i <- 0..degree, do: Enum.sum(for {z, y} <- points, do: y * z ** i)
solve(a, b)
end
def rmse(coefs, points) do
predict = fn z ->
coefs |> Enum.with_index() |> Enum.map(fn {c, i} -> c * z ** i end) |> Enum.sum()
end
total = Enum.sum(for {z, y} <- points, do: (predict.(z) - y) ** 2)
:math.sqrt(total / length(points))
end
def say(name, value), do: IO.puts("#{name} #{:erlang.float_to_binary(value, decimals: 6)}")
end
[path] = System.argv()
drives = TrainAndTest.read_drives(path)
train = for {"train", x, y} <- drives, do: {x, y}
test = for {"test", x, y} <- drives, do: {x, y}
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately.
{low, high} = train |> Enum.map(fn {x, _} -> x end) |> Enum.min_max()
mid = (low + high) / 2
half = (high - low) / 2
scale = fn points -> for {x, y} <- points, do: {(x - mid) / half, y} end
train_z = scale.(train)
test_z = scale.(test)
IO.puts("train_count #{length(train)}")
IO.puts("test_count #{length(test)}")
fits = for d <- 1..TrainAndTest.max_degree(), do: {d, TrainAndTest.fit(train_z, d)}
for {d, coefs} <- fits, do: TrainAndTest.say("train_rmse_d#{d}", TrainAndTest.rmse(coefs, train_z))
for {d, coefs} <- fits, do: TrainAndTest.say("test_rmse_d#{d}", TrainAndTest.rmse(coefs, test_z))A matrix is a list of row lists, and nothing is changed in place, so the
elimination is recursive. Enum.max_by picks the row with the largest
first entry (partial pivoting), List.delete removes it, and a
comprehension subtracts a multiple of it from every other row. The smaller
system that’s left goes back into eliminate.
Back substitution splits off each row’s right-hand side with
Enum.split(rest, -1). Powers use **, and each normal-equation entry
is an Enum.sum over the training drives. No solver or fitting function
is called.
! Train and test error for polynomials of degree 1 to 8, using only standard
! Fortran.
!
! Each polynomial is fitted by least squares to the training drives only: we
! build the normal equations and solve them by Gaussian elimination. Then we
! measure its root-mean-square error on the training drives and on the test
! drives it never saw.
program train_and_test
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
integer, parameter :: max_degree = 8
real(dp), allocatable :: speed(:), litres(:), train_z(:), train_y(:), test_z(:), test_y(:)
logical, allocatable :: is_train(:)
real(dp) :: mid, half, coefs(0:max_degree), train_rmse(max_degree), test_rmse(max_degree)
character(len=20) :: name
integer :: d
call read_drives(speed, litres, is_train)
! Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
! largest entry at degree 8) would swamp the small ones and make the
! equations hard to solve accurately. pack keeps the elements where the
! mask is .true.
train_y = pack(litres, is_train)
test_y = pack(litres, .not. is_train)
train_z = pack(speed, is_train)
mid = (minval(train_z) + maxval(train_z)) / 2
half = (maxval(train_z) - minval(train_z)) / 2
train_z = (train_z - mid) / half
test_z = (pack(speed, .not. is_train) - mid) / half
do d = 1, max_degree
coefs(:d) = fit(train_z, train_y, d)
train_rmse(d) = rmse(coefs(:d), train_z, train_y)
test_rmse(d) = rmse(coefs(:d), test_z, test_y)
end do
print "(a, 1x, i0)", "train_count", size(train_y)
print "(a, 1x, i0)", "test_count", size(test_y)
do d = 1, max_degree
write(name, "(a, i0)") "train_rmse_d", d
call say(trim(name), train_rmse(d))
end do
do d = 1, max_degree
write(name, "(a, i0)") "test_rmse_d", d
call say(trim(name), test_rmse(d))
end do
contains
! The design matrix: row i is 1, z_i, z_i^2, ..., z_i^degree. Its columns
! count from 0, so column j holds z to the power j.
function powers(z, degree) result(m)
real(dp), intent(in) :: z(:)
integer, intent(in) :: degree
real(dp) :: m(size(z), 0:degree)
integer :: j
do j = 0, degree
m(:, j) = z**j
end do
end function powers
! Least squares: the coefficients c solve (Z^T Z) c = Z^T y.
function fit(z, y, degree) result(c)
real(dp), intent(in) :: z(:), y(:)
integer, intent(in) :: degree
real(dp) :: c(0:degree), zm(size(z), 0:degree)
zm = powers(z, degree)
c = solve(matmul(transpose(zm), zm), matmul(transpose(zm), y))
end function fit
real(dp) function rmse(c, z, y)
real(dp), intent(in) :: c(0:), z(:), y(:)
rmse = sqrt(sum((matmul(powers(z, size(c) - 1), c) - y)**2) / size(y))
end function rmse
! Solve a x = b by Gaussian elimination with partial pivoting.
function solve(a, b) result(x)
real(dp), intent(in) :: a(:, :), b(:)
real(dp) :: x(size(b)), m(size(b), size(b) + 1), row(size(b) + 1), factor
integer :: n, col, pivot, r
n = size(b)
! The augmented matrix: a with b as an extra last column.
m(:, :n) = a
m(:, n + 1) = b
do col = 1, n
! Swap up the row with the largest entry in this column, which keeps
! the rounding error small. maxloc gives the position within the
! section col:n, so we shift it back to a row number.
pivot = maxloc(abs(m(col:, col)), 1) + col - 1
row = m(col, :)
m(col, :) = m(pivot, :)
m(pivot, :) = row
do r = col + 1, n
factor = m(r, col) / m(col, col)
! One statement updates the whole rest of the row.
m(r, col:) = m(r, col:) - factor * m(col, col:)
end do
end do
! Back substitution, from the last unknown up.
do r = n, 1, -1
x(r) = (m(r, n + 1) - dot_product(m(r, r + 1:n), x(r + 1:n))) / m(r, r)
end do
end function solve
subroutine read_drives(speed, litres, is_train)
real(dp), allocatable, intent(out) :: speed(:), litres(:)
logical, allocatable, intent(out) :: is_train(:)
character(len=256) :: path, line, text
integer :: unit, status, speed_col, litres_col, split_col
real(dp) :: x, y
call get_command_argument(1, path)
open(newunit=unit, file=path, status="old", action="read")
read(unit, "(a)") line
speed_col = field_index(line, "speed")
litres_col = field_index(line, "litres")
split_col = field_index(line, "split")
allocate(speed(0), litres(0), is_train(0))
do
read(unit, "(a)", iostat=status) line
if (status /= 0) exit
if (len_trim(line) == 0) cycle
! An internal read parses a number. It must read from a variable,
! not straight from the result of field(), so we copy into text.
text = field(line, speed_col)
read(text, *) x
text = field(line, litres_col)
read(text, *) y
speed = [speed, x] ! grows the array by one
litres = [litres, y]
is_train = [is_train, field(line, split_col) == "train"]
end do
close(unit)
end subroutine read_drives
! The position of a name in a comma-separated header line.
integer function field_index(header, name)
character(len=*), intent(in) :: header, name
integer :: i, fields
fields = 1
do i = 1, len_trim(header)
if (header(i:i) == ",") fields = fields + 1
end do
do field_index = 1, fields
if (field(header, field_index) == name) return
end do
error stop "column not found: " // name
end function field_index
! The k-th comma-separated field of a line.
function field(line, k) result(text)
character(len=*), intent(in) :: line
integer, intent(in) :: k
character(len=:), allocatable :: text
integer :: start, i, comma
start = 1
do i = 1, k - 1
start = start + index(line(start:), ",")
end do
comma = index(line(start:), ",")
if (comma == 0) then
text = trim(line(start:))
else
text = line(start:start + comma - 2)
end if
end function field
! Print "name value" with six decimals. The f0.6 format drops the zero
! before the point (0.5 prints as .500000), so we put it back.
subroutine say(name, x)
character(len=*), intent(in) :: name
real(dp), intent(in) :: x
character(len=40) :: text
write(text, "(f0.6)") x
if (text(1:1) == ".") text = "0" // text
if (text(1:2) == "-.") text = "-0" // text(2:)
print "(a, 1x, a)", name, trim(text)
end subroutine say
end program train_and_testFortran’s whole-array operations carry most of the work. The design matrix
has a column of z to each power, and matmul(transpose(zm), zm) builds the
normal equations; multiplying matrices is a general tool, not the fit
itself. The Gaussian elimination is written out: maxloc finds the pivot
row, a temporary array swaps two rows, and one statement updates the rest
of a row at once.
Fortran’s usual solver library, LAPACK, isn’t part of the language standard, and it isn’t installed in this image. Coefficient arrays count from 0, matching the powers.
// Train and test error for polynomials of degree 1 to 8, using only the
// standard library.
//
// Each polynomial is fitted by least squares to the training drives only:
// we build the normal equations and solve them by Gaussian elimination. Then
// we measure its root-mean-square error on the training drives and on the
// test drives it never saw.
package main
import (
"encoding/csv"
"fmt"
"math"
"os"
"strconv"
)
const maxDegree = 8
type point struct{ speed, litres float64 }
func read(path string) (train, test []point, err error) {
f, err := os.Open(path)
if err != nil {
return nil, nil, err
}
defer f.Close()
records, err := csv.NewReader(f).ReadAll()
if err != nil {
return nil, nil, err
}
// Find each column by name in the header.
iSpeed, iLitres, iSplit := -1, -1, -1
for i, name := range records[0] {
switch name {
case "speed":
iSpeed = i
case "litres":
iLitres = i
case "split":
iSplit = i
}
}
if iSpeed < 0 || iLitres < 0 || iSplit < 0 {
return nil, nil, fmt.Errorf("missing a column (speed, litres, split)")
}
for _, record := range records[1:] {
speed, err := strconv.ParseFloat(record[iSpeed], 64)
if err != nil {
return nil, nil, err
}
litres, err := strconv.ParseFloat(record[iLitres], 64)
if err != nil {
return nil, nil, err
}
p := point{speed, litres}
if record[iSplit] == "train" {
train = append(train, p)
} else {
test = append(test, p)
}
}
return train, test, nil
}
// solve solves a x = b by Gaussian elimination with partial pivoting.
func solve(a [][]float64, b []float64) []float64 {
n := len(b)
// The augmented matrix, a fresh copy so a and b are untouched.
m := make([][]float64, n)
for i := range a {
m[i] = append(append([]float64{}, a[i]...), b[i])
}
for col := 0; col < n; col++ {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small. Rows are slices, so this swaps headers,
// not the numbers inside.
pivot := col
for r := col + 1; r < n; r++ {
if math.Abs(m[r][col]) > math.Abs(m[pivot][col]) {
pivot = r
}
}
m[col], m[pivot] = m[pivot], m[col]
for r := col + 1; r < n; r++ {
factor := m[r][col] / m[col][col]
for c := col; c <= n; c++ {
m[r][c] -= factor * m[col][c]
}
}
}
// Back substitution, from the last unknown up.
x := make([]float64, n)
for r := n - 1; r >= 0; r-- {
total := m[r][n]
for c := r + 1; c < n; c++ {
total -= m[r][c] * x[c]
}
x[r] = total / m[r][r]
}
return x
}
// fit finds least-squares coefficients: c solves (Z^T Z) c = Z^T y, where
// row i of Z is 1, z_i, z_i^2, ..., z_i^degree.
func fit(points []point, degree int, scale func(float64) float64) []float64 {
size := degree + 1
a := make([][]float64, size)
for i := range a {
a[i] = make([]float64, size)
}
b := make([]float64, size)
for _, p := range points {
z := scale(p.speed)
for i := 0; i < size; i++ {
b[i] += p.litres * math.Pow(z, float64(i))
for j := 0; j < size; j++ {
a[i][j] += math.Pow(z, float64(i+j))
}
}
}
return solve(a, b)
}
func rmse(coefs []float64, points []point, scale func(float64) float64) float64 {
total := 0.0
for _, p := range points {
z := scale(p.speed)
predicted := 0.0
for i, c := range coefs {
predicted += c * math.Pow(z, float64(i))
}
total += (predicted - p.litres) * (predicted - p.litres)
}
return math.Sqrt(total / float64(len(points)))
}
func main() {
train, test, err := read(os.Args[1])
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power
// (the largest entry at degree 8) would swamp the small ones and make
// the equations hard to solve accurately.
low, high := train[0].speed, train[0].speed
for _, p := range train {
low = math.Min(low, p.speed)
high = math.Max(high, p.speed)
}
mid, half := (low+high)/2, (high-low)/2
scale := func(x float64) float64 { return (x - mid) / half }
fmt.Printf("train_count %d\n", len(train))
fmt.Printf("test_count %d\n", len(test))
var trainErrors, testErrors [maxDegree + 1]float64 // index 0 unused
for degree := 1; degree <= maxDegree; degree++ {
coefs := fit(train, degree, scale)
trainErrors[degree] = rmse(coefs, train, scale)
testErrors[degree] = rmse(coefs, test, scale)
}
for degree := 1; degree <= maxDegree; degree++ {
fmt.Printf("train_rmse_d%d %.6f\n", degree, trainErrors[degree])
}
for degree := 1; degree <= maxDegree; degree++ {
fmt.Printf("test_rmse_d%d %.6f\n", degree, testErrors[degree])
}
}encoding/csv reads the file, and each drive becomes a point struct
sorted into train or test. The matrix is a slice of row slices, so a
parallel assignment, m[a], m[b] = m[b], m[a], swaps two rows by swapping
slice headers, without copying the numbers inside. Elimination and back
substitution are plain nested for loops.
Powers come from math.Pow, which takes and returns float64, so each
exponent is wrapped in float64(...). Go’s standard library has no
linear algebra package, so there was no solver to skip; gonum has one, but
it isn’t installed.
-- Train and test error for polynomials of degree 1 to 8, using only
-- Haskell's base library.
--
-- Each polynomial is fitted by least squares to the training drives only:
-- we build the normal equations and solve them by Gaussian elimination. Then
-- we measure its root-mean-square error on the training drives and on the
-- test drives it never saw. A matrix is a list of rows.
import Control.Monad (forM_)
import Data.List (elemIndex, maximumBy)
import Data.Ord (comparing)
import System.Environment (getArgs)
import Text.Printf (printf)
maxDegree :: Int
maxDegree = 8
-- Split a line at every comma. break returns the text before the first comma
-- and the rest, which still starts with the comma.
splitOnComma :: String -> [String]
splitOnComma line = case break (== ',') line of
(field, ',' : rest) -> field : splitOnComma rest
(field, _) -> [field]
-- Every drive as (split, (speed, litres)), in file order.
readDrives :: String -> [(String, (Double, Double))]
readDrives text = case lines text of
[] -> error "empty file"
headerLine : rows ->
[ (field "split", (read (field "speed"), read (field "litres")))
| row <- rows
, not (null row)
, let fields = splitOnComma row
field name = case elemIndex name header of
Just index -> fields !! index
Nothing -> error ("no column " ++ name)
]
where
header = splitOnComma headerLine
-- Solve a x = b by Gaussian elimination with partial pivoting. Each row of
-- the augmented matrix is a row of a with its entry of b on the end.
solve :: [[Double]] -> [Double] -> [Double]
solve a b = backSubstitute (eliminate (zipWith (\row v -> row ++ [v]) a b))
-- Move the row with the largest first entry to the top, which keeps the
-- rounding error small, and subtract a multiple of it from every other row
-- so their first entries become 0. Those rows, without their first column,
-- are a smaller system of the same kind, so we recurse. The result is the
-- upper triangle, one row shorter each time.
eliminate :: [[Double]] -> [[Double]]
eliminate [] = []
eliminate rows = pivot : eliminate (map reduce others)
where
leading row = case row of
x : _ -> abs x
[] -> 0
pivotIndex = fst (maximumBy (comparing (leading . snd)) (zip [0 :: Int ..] rows))
pivot = rows !! pivotIndex
others = [row | (i, row) <- zip [0 ..] rows, i /= pivotIndex]
reduce row = case (pivot, row) of
(p : ps, r : rs) -> zipWith (\pc rc -> rc - (r / p) * pc) ps rs
_ -> error "empty row"
-- Back substitution, from the last unknown up: each triangle row is
-- [coefficient, entries for the later unknowns ..., right-hand side].
backSubstitute :: [[Double]] -> [Double]
backSubstitute [] = []
backSubstitute (row : rows) = case row of
coefficient : rest ->
let later = backSubstitute rows
rhs = last rest
known = sum (zipWith (*) (init rest) later)
in (rhs - known) / coefficient : later
[] -> error "empty row"
-- Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
-- Z is 1, z_i, z_i^2, ..., z_i^degree.
fit :: [(Double, Double)] -> Int -> [Double]
fit points degree = solve a b
where
powers = [0 .. degree]
a = [[sum [z ^ (i + j) | (z, _) <- points] | j <- powers] | i <- powers]
b = [sum [y * z ^ i | (z, y) <- points] | i <- powers]
rmse :: [Double] -> [(Double, Double)] -> Double
rmse coefs points = sqrt (total / fromIntegral (length points))
where
predict z = sum [c * z ^ i | (c, i) <- zip coefs [0 :: Int ..]]
total = sum [(predict z - y) ^ (2 :: Int) | (z, y) <- points]
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let drives = readDrives text
train = [p | ("train", p) <- drives]
test = [p | ("test", p) <- drives]
-- Map the training speeds onto -1 to 1. Raw speeds to the 16th power
-- (the largest entry at degree 8) would swamp the small ones and make
-- the equations hard to solve accurately.
low = minimum (map fst train)
high = maximum (map fst train)
mid = (low + high) / 2
half = (high - low) / 2
scaled points = [((x - mid) / half, y) | (x, y) <- points]
fits = [fit (scaled train) d | d <- [1 .. maxDegree]]
degrees = zip [1 :: Int ..] fits
printf "train_count %d\n" (length train)
printf "test_count %d\n" (length test)
forM_ degrees $ \(d, coefs) ->
printf "train_rmse_d%d %.6f\n" d (rmse coefs (scaled train))
forM_ degrees $ \(d, coefs) ->
printf "test_rmse_d%d %.6f\n" d (rmse coefs (scaled test))A matrix is a list of rows, and the elimination is recursive rather than a
loop. eliminate takes out the row with the largest first entry (partial
pivoting), subtracts a multiple of it from every other row, and calls
itself on the smaller system that’s left. Back substitution then works up
the triangle, one unknown per row.
Each normal-equation entry is a list comprehension summing z ^ (i + j)
over the training drives. Haskell’s ^ takes a whole-number exponent and
multiplies, so no ** is needed. Nothing here calls a fitting or solving
function: every step of the least-squares fit is written out.
// Train and test error for polynomials of degree 1 to 8, using only the Java
// standard library.
//
// Each polynomial is fitted by least squares to the training drives only: we
// build the normal equations and solve them by Gaussian elimination. Then we
// measure its root-mean-square error on the training drives and on the test
// drives it never saw.
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
import java.util.function.DoubleUnaryOperator;
class Main {
static final int MAX_DEGREE = 8;
// A record is a small immutable class; Java writes its constructor and
// accessors, so a point's speed is p.x().
record Point(double x, double y) {}
// Solve a x = b by Gaussian elimination with partial pivoting.
static double[] solve(double[][] a, double[] b) {
int n = b.length;
double[][] m = new double[n][];
for (int i = 0; i < n; i++) {
m[i] = Arrays.copyOf(a[i], n + 1); // the row, plus room for b
m[i][n] = b[i];
}
for (int col = 0; col < n; col++) {
// Swap up the row with the largest entry in this column, which
// keeps the rounding error small. A double[][] is an array of row
// arrays, so the swap moves two references, not the entries.
int pivot = col;
for (int r = col + 1; r < n; r++) {
if (Math.abs(m[r][col]) > Math.abs(m[pivot][col])) {
pivot = r;
}
}
double[] row = m[col];
m[col] = m[pivot];
m[pivot] = row;
for (int r = col + 1; r < n; r++) {
double factor = m[r][col] / m[col][col];
for (int c = col; c <= n; c++) {
m[r][c] -= factor * m[col][c];
}
}
}
// Back substitution, from the last unknown up.
double[] x = new double[n];
for (int r = n - 1; r >= 0; r--) {
double s = m[r][n];
for (int c = r + 1; c < n; c++) {
s -= m[r][c] * x[c];
}
x[r] = s / m[r][r];
}
return x;
}
// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
// of Z is 1, z_i, z_i^2, ..., z_i^degree.
static double[] fit(List<Point> points, int degree, DoubleUnaryOperator scale) {
int size = degree + 1;
double[][] a = new double[size][size]; // new arrays start at 0.0
double[] b = new double[size];
for (Point p : points) {
double z = scale.applyAsDouble(p.x());
for (int i = 0; i < size; i++) {
b[i] += p.y() * Math.pow(z, i);
for (int j = 0; j < size; j++) {
a[i][j] += Math.pow(z, i + j);
}
}
}
return solve(a, b);
}
static double rmse(double[] coefs, List<Point> points, DoubleUnaryOperator scale) {
double total = 0;
for (Point p : points) {
double z = scale.applyAsDouble(p.x());
double predicted = 0;
for (int i = 0; i < coefs.length; i++) {
predicted += coefs[i] * Math.pow(z, i);
}
total += (predicted - p.y()) * (predicted - p.y());
}
return Math.sqrt(total / points.size());
}
// Reading a file can fail with an IOException, and Java makes us say so
// with "throws".
public static void main(String[] args) throws IOException {
List<String> lines = Files.readAllLines(Path.of(args[0]));
List<String> header = Arrays.asList(lines.get(0).split(","));
int speed = header.indexOf("speed");
int litres = header.indexOf("litres");
int split = header.indexOf("split");
List<Point> train = new ArrayList<>();
List<Point> test = new ArrayList<>();
for (String line : lines.subList(1, lines.size())) {
if (line.isBlank()) {
continue;
}
String[] fields = line.split(",");
Point p = new Point(Double.parseDouble(fields[speed]), Double.parseDouble(fields[litres]));
(fields[split].equals("train") ? train : test).add(p);
}
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power
// (the largest entry at degree 8) would swamp the small ones and make
// the equations hard to solve accurately.
double low = train.stream().mapToDouble(Point::x).min().orElseThrow();
double high = train.stream().mapToDouble(Point::x).max().orElseThrow();
double mid = (low + high) / 2;
double half = (high - low) / 2;
DoubleUnaryOperator scale = x -> (x - mid) / half;
double[] trainError = new double[MAX_DEGREE + 1];
double[] testError = new double[MAX_DEGREE + 1];
for (int d = 1; d <= MAX_DEGREE; d++) {
double[] coefs = fit(train, d, scale);
trainError[d] = rmse(coefs, train, scale);
testError[d] = rmse(coefs, test, scale);
}
// printf follows the computer's language settings unless told
// otherwise; Locale.ROOT keeps the decimal point a point.
Locale root = Locale.ROOT;
System.out.printf(root, "train_count %d%n", train.size());
System.out.printf(root, "test_count %d%n", test.size());
for (int d = 1; d <= MAX_DEGREE; d++) {
System.out.printf(root, "train_rmse_d%d %.6f%n", d, trainError[d]);
}
for (int d = 1; d <= MAX_DEGREE; d++) {
System.out.printf(root, "test_rmse_d%d %.6f%n", d, testError[d]);
}
}
}We read each drive into a Point record and sort it into a training or
test list by its split column. fit adds up the normal equations in a
double[][], and solve copies each row one entry longer to hold b,
then runs Gaussian elimination with partial pivoting.
Because a double[][] is an array of row arrays, the pivot swap trades
two references in three lines and copies no entries. Beyond Math.pow, Math.abs, and
Math.sqrt, every step is a loop we wrote.
// Train and test error for polynomials of degree 1 to 8, using only Node's
// built-in modules.
//
// Each polynomial is fitted by least squares to the training drives only:
// we build the normal equations and solve them by Gaussian elimination. Then
// we measure its root-mean-square error on the training drives and on the
// test drives it never saw.
import { readFileSync } from "node:fs";
const MAX_DEGREE = 8;
// No field here is quoted, so we split on commas and find columns by name.
// Each drive becomes [speed, litres], sorted into train or test.
function read(path) {
const [header, ...lines] = readFileSync(path, "utf8").trim().split("\n");
const columns = header.split(",");
const iSpeed = columns.indexOf("speed");
const iLitres = columns.indexOf("litres");
const iSplit = columns.indexOf("split");
const train = [], test = [];
for (const line of lines) {
const fields = line.split(",");
const point = [Number(fields[iSpeed]), Number(fields[iLitres])];
(fields[iSplit] === "train" ? train : test).push(point);
}
return { train, test };
}
// Solve a x = b by Gaussian elimination with partial pivoting.
function solve(a, b) {
const n = b.length;
const m = a.map((row, i) => [...row, b[i]]); // augmented matrix, copied
for (let col = 0; col < n; col++) {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small.
let pivot = col;
for (let r = col + 1; r < n; r++) {
if (Math.abs(m[r][col]) > Math.abs(m[pivot][col])) pivot = r;
}
[m[col], m[pivot]] = [m[pivot], m[col]];
for (let r = col + 1; r < n; r++) {
const factor = m[r][col] / m[col][col];
for (let c = col; c <= n; c++) m[r][c] -= factor * m[col][c];
}
}
// Back substitution, from the last unknown up.
const x = new Array(n).fill(0);
for (let r = n - 1; r >= 0; r--) {
let total = m[r][n];
for (let c = r + 1; c < n; c++) total -= m[r][c] * x[c];
x[r] = total / m[r][r];
}
return x;
}
// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
// of Z is 1, z_i, z_i^2, ..., z_i^degree.
function fit(points, degree, scale) {
const size = degree + 1;
const a = Array.from({ length: size }, () => new Array(size).fill(0));
const b = new Array(size).fill(0);
for (const [x, y] of points) {
const z = scale(x);
for (let i = 0; i < size; i++) {
b[i] += y * z ** i;
for (let j = 0; j < size; j++) a[i][j] += z ** (i + j);
}
}
return solve(a, b);
}
function rmse(coefs, points, scale) {
let total = 0;
for (const [x, y] of points) {
const z = scale(x);
const predicted = coefs.reduce((acc, c, i) => acc + c * z ** i, 0);
total += (predicted - y) ** 2;
}
return Math.sqrt(total / points.length);
}
const { train, test } = read(process.argv[2]);
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
// largest entry at degree 8) would swamp the small ones and make the
// equations hard to solve accurately.
const speeds = train.map(([x]) => x);
const low = Math.min(...speeds);
const high = Math.max(...speeds);
const mid = (low + high) / 2, half = (high - low) / 2;
const scale = (x) => (x - mid) / half;
console.log(`train_count ${train.length}`);
console.log(`test_count ${test.length}`);
const errors = [];
for (let degree = 1; degree <= MAX_DEGREE; degree++) {
const coefs = fit(train, degree, scale);
errors.push([degree, rmse(coefs, train, scale), rmse(coefs, test, scale)]);
}
for (const [degree, trainError] of errors) {
console.log(`train_rmse_d${degree} ${trainError.toFixed(6)}`);
}
for (const [degree, , testError] of errors) {
console.log(`test_rmse_d${degree} ${testError.toFixed(6)}`);
}Each drive becomes a [speed, litres] pair, pushed onto the train or test
array by its split column. For each degree we fill the normal equations
with nested loops, then solve runs Gaussian elimination on an augmented
copy, swapping up the row with the largest entry in each column. A
destructuring assignment, [m[a], m[b]] = [m[b], m[a]], swaps two rows in
one line.
Speeds are mapped onto −1 to 1 from the training range, which
Math.min(...speeds) and Math.max(...speeds) find by spreading the array
into arguments. JavaScript’s ** operator raises a negative z to a whole
power without trouble.
# Train and test error for polynomials of degree 1 to 8, using only Base Julia.
#
# Each polynomial is fitted by least squares to the training drives only:
# we build the normal equations and solve them by Gaussian elimination,
# written out instead of the backslash operator. Then we measure its
# root-mean-square error on the training drives and on the test drives it
# never saw.
using Printf
const MAX_DEGREE = 8
# Each drive becomes a (speed, litres) tuple; columns are found by name.
function read_drives(path)
lines = readlines(path)
columns = split(lines[1], ",")
i_speed = findfirst(==("speed"), columns)
i_litres = findfirst(==("litres"), columns)
i_split = findfirst(==("split"), columns)
train, test = Tuple{Float64,Float64}[], Tuple{Float64,Float64}[]
for line in lines[2:end]
fields = split(line, ",")
point = (parse(Float64, fields[i_speed]), parse(Float64, fields[i_litres]))
push!(fields[i_split] == "train" ? train : test, point)
end
return train, test
end
# Solve A x = b by Gaussian elimination with partial pivoting.
function gauss(A, b)
n = length(b)
M = [A b] # augmented matrix: a new copy, with b as an extra column
for col in 1:n
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small. argmax gives a position within col:n.
pivot = col - 1 + argmax(abs.(M[col:n, col]))
M[[col, pivot], :] = M[[pivot, col], :]
for r in col+1:n
factor = M[r, col] / M[col, col]
M[r, :] .-= factor .* M[col, :] # a whole row at once
end
end
# Back substitution, from the last unknown up.
x = zeros(n)
for r in n:-1:1
x[r] = (M[r, n+1] - sum(M[r, r+1:n] .* x[r+1:n])) / M[r, r]
end
return x
end
# Least squares: the coefficients c solve (Z' Z) c = Z' y, where row i of Z
# is 1, z_i, z_i^2, ..., z_i^degree. Broadcasting a column of z against a
# row of powers builds Z in one line.
powers(z, degree) = z .^ (0:degree)'
function fit(points, degree, scale)
Z = powers(scale.(first.(points)), degree)
y = last.(points)
return gauss(Z' * Z, Z' * y)
end
function rmse(coefs, points, scale)
predicted = powers(scale.(first.(points)), length(coefs) - 1) * coefs
return sqrt(sum((predicted .- last.(points)) .^ 2) / length(points))
end
function main(path)
train, test = read_drives(path)
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately.
low, high = extrema(first.(train))
mid, half = (low + high) / 2, (high - low) / 2
scale(x) = (x - mid) / half
println("train_count $(length(train))")
println("test_count $(length(test))")
errors = map(1:MAX_DEGREE) do degree
coefs = fit(train, degree, scale)
(rmse(coefs, train, scale), rmse(coefs, test, scale))
end
for (degree, (train_error, _)) in enumerate(errors)
@printf("train_rmse_d%d %.6f\n", degree, train_error)
end
for (degree, (_, test_error)) in enumerate(errors)
@printf("test_rmse_d%d %.6f\n", degree, test_error)
end
end
main(ARGS[1])z .^ (0:degree)' broadcasts a column of scaled speeds against a row of
powers, building the whole matrix at once, and Z' * Z gives the normal
equations. Our gauss function runs the elimination: argmax finds the
pivot, M[[col, pivot], :] swaps two rows, and .-= updates a row in
place.
Base Julia’s backslash, A \ b, would solve the system in one step; we
avoided it. Scaling matters here: on raw speeds, backslash on the degree-8
normal equations gives a training error of 0.073996, not 0.073821.
-- Train and test error for polynomials of degree 1 to 8, using only Lua's
-- standard library.
--
-- Each polynomial is fitted by least squares to the training drives only:
-- we build the normal equations and solve them by Gaussian elimination. Then
-- we measure its root-mean-square error on the training drives and on the
-- test drives it never saw. A matrix is a table of row tables, indexed
-- from 1.
local MAX_DEGREE = 8
-- Returns two arrays of { x = speed, y = litres }: the training drives and
-- the test drives.
local function read_drives(path)
local f = assert(io.open(path, "r"))
local col = {} -- column number for each name
local i = 0
for field in f:read("l"):gmatch("[^,]+") do
i = i + 1
col[field] = i
end
local train, test = {}, {}
for line in f:lines() do
local fields = {}
for field in line:gmatch("[^,]+") do
fields[#fields + 1] = field
end
local point = { x = tonumber(fields[col.speed]), y = tonumber(fields[col.litres]) }
local into = fields[col.split] == "train" and train or test
into[#into + 1] = point
end
f:close()
return train, test
end
-- Solve a x = b by Gaussian elimination with partial pivoting.
local function solve(a, b)
local n = #b
local m = {} -- the augmented matrix: b is column n + 1
for i = 1, n do
m[i] = table.move(a[i], 1, n, 1, {}) -- a copy of row i
m[i][n + 1] = b[i]
end
for col = 1, n do
-- Swap up the row with the largest entry in this column, which keeps the
-- rounding error small. Swapping two rows swaps two table references.
local pivot = col
for r = col + 1, n do
if math.abs(m[r][col]) > math.abs(m[pivot][col]) then
pivot = r
end
end
m[col], m[pivot] = m[pivot], m[col]
for r = col + 1, n do
local factor = m[r][col] / m[col][col]
for c = col, n + 1 do
m[r][c] = m[r][c] - factor * m[col][c]
end
end
end
-- Back substitution, from the last unknown up.
local x = {}
for r = n, 1, -1 do
local total = m[r][n + 1]
for c = r + 1, n do
total = total - m[r][c] * x[c]
end
x[r] = total / m[r][r]
end
return x
end
-- Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
-- Z is 1, z_i, z_i^2, ..., z_i^degree. Lua counts from 1, so entry i of a
-- row holds power i - 1.
local function fit(points, degree, scale)
local size = degree + 1
local a, b = {}, {}
for i = 1, size do
a[i] = {}
for j = 1, size do
a[i][j] = 0.0
end
b[i] = 0.0
end
for _, p in ipairs(points) do
local z = scale(p.x)
for i = 1, size do
b[i] = b[i] + p.y * z ^ (i - 1)
for j = 1, size do
a[i][j] = a[i][j] + z ^ (i + j - 2)
end
end
end
return solve(a, b)
end
local function rmse(coefs, points, scale)
local total = 0.0
for _, p in ipairs(points) do
local z = scale(p.x)
local predicted = 0.0
for i, c in ipairs(coefs) do
predicted = predicted + c * z ^ (i - 1)
end
total = total + (predicted - p.y) ^ 2
end
return math.sqrt(total / #points)
end
local train, test = read_drives(arg[1])
-- Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
-- largest entry at degree 8) would swamp the small ones and make the
-- equations hard to solve accurately.
local low, high = math.huge, -math.huge
for _, p in ipairs(train) do
low = math.min(low, p.x)
high = math.max(high, p.x)
end
local mid, half = (low + high) / 2, (high - low) / 2
local function scale(x)
return (x - mid) / half
end
print("train_count " .. #train)
print("test_count " .. #test)
local train_errors, test_errors = {}, {}
for degree = 1, MAX_DEGREE do
local coefs = fit(train, degree, scale)
train_errors[degree] = rmse(coefs, train, scale)
test_errors[degree] = rmse(coefs, test, scale)
end
for degree, e in ipairs(train_errors) do
print(string.format("train_rmse_d%d %.6f", degree, e))
end
for degree, e in ipairs(test_errors) do
print(string.format("test_rmse_d%d %.6f", degree, e))
endThe augmented matrix is a table of row tables, indexed from 1, so entry i
of a row holds power i - 1. Swapping two rows is one multiple assignment:
m[col], m[pivot] on the left takes m[pivot], m[col], which trades the two row references
without copying a number. table.move copies each row of the equations
into the matrix first, so the caller’s stays as it was.
Elimination and back substitution are plain numeric for loops, and ^
raises the scaled speeds to each power.
(* Train and test error for polynomials of degree 1 to 8, using only
OCaml's standard library.
Each polynomial is fitted by least squares to the training drives only:
we build the normal equations and solve them by Gaussian elimination.
Then we measure its root-mean-square error on the training drives and on
the test drives it never saw. A matrix is an array of row arrays, so
m.(i).(j) is row i, column j, counting from 0. *)
let max_degree = 8
(* Every drive as (split, speed, litres), with each column found by its
header name. *)
let read_drives path =
In_channel.with_open_text path (fun ic ->
let lines = In_channel.input_lines ic in
let header = String.split_on_char ',' (List.hd lines) in
let index name =
match List.find_index (String.equal name) header with
| Some i -> i
| None -> failwith ("column not found: " ^ name)
in
List.filter_map
(fun line ->
if line = "" then None
else
let fields = String.split_on_char ',' line in
let field name = List.nth fields (index name) in
Some
( field "split",
float_of_string (field "speed"),
float_of_string (field "litres") ))
(List.tl lines))
(* Solve a x = b by Gaussian elimination with partial pivoting. *)
let solve a b =
let n = Array.length b in
(* The augmented matrix: each row of a with its entry of b on the end. *)
let m = Array.init n (fun i -> Array.append a.(i) [| b.(i) |]) in
for col = 0 to n - 1 do
(* Swap up the row with the largest entry in this column, which keeps
the rounding error small. Swapping two rows just swaps the two row
arrays in m. *)
let pivot = ref col in
for r = col + 1 to n - 1 do
if Float.abs m.(r).(col) > Float.abs m.(!pivot).(col) then pivot := r
done;
let top = m.(col) in
m.(col) <- m.(!pivot);
m.(!pivot) <- top;
for r = col + 1 to n - 1 do
let factor = m.(r).(col) /. m.(col).(col) in
for c = col to n do
m.(r).(c) <- m.(r).(c) -. (factor *. m.(col).(c))
done
done
done;
(* Back substitution, from the last unknown up. *)
let x = Array.make n 0. in
for r = n - 1 downto 0 do
let known = ref 0. in
for c = r + 1 to n - 1 do
known := !known +. (m.(r).(c) *. x.(c))
done;
x.(r) <- (m.(r).(n) -. !known) /. m.(r).(r)
done;
x
(* z to a whole-number power. ** takes two floats, so the power converts. *)
let power z i = z ** float_of_int i
(* Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
Z is 1, z_i, z_i^2, ..., z_i^degree. points holds (z, y) pairs. *)
let fit points degree =
let size = degree + 1 in
let a = Array.make_matrix size size 0. and b = Array.make size 0. in
List.iter
(fun (z, y) ->
for i = 0 to degree do
b.(i) <- b.(i) +. (y *. power z i);
for j = 0 to degree do
a.(i).(j) <- a.(i).(j) +. power z (i + j)
done
done)
points;
solve a b
let rmse coefs points =
let predict z = Array.fold_left ( +. ) 0. (Array.mapi (fun i c -> c *. power z i) coefs) in
let total = List.fold_left (fun acc (z, y) -> acc +. ((predict z -. y) ** 2.)) 0. points in
sqrt (total /. float_of_int (List.length points))
let () =
let drives = read_drives Sys.argv.(1) in
let part name =
List.filter_map (fun (s, x, y) -> if s = name then Some (x, y) else None) drives
in
let train = part "train" and test = part "test" in
(* Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
largest entry at degree 8) would swamp the small ones and make the
equations hard to solve accurately. *)
let speeds = List.map fst train in
let low = List.fold_left Float.min infinity speeds in
let high = List.fold_left Float.max neg_infinity speeds in
let mid = (low +. high) /. 2. and half = (high -. low) /. 2. in
let scaled points = List.map (fun (x, y) -> ((x -. mid) /. half, y)) points in
let train_z = scaled train and test_z = scaled test in
Printf.printf "train_count %d\n" (List.length train);
Printf.printf "test_count %d\n" (List.length test);
let fits = List.init max_degree (fun i -> (i + 1, fit train_z (i + 1))) in
List.iter (fun (d, c) -> Printf.printf "train_rmse_d%d %.6f\n" d (rmse c train_z)) fits;
List.iter (fun (d, c) -> Printf.printf "test_rmse_d%d %.6f\n" d (rmse c test_z)) fitsThe augmented matrix is an array of row arrays, so partial pivoting swaps
two rows by swapping the two array references in m, without copying any
numbers. The elimination and back substitution are for loops that update
the entries in place with <-.
OCaml’s ** takes two floats; writing 2.0 ** 3 is a compile error. A
small power helper converts the exponent with float_of_int. Speeds are
mapped onto −1 to 1 using the training range only, and every fit and error
is computed from the drives’ (z, litres) pairs. No solver or fitting
function is called.
% Train and test error for polynomials of degree 1 to 8, using only core GNU
% Octave.
%
% Each polynomial is fitted by least squares to the training drives only:
% we build the normal equations and solve them by Gaussian elimination. Then
% we measure its root-mean-square error on the training drives and on the
% test drives it never saw. Octave could fit each one with polyfit, or solve
% the equations with backslash, A \ b; the elimination is written out instead.
MAX_DEGREE = 8;
% Solve A x = b by Gaussian elimination with partial pivoting.
function x = solve(A, b)
n = numel(b);
M = [A, b]; % the augmented matrix: b becomes column n + 1
for col = 1:n
% Swap up the row with the largest entry in this column, which keeps the
% rounding error small. max returns the value and its position, counted
% from the top of the slice, so add col - 1 to get the row.
[~, offset] = max(abs(M(col:n, col)));
pivot = col + offset - 1;
M([col, pivot], :) = M([pivot, col], :);
for r = col + 1:n
factor = M(r, col) / M(col, col);
M(r, col:end) -= factor * M(col, col:end);
end
end
% Back substitution, from the last unknown up. M(r, r + 1:n) is a row and
% x(r + 1:n) a column, so * is their dot product.
x = zeros(n, 1);
for r = n:-1:1
x(r) = (M(r, n + 1) - M(r, r + 1:n) * x(r + 1:n)) / M(r, r);
end
end
% Row i of the result is 1, z_i, z_i^2, ..., z_i^degree. z is a column and
% 0:degree a row, so .^ broadcasts them into a matrix.
function P = powers(z, degree)
P = z .^ (0:degree);
end
function e = rmse(coefs, z, y)
predicted = powers(z, numel(coefs) - 1) * coefs;
e = sqrt(sum((predicted - y) .^ 2) / numel(y));
end
args = argv();
lines = strsplit(strtrim(fileread(args{1})), "\n");
header = strsplit(strtrim(lines{1}), ",");
speed_col = find(strcmp(header, "speed"));
litres_col = find(strcmp(header, "litres"));
split_col = find(strcmp(header, "split"));
n = numel(lines) - 1;
speed = zeros(n, 1);
litres = zeros(n, 1);
is_train = false(n, 1);
for row = 1:n
fields = strsplit(strtrim(lines{row + 1}), ",");
speed(row) = str2double(fields{speed_col});
litres(row) = str2double(fields{litres_col});
is_train(row) = strcmp(fields{split_col}, "train");
end
% A logical vector picks out the rows where it is true.
x_train = speed(is_train);
y_train = litres(is_train);
x_test = speed(!is_train);
y_test = litres(!is_train);
% Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
% largest entry at degree 8) would swamp the small ones and make the
% equations hard to solve accurately.
low = min(x_train);
high = max(x_train);
mid = (low + high) / 2;
half = (high - low) / 2;
z_train = (x_train - mid) / half;
z_test = (x_test - mid) / half;
printf("train_count %d\n", numel(y_train));
printf("test_count %d\n", numel(y_test));
train_errors = zeros(MAX_DEGREE, 1);
test_errors = zeros(MAX_DEGREE, 1);
for degree = 1:MAX_DEGREE
% Least squares: the coefficients c solve (P' P) c = P' y. Entry (i, j) of
% P' P is the sum of z^(i + j), and entry i of P' y the sum of y z^i.
P = powers(z_train, degree);
coefs = solve(P' * P, P' * y_train);
train_errors(degree) = rmse(coefs, z_train, y_train);
test_errors(degree) = rmse(coefs, z_test, y_test);
end
for degree = 1:MAX_DEGREE
printf("train_rmse_d%d %.6f\n", degree, train_errors(degree));
end
for degree = 1:MAX_DEGREE
printf("test_rmse_d%d %.6f\n", degree, test_errors(degree));
endOctave can fit each curve in one call: polyfit(x, y, d) returns the
coefficients, and the backslash solve A \ b would answer the normal
equations directly. We skip both and write the elimination out in solve.
Octave still keeps the rest short. z .^ (0:degree) broadcasts a column of
scaled speeds against a row of powers, so it builds every row of
1, z, z², and so on at once, and P' * P and P' * y are the normal
equations. Swapping two rows is one indexed assignment:
M([col, pivot], :) on the left takes M([pivot, col], :). To find the
pivot, max returns the position of the largest entry along with its value.
{ Train and test error for polynomials of degree 1 to 8, using only Free
Pascal's run-time library.
Each polynomial is fitted by least squares to the training drives only:
we build the normal equations and solve them by Gaussian elimination.
Then we measure its root-mean-square error on the training drives and on
the test drives it never saw. }
program TrainAndTest;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
const
MaxDegree = 8;
MaxSize = MaxDegree + 1; // unknowns in the largest system
type
TPoint = record
X, Y: Double;
end;
TPoints = array of TPoint; // a dynamic array, counted from 0
// Fixed arrays big enough for the largest system; a smaller degree uses
// the top-left corner. Column MaxSize + 1 holds the right-hand side.
TMatrix = array[0..MaxSize - 1, 0..MaxSize] of Double;
TVector = array[0..MaxSize - 1] of Double;
var
Mid, Half: Double; // the scale, set from the training speeds
procedure Append(var Points: TPoints; X, Y: Double);
begin
SetLength(Points, Length(Points) + 1); // grow by one
Points[High(Points)].X := X;
Points[High(Points)].Y := Y;
end;
{ Read every drive into the training or the test points. }
procedure ReadDrives(const Path: String; out Train, Test: TPoints);
var
F: TextFile;
Line: String;
Fields: TStringArray;
I, SpeedCol, LitresCol, SplitCol: Integer;
X, Y: Double;
begin
Train := nil;
Test := nil;
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
SpeedCol := -1;
LitresCol := -1;
SplitCol := -1;
for I := 0 to High(Fields) do
if Fields[I] = 'speed' then SpeedCol := I
else if Fields[I] = 'litres' then LitresCol := I
else if Fields[I] = 'split' then SplitCol := I;
if (SpeedCol < 0) or (LitresCol < 0) or (SplitCol < 0) then
raise Exception.Create('missing column');
while not Eof(F) do
begin
ReadLn(F, Line);
Fields := Line.Split(',');
X := StrToFloat(Fields[SpeedCol]);
Y := StrToFloat(Fields[LitresCol]);
if Fields[SplitCol] = 'train' then
Append(Train, X, Y)
else
Append(Test, X, Y);
end;
CloseFile(F);
end;
function Scale(X: Double): Double;
begin
Result := (X - Mid) / Half;
end;
{ Solve the N equations in the augmented matrix M by Gaussian elimination
with partial pivoting. M is a copy (no var), so the caller's is kept. }
function Solve(M: TMatrix; N: Integer): TVector;
var
Col, R, C, Pivot: Integer;
Factor, Temp, Total: Double;
X: TVector;
begin
// With warnings as errors, X must be set before the compiler sees it read.
FillChar(X, SizeOf(X), 0);
for Col := 0 to N - 1 do
begin
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small.
Pivot := Col;
for R := Col + 1 to N - 1 do
if Abs(M[R, Col]) > Abs(M[Pivot, Col]) then
Pivot := R;
for C := 0 to N do
begin
Temp := M[Col, C];
M[Col, C] := M[Pivot, C];
M[Pivot, C] := Temp;
end;
for R := Col + 1 to N - 1 do
begin
Factor := M[R, Col] / M[Col, Col];
for C := Col to N do
M[R, C] := M[R, C] - Factor * M[Col, C];
end;
end;
// Back substitution, from the last unknown up.
for R := N - 1 downto 0 do
begin
Total := M[R, N];
for C := R + 1 to N - 1 do
Total := Total - M[R, C] * X[C];
X[R] := Total / M[R, R];
end;
Result := X;
end;
{ Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
Z is 1, z_i, z_i^2, ..., z_i^Degree. }
function Fit(const Points: TPoints; Degree: Integer): TVector;
var
M: TMatrix;
Powers: array[0..2 * MaxDegree] of Double;
P, I, J: Integer;
begin
FillChar(M, SizeOf(M), 0); // set every byte to 0, which is 0.0
for P := 0 to High(Points) do
begin
Powers[0] := 1;
for I := 1 to 2 * Degree do
Powers[I] := Powers[I - 1] * Scale(Points[P].X);
for I := 0 to Degree do
begin
M[I, Degree + 1] := M[I, Degree + 1] + Points[P].Y * Powers[I];
for J := 0 to Degree do
M[I, J] := M[I, J] + Powers[I + J];
end;
end;
Result := Solve(M, Degree + 1);
end;
function Rmse(const Coefs: TVector; Degree: Integer; const Points: TPoints): Double;
var
P, I: Integer;
Predicted, Total: Double;
begin
Total := 0;
for P := 0 to High(Points) do
begin
// Horner's rule: c0 + z (c1 + z (c2 + ...)).
Predicted := 0;
for I := Degree downto 0 do
Predicted := Predicted * Scale(Points[P].X) + Coefs[I];
Total := Total + Sqr(Predicted - Points[P].Y);
end;
Result := Sqrt(Total / Length(Points));
end;
var
Train, Test: TPoints;
Coefs: TVector;
TrainErrors, TestErrors: array[1..MaxDegree] of Double;
LowX, HighX: Double;
Degree, I: Integer;
begin
ReadDrives(ParamStr(1), Train, Test);
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power
// (the largest entry at degree 8) would swamp the small ones and make
// the equations hard to solve accurately.
LowX := Train[0].X;
HighX := LowX;
for I := 1 to High(Train) do
begin
if Train[I].X < LowX then LowX := Train[I].X;
if Train[I].X > HighX then HighX := Train[I].X;
end;
Mid := (LowX + HighX) / 2;
Half := (HighX - LowX) / 2;
WriteLn('train_count ', Length(Train));
WriteLn('test_count ', Length(Test));
for Degree := 1 to MaxDegree do
begin
Coefs := Fit(Train, Degree);
TrainErrors[Degree] := Rmse(Coefs, Degree, Train);
TestErrors[Degree] := Rmse(Coefs, Degree, Test);
end;
for Degree := 1 to MaxDegree do
WriteLn('train_rmse_d', Degree, ' ', TrainErrors[Degree]:0:6); // six decimals
for Degree := 1 to MaxDegree do
WriteLn('test_rmse_d', Degree, ' ', TestErrors[Degree]:0:6);
end.The equations live in a fixed TMatrix, sized for degree 8, with the
right-hand side in its last column; lower degrees use its top-left corner.
Solve takes the matrix by value, so it works on its own copy, and swaps
two rows one entry at a time through Temp.
No installed unit has a linear solver to skip. Free Pascal’s NumLib has
one, but it isn’t installed in this image. With warnings treated as errors,
the compiler rejected a function result filled in a loop as possibly
uninitialized, so Solve clears its answer with FillChar first.
# Train and test error for polynomials of degree 1 to 8, using only Perl and
# its core modules.
#
# Each polynomial is fitted by least squares to the training drives only: we
# build the normal equations and solve them by Gaussian elimination. Then we
# measure its root-mean-square error on the training drives and on the test
# drives it never saw. A matrix is an array of array refs: $m->[$i][$j] is
# row $i, column $j.
use strict;
use warnings;
use List::Util qw(sum min max);
my $MAX_DEGREE = 8;
# Returns two array refs of [speed, litres] pairs: the training drives and
# the test drives.
sub read_drives {
my ($path) = @_;
open my $fh, '<', $path or die "can't open $path: $!";
chomp(my $header = <$fh>);
my @names = split /,/, $header;
my (@train, @test);
while (my $line = <$fh>) {
chomp $line;
my %row;
@row{@names} = split /,/, $line; # a hash slice: names to fields
my $point = [ $row{speed}, $row{litres} ];
if ($row{split} eq 'train') { push @train, $point } else { push @test, $point }
}
close $fh;
return (\@train, \@test);
}
# Solve a x = b by Gaussian elimination with partial pivoting. ($a and $b
# are reserved for sort, so the arguments are $lhs and $rhs.)
sub solve {
my ($lhs, $rhs) = @_;
my $n = @$rhs; # an array in scalar context is its length
# The augmented matrix: a copy of each row, with b's entry on the end.
my @m = map { [ @{ $lhs->[$_] }, $rhs->[$_] ] } 0 .. $n - 1;
for my $col (0 .. $n - 1) {
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small. The rows are refs, so a slice assignment
# swaps them without copying any numbers.
my $pivot = $col;
for my $r ($col + 1 .. $n - 1) {
$pivot = $r if abs $m[$r][$col] > abs $m[$pivot][$col];
}
@m[$col, $pivot] = @m[$pivot, $col];
for my $r ($col + 1 .. $n - 1) {
my $factor = $m[$r][$col] / $m[$col][$col];
$m[$r][$_] -= $factor * $m[$col][$_] for $col .. $n;
}
}
# Back substitution, from the last unknown up.
my @x = (0) x $n;
for my $r (reverse 0 .. $n - 1) {
my $total = $m[$r][$n];
$total -= $m[$r][$_] * $x[$_] for $r + 1 .. $n - 1;
$x[$r] = $total / $m[$r][$r];
}
return \@x;
}
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
# Z is 1, z_i, z_i^2, ..., z_i^degree.
sub fit {
my ($points, $degree, $scale) = @_;
my @a = map { [ (0) x ($degree + 1) ] } 0 .. $degree;
my @b = (0) x ($degree + 1);
for my $p (@$points) {
my ($x, $y) = @$p;
my $z = $scale->($x);
for my $i (0 .. $degree) {
$b[$i] += $y * $z ** $i;
$a[$i][$_] += $z ** ($i + $_) for 0 .. $degree;
}
}
return solve(\@a, \@b);
}
sub rmse {
my ($coefs, $points, $scale) = @_;
my $total = 0;
for my $p (@$points) {
my ($x, $y) = @$p;
my $z = $scale->($x);
my $predicted = sum(map { $coefs->[$_] * $z ** $_ } 0 .. $#$coefs);
$total += ($predicted - $y) ** 2;
}
return sqrt($total / @$points);
}
my ($train, $test) = read_drives($ARGV[0]);
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately.
my @speeds = map { $_->[0] } @$train;
my ($low, $high) = (min(@speeds), max(@speeds));
my ($mid, $half) = (($low + $high) / 2, ($high - $low) / 2);
my $scale = sub { ($_[0] - $mid) / $half }; # an anonymous sub, called as $scale->(x)
printf "train_count %d\n", scalar @$train;
printf "test_count %d\n", scalar @$test;
my (@train_errors, @test_errors);
for my $degree (1 .. $MAX_DEGREE) {
my $coefs = fit($train, $degree, $scale);
push @train_errors, rmse($coefs, $train, $scale);
push @test_errors, rmse($coefs, $test, $scale);
}
printf "train_rmse_d%d %.6f\n", $_ + 1, $train_errors[$_] for 0 .. $#train_errors;
printf "test_rmse_d%d %.6f\n", $_ + 1, $test_errors[$_] for 0 .. $#test_errors;Each row of the augmented matrix is an array ref, so the pivot swap is one
slice assignment: @m[$col, $pivot] takes @m[$pivot, $col], which trades two
refs without copying a number. The rest are loops written with statement
modifiers, such as ... for $col .. $n.
solve names its arguments $lhs and $rhs, because $a and $b are
the variables sort uses.
PDL and Math::Matrix would solve the system, but neither is installed.
<?php
// Train and test error for polynomials of degree 1 to 8, using only PHP and
// the extensions built into it.
//
// Each polynomial is fitted by least squares to the training drives only: we
// build the normal equations and solve them by Gaussian elimination. Then we
// measure its root-mean-square error on the training drives and on the test
// drives it never saw. A matrix is an array of row arrays.
declare(strict_types=1);
const MAX_DEGREE = 8;
// No field is quoted, so explode splits each line on commas. Returns the
// training drives and the test drives, each a list of [speed, litres].
function read_drives(string $path): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
$header = explode(',', array_shift($lines));
$train = [];
$test = [];
foreach ($lines as $line) {
// array_combine pairs each header name with its field.
$row = array_combine($header, explode(',', $line));
$point = [(float) $row['speed'], (float) $row['litres']];
if ($row['split'] === 'train') {
$train[] = $point;
} else {
$test[] = $point;
}
}
return [$train, $test];
}
// Solve a x = b by Gaussian elimination with partial pivoting. PHP passes
// arrays by value, so changing $m here leaves the caller's copy alone.
function solve(array $a, array $b): array
{
$n = count($b);
$m = []; // the augmented matrix: b is column n
foreach ($a as $i => $row) {
$m[] = [...$row, $b[$i]];
}
for ($col = 0; $col < $n; $col++) {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small.
$pivot = $col;
for ($r = $col + 1; $r < $n; $r++) {
if (abs($m[$r][$col]) > abs($m[$pivot][$col])) {
$pivot = $r;
}
}
[$m[$col], $m[$pivot]] = [$m[$pivot], $m[$col]];
for ($r = $col + 1; $r < $n; $r++) {
$factor = $m[$r][$col] / $m[$col][$col];
for ($c = $col; $c <= $n; $c++) {
$m[$r][$c] -= $factor * $m[$col][$c];
}
}
}
// Back substitution, from the last unknown up.
$x = array_fill(0, $n, 0.0);
for ($r = $n - 1; $r >= 0; $r--) {
$total = $m[$r][$n];
for ($c = $r + 1; $c < $n; $c++) {
$total -= $m[$r][$c] * $x[$c];
}
$x[$r] = $total / $m[$r][$r];
}
return $x;
}
// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
// Z is 1, z_i, z_i^2, ..., z_i^degree.
function fit(array $points, int $degree, callable $scale): array
{
$size = $degree + 1;
$a = array_fill(0, $size, array_fill(0, $size, 0.0));
$b = array_fill(0, $size, 0.0);
foreach ($points as [$x, $y]) {
$z = $scale($x);
for ($i = 0; $i < $size; $i++) {
$b[$i] += $y * $z ** $i;
for ($j = 0; $j < $size; $j++) {
$a[$i][$j] += $z ** ($i + $j);
}
}
}
return solve($a, $b);
}
function rmse(array $coefs, array $points, callable $scale): float
{
$total = 0.0;
foreach ($points as [$x, $y]) {
$z = $scale($x);
$predicted = 0.0;
foreach ($coefs as $i => $c) {
$predicted += $c * $z ** $i;
}
$total += ($predicted - $y) ** 2;
}
return sqrt($total / count($points));
}
[$train, $test] = read_drives($argv[1]);
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
// largest entry at degree 8) would swamp the small ones and make the
// equations hard to solve accurately.
$speeds = array_column($train, 0);
$low = min($speeds);
$high = max($speeds);
$mid = ($low + $high) / 2;
$half = ($high - $low) / 2;
$scale = fn(float $x): float => ($x - $mid) / $half; // fn captures $mid and $half
printf("train_count %d\n", count($train));
printf("test_count %d\n", count($test));
$train_errors = [];
$test_errors = [];
for ($degree = 1; $degree <= MAX_DEGREE; $degree++) {
$coefs = fit($train, $degree, $scale);
$train_errors[$degree] = rmse($coefs, $train, $scale);
$test_errors[$degree] = rmse($coefs, $test, $scale);
}
foreach ($train_errors as $degree => $e) {
printf("train_rmse_d%d %.6f\n", $degree, $e);
}
foreach ($test_errors as $degree => $e) {
printf("test_rmse_d%d %.6f\n", $degree, $e);
}The augmented matrix is an array of row arrays, built with the spread
operator: [...$row, $b[$i]]. Swapping two rows is one destructuring
assignment, [$m[$col], $m[$pivot]] = ... with the pair reversed on the
right. PHP passes arrays by value, so solve works on its own copy.
There’s no matrix type or solver to skip, so elimination and back
substitution are for loops. $scale is an arrow function, fn, which
picks up $mid and $half from the code around it without a use list.
-- Train and test error for polynomials of degree 1 to 8, in PostgreSQL.
--
-- Each polynomial is fitted by least squares to the training drives only: we
-- build the normal equations and solve them by Gaussian elimination. Then we
-- measure its root-mean-square error on the training drives and on the test
-- drives it never saw. Each system is a two-dimensional double precision[][]
-- array, and elimination is a WITH RECURSIVE query whose rows are one step
-- each: (degree, column, matrix).
WITH RECURSIVE range AS (
-- Map the training speeds onto -1 to 1. Raw speeds to the 16th power
-- (the largest entry at degree 8) would swamp the small ones and make
-- the equations hard to solve accurately.
SELECT (min(speed) + max(speed)) / 2 AS mid, (max(speed) - min(speed)) / 2 AS half
FROM data
WHERE split = 'train'
),
drives AS (
SELECT split, (speed - mid) / half AS z, litres AS y
FROM data, range
),
degrees AS (
SELECT d FROM generate_series(1, 8) AS d
),
entries AS (
-- The augmented normal equations, one row per entry (array positions
-- start at 1). Columns 1 to d + 1 hold sum z^(i+j); the last column,
-- d + 2, holds b, sum y z^i.
SELECT d, i, j,
CASE WHEN j <= d + 1 THEN sum(z ^ (i + j - 2))
ELSE sum(y * z ^ (i - 1)) END AS v
FROM degrees,
generate_series(1, d + 1) AS i,
generate_series(1, d + 2) AS j,
drives
WHERE split = 'train'
GROUP BY d, i, j
),
systems AS (
-- array_agg of one-dimensional arrays stacks them into a matrix.
SELECT d, array_agg(row ORDER BY i) AS m
FROM (SELECT d, i, array_agg(v ORDER BY j) AS row FROM entries GROUP BY d, i) AS rows
GROUP BY d
),
elimination (d, k, m) AS (
SELECT d, 0, m FROM systems
UNION ALL
-- k counts the columns already cleared, and this step clears column
-- k + 1 below the diagonal. First the row with the largest entry in that
-- column (at or below row k + 1) swaps up to row k + 1, which keeps the
-- rounding error small; then each lower row subtracts a multiple of it. Aggregates aren't allowed directly in the
-- recursive step, so they sit in LATERAL.
SELECT e.d, e.k + 1, next.m
FROM elimination e,
LATERAL (
SELECT r AS pivot
FROM generate_series(e.k + 1, e.d + 1) AS r
ORDER BY abs(e.m[r][e.k + 1]) DESC, r
LIMIT 1
) AS p,
LATERAL (
-- s is the row of the old matrix that lands at row r after the swap,
-- so swapped row r, column c is e.m[s][c], and the new row k is
-- e.m[p.pivot].
SELECT array_agg(row ORDER BY r) AS m
FROM (
SELECT r, array_agg(
CASE WHEN r <= e.k + 1 OR c < e.k + 1 THEN e.m[s][c]
ELSE e.m[s][c] - e.m[s][e.k + 1] / e.m[p.pivot][e.k + 1]
* e.m[p.pivot][c] END
ORDER BY c) AS row
FROM generate_series(1, e.d + 1) AS r,
LATERAL (SELECT CASE r WHEN e.k + 1 THEN p.pivot
WHEN p.pivot THEN e.k + 1
ELSE r END AS s) AS swap,
generate_series(1, e.d + 2) AS c
GROUP BY r
) AS rows
) AS next
WHERE e.k < e.d -- the last column has no rows below it
),
triangular AS (
SELECT d, m FROM elimination WHERE k = d
),
back (d, r, x) AS (
-- Back substitution, from the last unknown up. x starts as zeros, and
-- each step fills in x[r].
SELECT d, d + 2, array_fill(0::double precision, ARRAY[d + 1])
FROM triangular
UNION ALL
SELECT b.d, b.r - 1, next.x
FROM back b JOIN triangular t USING (d),
LATERAL (
SELECT b.x[:b.r - 2]
|| (t.m[b.r - 1][b.d + 2] - coalesce(sum(t.m[b.r - 1][c] * b.x[c]), 0))
/ t.m[b.r - 1][b.r - 1]
|| b.x[b.r:] AS x
FROM generate_series(b.r, b.d + 1) AS c
) AS next
WHERE b.r > 1
),
coefficients AS (
SELECT d, x FROM back WHERE r = 1
),
errors AS (
-- Root-mean-square error of each fit, on each split.
SELECT c.d, s.split, sqrt(sum((s.predicted - s.y) ^ 2) / count(*)) AS rmse
FROM coefficients c,
LATERAL (
SELECT dr.split, dr.y, (SELECT sum(c.x[i] * dr.z ^ (i - 1))
FROM generate_series(1, c.d + 1) AS i) AS predicted
FROM drives dr
) AS s
GROUP BY c.d, s.split
),
results (position, name, value) AS (
SELECT 1, 'train_count', count(*)::text FROM data WHERE split = 'train'
UNION ALL
SELECT 2, 'test_count', count(*)::text FROM data WHERE split = 'test'
UNION ALL
SELECT CASE split WHEN 'train' THEN 10 ELSE 20 END + d,
split || '_rmse_d' || d, round(rmse::numeric, 6)::text
FROM errors
)
SELECT name, value FROM results ORDER BY position;We store each degree’s augmented normal equations as a two-dimensional
double precision[][], stacked from one-row arrays by array_agg.
Elimination is WITH RECURSIVE, one row per cleared column. A query
can’t assign to one array entry, so each step builds a whole new
matrix: a CASE maps every row to its source row, which performs the
pivot swap, and the lower rows subtract a multiple of the pivot row. A
second recursive query does back substitution. Note that m[2] on a 2D
array is NULL, not row 2. We skip regr_slope, which fits only a line;
core PostgreSQL has no polynomial fit or linear solver.
% Train and test error for polynomials of degree 1 to 8, using SWI-Prolog
% and the libraries that come with it.
%
% Each polynomial is fitted by least squares to the training drives only:
% we build the normal equations and solve them by Gaussian elimination.
% Then we measure its root-mean-square error on the training drives and on
% the test drives it never saw. A matrix is a list of rows.
:- initialization(main, main).
max_degree(8).
% read_drives(+Path, -Drives): every drive as Split-Speed-Litres, in file
% order, with each column found by its header name.
read_drives(Path, Drives) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
nth0(S, Header, "speed"),
nth0(L, Header, "litres"),
nth0(P, Header, "split"),
findall(Split-Speed-Litres,
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
nth0(S, Fields, SpeedText),
nth0(L, Fields, LitresText),
nth0(P, Fields, Split),
number_string(Speed, SpeedText),
number_string(Litres, LitresText) ),
Drives).
dot(U, V, Dot) :-
foldl([A, B, Acc0, Acc]>>(Acc is Acc0 + A * B), U, V, 0, Dot).
% solve(+A, +B, -X): solve A x = B by Gaussian elimination with partial
% pivoting. Each row of the augmented matrix is a row of A with its entry
% of B on the end.
solve(A, B, X) :-
maplist([Row, V, Aug]>>append(Row, [V], Aug), A, B, Augmented),
eliminate(Augmented, Triangle),
back_substitute(Triangle, X).
% eliminate(+Rows, -Triangle): take out the row with the largest first
% entry, which keeps the rounding error small, and subtract a multiple of
% it from every other row so their first entries become 0. Those rows,
% without their first column, are a smaller system of the same kind, so we
% recurse. Each row of the triangle is one entry shorter than the last.
eliminate([], []).
eliminate([First|Rest], [Pivot|Triangle]) :-
foldl(larger_leading, Rest, First, Pivot),
selectchk(Pivot, [First|Rest], Others),
Pivot = [P|PivotRest],
maplist(reduce(P, PivotRest), Others, Reduced),
eliminate(Reduced, Triangle).
% larger_leading(+Row, +Best0, -Best): keep whichever row has the larger
% first entry, ignoring sign.
larger_leading(Row, Best0, Best) :-
Row = [A|_],
Best0 = [B|_],
( abs(A) > abs(B)
-> Best = Row
; Best = Best0
).
% reduce(+P, +PivotRest, +Row, -Reduced): Row minus Factor times the pivot
% row, dropping the first column, which is now 0.
reduce(P, PivotRest, [R|RowRest], Reduced) :-
Factor is R / P,
maplist([PC, RC, NC]>>(NC is RC - Factor * PC), PivotRest, RowRest, Reduced).
% back_substitute(+Triangle, -X): from the last unknown up. Each row is
% [Coefficient, entries for the later unknowns ..., right-hand side].
back_substitute([], []).
back_substitute([[C|Rest]|Rows], [X|Later]) :-
back_substitute(Rows, Later),
append(Entries, [Rhs], Rest),
dot(Entries, Later, Known),
X is (Rhs - Known) / C.
% fit(+Points, +Degree, -Coefs): least squares. The coefficients solve
% (Z^T Z) c = Z^T y, where row i of Z is 1, z_i, z_i^2, ..., z_i^Degree.
% Points holds Z-Y pairs.
fit(Points, Degree, Coefs) :-
numlist(0, Degree, Powers),
findall(Row,
( member(I, Powers),
findall(Sum,
( member(J, Powers),
aggregate_all(sum(Z ** (I + J)), member(Z-_, Points), Sum) ),
Row) ),
A),
findall(Sum,
( member(I, Powers),
aggregate_all(sum(Y * Z ** I), member(Z-Y, Points), Sum) ),
B),
solve(A, B, Coefs).
% predict(+Coefs, +Z, -Y): c0 + c1 z + c2 z^2 + ...
predict(Coefs, Z, Y) :-
length(Coefs, N),
N1 is N - 1,
numlist(0, N1, Powers),
foldl([C, I, Acc0, Acc]>>(Acc is Acc0 + C * Z ** I), Coefs, Powers, 0, Y).
rmse(Coefs, Points, RMSE) :-
aggregate_all(sum((P - Y) ** 2),
( member(Z-Y, Points), predict(Coefs, Z, P) ),
Total),
length(Points, N),
RMSE is sqrt(Total / N).
main :-
current_prolog_flag(argv, [Path|_]),
read_drives(Path, Drives),
findall(X-Y, member("train"-X-Y, Drives), Train),
findall(X-Y, member("test"-X-Y, Drives), Test),
% Map the training speeds onto -1 to 1. Raw speeds to the 16th power
% (the largest entry at degree 8) would swamp the small ones and make
% the equations hard to solve accurately. The speeds are integers, so
% float/1 makes each z a float rather than an integer or a fraction.
findall(X, member(X-_, Train), Speeds),
min_list(Speeds, Low),
max_list(Speeds, High),
Mid is (Low + High) / 2,
Half is (High - Low) / 2,
Scale = [X-Y, Z-Y]>>(Z is float((X - Mid) / Half)),
maplist(Scale, Train, TrainZ),
maplist(Scale, Test, TestZ),
length(Train, NTrain),
length(Test, NTest),
format("train_count ~d~n", [NTrain]),
format("test_count ~d~n", [NTest]),
max_degree(Max),
numlist(1, Max, Degrees),
maplist(fit(TrainZ), Degrees, Fits),
forall(nth1(D, Fits, Coefs),
( rmse(Coefs, TrainZ, E), format("train_rmse_d~d ~6f~n", [D, E]) )),
forall(nth1(D, Fits, Coefs),
( rmse(Coefs, TestZ, E), format("test_rmse_d~d ~6f~n", [D, E]) )).A matrix is a list of rows, and the elimination is recursive. foldl
walks the rows to find the one with the largest first entry (partial
pivoting), selectchk takes it out, and every other row has a multiple
of it subtracted. The smaller system that’s left goes back into
eliminate.
The speeds are read as integers, and SWI-Prolog’s / returns an integer
when it divides exactly: (75-42)/33 is 1, but 7/2 is 3.5. So
float/1 makes every z a float. The sums use aggregate_all(sum(...)),
and no solver or fitting predicate is called.
# Train and test error for polynomials of degree 1 to 8, using only base R.
#
# Each polynomial is fitted by least squares to the training drives only:
# we build the normal equations and solve them by Gaussian elimination,
# written out instead of lm() or solve(). Then we measure its root-mean-square
# error on the training drives and on the test drives it never saw.
max_degree <- 8
args <- commandArgs(trailingOnly = TRUE)
data <- read.csv(args[1])
is_train <- data$split == "train"
train_x <- data$speed[is_train]
train_y <- data$litres[is_train]
test_x <- data$speed[!is_train]
test_y <- data$litres[!is_train]
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately.
low <- min(train_x)
high <- max(train_x)
mid <- (low + high) / 2
half <- (high - low) / 2
scale <- function(x) (x - mid) / half
# Solve a x = b by Gaussian elimination with partial pivoting.
gauss <- function(a, b) {
n <- length(b)
m <- cbind(a, b) # augmented matrix; R copies on change, so a is untouched
for (col in 1:n) {
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small. which.max() finds the first largest.
pivot <- col - 1 + which.max(abs(m[col:n, col]))
m[c(col, pivot), ] <- m[c(pivot, col), ]
if (col < n) {
for (r in (col + 1):n) {
factor <- m[r, col] / m[col, col]
m[r, ] <- m[r, ] - factor * m[col, ] # a whole row at once
}
}
}
# Back substitution, from the last unknown up.
x <- numeric(n)
for (r in n:1) {
later <- seq_len(n - r) + r # the columns after r (none for the last)
x[r] <- (m[r, n + 1] - sum(m[r, later] * x[later])) / m[r, r]
}
x
}
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
# of Z is 1, z_i, z_i^2, ..., z_i^degree. outer() builds Z, one column per
# power, and %*% is plain matrix multiplication.
fit <- function(degree) {
z <- outer(scale(train_x), 0:degree, "^")
gauss(t(z) %*% z, drop(t(z) %*% train_y))
}
rmse <- function(coefs, x, y) {
z <- outer(scale(x), 0:(length(coefs) - 1), "^")
predicted <- drop(z %*% coefs)
sqrt(sum((predicted - y)^2) / length(y))
}
cat(sprintf("train_count %d\n", length(train_x)))
cat(sprintf("test_count %d\n", length(test_x)))
degrees <- 1:max_degree
coefs <- lapply(degrees, fit)
train_rmse <- sapply(coefs, rmse, train_x, train_y)
test_rmse <- sapply(coefs, rmse, test_x, test_y)
cat(sprintf("train_rmse_d%d %.6f\n", degrees, train_rmse), sep = "")
cat(sprintf("test_rmse_d%d %.6f\n", degrees, test_rmse), sep = "")outer(z, 0:degree, "^") builds the matrix of powers in one call, and
t(z) %*% z gives the normal equations. Our gauss function then runs
Gaussian elimination: which.max finds the pivot row, and
m[c(col, pivot), ] swaps two rows by indexing. Each elimination step
subtracts a whole row at once.
Base R could do all of this with solve(), or with lm() on
poly(speed, 8, raw = TRUE); we avoided both. On this data,
lm gives the same train and test errors to six decimals, even on raw
speeds, because it uses a QR decomposition rather than the normal
equations.
# Train and test error for polynomials of degree 1 to 8, using only the
# standard library.
#
# Each polynomial is fitted by least squares to the training drives only:
# we build the normal equations and solve them by Gaussian elimination. Then
# we measure its root-mean-square error on the training drives and on the
# test drives it never saw.
MAX_DEGREE = 8
# The csv library is a bundled gem since Ruby 3.4, so split lines by hand
# (no quoted fields here). Each drive becomes [speed, litres].
def read(path)
header, *lines = File.readlines(path, chomp: true)
columns = header.split(",")
i_speed = columns.index("speed")
i_litres = columns.index("litres")
i_split = columns.index("split")
train = []
test = []
lines.each do |line|
fields = line.split(",")
point = [Float(fields[i_speed]), Float(fields[i_litres])]
(fields[i_split] == "train" ? train : test) << point
end
[train, test]
end
# Solve a x = b by Gaussian elimination with partial pivoting.
def solve(a, b)
n = b.length
m = a.each_with_index.map { |row, i| row + [b[i]] } # augmented matrix, copied
(0...n).each do |col|
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small.
pivot = (col...n).max_by { |r| m[r][col].abs }
m[col], m[pivot] = m[pivot], m[col]
(col + 1...n).each do |r|
factor = m[r][col] / m[col][col]
(col..n).each { |c| m[r][c] -= factor * m[col][c] }
end
end
# Back substitution, from the last unknown up.
x = Array.new(n, 0.0)
(n - 1).downto(0) do |r|
rest = (r + 1...n).sum { |c| m[r][c] * x[c] }
x[r] = (m[r][n] - rest) / m[r][r]
end
x
end
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
# of Z is 1, z_i, z_i^2, ..., z_i^degree.
def fit(points, degree, scale)
size = degree + 1
a = Array.new(size) { Array.new(size, 0.0) }
b = Array.new(size, 0.0)
points.each do |x, y|
z = scale.call(x)
size.times do |i|
b[i] += y * z**i
size.times { |j| a[i][j] += z**(i + j) }
end
end
solve(a, b)
end
def rmse(coefs, points, scale)
total = points.sum do |x, y|
z = scale.call(x)
predicted = coefs.each_with_index.sum { |c, i| c * z**i }
(predicted - y)**2
end
Math.sqrt(total / points.length)
end
train, test = read(ARGV[0])
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately.
low, high = train.map(&:first).minmax
mid = (low + high) / 2
half = (high - low) / 2
scale = ->(x) { (x - mid) / half }
puts "train_count #{train.length}"
puts "test_count #{test.length}"
errors = (1..MAX_DEGREE).map do |degree|
coefs = fit(train, degree, scale)
[degree, rmse(coefs, train, scale), rmse(coefs, test, scale)]
end
errors.each { |degree, train_error, _| puts format("train_rmse_d%d %.6f", degree, train_error) }
errors.each { |degree, _, test_error| puts format("test_rmse_d%d %.6f", degree, test_error) }Each drive becomes a [speed, litres] pair, sorted into train and test as
the lines are split. For each degree we fill the normal equations, then
solve runs Gaussian elimination on an augmented copy. At each column,
max_by finds the row with the largest entry, and a parallel assignment,
m[a], m[b] = m[b], m[a], swaps the two rows.
minmax gives the training range in one call, and a lambda maps each speed
onto −1 to 1. The matrix gem’s Matrix#lup could solve the system
directly, but it’s a bundled gem, not part of the standard library, so we
leave it out.
//! Train and test error for polynomials of degree 1 to 8, using only the
//! standard library.
//!
//! Each polynomial is fitted by least squares to the training drives only:
//! we build the normal equations and solve them by Gaussian elimination.
//! Then we measure its root-mean-square error on the training drives and on
//! the test drives it never saw.
use std::env;
use std::fs;
const MAX_DEGREE: usize = 8;
/// Each drive becomes (speed, litres), sorted into train and test. The file
/// has no quoted fields, so commas suffice; columns are found by name.
fn read(path: &str) -> (Vec<(f64, f64)>, Vec<(f64, f64)>) {
let text = fs::read_to_string(path).expect("cannot read data file");
let mut lines = text.lines();
let columns: Vec<&str> = lines.next().expect("empty file").split(',').collect();
let column = |name| columns.iter().position(|&c| c == name).expect("missing column");
let (i_speed, i_litres, i_split) = (column("speed"), column("litres"), column("split"));
let (mut train, mut test) = (Vec::new(), Vec::new());
for line in lines {
let fields: Vec<&str> = line.split(',').collect();
let speed: f64 = fields[i_speed].parse().expect("speed must be a number");
let litres: f64 = fields[i_litres].parse().expect("litres must be a number");
if fields[i_split] == "train" {
train.push((speed, litres));
} else {
test.push((speed, litres));
}
}
(train, test)
}
/// Solve a x = b by Gaussian elimination with partial pivoting.
fn solve(a: &[Vec<f64>], b: &[f64]) -> Vec<f64> {
let n = b.len();
// The augmented matrix, a fresh copy so a and b are untouched.
let mut m: Vec<Vec<f64>> = a
.iter()
.zip(b)
.map(|(row, &bi)| row.iter().copied().chain([bi]).collect())
.collect();
for col in 0..n {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small. f64 isn't Ord, so compare with total_cmp.
let pivot = (col..n)
.max_by(|&r, &s| m[r][col].abs().total_cmp(&m[s][col].abs()))
.unwrap();
m.swap(col, pivot);
for r in col + 1..n {
let factor = m[r][col] / m[col][col];
for c in col..=n {
m[r][c] -= factor * m[col][c];
}
}
}
// Back substitution, from the last unknown up.
let mut x = vec![0.0; n];
for r in (0..n).rev() {
let rest: f64 = (r + 1..n).map(|c| m[r][c] * x[c]).sum();
x[r] = (m[r][n] - rest) / m[r][r];
}
x
}
/// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i
/// of Z is 1, z_i, z_i^2, ..., z_i^degree.
fn fit(points: &[(f64, f64)], degree: usize, scale: impl Fn(f64) -> f64) -> Vec<f64> {
let size = degree + 1;
let mut a = vec![vec![0.0; size]; size];
let mut b = vec![0.0; size];
for &(x, y) in points {
let z = scale(x);
for i in 0..size {
b[i] += y * z.powi(i as i32);
for j in 0..size {
a[i][j] += z.powi((i + j) as i32);
}
}
}
solve(&a, &b)
}
fn rmse(coefs: &[f64], points: &[(f64, f64)], scale: impl Fn(f64) -> f64) -> f64 {
let total: f64 = points
.iter()
.map(|&(x, y)| {
let z = scale(x);
let predicted: f64 = coefs.iter().enumerate().map(|(i, c)| c * z.powi(i as i32)).sum();
(predicted - y).powi(2)
})
.sum();
(total / points.len() as f64).sqrt()
}
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let (train, test) = read(&path);
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power
// (the largest entry at degree 8) would swamp the small ones and make
// the equations hard to solve accurately.
let low = train.iter().map(|p| p.0).fold(f64::INFINITY, f64::min);
let high = train.iter().map(|p| p.0).fold(f64::NEG_INFINITY, f64::max);
let (mid, half) = ((low + high) / 2.0, (high - low) / 2.0);
let scale = |x: f64| (x - mid) / half;
println!("train_count {}", train.len());
println!("test_count {}", test.len());
let errors: Vec<(f64, f64)> = (1..=MAX_DEGREE)
.map(|degree| {
let coefs = fit(&train, degree, scale);
(rmse(&coefs, &train, scale), rmse(&coefs, &test, scale))
})
.collect();
for (degree, (train_error, _)) in (1..).zip(&errors) {
println!("train_rmse_d{degree} {train_error:.6}");
}
for (degree, (_, test_error)) in (1..).zip(&errors) {
println!("test_rmse_d{degree} {test_error:.6}");
}
}Each drive becomes an (f64, f64) tuple, and the matrix is a
Vec<Vec<f64>>, so m.swap(col, pivot) exchanges two rows in one call.
Finding the pivot needs care: f64 doesn’t implement Ord, because NaN
has no place in a total order, so max_by_key on an f64 doesn’t
compile. We use max_by with total_cmp instead.
Powers come from powi, which raises to an i32 exponent. The scaling
closure is passed to fit and rmse as impl Fn(f64) -> f64. No crates
are installed in this image, and the run uses rustc directly.
// Train and test error for polynomials of degree 1 to 8, using only Swift's
// standard library and Foundation, which ships with it.
//
// Each polynomial is fitted by least squares to the training drives only: we
// build the normal equations and solve them by Gaussian elimination. Then we
// measure its root-mean-square error on the training drives and on the test
// drives it never saw.
import Foundation
let maxDegree = 8
func six(_ x: Double) -> String { String(format: "%.6f", x) }
// Read speed,litres,split lines into two lists of (x, y) tuples, finding each
// column by name. split skips empty pieces, so a blank last line doesn't
// become a row. Double(_:) gives an optional, which ! unwraps.
typealias Point = (x: Double, y: Double)
let text = try String(contentsOfFile: CommandLine.arguments[1], encoding: .utf8)
var lines = text.split(separator: "\n")
let header = lines.removeFirst().split(separator: ",").map(String.init)
let speedCol = header.firstIndex(of: "speed")!
let litresCol = header.firstIndex(of: "litres")!
let splitCol = header.firstIndex(of: "split")!
var train: [Point] = []
var test: [Point] = []
for line in lines {
let fields = line.split(separator: ",")
let point = (x: Double(fields[speedCol])!, y: Double(fields[litresCol])!)
if fields[splitCol] == "train" {
train.append(point)
} else {
test.append(point)
}
}
// Solve a x = b by Gaussian elimination with partial pivoting. Arrays are
// values in Swift, so m is our own copy and changing it leaves a alone.
func solve(_ a: [[Double]], _ b: [Double]) -> [Double] {
let n = b.count
var m = (0..<n).map { a[$0] + [b[$0]] } // augmented matrix
for col in 0..<n {
// Swap up the row with the largest entry in this column, which keeps
// the rounding error small.
let pivot = (col..<n).max { abs(m[$0][col]) < abs(m[$1][col]) }!
m.swapAt(col, pivot)
for r in (col + 1)..<n {
let factor = m[r][col] / m[col][col]
for c in col...n {
m[r][c] -= factor * m[col][c]
}
}
}
// Back substitution, from the last unknown up.
var x = [Double](repeating: 0, count: n)
for r in (0..<n).reversed() {
var s = m[r][n]
for c in (r + 1)..<n {
s -= m[r][c] * x[c]
}
x[r] = s / m[r][r]
}
return x
}
// Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
// Z is 1, z_i, z_i^2, ..., z_i^degree.
func fit(_ points: [Point], _ degree: Int, _ scale: (Double) -> Double) -> [Double] {
let size = degree + 1
var a = [[Double]](repeating: [Double](repeating: 0, count: size), count: size)
var b = [Double](repeating: 0, count: size)
for p in points {
let z = scale(p.x)
for i in 0..<size {
b[i] += p.y * pow(z, Double(i))
for j in 0..<size {
a[i][j] += pow(z, Double(i + j))
}
}
}
return solve(a, b)
}
func rmse(_ coefs: [Double], _ points: [Point], _ scale: (Double) -> Double) -> Double {
var total = 0.0
for p in points {
let z = scale(p.x)
var predicted = 0.0
for (i, c) in coefs.enumerated() {
predicted += c * pow(z, Double(i))
}
total += (predicted - p.y) * (predicted - p.y)
}
return (total / Double(points.count)).squareRoot()
}
// Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
// largest entry at degree 8) would swamp the small ones and make the
// equations hard to solve accurately.
let low = train.map { $0.x }.min()!
let high = train.map { $0.x }.max()!
let mid = (low + high) / 2
let half = (high - low) / 2
func scale(_ x: Double) -> Double { (x - mid) / half }
print("train_count", train.count)
print("test_count", test.count)
var errors: [(train: Double, test: Double)] = []
for degree in 1...maxDegree {
let coefs = fit(train, degree, scale)
errors.append((train: rmse(coefs, train, scale), test: rmse(coefs, test, scale)))
}
for (i, e) in errors.enumerated() {
print("train_rmse_d\(i + 1)", six(e.train))
}
for (i, e) in errors.enumerated() {
print("test_rmse_d\(i + 1)", six(e.test))
}We read each drive into an (x, y) tuple and sort it into a training or
test list. fit adds up the normal equations in a [[Double]], and
solve appends b to each row, then runs Gaussian elimination with
partial pivoting. The pivot row is the largest entry found by max with
a comparison closure, and the swap is one call to swapAt.
Swift arrays are values, so the matrix solve changes is its own copy
and the caller’s equations stay untouched. Apple’s Accelerate framework
has a solver, but it isn’t available in our Linux image.
# Train and test error for polynomials of degree 1 to 8, using only Tcl's
# built-in commands.
#
# Each polynomial is fitted by least squares to the training drives only: we
# build the normal equations and solve them by Gaussian elimination. Then we
# measure its root-mean-square error on the training drives and on the test
# drives it never saw. A matrix is a list of row lists; lindex reads an
# entry and lset writes one in place.
set max_degree 8
# Returns two lists of {speed litres} pairs: the training drives and the
# test drives.
proc read_drives {path} {
set f [open $path r]
set header [split [gets $f] ,]
set speed_col [lsearch -exact $header speed]
set litres_col [lsearch -exact $header litres]
set split_col [lsearch -exact $header split]
set train {}
set test {}
while {[gets $f line] >= 0} {
set fields [split $line ,]
set point [list [lindex $fields $speed_col] [lindex $fields $litres_col]]
if {[lindex $fields $split_col] eq "train"} {
lappend train $point
} else {
lappend test $point
}
}
close $f
return [list $train $test]
}
# Solve a x = b by Gaussian elimination with partial pivoting.
proc solve {a b} {
set n [llength $b]
# The augmented matrix: b becomes each row's last entry, column n.
set m {}
foreach row $a value $b {
lappend m [concat $row $value]
}
for {set col 0} {$col < $n} {incr col} {
# Swap up the row with the largest entry in this column, which keeps
# the rounding error small.
set pivot $col
for {set r [expr {$col + 1}]} {$r < $n} {incr r} {
if {abs([lindex $m $r $col]) > abs([lindex $m $pivot $col])} {
set pivot $r
}
}
set top [lindex $m $col]
lset m $col [lindex $m $pivot]
lset m $pivot $top
for {set r [expr {$col + 1}]} {$r < $n} {incr r} {
set factor [expr {[lindex $m $r $col] / [lindex $m $col $col]}]
for {set c $col} {$c <= $n} {incr c} {
lset m $r $c [expr {[lindex $m $r $c] - $factor * [lindex $m $col $c]}]
}
}
}
# Back substitution, from the last unknown up.
set x [lrepeat $n 0.0]
for {set r [expr {$n - 1}]} {$r >= 0} {incr r -1} {
set total [lindex $m $r $n]
for {set c [expr {$r + 1}]} {$c < $n} {incr c} {
set total [expr {$total - [lindex $m $r $c] * [lindex $x $c]}]
}
lset x $r [expr {$total / [lindex $m $r $r]}]
}
return $x
}
# Least squares: the coefficients c solve (Z^T Z) c = Z^T y, where row i of
# Z is 1, z_i, z_i^2, ..., z_i^degree. zs and ys are the scaled speeds and
# the litres.
proc fit {zs ys degree} {
set size [expr {$degree + 1}]
set a [lrepeat $size [lrepeat $size 0.0]]
set b [lrepeat $size 0.0]
foreach z $zs y $ys {
for {set i 0} {$i < $size} {incr i} {
lset b $i [expr {[lindex $b $i] + $y * $z ** $i}]
for {set j 0} {$j < $size} {incr j} {
lset a $i $j [expr {[lindex $a $i $j] + $z ** ($i + $j)}]
}
}
}
return [solve $a $b]
}
proc rmse {coefs zs ys} {
set total 0.0
foreach z $zs y $ys {
set predicted 0.0
set i 0
foreach c $coefs {
set predicted [expr {$predicted + $c * $z ** $i}]
incr i
}
set total [expr {$total + ($predicted - $y) ** 2}]
}
return [expr {sqrt($total / [llength $ys])}]
}
lassign [read_drives [lindex $argv 0]] train test
# Map the training speeds onto -1 to 1. Raw speeds to the 16th power (the
# largest entry at degree 8) would swamp the small ones and make the
# equations hard to solve accurately. The speeds are whole numbers, and
# dividing two integers would round down, so we divide by 2.0.
set train_speeds [lmap point $train { lindex $point 0 }]
set low [tcl::mathfunc::min {*}$train_speeds]
set high [tcl::mathfunc::max {*}$train_speeds]
set mid [expr {($low + $high) / 2.0}]
set half [expr {($high - $low) / 2.0}]
proc scaled_speeds {points mid half} {
return [lmap point $points { expr {([lindex $point 0] - $mid) / $half} }]
}
proc litres {points} {
return [lmap point $points { lindex $point 1 }]
}
set z_train [scaled_speeds $train $mid $half]
set z_test [scaled_speeds $test $mid $half]
set y_train [litres $train]
set y_test [litres $test]
puts "train_count [llength $train]"
puts "test_count [llength $test]"
set train_errors {}
set test_errors {}
for {set degree 1} {$degree <= $max_degree} {incr degree} {
set coefs [fit $z_train $y_train $degree]
lappend train_errors [rmse $coefs $z_train $y_train]
lappend test_errors [rmse $coefs $z_test $y_test]
}
set degree 1
foreach e $train_errors {
puts "train_rmse_d$degree [format %.6f $e]"
incr degree
}
set degree 1
foreach e $test_errors {
puts "test_rmse_d$degree [format %.6f $e]"
incr degree
}The augmented matrix is a list of row lists. lindex $m $r $c reads an
entry and lset m $r $c writes one in place, and swapping two rows is two
lset calls on whole rows. Lists are values, so solve changes only its
own copy, never the caller’s.
The speeds are whole numbers, and Tcl’s / rounds down when both sides are
integers, so the scale divides by 2.0. There’s no solver to skip: Tcl’s
own commands have no matrix type, and Tcllib’s math::linearalgebra, which
does, isn’t installed.