Concept · advanced
Principal component analysis
Find the few directions that capture most of the variation in many columns.
Principal component analysis (PCA) finds the directions in which your data varies most. When several columns move together, a few of these directions, the principal components, can stand in for all of them. PCA is used to compress data, to plot many columns in two dimensions, and to find the hidden factors behind a set of measurements.
The dataset
data.csv holds 100 synthetic students’ scores, from 0 to 100, in math,
physics, reading, and writing. We drew it once with make_data.py
(fixed seed 20261103) from two hidden traits: a general ability that lifts
all four scores, and a lean that raises math and physics while lowering
reading and writing, or the reverse. Then we added noise. PCA never sees
those traits; it has only the scores.
The method
- Standardize each column to mean 0 and (sample) standard deviation 1, so that no subject counts more just because its scores spread wider.
- Build the correlation matrix, a 4×4 table of how strongly each pair of subjects moves together. Math and physics correlate at 0.824, and reading and writing at 0.767, but math and reading only at 0.278.
- Find its top eigenvector by power iteration. Start with , then repeat . Each multiplication stretches most along the direction the matrix stretches most, so turns toward it. The eigenvalue, how much variance lies along that direction, is .
- Deflate: subtract from the matrix, which removes the first component, then repeat step 3 for the second.
- Fix the sign. If is an eigenvector, so is , so we agree to make each one’s first entry positive. Starting from already lands on that version, but every implementation checks anyway, so the rule doesn’t depend on the starting point.
We run 1,000 steps, which is generous. After one step, is , still leaning toward where it started. After 5 steps it’s , and from step 18 on, all six decimals stay fixed. The second component settles even faster, in 9 steps, because its eigenvalue stands further above the next one.
Reading the results
| Component 1 | Component 2 | |
|---|---|---|
| Eigenvalue | 2.496154 | 1.110197 |
| Share of variance | 0.624 | 0.278 |
| math | 0.504 | 0.498 |
| physics | 0.518 | 0.469 |
| reading | 0.457 | −0.585 |
| writing | 0.518 | −0.437 |
The four standardized columns have a total variance of 4, so the first component carries of it, and the two together carry 90.2%. Two numbers per student keep nine-tenths of what four did.
The components are the hidden traits. The first weights every subject about equally: that’s general ability, and the true direction was 0.5 for each. The second gives math and physics positive weights and reading and writing negative ones: that’s the lean, and the true direction was .
Take the first student, who scored 82, 91, 74, and 72. Standardized, that’s . Multiplying by each component gives 1.712 on the first, well above average overall, and 1.114 on the second, leaning toward math and physics.
What to watch for across languages
None of the from-scratch implementations calls a ready-made PCA or
eigenvalue routine. R, Julia, Octave, and Ada have them, and their tabs
have a Built-in switch that shows R’s prcomp(), Julia’s eigen,
Octave’s eig, and Ada’s Eigensystem. All four find the same
eigenvalues as power iteration, but R returns the first component with
every sign flipped, Julia flips both, and Ada flips the second. That’s why
the sign rule matters: with it, the built-ins print the same results.
Octave’s eig lists the eigenvalues smallest first, so it sorts them.
Ruby’s
Matrix#eigensystem would work too, but matrix is a bundled gem rather
than part of the standard library, so it’s left out.
Watch how each holds a matrix. R, Julia, Octave, Fortran, and APL have
matrices built in, and Ada has them in its standard library, so a
power-iteration step is one product: %*%, * in Julia, Octave, and Ada,
matmul, or +.×. APL also needs no loop, since ⍣ repeats the step 1,000 times,
and neither does Haskell, whose iterate builds the endless list of steps.
Python, Ruby, Java, JavaScript, Swift, Haskell, OCaml, Elixir, Common Lisp, Prolog, Tcl, Lua, Perl, and PHP use arrays or lists of rows with
small helper functions, and Go, Rust, C, C++, and Pascal use fixed-size arrays whose length is part of the type. BASIC uses 2-D arrays sized by DIM, with a loop for every product. COBOL nests one
table inside another, awk keys a single array by m[i, j], and PostgreSQL
stores one row per entry, so a matrix-vector product becomes a GROUP BY.
And watch COBOL’s precision. It works in fixed-point decimals, and its results match the floating-point languages to 14 decimal places. Power iteration corrects its own rounding at every step, so even 9 decimal places would pass.
Data preview
| math | physics | reading | writing |
|---|---|---|---|
| 82 | 91 | 74 | 72 |
| 79 | 75 | 57 | 62 |
| 60 | 71 | 64 | 77 |
| 76 | 74 | 65 | 68 |
| 75 | 75 | 53 | 58 |
| 77 | 83 | 74 | 80 |
| 99 | 88 | 100 | 100 |
| 62 | 70 | 80 | 79 |
First 8 of 100 rows · Download data.csv
Expected results
Every implementation must print these, to within a relative tolerance of 10⁻⁴.
count | 100 |
eigenvalue_1 | 2.496154 |
explained_1 | 0.624038 |
pc1_math | 0.504100 |
pc1_physics | 0.517770 |
pc1_reading | 0.457228 |
pc1_writing | 0.518401 |
eigenvalue_2 | 1.110197 |
explained_2 | 0.277549 |
pc2_math | 0.497992 |
pc2_physics | 0.468590 |
pc2_reading | -0.584599 |
pc2_writing | -0.436659 |
student_1_pc1 | 1.712316 |
student_1_pc2 | 1.113676 |
Implementations
"""Principal component analysis by power iteration, using only the standard library.
Standardizes the four scores, builds their correlation matrix, and finds the
top two eigenvectors (the principal components): power iteration finds the
first, then deflation removes it from the matrix so the same steps find the
second.
"""
import csv
import math
import sys
SUBJECTS = ("math", "physics", "reading", "writing")
STEPS = 1000
def mat_vec(m, v):
return [sum(a * b for a, b in zip(row, v)) for row in m]
def dot(u, v):
return sum(a * b for a, b in zip(u, v))
def normalize(v):
length = math.sqrt(dot(v, v))
return [x / length for x in v]
def top_eigenvector(m):
"""Power iteration: multiply by the matrix and rescale, STEPS times."""
v = [1.0] + [0.0] * (len(m) - 1)
for _ in range(STEPS):
v = normalize(mat_vec(m, v))
# An eigenvector can point either way; make its first entry positive.
if v[0] < 0:
v = [-x for x in v]
return dot(v, mat_vec(m, v)), v # eigenvalue, eigenvector
def main():
with open(sys.argv[1], newline="") as f:
rows = [[float(row[s]) for s in SUBJECTS] for row in csv.DictReader(f)]
n, k = len(rows), len(SUBJECTS)
# Standardize each column to mean 0 and (sample) standard deviation 1.
columns = list(zip(*rows))
means = [sum(c) / n for c in columns]
sds = [math.sqrt(sum((x - m) ** 2 for x in c) / (n - 1)) for c, m in zip(columns, means)]
z = [[(x - m) / s for x, m, s in zip(row, means, sds)] for row in rows]
# The correlation matrix: the covariance of the standardized columns.
corr = [[sum(r[i] * r[j] for r in z) / (n - 1) for j in range(k)] for i in range(k)]
value_1, pc_1 = top_eigenvector(corr)
# Deflation: subtract the first component's share, leaving the rest.
deflated = [[corr[i][j] - value_1 * pc_1[i] * pc_1[j] for j in range(k)] for i in range(k)]
value_2, pc_2 = top_eigenvector(deflated)
print(f"count {n}")
for c, (value, pc) in enumerate([(value_1, pc_1), (value_2, pc_2)], start=1):
print(f"eigenvalue_{c} {value:.6f}")
# Standardized columns each have variance 1, so the total is k.
print(f"explained_{c} {value / k:.6f}")
for subject, loading in zip(SUBJECTS, pc):
print(f"pc{c}_{subject} {loading:.6f}")
print(f"student_1_pc1 {dot(z[0], pc_1):.6f}")
print(f"student_1_pc2 {dot(z[0], pc_2):.6f}")
if __name__ == "__main__":
main()Python’s standard library has no matrix type and no eigenvalue routine, and
NumPy, where numpy.linalg.eigh would do this in one call, isn’t
installed. So the matrices are lists of lists, and two small helpers do
the linear algebra: mat_vec multiplies a matrix by a vector with one
sum per row, and dot multiplies two vectors. Everything else is built
from those two.
The correlation matrix is one nested list comprehension, and deflation is
another. top_eigenvector runs for both components unchanged: the deflated
matrix just has a different largest eigenvalue.
-- Principal component analysis by power iteration, using only Ada's
-- standard library.
--
-- Standardizes the four scores, builds their correlation matrix, and finds
-- the top two eigenvectors (the principal components): power iteration
-- finds the first, then deflation removes it from the matrix so the same
-- steps find the second.
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
K : constant := 4;
Steps : constant := 1000;
Subjects : constant array (1 .. K) of String (1 .. 7) :=
["math ", "physics", "reading", "writing"]; -- padded to one length
-- Real_Vector and Real_Matrix come from Annex G of the standard, with
-- "*" for matrix and vector products and "abs" for a vector's length.
subtype Row is Real_Vector (1 .. K);
package Row_Vectors is new Ada.Containers.Vectors (Positive, Row);
-- 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;
-- Read the four subject columns into an n by k matrix.
function Read_Scores (Path : String) return Real_Matrix is
File : File_Type;
Rows : Row_Vectors.Vector;
Cols : array (1 .. K) of Natural := [others => 0];
begin
Open (File, In_File, Path);
declare
Header : constant String := Get_Line (File);
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
for J in 1 .. K loop
if Field (Header, C) = Trim (Subjects (J), Ada.Strings.Right)
then
Cols (J) := C;
end if;
end loop;
end loop;
end;
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
Scores : Row;
begin
for J in 1 .. K loop
Scores (J) := Long_Float'Value (Field (Line, Cols (J)));
end loop;
Rows.Append (Scores);
end;
end loop;
Close (File);
return Result : Real_Matrix (1 .. Natural (Rows.Length), 1 .. K) do
for I in Result'Range (1) loop
for J in 1 .. K loop
Result (I, J) := Rows (I) (J);
end loop;
end loop;
end return;
end Read_Scores;
-- Power iteration: multiply by the matrix and rescale, Steps times.
procedure Top_Eigenvector (M : Real_Matrix; Value : out Long_Float;
V : out Row) is
begin
V := [1.0, 0.0, 0.0, 0.0];
for Step in 1 .. Steps loop
V := M * V;
V := V / abs V;
end loop;
-- An eigenvector can point either way; make its first entry positive.
if V (1) < 0.0 then
V := -V;
end if;
Value := V * (M * V); -- the Rayleigh quotient, an inner product
end Top_Eigenvector;
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;
procedure Report (C : Positive; Value : Long_Float; PC : Row) is
-- C'Image is " 1", so the slice drops the leading space.
Tag : constant String := C'Image (2 .. C'Image'Last);
begin
Say ("eigenvalue_" & Tag, Value);
-- Standardized columns each have variance 1, so the total is k.
Say ("explained_" & Tag, Value / Long_Float (K));
for J in 1 .. K loop
Say ("pc" & Tag & "_" & Trim (Subjects (J), Ada.Strings.Right),
PC (J));
end loop;
end Report;
Z : Real_Matrix := Read_Scores (Ada.Command_Line.Argument (1));
N : constant Positive := Z'Length (1);
Value_1, Value_2 : Long_Float;
PC_1, PC_2 : Row;
begin
-- Standardize each column to mean 0 and (sample) standard deviation 1.
for J in 1 .. K loop
declare
Total, Squares : Long_Float := 0.0;
Mean, Sd : Long_Float;
begin
for I in 1 .. N loop
Total := Total + Z (I, J);
end loop;
Mean := Total / Long_Float (N);
for I in 1 .. N loop
Squares := Squares + (Z (I, J) - Mean) ** 2;
end loop;
Sd := Sqrt (Squares / Long_Float (N - 1));
for I in 1 .. N loop
Z (I, J) := (Z (I, J) - Mean) / Sd;
end loop;
end;
end loop;
declare
-- The correlation matrix: the covariance of the standardized columns.
-- Transpose (Z) is k by n, so the product is k by k.
Corr : constant Real_Matrix := Transpose (Z) * Z / Long_Float (N - 1);
begin
Top_Eigenvector (Corr, Value_1, PC_1);
-- Deflation: subtract the first component's share, leaving the rest.
-- "*" of two vectors gives a matrix (the outer product) when a matrix
-- is wanted, and a number (the inner product) when a number is.
Top_Eigenvector (Corr - Value_1 * (PC_1 * PC_1), Value_2, PC_2);
end;
Put_Line ("count" & N'Image);
Report (1, Value_1, PC_1);
Report (2, Value_2, PC_2);
declare
First : Row; -- the first student's standardized scores
begin
for J in 1 .. K loop
First (J) := Z (1, J);
end loop;
Say ("student_1_pc1", First * PC_1);
Say ("student_1_pc2", First * PC_2);
end;
end Main;Ada’s Annex G gives us Real_Vector and Real_Matrix, and we use their
* as a general tool, the way Fortran uses matmul: Transpose (Z) * Z
builds the correlation matrix, M * V is one power-iteration step, and
abs V is a vector’s length. The same * on two vectors returns a number
or a matrix, whichever the context asks for, so PC_1 * PC_1 is the
outer product that deflation subtracts. Annex G also has Eigensystem,
which finds the same components up to sign, but it performs the method
itself, so we leave it out. Alire crates aren’t installed.
-- Principal component analysis with Eigensystem, from Ada's standard
-- library.
--
-- The built-in version of ../ada/main.adb. Eigensystem finds every
-- eigenvalue and eigenvector of the correlation matrix at once, in place of
-- power iteration and deflation.
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
K : constant := 4;
Subjects : constant array (1 .. K) of String (1 .. 7) :=
["math ", "physics", "reading", "writing"]; -- padded to one length
-- Real_Vector and Real_Matrix come from Annex G of the standard, with
-- "*" for matrix and vector products and "abs" for a vector's length.
subtype Row is Real_Vector (1 .. K);
package Row_Vectors is new Ada.Containers.Vectors (Positive, Row);
-- 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;
-- Read the four subject columns into an n by k matrix.
function Read_Scores (Path : String) return Real_Matrix is
File : File_Type;
Rows : Row_Vectors.Vector;
Cols : array (1 .. K) of Natural := [others => 0];
begin
Open (File, In_File, Path);
declare
Header : constant String := Get_Line (File);
begin
for C in 1 .. Ada.Strings.Fixed.Count (Header, ",") + 1 loop
for J in 1 .. K loop
if Field (Header, C) = Trim (Subjects (J), Ada.Strings.Right)
then
Cols (J) := C;
end if;
end loop;
end loop;
end;
while not End_Of_File (File) loop
declare
Line : constant String := Get_Line (File);
Scores : Row;
begin
for J in 1 .. K loop
Scores (J) := Long_Float'Value (Field (Line, Cols (J)));
end loop;
Rows.Append (Scores);
end;
end loop;
Close (File);
return Result : Real_Matrix (1 .. Natural (Rows.Length), 1 .. K) do
for I in Result'Range (1) loop
for J in 1 .. K loop
Result (I, J) := Rows (I) (J);
end loop;
end loop;
end return;
end Read_Scores;
-- Column C of Vectors, with the sign rule applied: an eigenvector can
-- point either way, so make its first entry positive.
function Component (Vectors : Real_Matrix; C : Positive) return Row is
V : Row;
begin
for J in 1 .. K loop
V (J) := Vectors (J, C);
end loop;
return (if V (1) < 0.0 then -V else V);
end Component;
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;
procedure Report (C : Positive; Value : Long_Float; PC : Row) is
-- C'Image is " 1", so the slice drops the leading space.
Tag : constant String := C'Image (2 .. C'Image'Last);
begin
Say ("eigenvalue_" & Tag, Value);
-- Standardized columns each have variance 1, so the total is k.
Say ("explained_" & Tag, Value / Long_Float (K));
for J in 1 .. K loop
Say ("pc" & Tag & "_" & Trim (Subjects (J), Ada.Strings.Right),
PC (J));
end loop;
end Report;
Z : Real_Matrix := Read_Scores (Ada.Command_Line.Argument (1));
N : constant Positive := Z'Length (1);
Value_1, Value_2 : Long_Float;
PC_1, PC_2 : Row;
begin
-- Standardize each column to mean 0 and (sample) standard deviation 1.
for J in 1 .. K loop
declare
Total, Squares : Long_Float := 0.0;
Mean, Sd : Long_Float;
begin
for I in 1 .. N loop
Total := Total + Z (I, J);
end loop;
Mean := Total / Long_Float (N);
for I in 1 .. N loop
Squares := Squares + (Z (I, J) - Mean) ** 2;
end loop;
Sd := Sqrt (Squares / Long_Float (N - 1));
for I in 1 .. N loop
Z (I, J) := (Z (I, J) - Mean) / Sd;
end loop;
end;
end loop;
declare
-- The correlation matrix: the covariance of the standardized columns.
-- Transpose (Z) is k by n, so the product is k by k.
Corr : constant Real_Matrix := Transpose (Z) * Z / Long_Float (N - 1);
Values : Row;
Vectors : Real_Matrix (1 .. K, 1 .. K);
begin
-- Eigensystem, for a symmetric matrix, lists the eigenvalues largest
-- first, with each one's eigenvector in the matching column.
Eigensystem (Corr, Values, Vectors);
Value_1 := Values (1);
Value_2 := Values (2);
PC_1 := Component (Vectors, 1);
PC_2 := Component (Vectors, 2);
end;
Put_Line ("count" & N'Image);
Report (1, Value_1, PC_1);
Report (2, Value_2, PC_2);
declare
First : Row; -- the first student's standardized scores
begin
for J in 1 .. K loop
First (J) := Z (1, J);
end loop;
Say ("student_1_pc1", First * PC_1);
Say ("student_1_pc2", First * PC_2);
end;
end Main;Eigensystem comes with the matrix types in Annex G of the Ada standard,
so it’s part of the standard library. Given a symmetric matrix, it fills
Values with every eigenvalue, largest first, and Vectors with the
matching eigenvectors as columns, replacing both power iteration and
deflation.
The rest is the from-scratch version: the same loops standardize the
scores, and Transpose (Z) * Z builds the correlation matrix. Eigensystem
returns the second component with every sign flipped, so Component
applies the sign rule, making each one’s first entry positive.
⍝ Principal component analysis by power iteration, from scratch in GNU APL.
⍝ Standardizes the four scores, builds their correlation matrix, and finds
⍝ the top two eigenvectors (the principal components): power iteration finds
⍝ the first, then deflation removes it from the matrix so the same steps find
⍝ the second. APL reads right to left, with no operator precedence: 2×3+4 is
⍝ 2×(3+4). ⍝ starts a comment; ← assigns.
∇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 ¯
∇
STEPS←1000
SUBJECTS←'math' 'physics' 'reading' 'writing'
⍝ 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↓); ⍎¨ makes numbers.
COLUMN←{⍎¨1↓ROWS[;ROWS[1;]⍳⊂⍵]}
X←⍉⊃COLUMN¨SUBJECTS ⍝ one row per student, one column per subject
N←≢X ⋄ K←≢SUBJECTS ⍝ ≢ counts the items
⍝ Standardize each column to mean 0 and sample standard deviation 1. ⌿ works
⍝ down the columns: +⌿ sums each one. (⍴X)⍴ repeats a row of column values
⍝ on every row, and *0.5 is the square root.
C←X-(⍴X)⍴(+⌿X)÷N
Z←C÷(⍴X)⍴((+⌿C*2)÷N-1)*0.5
⍝ The correlation matrix. ⍉ transposes, and the matrix product (⍉Z)+.×Z
⍝ sums z_i × z_j over the students for every pair of columns at once.
CORR←((⍉Z)+.×Z)÷N-1
NORMALIZE←{⍵÷(+/⍵×⍵)*0.5} ⍝ divide by the length
∇V←TOP M;START
⍝ Power iteration: multiply by M and rescale, STEPS times. The power
⍝ operator ⍣ repeats the step, starting from 1 0 0 0: (≢M)↑1 pads 1 with
⍝ zeros to M's size.
START←(≢M)↑1
V←({NORMALIZE M+.×⍵}⍣STEPS)START
⍝ An eigenvector can point either way; make its first entry positive.
⍝ ¯1*1 is ¯1 and ¯1*0 is 1, so this flips V only when V[1]<0.
V←Vׯ1*V[1]<0
∇
⍝ The eigenvalue is the Rayleigh quotient v·(M v).
EIGENVALUE←{⍵+.×⍺+.×⍵}
PC1←TOP CORR
VALUE1←CORR EIGENVALUE PC1
⍝ Deflation: subtract the first component's share, leaving the rest.
⍝ The outer product PC1∘.×PC1 is the K-by-K table of v_i × v_j.
DEFLATED←CORR-VALUE1×PC1∘.×PC1
PC2←TOP DEFLATED
VALUE2←DEFLATED EIGENVALUE PC2
'count' SAY N
⍝ Standardized columns each have variance 1, so the total is K.
'eigenvalue_1' SAY6 VALUE1
'explained_1' SAY6 VALUE1÷K
⍝ Each label is 'pc1_' joined (,) to a subject name, and SAY6¨ pairs each
⍝ label with its loading.
({'pc1_',⍵}¨SUBJECTS) SAY6¨PC1
'eigenvalue_2' SAY6 VALUE2
'explained_2' SAY6 VALUE2÷K
({'pc2_',⍵}¨SUBJECTS) SAY6¨PC2
⍝ A student's score on a component is the dot product z·v.
'student_1_pc1' SAY6 Z[1;]+.×PC1
'student_1_pc2' SAY6 Z[1;]+.×PC2
)OFFEvery matrix step is one APL product. The inner product (⍉Z)+.×Z builds
the whole correlation matrix, M+.×⍵ multiplies it by the current vector,
and the outer product PC1∘.×PC1 gives the table of vᵢ × vⱼ that
deflation subtracts. Power iteration needs no loop: ⍣ applies the
multiply-and-rescale step 1000 times from 1 0 0 0.
One parsing surprise: in GNU APL, (f⍣STEPS)(≢M)↑1 raises a DOMAIN ERROR,
though the same line with a literal 1000 runs. We name the start vector
START first. The sign flip has no if: Vׯ1*V[1]<0 multiplies by ¯1
to the power 1 or 0.
# Principal component analysis by power iteration, in POSIX awk.
#
# Standardizes the four scores, builds their correlation matrix, and finds the
# top two eigenvectors (the principal components): power iteration finds the
# first, then deflation removes it from the matrix so the same steps find the
# second. Rows are stored as they are read; the analysis happens in END.
BEGIN { FS = ","; K = 4; STEPS = 1000; split("math physics reading writing", subject, " ") }
# The header line: find each subject's column by name.
NR == 1 {
for (i = 1; i <= NF; i++)
for (j = 1; j <= K; j++)
if ($i == subject[j]) col[j] = i
for (j = 1; j <= K; j++) {
if (!col[j]) {
print "missing a column (math, physics, reading, writing)" > "/dev/stderr"
failed = 1
exit 1
}
}
next
}
# x[i, j] is subject j of row i; awk fakes 2-D arrays with a joined key.
{
n++
for (j = 1; j <= K; j++) x[n, j] = $col[j] + 0
}
# Power iteration on the K-by-K matrix m. Functions can't return arrays, so
# the eigenvector is written into v and the eigenvalue is returned. The extra
# parameters are awk's way of declaring local variables.
function top_eigenvector(m, v, step, i, j, mv, len, value) {
v[1] = 1
for (i = 2; i <= K; i++) v[i] = 0
for (step = 1; step <= STEPS; step++) {
# v <- m v, then rescale to length 1.
len = 0
for (i = 1; i <= K; i++) {
mv[i] = 0
for (j = 1; j <= K; j++) mv[i] += m[i, j] * v[j]
len += mv[i] ^ 2
}
len = sqrt(len)
for (i = 1; i <= K; i++) v[i] = mv[i] / len
}
# An eigenvector can point either way; make its first entry positive.
if (v[1] < 0)
for (i = 1; i <= K; i++) v[i] = -v[i]
# The eigenvalue is v . (m v).
value = 0
for (i = 1; i <= K; i++)
for (j = 1; j <= K; j++) value += v[i] * m[i, j] * v[j]
return value
}
END {
if (failed) exit 1
# Standardize each column to mean 0 and sample standard deviation 1
# (divide by n - 1), in place.
for (j = 1; j <= K; j++) {
total = 0
for (i = 1; i <= n; i++) total += x[i, j]
mean = total / n
ss = 0
for (i = 1; i <= n; i++) ss += (x[i, j] - mean) ^ 2
sd = sqrt(ss / (n - 1))
for (i = 1; i <= n; i++) x[i, j] = (x[i, j] - mean) / sd
}
# The correlation matrix: the covariance of the standardized columns.
for (a = 1; a <= K; a++) {
for (b = 1; b <= K; b++) {
total = 0
for (i = 1; i <= n; i++) total += x[i, a] * x[i, b]
corr[a, b] = total / (n - 1)
}
}
value[1] = top_eigenvector(corr, pc1)
# Deflation: subtract the first component's share, leaving the rest.
for (a = 1; a <= K; a++)
for (b = 1; b <= K; b++) deflated[a, b] = corr[a, b] - value[1] * pc1[a] * pc1[b]
value[2] = top_eigenvector(deflated, pc2)
printf "count %d\n", n
for (c = 1; c <= 2; c++) {
printf "eigenvalue_%d %.6f\n", c, value[c]
# Standardized columns each have variance 1, so the total is K.
printf "explained_%d %.6f\n", c, value[c] / K
for (j = 1; j <= K; j++)
printf "pc%d_%s %.6f\n", c, subject[j], (c == 1 ? pc1[j] : pc2[j])
}
# The first student's scores in the new coordinates.
s1 = s2 = 0
for (j = 1; j <= K; j++) { s1 += x[1, j] * pc1[j]; s2 += x[1, j] * pc2[j] }
printf "student_1_pc1 %.6f\n", s1
printf "student_1_pc2 %.6f\n", s2
}awk has no matrices, so the data, the correlation matrix, and its deflated
copy all live in arrays keyed like m[i, j], which awk joins into one
string key with SUBSEP. Every product is a written-out loop. Power
iteration is one function used twice. An awk function can’t return an
array, but arrays are passed by reference, so top_eigenvector fills in
the vector it’s given and returns the eigenvalue. Its extra parameters are
how awk declares local variables. awk has no linear algebra built in, so
there was nothing to skip. It sticks to POSIX awk and runs on mawk and on
macOS awk.
10 REM Principal component analysis by power iteration, in classic BASIC.
20 REM
30 REM Standardizes the four scores, builds their correlation matrix, and
40 REM finds the top two eigenvectors (the principal components): power
50 REM iteration finds the first, then deflation removes it from the
60 REM matrix so the same steps find the second.
70 REM
80 REM run.sh passes the data path in the DATA environment variable,
90 REM because bwbasic can't read command-line arguments.
100 K = 4 : STEPS = 1000
110 REM PC(1, J) and PC(2, J) hold the two components; EV(1), EV(2) their
120 REM eigenvalues. M is the matrix that power iteration works on.
130 DIM SUBJ$(4), COL(4), Z(200, 4), C(4, 4), M(4, 4), V(4), MV(4)
140 DIM PC(2, 4), EV(2), FL$(20)
150 SUBJ$(1) = "math" : SUBJ$(2) = "physics"
160 SUBJ$(3) = "reading" : SUBJ$(4) = "writing"
170 GOSUB 1000 : REM read the four columns into Z(1..N, 1..4)
180 GOSUB 2000 : REM standardize them and fill the correlation matrix C
190 FOR A = 1 TO K
200 FOR B = 1 TO K
210 M(A, B) = C(A, B)
220 NEXT B
230 NEXT A
240 CI = 1 : GOSUB 5000 : REM power iteration on M fills PC(1, ...) and EV(1)
250 REM Deflation: subtract the first component's share, leaving the rest.
260 FOR A = 1 TO K
270 FOR B = 1 TO K
280 M(A, B) = C(A, B) - EV(1) * PC(1, A) * PC(1, B)
290 NEXT B
300 NEXT A
310 CI = 2 : GOSUB 5000
320 PRINT "count"; N
330 CI = 1 : GOSUB 6000 : REM print the first component
340 CI = 2 : GOSUB 6000
350 REM The first student's scores, projected onto each component.
360 S1 = 0 : S2 = 0
370 FOR A = 1 TO K
380 S1 = S1 + Z(1, A) * PC(1, A) : S2 = S2 + Z(1, A) * PC(2, A)
390 NEXT A
400 L$ = "student_1_pc1" : Y = S1 : GOSUB 3000
410 L$ = "student_1_pc2" : Y = S2 : GOSUB 3000
420 SYSTEM : REM without SYSTEM, bwbasic waits at its prompt
1000 REM Read the columns named in SUBJ$ into Z(1..N, 1..4).
1010 F$ = ENVIRON$("DATA")
1020 OPEN "I", #1, F$
1030 LINE INPUT #1, L$
1040 GOSUB 4000 : REM find each subject's position in the header
1050 FOR J = 1 TO K
1060 FOR I = 1 TO NF
1070 IF FL$(I) = SUBJ$(J) THEN COL(J) = I
1080 NEXT I
1090 NEXT J
1100 N = 0
1110 IF EOF(1) THEN 1190
1120 LINE INPUT #1, L$
1130 GOSUB 4000
1140 N = N + 1
1150 FOR J = 1 TO K
1160 Z(N, J) = VAL(FL$(COL(J)))
1170 NEXT J
1180 GOTO 1110
1190 CLOSE #1
1200 RETURN
2000 REM Standardize each column to mean 0 and (sample) standard deviation
2010 REM 1, then fill C, the correlation matrix: their covariance.
2020 FOR J = 1 TO K
2030 T = 0
2040 FOR I = 1 TO N
2050 T = T + Z(I, J)
2060 NEXT I
2070 MU = T / N : T = 0
2080 FOR I = 1 TO N
2090 T = T + (Z(I, J) - MU) ^ 2
2100 NEXT I
2110 SD = SQR(T / (N - 1))
2120 FOR I = 1 TO N
2130 Z(I, J) = (Z(I, J) - MU) / SD
2140 NEXT I
2150 NEXT J
2160 FOR A = 1 TO K
2170 FOR B = 1 TO K
2180 T = 0
2190 FOR I = 1 TO N
2200 T = T + Z(I, A) * Z(I, B)
2210 NEXT I
2220 C(A, B) = T / (N - 1)
2230 NEXT B
2240 NEXT A
2250 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 Power iteration: multiply V by M and rescale it, STEPS times.
5010 REM Stores V in PC(CI, ...) and its eigenvalue in EV(CI).
5020 V(1) = 1 : V(2) = 0 : V(3) = 0 : V(4) = 0
5030 FOR S = 1 TO STEPS
5040 GOSUB 7000 : REM MV = M times V
5050 T = 0
5060 FOR A = 1 TO K
5070 T = T + MV(A) ^ 2
5080 NEXT A
5090 T = SQR(T)
5100 FOR A = 1 TO K
5110 V(A) = MV(A) / T
5120 NEXT A
5130 NEXT S
5140 REM An eigenvector can point either way; make its first entry positive.
5150 IF V(1) >= 0 THEN 5190
5160 FOR A = 1 TO K
5170 V(A) = -V(A)
5180 NEXT A
5190 GOSUB 7000
5200 EV(CI) = 0 : REM the Rayleigh quotient V . (M V)
5210 FOR A = 1 TO K
5220 EV(CI) = EV(CI) + V(A) * MV(A) : PC(CI, A) = V(A)
5230 NEXT A
5240 RETURN
6000 REM Print component CI: its eigenvalue, share, and loadings.
6010 REM STR$ puts a space before the digits, so MID$ starts at the 2nd.
6020 D$ = MID$(STR$(CI), 2)
6030 L$ = "eigenvalue_" + D$ : Y = EV(CI) : GOSUB 3000
6040 REM Standardized columns each have variance 1, so the total is K.
6050 L$ = "explained_" + D$ : Y = EV(CI) / K : GOSUB 3000
6060 FOR J = 1 TO K
6070 L$ = "pc" + D$ + "_" + SUBJ$(J) : Y = PC(CI, J) : GOSUB 3000
6080 NEXT J
6090 RETURN
7000 REM MV = M times V.
7010 FOR A = 1 TO K
7020 MV(A) = 0
7030 FOR B = 1 TO K
7040 MV(A) = MV(A) + M(A, B) * V(B)
7050 NEXT B
7060 NEXT A
7070 RETURNThe 4 × 4 correlation matrix is a 2-D array, C(4, 4), and the two
components fill the rows of PC(2, 4). Power iteration is the subroutine
at line 5000. It can’t take arguments, so the caller copies the matrix into
M and sets CI to 1 or 2 to say which row to fill. Multiplying M by the
vector is another subroutine, at line 7000. There are no matrix or
eigenvector functions to avoid; even a vector’s length is a loop and SQR.
STR$(1) gives " 1", with a leading space, so the labels use MID$ to
drop it. The whole run takes under a second.
/* Principal component analysis by power iteration, using only the C
* standard library.
*
* Standardizes the four scores, builds their correlation matrix, and finds
* the top two eigenvectors (the principal components): power iteration finds
* the first, then deflation removes it from the matrix so the same steps
* find the second. */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAX_LINE 1024
#define MAX_FIELDS 16
#define K 4 /* subjects */
#define STEPS 1000
static const char *SUBJECTS[K] = {"math", "physics", "reading", "writing"};
/* Cut a line at its commas into fields[], returning how many there are.
* strtok writes '\0' over each comma, so the fields point into line. */
static int split(char *line, char *fields[]) {
int count = 0;
for (char *f = strtok(line, ",\r\n"); f && count < MAX_FIELDS; f = strtok(NULL, ",\r\n"))
fields[count++] = f;
return count;
}
/* Read the four subject columns, found by name, into rows of K scores.
* double (*)[K] is a pointer to rows of K doubles, so scores[i][j] works
* on one block of memory that realloc can grow. */
static double (*read_scores(const char *path, size_t *count))[K] {
FILE *file = fopen(path, "r");
if (!file) { perror(path); exit(1); }
char line[MAX_LINE];
char *fields[MAX_FIELDS];
if (!fgets(line, sizeof line, file)) { fprintf(stderr, "empty file\n"); exit(1); }
int fields_in_header = split(line, fields), cols[K];
for (int j = 0; j < K; j++) {
cols[j] = -1;
for (int i = 0; i < fields_in_header; i++)
if (strcmp(fields[i], SUBJECTS[j]) == 0) cols[j] = i;
if (cols[j] < 0) { fprintf(stderr, "no column %s\n", SUBJECTS[j]); exit(1); }
}
size_t n = 0, capacity = 64;
double (*scores)[K] = malloc(capacity * sizeof *scores);
while (fgets(line, sizeof line, file)) {
if (split(line, fields) < fields_in_header) continue;
if (n == capacity) {
capacity *= 2;
scores = realloc(scores, capacity * sizeof *scores);
}
for (int j = 0; j < K; j++) scores[n][j] = strtod(fields[cols[j]], NULL);
n++;
}
fclose(file);
*count = n;
return scores;
}
/* out = m v. An array parameter like m[K][K] is really a pointer, so the
* result goes into an array the caller provides. (m can't be marked const:
* before C23, passing a plain matrix to a const one is an error.) */
static void mat_vec(double m[K][K], const double v[K], double out[K]) {
for (int i = 0; i < K; i++) {
out[i] = 0;
for (int j = 0; j < K; j++) out[i] += m[i][j] * v[j];
}
}
static double dot(const double u[K], const double v[K]) {
double total = 0;
for (int i = 0; i < K; i++) total += u[i] * v[i];
return total;
}
/* Power iteration: multiply by the matrix and rescale, STEPS times. Fills
* v with the eigenvector and returns its eigenvalue. */
static double top_eigenvector(double m[K][K], double v[K]) {
double mv[K];
for (int i = 0; i < K; i++) v[i] = i == 0 ? 1 : 0;
for (int step = 0; step < STEPS; step++) {
mat_vec(m, v, mv);
double length = sqrt(dot(mv, mv));
for (int i = 0; i < K; i++) v[i] = mv[i] / length;
}
/* An eigenvector can point either way; make its first entry positive. */
if (v[0] < 0)
for (int i = 0; i < K; i++) v[i] = -v[i];
mat_vec(m, v, mv);
return dot(v, mv); /* the Rayleigh quotient v . (m v) */
}
static void report(int c, double value, const double pc[K]) {
printf("eigenvalue_%d %.6f\n", c, value);
/* Standardized columns each have variance 1, so the total is K. */
printf("explained_%d %.6f\n", c, value / K);
for (int j = 0; j < K; j++) printf("pc%d_%s %.6f\n", c, SUBJECTS[j], pc[j]);
}
int main(int argc, char **argv) {
if (argc < 2) { fprintf(stderr, "usage: %s data.csv\n", argv[0]); return 1; }
size_t n;
double (*z)[K] = read_scores(argv[1], &n);
/* Standardize each column to mean 0 and (sample) standard deviation 1. */
for (int j = 0; j < K; j++) {
double total = 0, squares = 0;
for (size_t i = 0; i < n; i++) total += z[i][j];
double mean = total / n;
for (size_t i = 0; i < n; i++) squares += (z[i][j] - mean) * (z[i][j] - mean);
double sd = sqrt(squares / (n - 1));
for (size_t i = 0; i < n; i++) z[i][j] = (z[i][j] - mean) / sd;
}
/* The correlation matrix: the covariance of the standardized columns. */
double corr[K][K];
for (int a = 0; a < K; a++)
for (int b = 0; b < K; b++) {
double total = 0;
for (size_t i = 0; i < n; i++) total += z[i][a] * z[i][b];
corr[a][b] = total / (n - 1);
}
double pc_1[K], pc_2[K], deflated[K][K];
double value_1 = top_eigenvector(corr, pc_1);
/* Deflation: subtract the first component's share, leaving the rest. */
for (int a = 0; a < K; a++)
for (int b = 0; b < K; b++) deflated[a][b] = corr[a][b] - value_1 * pc_1[a] * pc_1[b];
double value_2 = top_eigenvector(deflated, pc_2);
printf("count %zu\n", n);
report(1, value_1, pc_1);
report(2, value_2, pc_2);
/* Arrays count from 0, so the first student is z[0]. */
printf("student_1_pc1 %.6f\n", dot(z[0], pc_1));
printf("student_1_pc2 %.6f\n", dot(z[0], pc_2));
free(z);
return 0;
}The correlation matrix is a plain double corr[4][4], since the four
subjects fix its size, and the students’ scores live in one block of
memory typed double (*)[4], a pointer to rows of four, so z[i][j]
works while realloc grows it. The CSV columns are found by name after
strtok cuts each line at its commas.
C has no matrix product, so mat_vec writes one out as two nested loops,
and because a function can’t return an array, it fills one the caller
passes in. Power iteration, deflation, and the sign rule are a few more
short loops.
// Principal component analysis by power iteration, using only the C++
// standard library.
//
// Standardizes the four scores, builds their correlation matrix, and finds
// the top two eigenvectors (the principal components): power iteration finds
// the first, then deflation removes it from the matrix so the same steps
// find the second. C++ has no matrix type and no eigenvalue routine.
#include <algorithm>
#include <array>
#include <cmath>
#include <fstream>
#include <numeric>
#include <print>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
constexpr std::size_t K = 4; // subjects
constexpr int STEPS = 1000;
const std::array<std::string, K> SUBJECTS = {"math", "physics", "reading", "writing"};
// The size of a std::array is part of its type, so the compiler knows a
// Matrix is 4 by 4. Unlike a C array, a std::array can be returned from a
// function, and assigning one copies every entry.
using Vector = std::array<double, K>;
using Matrix = std::array<Vector, K>;
// 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 four subject columns, found by name, as rows of four scores.
std::vector<Vector> read_scores(const std::string& path) {
std::ifstream file(path);
std::string line;
std::getline(file, line);
std::vector<std::string> header = split(line);
std::array<std::size_t, K> cols;
for (std::size_t j = 0; j < K; j++) {
auto found = std::find(header.begin(), header.end(), SUBJECTS[j]);
if (found == header.end()) {
throw std::runtime_error("no column " + SUBJECTS[j]);
}
cols[j] = found - header.begin();
}
std::vector<Vector> rows;
while (std::getline(file, line)) {
if (line.empty()) {
continue;
}
std::vector<std::string> fields = split(line);
Vector row;
for (std::size_t j = 0; j < K; j++) row[j] = std::stod(fields.at(cols[j]));
rows.push_back(row);
}
return rows;
}
// inner_product multiplies two ranges entry by entry and adds; the 0.0
// makes it add doubles.
double dot(const Vector& u, const Vector& v) {
return std::inner_product(u.begin(), u.end(), v.begin(), 0.0);
}
// m v: each entry is one row of m dotted with v.
Vector mat_vec(const Matrix& m, const Vector& v) {
Vector out;
for (std::size_t i = 0; i < K; i++) out[i] = dot(m[i], v);
return out;
}
// Power iteration: multiply by the matrix and rescale, STEPS times. Fills
// v with the eigenvector and returns its eigenvalue.
double top_eigenvector(const Matrix& m, Vector& v) {
v = {1, 0, 0, 0};
for (int step = 0; step < STEPS; step++) {
Vector mv = mat_vec(m, v);
double length = std::sqrt(dot(mv, mv));
for (std::size_t i = 0; i < K; i++) v[i] = mv[i] / length;
}
// An eigenvector can point either way; make its first entry positive.
if (v[0] < 0) {
for (double& x : v) x = -x;
}
return dot(v, mat_vec(m, v)); // the Rayleigh quotient v . (m v)
}
void report(int c, double value, const Vector& pc) {
std::println("eigenvalue_{} {:.6f}", c, value);
// Standardized columns each have variance 1, so the total is K.
std::println("explained_{} {:.6f}", c, value / K);
for (std::size_t j = 0; j < K; j++) {
std::println("pc{}_{} {:.6f}", c, SUBJECTS[j], pc[j]);
}
}
int main(int argc, char* argv[]) {
if (argc != 2) {
std::println(stderr, "usage: program data.csv");
return 1;
}
std::vector<Vector> z = read_scores(argv[1]);
std::size_t n = z.size();
// Standardize each column to mean 0 and (sample) standard deviation 1.
// The & makes row a reference, so the last loop changes z itself.
for (std::size_t j = 0; j < K; j++) {
double total = 0.0;
for (const Vector& row : z) total += row[j];
double mean = total / n;
double squares = 0.0;
for (const Vector& row : z) squares += (row[j] - mean) * (row[j] - mean);
double sd = std::sqrt(squares / (n - 1));
for (Vector& row : z) row[j] = (row[j] - mean) / sd;
}
// The correlation matrix: the covariance of the standardized columns.
Matrix corr;
for (std::size_t a = 0; a < K; a++) {
for (std::size_t b = 0; b < K; b++) {
double total = 0.0;
for (const Vector& row : z) total += row[a] * row[b];
corr[a][b] = total / (n - 1);
}
}
Vector pc_1, pc_2;
double value_1 = top_eigenvector(corr, pc_1);
// Deflation: subtract the first component's share, leaving the rest.
// Assigning a std::array copies all 16 entries, so changing deflated
// leaves corr as it was.
Matrix deflated = corr;
for (std::size_t a = 0; a < K; a++) {
for (std::size_t b = 0; b < K; b++) deflated[a][b] -= value_1 * pc_1[a] * pc_1[b];
}
double value_2 = top_eigenvector(deflated, pc_2);
std::println("count {}", n);
report(1, value_1, pc_1);
report(2, value_2, pc_2);
// Vectors count from 0, so the first student is z[0].
std::println("student_1_pc1 {:.6f}", dot(z[0], pc_1));
std::println("student_1_pc2 {:.6f}", dot(z[0], pc_2));
}The correlation matrix is a std::array of four std::arrays of four
doubles, so its size is part of its type, as in C. Unlike a C array, though,
a std::array can be returned from a function, so mat_vec simply returns
the product, one std::inner_product per row.
Watch the copy. Matrix deflated = corr copies all 16 entries, so
subtracting the first component from deflated leaves corr untouched.
That’s the opposite of Java, where assignment shares the array, and in C,
assigning one array to another won’t compile. The standard library has no
eigenvalue routine to skip, and libraries such as Eigen aren’t installed.
*> Principal component analysis by power iteration, in COBOL.
*>
*> Standardizes the four scores, builds their correlation matrix,
*> and finds the top two eigenvectors (the principal components):
*> power iteration finds the first, then deflation removes it from
*> the matrix so the same steps find the second. Values are
*> fixed-point decimals carried to 18 decimal places.
IDENTIFICATION DIVISION.
PROGRAM-ID. PCA.
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.
*> The four subjects, and which CSV column holds each one.
01 WS-SUBJECT-NAMES.
05 FILLER PIC X(8) VALUE "math".
05 FILLER PIC X(8) VALUE "physics".
05 FILLER PIC X(8) VALUE "reading".
05 FILLER PIC X(8) VALUE "writing".
01 WS-SUBJECTS REDEFINES WS-SUBJECT-NAMES.
05 WS-SUBJECT PIC X(8) OCCURS 4 TIMES.
01 WS-COLUMNS.
05 WS-COLUMN PIC 99 VALUE 0 OCCURS 4 TIMES.
01 K PIC 9 VALUE 4.
*> Every student's four scores, standardized in place.
01 WS-N PIC 9(4) VALUE 0.
01 WS-ROWS.
05 WS-ROW OCCURS 5000 TIMES.
10 WS-Z PIC S9(6)V9(18) OCCURS 4 TIMES.
*> Each column's mean and sample standard deviation.
01 WS-STATS.
05 WS-MEAN PIC S9(6)V9(18) OCCURS 4 TIMES.
05 WS-SD PIC S9(6)V9(18) OCCURS 4 TIMES.
01 WS-SUM PIC S9(9)V9(18).
*> A 4 x 4 matrix is a table of rows, each a table of entries:
*> WS-CORR(I, J). WS-M is the matrix power iteration works on.
01 WS-CORR-TABLE.
05 WS-CORR-ROW OCCURS 4 TIMES.
10 WS-CORR PIC S9(3)V9(18) OCCURS 4 TIMES.
01 WS-M-TABLE.
05 WS-M-ROW OCCURS 4 TIMES.
10 WS-M PIC S9(3)V9(18) OCCURS 4 TIMES.
*> Power iteration's vector V, the product M V, and its length.
01 WS-V-TABLE.
05 WS-V PIC S9(3)V9(18) OCCURS 4 TIMES.
01 WS-MV-TABLE.
05 WS-MV PIC S9(3)V9(18) OCCURS 4 TIMES.
01 WS-LENGTH PIC S9(3)V9(18).
01 WS-EIGENVALUE PIC S9(3)V9(18).
01 WS-STEPS PIC 9(4) VALUE 1000.
*> The two components found: eigenvalues and eigenvectors.
01 WS-PC-TABLE.
05 WS-PC OCCURS 2 TIMES.
10 WS-PC-VALUE PIC S9(3)V9(18).
10 WS-PC-VECTOR PIC S9(3)V9(18) OCCURS 4 TIMES.
01 I PIC 9(4).
01 J PIC 9(4).
01 C PIC 9.
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
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
IF WS-FIELD(I) = WS-SUBJECT(J)
MOVE I TO WS-COLUMN(J)
END-IF
END-PERFORM
END-PERFORM
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
IF WS-COLUMN(J) = 0
DISPLAY "missing column" UPON SYSERR
MOVE 1 TO RETURN-CODE
STOP RUN
END-IF
END-PERFORM
PERFORM UNTIL END-OF-FILE
READ DATA-FILE
AT END SET END-OF-FILE TO TRUE
NOT AT END PERFORM ADD-ROW
END-READ
END-PERFORM
CLOSE DATA-FILE
PERFORM STANDARDIZE
*> The correlation matrix: the covariance of the standardized
*> columns.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
MOVE 0 TO WS-SUM
PERFORM VARYING S FROM 1 BY 1 UNTIL S > WS-N
COMPUTE WS-SUM ROUNDED =
WS-SUM + WS-Z(S, I) * WS-Z(S, J)
END-PERFORM
COMPUTE WS-CORR(I, J) ROUNDED = WS-SUM / (WS-N - 1)
END-PERFORM
END-PERFORM
MOVE WS-CORR-TABLE TO WS-M-TABLE
MOVE 1 TO C
PERFORM TOP-EIGENVECTOR
*> Deflation: subtract the first component's share, leaving the
*> rest.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
COMPUTE WS-M(I, J) ROUNDED = WS-CORR(I, J)
- WS-PC-VALUE(1) * WS-PC-VECTOR(1, I)
* WS-PC-VECTOR(1, J)
END-PERFORM
END-PERFORM
MOVE 2 TO C
PERFORM TOP-EIGENVECTOR
MOVE WS-N TO WS-COUNT-OUT
DISPLAY "count " FUNCTION TRIM(WS-COUNT-OUT)
PERFORM VARYING C FROM 1 BY 1 UNTIL C > 2
MOVE SPACES TO WS-NAME
STRING "eigenvalue_" C DELIMITED BY SIZE INTO WS-NAME
MOVE WS-PC-VALUE(C) TO WS-RESULT
PERFORM PRINT-RESULT
*> Standardized columns each have variance 1, so the total
*> is K.
MOVE SPACES TO WS-NAME
STRING "explained_" C DELIMITED BY SIZE INTO WS-NAME
COMPUTE WS-RESULT ROUNDED = WS-PC-VALUE(C) / K
PERFORM PRINT-RESULT
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
MOVE SPACES TO WS-NAME
STRING "pc" C "_" WS-SUBJECT(J)
DELIMITED BY SPACE INTO WS-NAME
MOVE WS-PC-VECTOR(C, J) TO WS-RESULT
PERFORM PRINT-RESULT
END-PERFORM
END-PERFORM
*> The first student's scores along each component.
PERFORM VARYING C FROM 1 BY 1 UNTIL C > 2
MOVE 0 TO WS-SUM
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
COMPUTE WS-SUM ROUNDED =
WS-SUM + WS-Z(1, J) * WS-PC-VECTOR(C, J)
END-PERFORM
MOVE SPACES TO WS-NAME
STRING "student_1_pc" C DELIMITED BY SIZE INTO WS-NAME
MOVE WS-SUM 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 student's four scores.
ADD-ROW.
PERFORM SPLIT-LINE
ADD 1 TO WS-N
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
COMPUTE WS-Z(WS-N, J) =
FUNCTION NUMVAL(WS-FIELD(WS-COLUMN(J)))
END-PERFORM.
*> Rescale each column to mean 0 and sample standard deviation 1
*> (dividing by n - 1).
STANDARDIZE.
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
MOVE 0 TO WS-SUM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
ADD WS-Z(I, J) TO WS-SUM
END-PERFORM
COMPUTE WS-MEAN(J) ROUNDED = WS-SUM / WS-N
MOVE 0 TO WS-SUM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
COMPUTE WS-SUM =
WS-SUM + (WS-Z(I, J) - WS-MEAN(J)) ** 2
END-PERFORM
COMPUTE WS-SD(J) ROUNDED =
FUNCTION SQRT(WS-SUM / (WS-N - 1))
PERFORM VARYING I FROM 1 BY 1 UNTIL I > WS-N
COMPUTE WS-Z(I, J) ROUNDED =
(WS-Z(I, J) - WS-MEAN(J)) / WS-SD(J)
END-PERFORM
END-PERFORM.
*> Power iteration on WS-M: multiply by the matrix and rescale,
*> WS-STEPS times. Saves the eigenvalue and eigenvector as
*> component C.
TOP-EIGENVECTOR.
MOVE 0 TO WS-V(2) WS-V(3) WS-V(4)
MOVE 1 TO WS-V(1)
PERFORM WS-STEPS TIMES
PERFORM MULTIPLY-M-V
MOVE 0 TO WS-SUM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
COMPUTE WS-SUM ROUNDED = WS-SUM + WS-MV(I) ** 2
END-PERFORM
COMPUTE WS-LENGTH ROUNDED = FUNCTION SQRT(WS-SUM)
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
COMPUTE WS-V(I) ROUNDED = WS-MV(I) / WS-LENGTH
END-PERFORM
END-PERFORM
*> An eigenvector can point either way; make its first entry
*> positive.
IF WS-V(1) < 0
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
COMPUTE WS-V(I) = - WS-V(I)
END-PERFORM
END-IF
*> The eigenvalue is the Rayleigh quotient V . (M V).
PERFORM MULTIPLY-M-V
MOVE 0 TO WS-SUM
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
COMPUTE WS-SUM ROUNDED = WS-SUM + WS-V(I) * WS-MV(I)
END-PERFORM
MOVE WS-SUM TO WS-PC-VALUE(C)
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
MOVE WS-V(I) TO WS-PC-VECTOR(C, I)
END-PERFORM.
*> WS-MV = WS-M times WS-V.
MULTIPLY-M-V.
PERFORM VARYING I FROM 1 BY 1 UNTIL I > K
MOVE 0 TO WS-MV(I)
PERFORM VARYING J FROM 1 BY 1 UNTIL J > K
COMPUTE WS-MV(I) ROUNDED =
WS-MV(I) + WS-M(I, J) * WS-V(J)
END-PERFORM
END-PERFORM.
*> 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).The 4 × 4 correlation matrix is a table of rows, each a table of entries,
so WS-CORR(I, J) reads like math, and a single MOVE copies the whole
matrix into the working copy that power iteration multiplies. Values are
fixed-point decimals with 18 places, and all 14 decimals we printed match
Python’s doubles. Rounding doesn’t build up over the 1000 normalizations,
because each step pulls the vector back toward the eigenvector: with only
9 places, the results still agree to about 1e-9. The program runs in about
13 milliseconds. GnuCOBOL has MEAN, VARIANCE, and STANDARD-DEVIATION,
which we skip, but nothing for matrices or eigenvectors. At 288 lines,
it’s four times the 72-line Python.
;;;; Principal component analysis by power iteration, using only standard
;;;; Common Lisp.
;;;;
;;;; Standardizes the four scores, builds their correlation matrix, and finds
;;;; the top two eigenvectors (the principal components): power iteration
;;;; finds the first, then deflation removes it from the matrix so the same
;;;; steps find the second. A vector is a list of numbers and a matrix is a
;;;; list of rows, so mapcar does most of the work. Common Lisp also has 2-D
;;;; arrays, read with (aref m i j), but lists need no index bookkeeping.
(defparameter *subjects* '("math" "physics" "reading" "writing"))
(defparameter *steps* 1000)
;; The scores are whole numbers, which parse-number turns into doubles. We
;; still read every decimal as a double, in the data and in this file, as
;; every program here does.
;; 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-rows (path names)
"The columns called NAMES in the CSV file at PATH, as a list of rows of
numbers, in the order NAMES gives."
(with-open-file (in path)
(let* ((header (split-on-comma (read-line in)))
(indexes (mapcar (lambda (name)
(or (position name header :test #'string=)
(error "no column ~a" name)))
names)))
(loop for line = (read-line in nil)
while line
unless (string= line "")
collect (let ((fields (split-on-comma line)))
(mapcar (lambda (i) (parse-number (nth i fields))) indexes))))))
(defun transpose (rows)
"The columns of ROWS. (mapcar #'list row1 row2 ...) collects the first
entry of every row, then the second, and so on; apply passes the rows as
those arguments."
(apply #'mapcar #'list rows))
(defun dot (u v)
"The sum of the products of U and V, entry by entry."
(reduce #'+ (mapcar #'* u v)))
(defun mat-vec (m v)
"A matrix times a vector: each row of M dotted with V."
(mapcar (lambda (row) (dot row v)) m))
(defun normalize (v)
(let ((length (sqrt (dot v v))))
(mapcar (lambda (x) (/ x length)) v)))
(defun standardize (column)
"COLUMN shifted to mean 0 and scaled to (sample) standard deviation 1."
(let* ((n (length column))
(mean (/ (reduce #'+ column) n))
(sd (sqrt (/ (reduce #'+ (mapcar (lambda (x) (expt (- x mean) 2)) column))
(1- n)))))
(mapcar (lambda (x) (/ (- x mean) sd)) column)))
(defun top-eigenvector (m)
"Power iteration: multiply by M and rescale, *steps* times. Returns two
values, the eigenvalue and the eigenvector."
(let ((v (cons 1.0 (make-list (1- (length m)) :initial-element 0.0))))
(loop repeat *steps*
do (setf v (normalize (mat-vec m v))))
;; An eigenvector can point either way; make its first entry positive.
(when (minusp (first v))
(setf v (mapcar #'- v)))
(values (dot v (mat-vec m v)) v)))
(let* ((rows (read-rows (second sb-ext:*posix-argv*) *subjects*))
(n (length rows))
(k (length *subjects*))
(columns (mapcar #'standardize (transpose rows)))
(z (transpose columns))
;; The correlation matrix: the covariance of the standardized columns.
(corr (mapcar (lambda (a)
(mapcar (lambda (b) (/ (dot a b) (1- n))) columns))
columns)))
;; multiple-value-bind catches both values top-eigenvector returns.
(multiple-value-bind (value-1 pc-1) (top-eigenvector corr)
;; Deflation: subtract the first component's share, leaving the rest.
(let ((deflated (mapcar (lambda (row a)
(mapcar (lambda (c b) (- c (* value-1 a b))) row pc-1))
corr pc-1)))
(multiple-value-bind (value-2 pc-2) (top-eigenvector deflated)
;; ~d prints an integer, ~a a string, ~,6f a number with six
;; decimals, and ~% a newline.
(format t "count ~d~%" n)
(loop for c from 1
for (value pc) in (list (list value-1 pc-1) (list value-2 pc-2))
do (format t "eigenvalue_~d ~,6f~%" c value)
;; Standardized columns each have variance 1, so the total is k.
(format t "explained_~d ~,6f~%" c (/ value k))
(loop for subject in *subjects*
for loading in pc
do (format t "pc~d_~a ~,6f~%" c subject loading)))
(format t "student_1_pc1 ~,6f~%" (dot (first z) pc-1))
(format t "student_1_pc2 ~,6f~%" (dot (first z) pc-2))))))A vector is a list of numbers and a matrix is a list of rows, so mat-vec
is one mapcar of dot over the rows, and (apply #'mapcar #'list rows)
turns rows into columns. Common Lisp has 2-D arrays too, read with
(aref m i j), but lists need no index bookkeeping here. There’s no
eigenvalue routine to skip.
Each power-iteration step builds a new list, and top-eigenvector returns
two values, the eigenvalue and the vector, which multiple-value-bind
catches by name. Flipping the sign is (mapcar #'- v), since - with one
argument negates it.
# Principal component analysis by power iteration, using only Elixir's
# standard library.
#
# Standardizes the four scores, builds their correlation matrix, and finds the
# top two eigenvectors (the principal components): power iteration finds the
# first, then deflation removes it from the matrix so the same steps find the
# second. A matrix is a list of row lists, and nothing is changed in place.
defmodule PCA do
@subjects ~w(math physics reading writing) # ~w splits words into a list
@steps 1000
def subjects, do: @subjects
# The subject columns of every row, read as floats, in the order above.
def read_rows(path) do
[header | lines] = path |> File.read!() |> String.split("\n", trim: true)
names = String.split(header, ",")
indexes = Enum.map(@subjects, fn s -> Enum.find_index(names, &(&1 == s)) || raise "column not found: #{s}" end)
Enum.map(lines, fn line ->
fields = String.split(line, ",")
# Float.parse returns {number, rest}; "82" parses as 82.0 too.
Enum.map(indexes, fn i -> {value, ""} = fields |> Enum.at(i) |> Float.parse(); value end)
end)
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)
# A matrix times a vector: each row dotted with the vector.
def mat_vec(m, v), do: Enum.map(m, &dot(&1, v))
def normalize(v) do
length = :math.sqrt(dot(v, v))
Enum.map(v, &(&1 / length))
end
# The columns of a list of rows: zip_with gathers the first of every row,
# then the second, and so on.
def columns(rows), do: Enum.zip_with(rows, & &1)
# Power iteration: multiply by the matrix and rescale, 1,000 times.
# Returns {eigenvalue, eigenvector}.
def top_eigenvector(m) do
start = [1.0 | List.duplicate(0.0, length(m) - 1)]
# reduce threads v through every step; _ ignores the step number.
v = Enum.reduce(1..@steps, start, fn _, v -> normalize(mat_vec(m, v)) end)
# An eigenvector can point either way; make its first entry positive.
v = if hd(v) < 0, do: Enum.map(v, &(-&1)), else: v
{dot(v, mat_vec(m, v)), v}
end
def say(name, value), do: IO.puts("#{name} #{:erlang.float_to_binary(value, decimals: 6)}")
end
[path] = System.argv()
rows = PCA.read_rows(path)
n = length(rows)
k = length(PCA.subjects())
# Standardize each column to mean 0 and (sample) standard deviation 1.
cols = PCA.columns(rows)
means = Enum.map(cols, fn c -> Enum.sum(c) / n end)
sds =
Enum.zip_with(cols, means, fn c, m ->
:math.sqrt((c |> Enum.map(&((&1 - m) ** 2)) |> Enum.sum()) / (n - 1))
end)
z = Enum.map(rows, fn row -> Enum.zip_with([row, means, sds], fn [x, m, s] -> (x - m) / s end) end)
# The correlation matrix: the covariance of the standardized columns.
zcols = PCA.columns(z)
corr = for a <- zcols, do: for(b <- zcols, do: PCA.dot(a, b) / (n - 1))
{value_1, pc_1} = PCA.top_eigenvector(corr)
# Deflation: subtract the first component's share, leaving the rest.
deflated =
Enum.zip_with(corr, pc_1, fn row, a ->
Enum.zip_with(row, pc_1, fn c, b -> c - value_1 * a * b end)
end)
{value_2, pc_2} = PCA.top_eigenvector(deflated)
IO.puts("count #{n}")
for {{value, pc}, c} <- Enum.with_index([{value_1, pc_1}, {value_2, pc_2}], 1) do
PCA.say("eigenvalue_#{c}", value)
# Standardized columns each have variance 1, so the total is k.
PCA.say("explained_#{c}", value / k)
for {subject, loading} <- Enum.zip(PCA.subjects(), pc), do: PCA.say("pc#{c}_#{subject}", loading)
end
PCA.say("student_1_pc1", PCA.dot(hd(z), pc_1))
PCA.say("student_1_pc2", PCA.dot(hd(z), pc_2))A matrix is a list of row lists. Elixir won’t index a list with
m[i][j]; it raises an ArgumentError that points you to Enum, so
mat_vec maps each row to its dot product with the vector, and columns
turns rows into columns with Enum.zip_with(rows, & &1). Power iteration
is Enum.reduce over 1..1000, which passes each new, normalized vector
into the next step, with nothing updated in place.
Neither Elixir nor Erlang ships a matrix type or an eigenvalue routine. Nx, the usual numerical library, is a Hex package, and it isn’t installed, so the correlation matrix, deflation, and every product are written out.
! Principal component analysis by power iteration, using only standard Fortran.
!
! Standardizes the four scores, builds their correlation matrix, and finds the
! top two eigenvectors (the principal components): power iteration finds the
! first, then deflation removes it from the matrix so the same steps find the
! second.
program pca
implicit none
integer, parameter :: dp = kind(1.0d0) ! double precision
integer, parameter :: k = 4, steps = 1000
character(len=*), parameter :: subjects(k) = &
[character(len=7) :: "math", "physics", "reading", "writing"]
real(dp), allocatable :: z(:, :) ! one row per student, one column per subject
real(dp) :: corr(k, k), deflated(k, k), pc_1(k), pc_2(k), value_1, value_2
real(dp) :: mean, sd
integer :: n, j
z = read_scores()
n = size(z, 1)
! Standardize each column to mean 0 and (sample) standard deviation 1.
do j = 1, k
mean = sum(z(:, j)) / n
sd = sqrt(sum((z(:, j) - mean)**2) / (n - 1))
z(:, j) = (z(:, j) - mean) / sd
end do
! The correlation matrix: the covariance of the standardized columns.
! transpose(z) is k by n, so the product is k by k.
corr = matmul(transpose(z), z) / (n - 1)
call top_eigenvector(corr, value_1, pc_1)
! Deflation: subtract the first component's share, leaving the rest.
! spread copies pc_1 into k columns (or rows), so their product is the
! outer product pc_1(i) * pc_1(j).
deflated = corr - value_1 * spread(pc_1, 2, k) * spread(pc_1, 1, k)
call top_eigenvector(deflated, value_2, pc_2)
print "(a, 1x, i0)", "count", n
call report(1, value_1, pc_1)
call report(2, value_2, pc_2)
! Fortran arrays count from 1, so the first student is row 1.
call say("student_1_pc1", dot_product(z(1, :), pc_1))
call say("student_1_pc2", dot_product(z(1, :), pc_2))
contains
! Power iteration: multiply by the matrix and rescale, steps times.
subroutine top_eigenvector(m, value, v)
real(dp), intent(in) :: m(k, k)
real(dp), intent(out) :: value, v(k)
integer :: step
v = [1, 0, 0, 0]
do step = 1, steps
v = matmul(m, v)
v = v / norm2(v) ! norm2 is the vector's length
end do
! An eigenvector can point either way; make its first entry positive.
if (v(1) < 0) v = -v
value = dot_product(v, matmul(m, v)) ! the Rayleigh quotient
end subroutine top_eigenvector
subroutine report(c, value, pc)
integer, intent(in) :: c
real(dp), intent(in) :: value, pc(k)
character(len=20) :: name
integer :: j
write(name, "(a, i0)") "eigenvalue_", c
call say(trim(name), value)
! Standardized columns each have variance 1, so the total is k.
write(name, "(a, i0)") "explained_", c
call say(trim(name), value / k)
do j = 1, k
write(name, "(a, i0, a, a)") "pc", c, "_", trim(subjects(j))
call say(trim(name), pc(j))
end do
end subroutine report
! Read the four subject columns into an n by k matrix.
function read_scores() result(scores)
real(dp), allocatable :: scores(:, :), grown(:, :)
character(len=256) :: path, header, line, text
integer :: unit, status, cols(k), j
call get_command_argument(1, path)
open(newunit=unit, file=path, status="old", action="read")
read(unit, "(a)") header
do j = 1, k
cols(j) = field_index(header, trim(subjects(j)))
end do
allocate(scores(0, k))
do
read(unit, "(a)", iostat=status) line
if (status /= 0) exit
! Fortran stores matrices column by column, so we rebuild the matrix
! one row longer; move_alloc hands grown's storage to scores.
allocate(grown(size(scores, 1) + 1, k))
grown(:size(scores, 1), :) = scores
do j = 1, k
text = field(line, cols(j))
read(text, *) grown(size(grown, 1), j) ! an internal read parses it
end do
call move_alloc(grown, scores)
end do
close(unit)
end function read_scores
! 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 pcaFortran was built for this kind of arithmetic. The correlation matrix is one
line, matmul(transpose(z), z) / (n - 1), and each power-iteration step is
v = matmul(m, v) followed by v = v / norm2(v), where norm2, the
vector’s length, is a Fortran 2008 intrinsic. Deflation needs the outer
product of pc_1 with itself, which we build with spread: it copies the
vector into the columns, or the rows, of a 4 by 4 matrix, and multiplying
the two element by element gives every pc_1(i) * pc_1(j). Fortran has no
eigenvalue intrinsic, and LAPACK, the usual Fortran library for it, isn’t
installed, so there was nothing to skip.
// Principal component analysis by power iteration, using only the standard
// library.
//
// Standardizes the four scores, builds their correlation matrix, and finds the
// top two eigenvectors (the principal components): power iteration finds the
// first, then deflation removes it from the matrix so the same steps find the
// second.
package main
import (
"encoding/csv"
"fmt"
"math"
"os"
"slices"
"strconv"
)
var subjects = [k]string{"math", "physics", "reading", "writing"}
const (
k = 4 // number of subjects
steps = 1000
)
// Arrays have a fixed length that is part of the type, and assigning one
// copies it, so a vector or matrix never shares storage by accident.
type (
vector [k]float64
matrix [k]vector
)
func read(path string) ([]vector, error) {
f, err := os.Open(path)
if err != nil {
return nil, err
}
defer f.Close()
records, err := csv.NewReader(f).ReadAll()
if err != nil {
return nil, err
}
var index [k]int
for j, name := range subjects {
index[j] = slices.Index(records[0], name)
if index[j] < 0 {
return nil, fmt.Errorf("no column %s", name)
}
}
rows := make([]vector, 0, len(records)-1)
for _, record := range records[1:] {
var row vector
for j, i := range index {
row[j], err = strconv.ParseFloat(record[i], 64)
if err != nil {
return nil, err
}
}
rows = append(rows, row)
}
return rows, nil
}
func dot(u, v vector) float64 {
total := 0.0
for i := range k {
total += u[i] * v[i]
}
return total
}
func (m matrix) times(v vector) vector {
var out vector
for i, row := range m {
out[i] = dot(row, v)
}
return out
}
func normalize(v vector) vector {
length := math.Sqrt(dot(v, v))
for i := range v {
v[i] /= length // v is a copy, so the caller's vector is untouched
}
return v
}
// topEigenvector runs power iteration: multiply by the matrix and rescale,
// steps times. It returns the eigenvalue and the eigenvector.
func topEigenvector(m matrix) (float64, vector) {
v := vector{1}
for range steps {
v = normalize(m.times(v))
}
// An eigenvector can point either way; make its first entry positive.
if v[0] < 0 {
for i := range v {
v[i] = -v[i]
}
}
return dot(v, m.times(v)), v
}
func main() {
rows, err := read(os.Args[1])
if err != nil {
fmt.Fprintln(os.Stderr, err)
os.Exit(1)
}
n := float64(len(rows))
// Standardize each column to mean 0 and (sample) standard deviation 1.
var means, sds vector
for _, row := range rows {
for j := range k {
means[j] += row[j]
}
}
for j := range k {
means[j] /= n
}
for _, row := range rows {
for j := range k {
sds[j] += (row[j] - means[j]) * (row[j] - means[j])
}
}
for j := range k {
sds[j] = math.Sqrt(sds[j] / (n - 1))
}
z := make([]vector, len(rows))
for r, row := range rows {
for j := range k {
z[r][j] = (row[j] - means[j]) / sds[j]
}
}
// The correlation matrix: the covariance of the standardized columns.
var corr matrix
for _, r := range z {
for i := range k {
for j := range k {
corr[i][j] += r[i] * r[j] / (n - 1)
}
}
}
value1, pc1 := topEigenvector(corr)
// Deflation: subtract the first component's share, leaving the rest.
deflated := corr // a full copy, since matrix is an array
for i := range k {
for j := range k {
deflated[i][j] -= value1 * pc1[i] * pc1[j]
}
}
value2, pc2 := topEigenvector(deflated)
fmt.Printf("count %d\n", len(rows))
for c, pc := range []struct {
value float64
v vector
}{{value1, pc1}, {value2, pc2}} {
fmt.Printf("eigenvalue_%d %.6f\n", c+1, pc.value)
// Standardized columns each have variance 1, so the total is k.
fmt.Printf("explained_%d %.6f\n", c+1, pc.value/k)
for j, subject := range subjects {
fmt.Printf("pc%d_%s %.6f\n", c+1, subject, pc.v[j])
}
}
fmt.Printf("student_1_pc1 %.6f\n", dot(z[0], pc1))
fmt.Printf("student_1_pc2 %.6f\n", dot(z[0], pc2))
}The vector and matrix are fixed-size arrays, type vector [k]float64 and
type matrix [k]vector. Assigning an array copies every element, so
deflated := corr is a fresh matrix we can change without touching corr.
A slice would share storage: after t := s; t[0][0] = 99, s changes too.
Arrays also fill unlisted entries with zero, so the starting vector
vector{1} is [1 0 0 0], and normalize can scale its argument in place
because it holds a copy. Go’s standard library has no linear algebra
package and no eigenvalue routine, so there was no built-in to skip.
-- Principal component analysis by power iteration, using only Haskell's base
-- library.
--
-- Standardizes the four scores, builds their correlation matrix, and finds
-- the top two eigenvectors (the principal components): power iteration finds
-- the first, then deflation removes it from the matrix so the same steps find
-- the second. A matrix is a list of rows, and a vector is a list of numbers.
import Control.Monad (forM_)
import Data.List (elemIndex, transpose)
import System.Environment (getArgs)
import Text.Printf (printf)
subjects :: [String]
subjects = ["math", "physics", "reading", "writing"]
steps :: Int
steps = 1000
-- Split a line at every comma, as in the descriptive statistics example.
splitOnComma :: String -> [String]
splitOnComma line = case break (== ',') line of
(field, ',' : rest) -> field : splitOnComma rest
(field, _) -> [field]
-- Where a column name sits in the header, counting from 0.
position :: [String] -> String -> Int
position header name = case elemIndex name header of
Just index -> index
Nothing -> error ("no column " ++ name)
-- The named columns of every row, read as Doubles, in the order asked for.
readRows :: [String] -> String -> [[Double]]
readRows names text = case lines text of
header : rows ->
[[read (fields !! i) | i <- indexes] | row <- rows, not (null row), let fields = splitOnComma row]
where
indexes = map (position (splitOnComma header)) names
[] -> error "empty file"
dot :: [Double] -> [Double] -> Double
dot u v = sum (zipWith (*) u v)
-- A matrix times a vector: each row of the matrix dotted with the vector.
matVec :: [[Double]] -> [Double] -> [Double]
matVec m v = map (`dot` v) m
normalize :: [Double] -> [Double]
normalize v = map (/ sqrt (dot v v)) v
-- Power iteration: multiply by the matrix and rescale, 1,000 times.
-- iterate builds the endless list v0, step v0, step (step v0), ..., and
-- !! picks the entry after 1,000 steps. Returns (eigenvalue, eigenvector).
topEigenvector :: [[Double]] -> (Double, [Double])
topEigenvector m = (dot v (matVec m v), v)
where
start = 1 : replicate (length m - 1) 0
v' = iterate (normalize . matVec m) start !! steps
-- An eigenvector can point either way; make its first entry positive.
v = if head v' < 0 then map negate v' else v'
main :: IO ()
main = do
[path] <- getArgs
text <- readFile path
let rows = readRows subjects text
n = length rows
k = length subjects
count = fromIntegral n :: Double
-- Standardize each column to mean 0 and (sample) standard deviation 1.
standardize column = map (\x -> (x - m) / s) column
where
m = sum column / count
s = sqrt (sum [(x - m) ^ (2 :: Int) | x <- column] / (count - 1))
columns = map standardize (transpose rows)
z = transpose columns
-- The correlation matrix: the covariance of the standardized columns.
corr = [[dot a b / (count - 1) | b <- columns] | a <- columns]
(value1, pc1) = topEigenvector corr
-- Deflation: subtract the first component's share, leaving the rest.
deflated = [[c - value1 * a * b | (c, b) <- zip row pc1] | (row, a) <- zip corr pc1]
(value2, pc2) = topEigenvector deflated
printf "count %d\n" n
forM_ (zip [1 :: Int ..] [(value1, pc1), (value2, pc2)]) $ \(c, (value, pc)) -> do
printf "eigenvalue_%d %.6f\n" c value
-- Standardized columns each have variance 1, so the total is k.
printf "explained_%d %.6f\n" c (value / fromIntegral k)
forM_ (zip subjects pc) $ \(subject, loading) ->
printf "pc%d_%s %.6f\n" c subject loading
printf "student_1_pc1 %.6f\n" (dot (head z) pc1)
printf "student_1_pc2 %.6f\n" (dot (head z) pc2)Here a matrix is a list of rows and a vector is a list of numbers, so a
matrix-vector product is map (`dot` v) m: each row dotted with the vector.
One power-iteration step is normalize . matVec m, the two functions
joined with ., and iterate repeats it, building the endless list of
every step’s vector. !! 1000 picks the one we want, and laziness means
no step past it is ever computed.
transpose turns the rows of scores into columns for standardizing and
back again, and the correlation matrix is a list comprehension over pairs
of columns. Deflation is another comprehension. Base has no matrix type
and no eigenvalue routine, so every step is written out here.
// Principal component analysis by power iteration, using only the Java
// standard library.
//
// Standardizes the four scores, builds their correlation matrix, and finds the
// top two eigenvectors (the principal components): power iteration finds the
// first, then deflation removes it from the matrix so the same steps find the
// second. The JDK has no matrix type or eigenvalue routine, so a matrix is a
// double[][] and the products are small helper methods.
import java.io.IOException;
import java.nio.file.Files;
import java.nio.file.Path;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
class Main {
static final String[] SUBJECTS = {"math", "physics", "reading", "writing"};
static final int STEPS = 1000;
// The named column of a CSV file, as an array of doubles. Reading a file
// can fail with an IOException, and Java makes us say so with "throws".
static double[] readColumn(Path path, String name) throws IOException {
List<String> lines = Files.readAllLines(path);
int index = Arrays.asList(lines.get(0).split(",")).indexOf(name);
if (index < 0) {
throw new IllegalArgumentException("no column " + name);
}
double[] values = new double[lines.size() - 1];
for (int row = 1; row < lines.size(); row++) {
values[row - 1] = Double.parseDouble(lines.get(row).split(",")[index]);
}
return values;
}
static double dot(double[] u, double[] v) {
double total = 0;
for (int i = 0; i < u.length; i++) {
total += u[i] * v[i];
}
return total;
}
// A double[][] is an array of row arrays, so m[i] is row i, ready for dot.
static double[] matVec(double[][] m, double[] v) {
double[] out = new double[m.length];
for (int i = 0; i < m.length; i++) {
out[i] = dot(m[i], v);
}
return out;
}
static double[] normalize(double[] v) {
double length = Math.sqrt(dot(v, v));
double[] out = new double[v.length];
for (int i = 0; i < v.length; i++) {
out[i] = v[i] / length;
}
return out;
}
// Power iteration: multiply by the matrix and rescale, STEPS times. The
// eigenvalue is v . (M v), which the caller computes from the vector.
static double[] topEigenvector(double[][] m) {
double[] v = new double[m.length]; // Java fills new arrays with zeros
v[0] = 1;
for (int step = 0; step < STEPS; step++) {
v = normalize(matVec(m, v));
}
// An eigenvector can point either way; make its first entry positive.
if (v[0] < 0) {
for (int i = 0; i < v.length; i++) {
v[i] = -v[i];
}
}
return v;
}
public static void main(String[] args) throws IOException {
int k = SUBJECTS.length;
// columns[j] holds every student's score in subject j.
double[][] columns = new double[k][];
for (int j = 0; j < k; j++) {
columns[j] = readColumn(Path.of(args[0]), SUBJECTS[j]);
}
int n = columns[0].length;
// Standardize each column to mean 0 and (sample) standard deviation
// 1. z[r][j] is student r's standardized score in subject j.
double[][] z = new double[n][k];
for (int j = 0; j < k; j++) {
double sum = 0;
for (double x : columns[j]) sum += x;
double mean = sum / n;
double squares = 0;
for (double x : columns[j]) squares += (x - mean) * (x - mean);
double sd = Math.sqrt(squares / (n - 1));
for (int r = 0; r < n; r++) {
z[r][j] = (columns[j][r] - mean) / sd;
}
}
// The correlation matrix: the covariance of the standardized columns.
double[][] corr = new double[k][k];
for (double[] row : z) {
for (int i = 0; i < k; i++) {
for (int j = 0; j < k; j++) {
corr[i][j] += row[i] * row[j] / (n - 1);
}
}
}
double[] pc1 = topEigenvector(corr);
double value1 = dot(pc1, matVec(corr, pc1));
// Deflation: subtract the first component's share, leaving the rest.
// We build a new matrix; "deflated = corr" would only copy the
// reference, and the loop below would change corr too.
double[][] deflated = new double[k][k];
for (int i = 0; i < k; i++) {
for (int j = 0; j < k; j++) {
deflated[i][j] = corr[i][j] - value1 * pc1[i] * pc1[j];
}
}
double[] pc2 = topEigenvector(deflated);
double value2 = dot(pc2, matVec(deflated, pc2));
// 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, "count %d%n", n);
double[] values = {value1, value2};
double[][] pcs = {pc1, pc2};
for (int c = 0; c < 2; c++) {
System.out.printf(root, "eigenvalue_%d %.6f%n", c + 1, values[c]);
// Standardized columns each have variance 1, so the total is k.
System.out.printf(root, "explained_%d %.6f%n", c + 1, values[c] / k);
for (int j = 0; j < k; j++) {
System.out.printf(root, "pc%d_%s %.6f%n", c + 1, SUBJECTS[j], pcs[c][j]);
}
}
System.out.printf(root, "student_1_pc1 %.6f%n", dot(z[0], pc1));
System.out.printf(root, "student_1_pc2 %.6f%n", dot(z[0], pc2));
}
}Java has no matrix type, so the correlation matrix is a double[][], an
array of row arrays. That makes m[i] a whole row, and matVec is one
call to dot per row. Power iteration, deflation, and the sign rule are
short loops around it.
Watch the copy. Writing deflated = corr copies only the reference, so
changing deflated would change corr too, and even corr.clone() shares
the rows. We build the deflated matrix as a new array instead. The JDK has
no eigenvalue routine to skip, and libraries such as Apache Commons Math
aren’t installed.
// Principal component analysis by power iteration, using only Node's built-in
// modules. Standardizes the four scores, builds their correlation matrix, and
// finds the top two eigenvectors (the principal components): power iteration
// finds the first, then deflation removes it from the matrix so the same steps
// find the second.
import { readFileSync } from "node:fs";
const SUBJECTS = ["math", "physics", "reading", "writing"];
const STEPS = 1000;
// Node has no CSV parser. No field here is quoted, so we split on commas and
// keep the four score columns of each row as numbers.
function readRows(path) {
const [header, ...lines] = readFileSync(path, "utf8").trim().split("\n");
const names = header.split(",");
const indexes = SUBJECTS.map((subject) => names.indexOf(subject));
return lines.map((line) => {
const fields = line.split(",");
return indexes.map((i) => Number(fields[i]));
});
}
const sum = (numbers) => numbers.reduce((acc, x) => acc + x, 0);
const dot = (u, v) => sum(u.map((a, i) => a * v[i]));
const matVec = (m, v) => m.map((row) => dot(row, v));
const normalize = (v) => {
const length = Math.sqrt(dot(v, v));
return v.map((x) => x / length);
};
// Power iteration: multiply by the matrix and rescale, STEPS times.
// Returns the eigenvalue and the eigenvector.
function topEigenvector(m) {
let v = m.map((_, i) => (i === 0 ? 1 : 0));
for (let step = 0; step < STEPS; step++) {
v = normalize(matVec(m, v));
}
// An eigenvector can point either way; make its first entry positive.
if (v[0] < 0) v = v.map((x) => -x);
return [dot(v, matVec(m, v)), v];
}
const rows = readRows(process.argv[2]);
const n = rows.length;
const k = SUBJECTS.length;
const indices = SUBJECTS.map((_, j) => j);
// Standardize each column to mean 0 and (sample) standard deviation 1.
const means = indices.map((j) => sum(rows.map((row) => row[j])) / n);
const sds = indices.map((j) =>
Math.sqrt(sum(rows.map((row) => (row[j] - means[j]) ** 2)) / (n - 1)),
);
const z = rows.map((row) => row.map((x, j) => (x - means[j]) / sds[j]));
// The correlation matrix: the covariance of the standardized columns.
const corr = indices.map((i) => indices.map((j) => sum(z.map((r) => r[i] * r[j])) / (n - 1)));
const [value1, pc1] = topEigenvector(corr);
// Deflation: subtract the first component's share, leaving the rest.
const deflated = corr.map((row, i) => row.map((x, j) => x - value1 * pc1[i] * pc1[j]));
const [value2, pc2] = topEigenvector(deflated);
// Prints one component: its eigenvalue, share of the variance, and loadings.
function report(c, value, pc) {
console.log(`eigenvalue_${c} ${value.toFixed(6)}`);
// Standardized columns each have variance 1, so the total is k.
console.log(`explained_${c} ${(value / k).toFixed(6)}`);
SUBJECTS.forEach((subject, j) => console.log(`pc${c}_${subject} ${pc[j].toFixed(6)}`));
}
console.log(`count ${n}`);
report(1, value1, pc1);
report(2, value2, pc2);
console.log(`student_1_pc1 ${dot(z[0], pc1).toFixed(6)}`);
console.log(`student_1_pc2 ${dot(z[0], pc2).toFixed(6)}`);Vectors are plain arrays and the matrix is an array of rows. Every step is a
map, which also passes each element’s index, so
corr.map((row, i) => row.map((x, j) => ...)) builds the deflated matrix
entry by entry. sum hands reduce a starting 0: without it, reduce on an
empty array throws a TypeError, “Reduce of empty array with no initial
value”. Node has no matrix type and no eigenvalue routine; none of its
built-in modules exports anything named for a matrix, an eigenvector, or PCA,
and DOMMatrix, which browsers provide, is undefined in Node.
# Principal component analysis by power iteration, using only Base Julia.
#
# Standardizes the four scores, builds their correlation matrix, and finds the
# top two eigenvectors (the principal components): power iteration finds the
# first, then deflation removes it from the matrix so the same steps find the
# second.
using Printf
const SUBJECTS = ["math", "physics", "reading", "writing"]
const STEPS = 1000
function top_eigenvector(M)
# Power iteration: multiply by the matrix and rescale, STEPS times.
v = [1.0; zeros(size(M, 1) - 1)]
for _ in 1:STEPS
w = M * v
v = w / sqrt(w' * w)
end
# An eigenvector can point either way; make its first entry positive.
v[1] < 0 && (v = -v)
return v' * M * v, v # eigenvalue, eigenvector
end
function main(path)
lines = readlines(path)
header = split(lines[1], ",")
order = [findfirst(==(s), header) for s in SUBJECTS]
# One row per student, one column per subject.
rows = [parse.(Float64, split(line, ","))[order] for line in lines[2:end]]
X = reduce(vcat, permutedims.(rows))
n, k = size(X)
# Standardize each column to mean 0 and (sample) standard deviation 1.
means = sum(X, dims = 1) / n
sds = sqrt.(sum((X .- means) .^ 2, dims = 1) / (n - 1))
Z = (X .- means) ./ sds
# The correlation matrix: the covariance of the standardized columns.
corr = Z' * Z / (n - 1)
value_1, pc_1 = top_eigenvector(corr)
# Deflation: subtract the first component's share, leaving the rest.
value_2, pc_2 = top_eigenvector(corr - value_1 * pc_1 * pc_1')
println("count $n")
for (c, (value, pc)) in enumerate([(value_1, pc_1), (value_2, pc_2)])
@printf("eigenvalue_%d %.6f\n", c, value)
# Standardized columns each have variance 1, so the total is k.
@printf("explained_%d %.6f\n", c, value / k)
for (subject, loading) in zip(SUBJECTS, pc)
@printf("pc%d_%s %.6f\n", c, subject, loading)
end
end
@printf("student_1_pc1 %.6f\n", Z[1, :]' * pc_1)
@printf("student_1_pc2 %.6f\n", Z[1, :]' * pc_2)
end
main(ARGS[1])Julia writes this the way the math does. corr = Z' * Z / (n - 1) is the
correlation matrix, M * v is a power-iteration step, and deflation is
corr - value_1 * pc_1 * pc_1', where pc_1 * pc_1' is an outer product
because ' turns a column into a row.
The standard library already has the answer: eigen from LinearAlgebra,
with cor from Statistics, gives the same eigenvalues and loadings. But
eigen returned both vectors pointing the other way, with every loading
negative on the first. That’s why every version here makes the first entry
positive.
# Principal component analysis with eigen, from the LinearAlgebra standard library.
#
# The built-in version of ../julia/main.jl. Statistics.cor builds the
# correlation matrix, and eigen finds all its eigenvalues and eigenvectors at
# once. eigen may return a vector with every sign flipped, so we apply the same
# rule as the from-scratch version: each one's first entry is positive.
using LinearAlgebra
using Printf
using Statistics
const SUBJECTS = ["math", "physics", "reading", "writing"]
function main(path)
lines = readlines(path)
header = split(lines[1], ",")
order = [findfirst(==(s), header) for s in SUBJECTS]
# One row per student, one column per subject.
rows = [parse.(Float64, split(line, ","))[order] for line in lines[2:end]]
X = reduce(vcat, permutedims.(rows))
n, k = size(X)
# Standardize with mean and std (std divides by n - 1).
Z = (X .- mean(X, dims = 1)) ./ std(X, dims = 1)
# For a symmetric matrix, eigen sorts eigenvalues from smallest to
# largest, so sortby = - puts the largest first.
values, vectors = eigen(Symmetric(cor(X)), sortby = -)
println("count $n")
pcs = Vector{Vector{Float64}}()
for c in 1:2
pc = vectors[:, c]
pc[1] < 0 && (pc = -pc)
push!(pcs, pc)
@printf("eigenvalue_%d %.6f\n", c, values[c])
# Standardized columns each have variance 1, so the total is k.
@printf("explained_%d %.6f\n", c, values[c] / k)
for (subject, loading) in zip(SUBJECTS, pc)
@printf("pc%d_%s %.6f\n", c, subject, loading)
end
end
@printf("student_1_pc1 %.6f\n", Z[1, :]' * pcs[1])
@printf("student_1_pc2 %.6f\n", Z[1, :]' * pcs[2])
end
main(ARGS[1])eigen comes from LinearAlgebra, a Julia standard library, and finds
every eigenvalue and eigenvector of a matrix at once. Statistics supplies
cor, which builds the correlation matrix directly from the raw scores.
Wrapping the matrix in Symmetric tells eigen it can use the faster
symmetric method, which also guarantees real results. That method returns
the eigenvalues from smallest to largest, so sortby = - reverses the
order. Here both components come back with every sign flipped, so the
first-entry-positive rule is what makes them match.
-- Principal component analysis by power iteration, using only Lua's
-- standard library.
--
-- Standardizes the four scores, builds their correlation matrix, and finds
-- the top two eigenvectors (the principal components): power iteration finds
-- the first, then deflation removes it from the matrix so the same steps find
-- the second. Lua has no matrix type or eigen solver, so a vector is a table
-- used as an array (indexed from 1) and a matrix is a table of row tables.
local SUBJECTS = { "math", "physics", "reading", "writing" }
local STEPS = 1000
-- Returns one array of numbers per row, in the order of the subjects.
local function read_rows(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 rows = {}
for line in f:lines() do
local fields = {}
for field in line:gmatch("[^,]+") do
fields[#fields + 1] = tonumber(field)
end
local row = {}
for j, name in ipairs(SUBJECTS) do
row[j] = fields[assert(col[name], "column not found: " .. name)]
end
rows[#rows + 1] = row
end
f:close()
return rows
end
local function dot(u, v)
local total = 0
for i = 1, #u do
total = total + u[i] * v[i]
end
return total
end
local function mat_vec(m, v)
local result = {}
for i, row in ipairs(m) do
result[i] = dot(row, v)
end
return result
end
local function normalize(v)
local length = math.sqrt(dot(v, v))
local result = {}
for i, x in ipairs(v) do
result[i] = x / length
end
return result
end
-- Power iteration: multiply by the matrix and rescale, STEPS times.
-- Returns two values, the eigenvalue and the eigenvector.
local function top_eigenvector(m)
local v = {}
for i = 1, #m do v[i] = 0.0 end
v[1] = 1.0
for _ = 1, STEPS do
v = normalize(mat_vec(m, v))
end
-- An eigenvector can point either way; make its first entry positive.
if v[1] < 0 then
for i = 1, #v do v[i] = -v[i] end
end
return dot(v, mat_vec(m, v)), v
end
local function say(name, value)
print(string.format("%s %.6f", name, value))
end
local rows = read_rows(arg[1])
local n, k = #rows, #SUBJECTS
-- Standardize each column to mean 0 and (sample) standard deviation 1.
local means, sds = {}, {}
for j = 1, k do
local total = 0
for _, row in ipairs(rows) do total = total + row[j] end
local mean = total / n
local sumsq = 0
for _, row in ipairs(rows) do sumsq = sumsq + (row[j] - mean) ^ 2 end
means[j], sds[j] = mean, math.sqrt(sumsq / (n - 1))
end
local z = {}
for i, row in ipairs(rows) do
z[i] = {}
for j = 1, k do
z[i][j] = (row[j] - means[j]) / sds[j]
end
end
-- The correlation matrix: the covariance of the standardized columns.
local corr = {}
for i = 1, k do
corr[i] = {}
for j = 1, k do
local total = 0
for _, r in ipairs(z) do total = total + r[i] * r[j] end
corr[i][j] = total / (n - 1)
end
end
local value_1, pc_1 = top_eigenvector(corr)
-- Deflation: subtract the first component's share, leaving the rest.
local deflated = {}
for i = 1, k do
deflated[i] = {}
for j = 1, k do
deflated[i][j] = corr[i][j] - value_1 * pc_1[i] * pc_1[j]
end
end
local value_2, pc_2 = top_eigenvector(deflated)
print("count " .. n)
for c, result in ipairs({ { value_1, pc_1 }, { value_2, pc_2 } }) do
local value, pc = result[1], result[2]
say("eigenvalue_" .. c, value)
-- Standardized columns each have variance 1, so the total is k.
say("explained_" .. c, value / k)
for j, subject in ipairs(SUBJECTS) do
say("pc" .. c .. "_" .. subject, pc[j])
end
end
say("student_1_pc1", dot(z[1], pc_1))
say("student_1_pc2", dot(z[1], pc_2))A vector is a table used as an array and the matrix is a table of row
tables, so corr[i][j] reads row i, column j, counting from 1. mat_vec
takes the dot product of each row with the vector, and power iteration is
a for loop that sets v = normalize(mat_vec(m, v)) 1000 times.
top_eigenvector returns two values at once, which
local value_1, pc_1 = ... unpacks.
The built-in we skip doesn’t exist here: Lua’s math library has square
roots and trigonometry but no matrices or eigenvalues. LuaRocks has matrix
modules, but neither LuaRocks nor any module is installed, so
require("matrix") fails with “module ‘matrix’ not found”.
(* Principal component analysis by power iteration, using only OCaml's
standard library.
Standardizes the four scores, builds their correlation matrix, and finds
the top two eigenvectors (the principal components): power iteration
finds the first, then deflation removes it from the matrix so the same
steps find the second. A matrix is an array of row arrays, so m.(i).(j)
is row i, column j, counting from 0. *)
let subjects = [| "math"; "physics"; "reading"; "writing" |]
let steps = 1000
(* The named columns of every row, read as floats, in the order asked for. *)
let read_rows 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 col name =
match List.find_index (String.equal name) header with
| Some i -> i
| None -> failwith ("column not found: " ^ name)
in
let indexes = Array.map col subjects in
Array.of_list
(List.map
(fun line ->
let f = Array.of_list (String.split_on_char ',' line) in
Array.map (fun i -> float_of_string f.(i)) indexes)
(List.tl lines)))
let dot u v = Array.fold_left ( +. ) 0. (Array.map2 ( *. ) u v)
(* A matrix times a vector: each row of the matrix dotted with the vector. *)
let mat_vec m v = Array.map (fun row -> dot row v) m
let normalize v =
let length = sqrt (dot v v) in
Array.map (fun x -> x /. length) v
(* Power iteration: multiply by the matrix and rescale, 1,000 times.
Returns the pair (eigenvalue, eigenvector). *)
let top_eigenvector m =
let v = ref (Array.init (Array.length m) (fun i -> if i = 0 then 1. else 0.)) in
for _ = 1 to steps do
v := normalize (mat_vec m !v) (* ! reads a ref, := replaces its value *)
done;
(* An eigenvector can point either way; make its first entry positive. *)
let v = if !v.(0) < 0. then Array.map Float.neg !v else !v in
(dot v (mat_vec m v), v)
let () =
let rows = read_rows Sys.argv.(1) in
let n = Array.length rows and k = Array.length subjects in
let count = float_of_int n in
(* Standardize each column to mean 0 and (sample) standard deviation 1. *)
let column j = Array.map (fun r -> r.(j)) rows in
let means = Array.init k (fun j -> Array.fold_left ( +. ) 0. (column j) /. count) in
let sds =
Array.init k (fun j ->
let sumsq = Array.fold_left (fun acc x -> acc +. ((x -. means.(j)) ** 2.)) 0. (column j) in
sqrt (sumsq /. (count -. 1.)))
in
let z = Array.map (Array.mapi (fun j x -> (x -. means.(j)) /. sds.(j))) rows in
(* The correlation matrix: the covariance of the standardized columns. *)
let zcol j = Array.map (fun r -> r.(j)) z in
let corr = Array.init k (fun i -> Array.init k (fun j -> dot (zcol i) (zcol j) /. (count -. 1.))) in
let value_1, pc_1 = top_eigenvector corr in
(* Deflation: subtract the first component's share, leaving the rest. *)
let deflated =
Array.mapi (fun i row -> Array.mapi (fun j c -> c -. (value_1 *. pc_1.(i) *. pc_1.(j))) row) corr
in
let value_2, pc_2 = top_eigenvector deflated in
Printf.printf "count %d\n" n;
List.iteri
(fun i (value, pc) ->
let c = i + 1 in
Printf.printf "eigenvalue_%d %.6f\n" c value;
(* Standardized columns each have variance 1, so the total is k. *)
Printf.printf "explained_%d %.6f\n" c (value /. float_of_int k);
Array.iteri (fun j subject -> Printf.printf "pc%d_%s %.6f\n" c subject pc.(j)) subjects)
[ (value_1, pc_1); (value_2, pc_2) ];
Printf.printf "student_1_pc1 %.6f\n" (dot z.(0) pc_1);
Printf.printf "student_1_pc2 %.6f\n" (dot z.(0) pc_2)A matrix is an array of row arrays, so corr.(i).(j) reads row i, column
j, counting from 0. mat_vec maps each row to its dot product with the
vector. Power iteration is a for loop that sets
v := normalize (mat_vec m !v) 1,000 times: v is a ref, a mutable
box, where !v reads it and := replaces it.
The built-in we skip doesn’t exist here. OCaml’s standard library has no
matrix module or eigenvalue routine, and Matrix.mul fails with
Unbound module "Matrix". Owl, the usual numerical library, comes through
opam, but neither opam nor Owl is installed, so the correlation matrix,
deflation, and every product are written out.
% Principal component analysis by power iteration, using only core GNU Octave.
%
% Standardizes the four scores, builds their correlation matrix, and finds the
% top two eigenvectors (the principal components): power iteration finds the
% first, then deflation removes it from the matrix so the same steps find the
% second. Octave has eig and svd, but those are the method itself, so we
% write the iteration out.
SUBJECTS = {"math", "physics", "reading", "writing"};
STEPS = 1000;
% The named column of a CSV file, as a column vector of numbers.
function values = read_column(path, name)
lines = strsplit(strtrim(fileread(path)), "\n");
header = strsplit(strtrim(lines{1}), ",");
index = find(strcmp(header, name));
if isempty(index)
error("no column %s", name);
end
values = zeros(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values(row - 1) = str2double(fields{index});
end
end
% Power iteration: multiply by the matrix and rescale, steps times.
function [value, v] = top_eigenvector(M, steps)
v = [1; zeros(rows(M) - 1, 1)];
for step = 1:steps
w = M * v;
v = w / norm(w); % norm is the vector's length, sqrt(w' * w)
end
% An eigenvector can point either way; make its first entry positive.
if v(1) < 0
v = -v;
end
value = v' * M * v;
end
args = argv();
% One row per student, one column per subject.
X = zeros(0, 0);
for s = 1:numel(SUBJECTS)
X(:, s) = read_column(args{1}, SUBJECTS{s});
end
[n, k] = size(X);
% Standardize each column to mean 0 and (sample) standard deviation 1.
% sum(X, 1) adds down each column, and subtracting the 1-by-k row of means
% subtracts it from every row.
means = sum(X, 1) / n;
sds = sqrt(sum((X - means) .^ 2, 1) / (n - 1));
Z = (X - means) ./ sds;
% The correlation matrix: the covariance of the standardized columns.
corr_matrix = Z' * Z / (n - 1);
[value_1, pc_1] = top_eigenvector(corr_matrix, STEPS);
% Deflation: subtract the first component's share, leaving the rest.
% pc_1 * pc_1' is a column times a row, the 4-by-4 outer product.
[value_2, pc_2] = top_eigenvector(corr_matrix - value_1 * (pc_1 * pc_1'), STEPS);
printf("count %d\n", n);
values = [value_1, value_2];
pcs = [pc_1, pc_2];
for c = 1:2
printf("eigenvalue_%d %.6f\n", c, values(c));
% Standardized columns each have variance 1, so the total is k.
printf("explained_%d %.6f\n", c, values(c) / k);
for s = 1:k
printf("pc%d_%s %.6f\n", c, SUBJECTS{s}, pcs(s, c));
end
end
printf("student_1_pc1 %.6f\n", Z(1, :) * pc_1);
printf("student_1_pc2 %.6f\n", Z(1, :) * pc_2);Octave has matrices built in, so most steps are one line. The correlation
matrix is Z' * Z / (n - 1), and each power-iteration step is
w = M * v followed by v = w / norm(w), where norm is the vector’s
length. Deflation needs the outer product of the first component with
itself, and here that’s just pc_1 * pc_1', a column times a row.
Octave ships eig, svd, cov, and corr, which would do the real work,
so we skip them. We call our matrix corr_matrix, because a variable named
corr would hide the built-in function. Call it after that, and the error
message notes that the variable “shadows function”.
% Principal component analysis with eig, Octave's eigenvalue function.
%
% The built-in version of ../octave/main.m. eig finds every eigenvalue and
% eigenvector of the correlation matrix at once, in place of power
% iteration and deflation.
SUBJECTS = {"math", "physics", "reading", "writing"};
% The named column of a CSV file, as a column vector of numbers.
function values = read_column(path, name)
lines = strsplit(strtrim(fileread(path)), "\n");
header = strsplit(strtrim(lines{1}), ",");
index = find(strcmp(header, name));
if isempty(index)
error("no column %s", name);
end
values = zeros(numel(lines) - 1, 1);
for row = 2:numel(lines)
fields = strsplit(strtrim(lines{row}), ",");
values(row - 1) = str2double(fields{index});
end
end
args = argv();
% One row per student, one column per subject.
X = zeros(0, 0);
for s = 1:numel(SUBJECTS)
X(:, s) = read_column(args{1}, SUBJECTS{s});
end
[n, k] = size(X);
% Standardize each column with mean and std (which divides by n - 1).
Z = (X - mean(X)) ./ std(X);
% eig returns the eigenvectors as the columns of V and the eigenvalues on
% the diagonal of D, smallest first, so sort puts the largest first.
[V, D] = eig(corr(X));
[values, order] = sort(diag(D), "descend");
pcs = V(:, order);
% An eigenvector can point either way; make each one's first entry positive.
pcs = pcs .* sign(pcs(1, :));
printf("count %d\n", n);
for c = 1:2
printf("eigenvalue_%d %.6f\n", c, values(c));
% Standardized columns each have variance 1, so the total is k.
printf("explained_%d %.6f\n", c, values(c) / k);
for s = 1:k
printf("pc%d_%s %.6f\n", c, SUBJECTS{s}, pcs(s, c));
end
end
printf("student_1_pc1 %.6f\n", Z(1, :) * pcs(:, 1));
printf("student_1_pc2 %.6f\n", Z(1, :) * pcs(:, 2));eig finds every eigenvalue and eigenvector of a matrix at once, so it
replaces both power iteration and deflation. Called as [V, D] = eig(M),
it returns the eigenvectors as the columns of V and the eigenvalues on
the diagonal of D, smallest first. sort with "descend" puts the
largest first, and its second result reorders the columns to match.
corr builds the correlation matrix, and mean and std standardize
the scores for the first student’s position. Here eig already returns
each component with a positive first entry, but we apply the sign rule
anyway, as the from-scratch version does.
{ Principal component analysis by power iteration, using only Free
Pascal's run-time library.
Standardizes the four scores, builds their correlation matrix, and finds
the top two eigenvectors (the principal components): power iteration
finds the first, then deflation removes it from the matrix so the same
steps find the second. The Math unit's Mean and StdDev are left out, and no
installed unit has an eigenvalue routine to skip. }
program PrincipalComponents;
{$mode objfpc}{$H+} // 32-bit Integer, and strings longer than 255 characters
uses
SysUtils;
const
K = 4; // subjects
Steps = 1000;
Subjects: array[1..K] of String = ('math', 'physics', 'reading', 'writing');
type
// Array bounds are part of the type, so these count from 1, not 0.
TVector = array[1..K] of Double;
TMatrix = array[1..K, 1..K] of Double;
TRows = array of TVector; // a dynamic array, counted from 0
{ Read the four subject columns, found by name, into rows of K scores. }
function ReadScores(const Path: String): TRows;
var
F: TextFile;
Line: String;
Fields: TStringArray;
Cols: array[1..K] of Integer;
I, J: Integer;
begin
Result := nil;
AssignFile(F, Path);
Reset(F);
ReadLn(F, Line);
Fields := Line.Split(',');
for J := 1 to K do
begin
Cols[J] := -1;
for I := 0 to High(Fields) do
if Fields[I] = Subjects[J] then
Cols[J] := I;
if Cols[J] < 0 then
raise Exception.Create('column not found: ' + Subjects[J]);
end;
while not Eof(F) do
begin
ReadLn(F, Line);
Fields := Line.Split(',');
SetLength(Result, Length(Result) + 1); // grow by one
for J := 1 to K do
Result[High(Result)][J] := StrToFloat(Fields[Cols[J]]);
end;
CloseFile(F);
end;
{ M V. Unlike C, a Pascal function can return a fixed-size array. }
function MatVec(const M: TMatrix; const V: TVector): TVector;
var
I, J: Integer;
begin
for I := 1 to K do
begin
Result[I] := 0;
for J := 1 to K do
Result[I] := Result[I] + M[I, J] * V[J];
end;
end;
function Dot(const U, V: TVector): Double;
var
I: Integer;
begin
Result := 0;
for I := 1 to K do
Result := Result + U[I] * V[I];
end;
{ Power iteration: multiply by the matrix and rescale, Steps times. Fills
V with the eigenvector and returns its eigenvalue. }
function TopEigenvector(const M: TMatrix; out V: TVector): Double;
var
MV: TVector;
Len: Double;
I, Step: Integer;
begin
for I := 1 to K do
V[I] := 0;
V[1] := 1;
for Step := 1 to Steps do
begin
MV := MatVec(M, V);
Len := Sqrt(Dot(MV, MV));
for I := 1 to K do
V[I] := MV[I] / Len;
end;
// An eigenvector can point either way; make its first entry positive.
if V[1] < 0 then
for I := 1 to K do
V[I] := -V[I];
Result := Dot(V, MatVec(M, V)); // the Rayleigh quotient v . (M v)
end;
procedure Say(const Name: String; Value: Double);
begin
WriteLn(Name, ' ', Value:0:6); // width 0, six decimals
end;
procedure Report(C: Integer; Value: Double; const PC: TVector);
var
J: Integer;
begin
Say('eigenvalue_' + IntToStr(C), Value);
// Standardized columns each have variance 1, so the total is K.
Say('explained_' + IntToStr(C), Value / K);
for J := 1 to K do
Say('pc' + IntToStr(C) + '_' + Subjects[J], PC[J]);
end;
var
Z: TRows;
Corr, Deflated: TMatrix;
PC1, PC2: TVector;
Value1, Value2, Total, Squares, Mean, Sd: Double;
N, I, A, B: Integer;
begin
Z := ReadScores(ParamStr(1));
N := Length(Z);
// Standardize each column to mean 0 and (sample) standard deviation 1.
for A := 1 to K do
begin
Total := 0;
for I := 0 to N - 1 do
Total := Total + Z[I][A];
Mean := Total / N;
Squares := 0;
for I := 0 to N - 1 do
Squares := Squares + Sqr(Z[I][A] - Mean);
Sd := Sqrt(Squares / (N - 1));
for I := 0 to N - 1 do
Z[I][A] := (Z[I][A] - Mean) / Sd;
end;
// The correlation matrix: the covariance of the standardized columns.
for A := 1 to K do
for B := 1 to K do
begin
Total := 0;
for I := 0 to N - 1 do
Total := Total + Z[I][A] * Z[I][B];
Corr[A, B] := Total / (N - 1);
end;
Value1 := TopEigenvector(Corr, PC1);
// Deflation: subtract the first component's share, leaving the rest.
for A := 1 to K do
for B := 1 to K do
Deflated[A, B] := Corr[A, B] - Value1 * PC1[A] * PC1[B];
Value2 := TopEigenvector(Deflated, PC2);
WriteLn('count ', N);
Report(1, Value1, PC1);
Report(2, Value2, PC2);
// Z is a dynamic array, counted from 0, so the first student is Z[0].
Say('student_1_pc1', Dot(Z[0], PC1));
Say('student_1_pc2', Dot(Z[0], PC2));
end.The correlation matrix is a fixed array[1..4, 1..4] of Double. In Pascal
the bounds are part of the type, so the subjects count from 1 while the
dynamic array of students counts from 0. No installed unit has an
eigenvalue routine, and we leave out the Math unit’s Mean and StdDev,
so standardizing, power iteration, deflation, and the sign rule are short
loops. Unlike a C function, MatVec returns a whole fixed-size array.
Watch the bounds: range checks are off by default. V[0] on an
array[1..4] fails to compile when the 0 is a constant, but with a
variable index it runs without an error.
# Principal component analysis by power iteration, using only Perl and its
# core modules.
#
# Standardizes the four scores, builds their correlation matrix, and finds
# the top two eigenvectors (the principal components): power iteration finds
# the first, then deflation removes it from the matrix so the same steps find
# the second. Perl has no matrix type, so a vector is an array ref and 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);
my @SUBJECTS = qw(math physics reading writing);
my $STEPS = 1000;
# Returns one array ref of numbers per row, in the order of the subjects.
sub read_rows {
my ($path) = @_;
open my $fh, '<', $path or die "can't open $path: $!";
chomp(my $header = <$fh>);
my @names = split /,/, $header;
my @rows;
while (my $line = <$fh>) {
chomp $line;
my %row;
@row{@names} = split /,/, $line; # a hash slice: names to fields
push @rows, [ @row{@SUBJECTS} ];
}
close $fh;
return @rows;
}
# $#$u is the last index of the array that $u refers to.
sub dot {
my ($u, $v) = @_;
return sum(map { $u->[$_] * $v->[$_] } 0 .. $#$u);
}
sub mat_vec {
my ($m, $v) = @_;
return [ map { dot($_, $v) } @$m ];
}
sub normalize {
my ($v) = @_;
my $length = sqrt dot($v, $v);
return [ map { $_ / $length } @$v ];
}
# Power iteration: multiply by the matrix and rescale, STEPS times.
# Returns two values, the eigenvalue and the eigenvector.
sub top_eigenvector {
my ($m) = @_;
my $v = [ 1.0, (0.0) x $#$m ]; # x repeats a list
$v = normalize(mat_vec($m, $v)) for 1 .. $STEPS;
# An eigenvector can point either way; make its first entry positive.
$v = [ map { -$_ } @$v ] if $v->[0] < 0;
return (dot($v, mat_vec($m, $v)), $v);
}
my @rows = read_rows($ARGV[0]);
my ($n, $k) = (scalar @rows, scalar @SUBJECTS);
# Standardize each column to mean 0 and (sample) standard deviation 1.
my (@means, @sds);
for my $j (0 .. $k - 1) {
my $mean = sum(map { $_->[$j] } @rows) / $n;
push @means, $mean;
push @sds, sqrt(sum(map { ($_->[$j] - $mean) ** 2 } @rows) / ($n - 1));
}
my @z = map {
my $row = $_;
[ map { ($row->[$_] - $means[$_]) / $sds[$_] } 0 .. $k - 1 ];
} @rows;
# The correlation matrix: the covariance of the standardized columns.
my @corr;
for my $i (0 .. $k - 1) {
for my $j (0 .. $k - 1) {
$corr[$i][$j] = sum(map { $_->[$i] * $_->[$j] } @z) / ($n - 1);
}
}
my ($value_1, $pc_1) = top_eigenvector(\@corr);
# Deflation: subtract the first component's share, leaving the rest.
my @deflated;
for my $i (0 .. $k - 1) {
for my $j (0 .. $k - 1) {
$deflated[$i][$j] = $corr[$i][$j] - $value_1 * $pc_1->[$i] * $pc_1->[$j];
}
}
my ($value_2, $pc_2) = top_eigenvector(\@deflated);
print "count $n\n";
my @results = ([ $value_1, $pc_1 ], [ $value_2, $pc_2 ]);
for my $c (1 .. 2) {
my ($value, $pc) = @{ $results[$c - 1] };
printf "eigenvalue_%d %.6f\n", $c, $value;
# Standardized columns each have variance 1, so the total is k.
printf "explained_%d %.6f\n", $c, $value / $k;
printf "pc%d_%s %.6f\n", $c, $SUBJECTS[$_], $pc->[$_] for 0 .. $k - 1;
}
printf "student_1_pc1 %.6f\n", dot($z[0], $pc_1);
printf "student_1_pc2 %.6f\n", dot($z[0], $pc_2);The matrix is an array of array refs, so $corr[$i][$j] reads row i,
column j, counting from 0, and assigning to it creates each row as needed.
Three small subs carry the linear algebra: dot sums products with map,
mat_vec takes the dot product of each row with the vector, and
normalize rescales to length 1. Power iteration is one line that
repeats normalize(mat_vec($m, $v)) 1,000 times.
Perl’s core has sqrt but no matrices or eigenvalues. CPAN’s PDL and
Math::MatrixReal have both, but neither is installed, so we write the
steps out, with sum from List::Util as the only import.
<?php
// Principal component analysis by power iteration, using only PHP and the
// extensions built into it.
//
// Standardizes the four scores, builds their correlation matrix, and finds
// the top two eigenvectors (the principal components): power iteration finds
// the first, then deflation removes it from the matrix so the same steps find
// the second.
declare(strict_types=1);
const SUBJECTS = ['math', 'physics', 'reading', 'writing'];
const STEPS = 1000;
// No field is quoted, so explode splits each line on commas. Each row keeps
// the four score columns, in SUBJECTS order.
function read_rows(string $path): array
{
$lines = file($path, FILE_IGNORE_NEW_LINES | FILE_SKIP_EMPTY_LINES);
$header = explode(',', array_shift($lines));
$indexes = [];
foreach (SUBJECTS as $s) {
$col = array_search($s, $header, true);
if ($col === false) {
throw new RuntimeException("column not found: $s");
}
$indexes[] = $col;
}
return array_map(function ($line) use ($indexes) {
$fields = explode(',', $line);
return array_map(fn($i) => (float) $fields[$i], $indexes);
}, $lines);
}
// Called with two arrays, array_map walks them side by side.
function dot(array $u, array $v): float
{
return array_sum(array_map(fn($a, $b) => $a * $b, $u, $v));
}
function mat_vec(array $m, array $v): array
{
return array_map(fn($row) => dot($row, $v), $m);
}
function normalize(array $v): array
{
$length = sqrt(dot($v, $v));
return array_map(fn($x) => $x / $length, $v);
}
// Power iteration: multiply by the matrix and rescale, STEPS times.
// Returns the eigenvalue and the eigenvector.
function top_eigenvector(array $m): array
{
$v = array_fill(0, count($m), 0.0);
$v[0] = 1.0;
for ($step = 0; $step < STEPS; $step++) {
$v = normalize(mat_vec($m, $v));
}
// An eigenvector can point either way; make its first entry positive.
if ($v[0] < 0) {
$v = array_map(fn($x) => -$x, $v);
}
return [dot($v, mat_vec($m, $v)), $v];
}
$rows = read_rows($argv[1]);
$n = count($rows);
$k = count(SUBJECTS);
// Standardize each column to mean 0 and (sample) standard deviation 1.
// array_column pulls out one subject's scores from every row.
$means = [];
$sds = [];
for ($j = 0; $j < $k; $j++) {
$column = array_column($rows, $j);
$m = array_sum($column) / $n;
$means[] = $m;
$sds[] = sqrt(array_sum(array_map(fn($x) => ($x - $m) ** 2, $column)) / ($n - 1));
}
$z = array_map(
fn($row) => array_map(fn($x, $m, $s) => ($x - $m) / $s, $row, $means, $sds),
$rows
);
// The correlation matrix, an array of row arrays: the covariance of the
// standardized columns. $corr[$i][$j] is row i, column j, from 0.
$corr = [];
for ($i = 0; $i < $k; $i++) {
for ($j = 0; $j < $k; $j++) {
$corr[$i][$j] = dot(array_column($z, $i), array_column($z, $j)) / ($n - 1);
}
}
[$value1, $pc1] = top_eigenvector($corr);
// Deflation: subtract the first component's share, leaving the rest.
$deflated = [];
for ($i = 0; $i < $k; $i++) {
for ($j = 0; $j < $k; $j++) {
$deflated[$i][$j] = $corr[$i][$j] - $value1 * $pc1[$i] * $pc1[$j];
}
}
[$value2, $pc2] = top_eigenvector($deflated);
echo "count $n\n";
foreach ([1 => [$value1, $pc1], 2 => [$value2, $pc2]] as $c => [$value, $pc]) {
printf("eigenvalue_%d %.6f\n", $c, $value);
// Standardized columns each have variance 1, so the total is k.
printf("explained_%d %.6f\n", $c, $value / $k);
foreach (SUBJECTS as $j => $subject) {
printf("pc%d_%s %.6f\n", $c, $subject, $pc[$j]);
}
}
printf("student_1_pc1 %.6f\n", dot($z[0], $pc1));
printf("student_1_pc2 %.6f\n", dot($z[0], $pc2));The matrix is an array of row arrays, so $corr[$i][$j] reads row i,
column j, counting from 0, and assigning to it creates each row as needed.
Three small functions carry the linear algebra: dot sums products with
array_map and array_sum, mat_vec takes the dot product of each row
with the vector, and normalize rescales to length 1. Power iteration
repeats normalize(mat_vec($m, $v)) 1,000 times.
PHP has sqrt but no matrices or eigenvalues, so we write the steps out.
array_column($z, $i) pulls one subject’s column out of every row, which
is how each correlation is built.
-- Principal component analysis by power iteration, in PostgreSQL.
--
-- Standardizes the four scores, builds their correlation matrix, and finds the
-- top two eigenvectors (the principal components): power iteration finds the
-- first, then deflation removes it from the matrix so the same steps find the
-- second. Matrices are stored long, one row per (i, j) entry, and each power
-- iteration is a WITH RECURSIVE query whose state is one row: the step number
-- and the vector, as an array of four numbers.
WITH RECURSIVE scores AS (
-- One row per (student, subject). j numbers the subjects 1 to 4, which is
-- also each subject's position in the vectors below.
SELECT d.row_id, s.j, s.subject, s.x
FROM data d,
LATERAL (VALUES (1, 'math', d.math), (2, 'physics', d.physics),
(3, 'reading', d.reading), (4, 'writing', d.writing)) AS s (j, subject, x)
),
means AS (
-- avg() would do this; sum and count keep the arithmetic visible.
SELECT j, count(*) AS n, sum(x) / count(*) AS mean
FROM scores
GROUP BY j
),
spread AS (
-- Sample standard deviation (divide by n - 1), instead of stddev_samp().
SELECT j, sqrt(sum((x - mean) ^ 2) / (n - 1)) AS sd
FROM scores JOIN means USING (j)
GROUP BY j, n
),
z AS (
-- Each score standardized to mean 0 and standard deviation 1.
SELECT row_id, j, subject, (x - mean) / sd AS z
FROM scores JOIN means USING (j) JOIN spread USING (j)
),
corr AS (
-- The correlation matrix: the covariance of the standardized columns,
-- pairing every student's subject i with the same student's subject j.
SELECT a.j AS i, b.j AS j, sum(a.z * b.z) / (count(*) - 1) AS r
FROM z a JOIN z b USING (row_id)
GROUP BY a.j, b.j
),
power_1 (step, v) AS (
SELECT 0, ARRAY[1, 0, 0, 0]::double precision[]
UNION ALL
-- One step: v <- M v, rescaled to length 1. Aggregates aren't allowed
-- directly in the recursive step, so they sit in LATERAL.
SELECT p.step + 1, next.v
FROM power_1 p,
LATERAL (
SELECT array_agg(mv / len ORDER BY i) AS v
FROM (
SELECT i, mv, sqrt(sum(mv ^ 2) OVER ()) AS len
FROM (SELECT i, sum(r * p.v[j]) AS mv FROM corr GROUP BY i) AS product
) AS lengths
) AS next
WHERE p.step < 1000
),
pc_1 AS (
-- unnest turns the array back into rows, numbered from 1. An eigenvector
-- can point either way; multiplying by sign(v[1]) makes v[1] positive.
SELECT u.j, u.v * sign(p.v[1]) AS v
FROM power_1 p, unnest(p.v) WITH ORDINALITY AS u (v, j)
WHERE p.step = 1000
),
value_1 AS (
-- The eigenvalue is v . (M v).
SELECT sum(a.v * c.r * b.v) AS value
FROM corr c JOIN pc_1 a ON a.j = c.i JOIN pc_1 b ON b.j = c.j
),
deflated AS (
-- Deflation: subtract the first component's share, leaving the rest.
SELECT c.i, c.j, c.r - value * a.v * b.v AS r
FROM corr c JOIN pc_1 a ON a.j = c.i JOIN pc_1 b ON b.j = c.j, value_1
),
power_2 (step, v) AS (
-- The same steps again, on the deflated matrix.
SELECT 0, ARRAY[1, 0, 0, 0]::double precision[]
UNION ALL
SELECT p.step + 1, next.v
FROM power_2 p,
LATERAL (
SELECT array_agg(mv / len ORDER BY i) AS v
FROM (
SELECT i, mv, sqrt(sum(mv ^ 2) OVER ()) AS len
FROM (SELECT i, sum(r * p.v[j]) AS mv FROM deflated GROUP BY i) AS product
) AS lengths
) AS next
WHERE p.step < 1000
),
pc_2 AS (
SELECT u.j, u.v * sign(p.v[1]) AS v
FROM power_2 p, unnest(p.v) WITH ORDINALITY AS u (v, j)
WHERE p.step = 1000
),
value_2 AS (
SELECT sum(a.v * c.r * b.v) AS value
FROM deflated c JOIN pc_2 a ON a.j = c.i JOIN pc_2 b ON b.j = c.j
),
components AS (
SELECT 1 AS c, value, j, v FROM value_1, pc_1
UNION ALL
SELECT 2, value, j, v FROM value_2, pc_2
),
student_1 AS (
-- The first student in file order, in the new coordinates.
SELECT p.c, sum(z.z * p.v) AS score
FROM z JOIN components p USING (j)
WHERE z.row_id = 1
GROUP BY p.c
),
results (position, name, value) AS (
SELECT 0, 'count', n::text FROM means WHERE j = 1
UNION ALL
SELECT c * 10, 'eigenvalue_' || c, round(value::numeric, 6)::text
FROM components WHERE j = 1
UNION ALL
-- Standardized columns each have variance 1, so the total is 4.
SELECT c * 10 + 1, 'explained_' || c, round((value / 4)::numeric, 6)::text
FROM components WHERE j = 1
UNION ALL
SELECT c * 10 + 1 + p.j, 'pc' || c || '_' || subject, round(v::numeric, 6)::text
FROM components p JOIN (SELECT DISTINCT j, subject FROM scores) AS s USING (j)
UNION ALL
SELECT 100 + c, 'student_1_pc' || c, round(score::numeric, 6)::text
FROM student_1
)
SELECT name, value FROM results ORDER BY position;SQL has no matrices, so we store them long: one row per entry, (i, j, r).
Then a matrix times a vector is a GROUP BY i with sum(r * v[j]), and
the correlation matrix itself is a self-join of the standardized scores on
row_id. Power iteration is WITH RECURSIVE twice, once per component,
each row holding the step number and the vector as a double precision[].
PostgreSQL rejects aggregates in the recursive step, so the product sits in
a LATERAL subquery. unnest ... WITH ORDINALITY turns the final array back
into numbered rows. We skip avg, stddev_samp, and corr, which would
give the same correlations; core PostgreSQL has no eigenvector function.
% Principal component analysis by power iteration, using SWI-Prolog and the
% libraries that come with it.
%
% Standardizes the four scores, builds their correlation matrix, and finds
% the top two eigenvectors (the principal components): power iteration
% finds the first, then deflation removes it from the matrix so the same
% steps find the second. A vector is a list of numbers and a matrix is a
% list of rows.
:- initialization(main, main).
% transpose/2, which turns rows into columns, comes from library(clpfd).
:- use_module(library(clpfd), [transpose/2]).
subjects(["math", "physics", "reading", "writing"]).
steps(1000).
% read_rows(+Path, +Names, -Rows): one list of numbers per line, taking the
% columns called Names in that order.
read_rows(Path, Names, Rows) :-
read_file_to_string(Path, Text, []),
split_string(Text, "\n", "", [HeaderLine|Lines]),
split_string(HeaderLine, ",", "", Header),
findall(Row,
( member(Line, Lines),
Line \== "",
split_string(Line, ",", "", Fields),
maplist([Name, Value]>>( nth0(I, Header, Name),
nth0(I, Fields, Field),
number_string(Value, Field) ),
Names, Row) ),
Rows).
dot(U, V, Dot) :-
foldl([A, B, Acc0, Acc]>>(Acc is Acc0 + A * B), U, V, 0, Dot).
% mat_vec(+M, +V, -MV): each entry is one row's dot product with V.
mat_vec(M, V, MV) :-
maplist([Row, X]>>dot(Row, V, X), M, MV).
normalize(V, Unit) :-
dot(V, V, Squared),
Length is sqrt(Squared),
maplist([X, Y]>>(Y is X / Length), V, Unit).
% power(+K, +M, +V, -Final): Final is V after K steps of v <- Mv / |Mv|.
% Each step builds a new list; the recursive call comes last, so Prolog
% reuses the same stack frame.
power(0, _, V, V) :- !.
power(K, M, V, Final) :-
mat_vec(M, V, MV),
normalize(MV, Next),
K1 is K - 1,
power(K1, M, Next, Final).
% top_eigenvector(+M, -Value, -Vector): start from (1, 0, 0, 0), iterate,
% and make the first entry positive, since -v is an eigenvector too.
top_eigenvector(M, Value, Vector) :-
length(M, K),
K0 is K - 1,
length(Zeros, K0),
maplist(=(0.0), Zeros),
steps(Steps),
power(Steps, M, [1.0|Zeros], V),
( V = [First|_], First < 0
-> maplist([X, Y]>>(Y is -X), V, Vector)
; Vector = V
),
mat_vec(M, Vector, MV),
dot(Vector, MV, Value).
% standardize(+Column, -Zs): mean 0 and (sample) standard deviation 1.
standardize(Column, Zs) :-
length(Column, N),
sum_list(Column, Sum),
Mean is Sum / N,
foldl([X, Acc0, Acc]>>(Acc is Acc0 + (X - Mean) ** 2), Column, 0, Squares),
SD is sqrt(Squares / (N - 1)),
maplist([X, Z]>>(Z is (X - Mean) / SD), Column, Zs).
print_component(C, Names, Value, Vector) :-
length(Names, K),
format("eigenvalue_~d ~6f~n", [C, Value]),
% Standardized columns each have variance 1, so the total is K.
Explained is Value / K,
format("explained_~d ~6f~n", [C, Explained]),
forall(nth1(I, Names, Name),
( nth1(I, Vector, Loading),
format("pc~d_~s ~6f~n", [C, Name, Loading]) )).
main :-
current_prolog_flag(argv, [Path|_]),
subjects(Names),
read_rows(Path, Names, Rows),
length(Rows, N),
transpose(Rows, Columns),
maplist(standardize, Columns, ZColumns),
% The correlation matrix: entry (i, j) is the dot product of
% standardized columns i and j, divided by n - 1.
findall(CorrRow,
( member(Zi, ZColumns),
maplist([Zj, R]>>(dot(Zi, Zj, D), R is D / (N - 1)),
ZColumns, CorrRow) ),
Corr),
top_eigenvector(Corr, Value1, PC1),
% Deflation: subtract Value1 * v v^T, leaving the rest.
maplist([Row, Vi, DRow]>>
maplist([Mij, Vj, D]>>(D is Mij - Value1 * Vi * Vj),
Row, PC1, DRow),
Corr, PC1, Deflated),
top_eigenvector(Deflated, Value2, PC2),
format("count ~d~n", [N]),
print_component(1, Names, Value1, PC1),
print_component(2, Names, Value2, PC2),
% The first student's row of standardized scores.
transpose(ZColumns, [Student1|_]),
dot(Student1, PC1, Score1),
dot(Student1, PC2, Score2),
format("student_1_pc1 ~6f~n", [Score1]),
format("student_1_pc2 ~6f~n", [Score2]).A vector is a list of numbers and a matrix is a list of rows, so
mat_vec is one maplist of dot over the rows, and dot is a
foldl that adds up products. transpose, from library(clpfd), turns
the rows of scores into columns, and library(yall) supplies the
[X, Y]>>Goal lambdas; both ship with SWI-Prolog. There’s no eigenvalue
or matrix library to skip.
Variables can’t be reassigned, so each power-iteration step builds a new
list, and power counts down from 1,000, passing it along. The recursive
call comes last, so Prolog reuses its stack frame: a million steps ran
inside a 2 MB stack.
# Principal component analysis by power iteration, using only base R.
#
# Standardizes the four scores, builds their correlation matrix, and finds the
# top two eigenvectors (the principal components). prcomp() or eigen(cor(x))
# would do it in one call; here power iteration finds the first component,
# then deflation removes it from the matrix so the same steps find the second.
subjects <- c("math", "physics", "reading", "writing")
steps <- 1000
args <- commandArgs(trailingOnly = TRUE)
x <- unname(as.matrix(read.csv(args[1])[subjects]))
n <- nrow(x)
k <- ncol(x)
# Standardize each column to mean 0 and sample standard deviation 1 (divide
# by n - 1), written out instead of scale() and sd(). sweep() applies the
# operation column by column (margin 2).
means <- colSums(x) / n
centered <- sweep(x, 2, means)
sds <- sqrt(colSums(centered^2) / (n - 1))
z <- sweep(centered, 2, sds, "/")
# The correlation matrix, instead of cor(): t(z) %*% z sums z_i * z_j over
# the students for every pair of columns at once.
corr <- t(z) %*% z / (n - 1)
top_eigenvector <- function(m) {
# Power iteration: multiply by the matrix and rescale, steps times.
v <- c(1, rep(0, nrow(m) - 1))
for (step in seq_len(steps)) {
v <- drop(m %*% v)
v <- v / sqrt(sum(v^2))
}
# An eigenvector can point either way; make its first entry positive.
if (v[1] < 0) v <- -v
# The Rayleigh quotient v . (M v) is the eigenvalue.
list(value = sum(v * (m %*% v)), vector = v)
}
pc1 <- top_eigenvector(corr)
# Deflation: subtract the first component's share, leaving the rest.
# outer(v, v) is the k-by-k matrix of v_i * v_j.
deflated <- corr - pc1$value * outer(pc1$vector, pc1$vector)
pc2 <- top_eigenvector(deflated)
cat(sprintf("count %d\n", n))
for (c in 1:2) {
pc <- list(pc1, pc2)[[c]]
cat(sprintf("eigenvalue_%d %.6f\n", c, pc$value))
# Standardized columns each have variance 1, so the total is k.
cat(sprintf("explained_%d %.6f\n", c, pc$value / k))
cat(sprintf("pc%d_%s %.6f\n", c, subjects, pc$vector), sep = "")
}
cat(sprintf("student_1_pc1 %.6f\n", sum(z[1, ] * pc1$vector)))
cat(sprintf("student_1_pc2 %.6f\n", sum(z[1, ] * pc2$vector)))R can do all of this in one call: prcomp(x, scale. = TRUE), or
eigen(cor(x)). We avoided both, along with scale(), sd(), and cor().
On this data, prcomp and eigen find the same two eigenvalues (2.496154
and 1.110197) and the same loadings, except that both return PC1 with every
sign flipped, so student 1 scores -1.712316 on it there. Either direction is
a valid eigenvector; we pick the one with a positive first entry.
Matrix arithmetic keeps the rest short. t(z) %*% z / (n - 1) builds the
whole correlation matrix at once, and outer(v, v) gives the table of
v[i] * v[j] that deflation subtracts.
# Principal component analysis with prcomp(), R's built-in PCA.
#
# The built-in version of ../r/main.R. scale. = TRUE standardizes each column
# first, so the components are the eigenvectors of the correlation matrix.
# prcomp() may return any component with every sign flipped, so we apply the
# same rule as the from-scratch version: each one's first entry is positive.
subjects <- c("math", "physics", "reading", "writing")
args <- commandArgs(trailingOnly = TRUE)
x <- read.csv(args[1])[subjects]
k <- ncol(x)
pca <- prcomp(x, scale. = TRUE)
# rotation holds one component per column; sdev holds each one's standard
# deviation, so its square is the eigenvalue.
signs <- sign(pca$rotation[1, ])
rotation <- sweep(pca$rotation, 2, signs, "*")
scores <- sweep(pca$x, 2, signs, "*") # each student's position on each component
eigenvalues <- pca$sdev^2
cat(sprintf("count %d\n", nrow(x)))
for (c in 1:2) {
cat(sprintf("eigenvalue_%d %.6f\n", c, eigenvalues[c]))
# Standardized columns each have variance 1, so the total is k.
cat(sprintf("explained_%d %.6f\n", c, eigenvalues[c] / k))
cat(sprintf("pc%d_%s %.6f\n", c, subjects, rotation[, c]), sep = "")
}
cat(sprintf("student_1_pc1 %.6f\n", scores[1, 1]))
cat(sprintf("student_1_pc2 %.6f\n", scores[1, 2]))prcomp() is base R’s PCA. With scale. = TRUE it standardizes each
column first (the trailing dot is part of the argument’s name), so its
components match the correlation matrix the from-scratch version builds.
It returns three pieces we use: rotation, with one component per column;
sdev, the standard deviation along each, whose square is the eigenvalue;
and x, every student’s score on every component. prcomp() returns the
first component with all four signs negative, so we flip any column whose
first entry is negative, and flip the scores to match.
# Principal component analysis by power iteration, using only the standard library.
# Standardizes the four scores, builds their correlation matrix, and finds the
# top two eigenvectors (the principal components): power iteration finds the
# first, then deflation removes it from the matrix so the same steps find the
# second.
SUBJECTS = %w[math physics reading writing].freeze
STEPS = 1000
# The csv library is a bundled gem since Ruby 3.4, so split lines by hand
# (no quoted fields here) and keep the four score columns of each row.
def read_rows(path)
header, *lines = File.readlines(path, chomp: true).map { |line| line.split(",") }
indexes = SUBJECTS.map { |s| header.index(s) or abort "no column #{s}" }
lines.map { |fields| fields.values_at(*indexes).map { |f| Float(f) } }
end
def dot(u, v)
u.zip(v).sum { |a, b| a * b }
end
def mat_vec(m, v)
m.map { |row| dot(row, v) }
end
def normalize(v)
length = Math.sqrt(dot(v, v))
v.map { |x| x / length }
end
# Power iteration: multiply by the matrix and rescale, STEPS times.
# Returns the eigenvalue and the eigenvector.
def top_eigenvector(m)
v = [1.0] + [0.0] * (m.length - 1)
STEPS.times { v = normalize(mat_vec(m, v)) }
# An eigenvector can point either way; make its first entry positive.
v = v.map { |x| -x } if v[0].negative?
[dot(v, mat_vec(m, v)), v]
end
rows = read_rows(ARGV[0])
n, k = rows.length, SUBJECTS.length
# Standardize each column to mean 0 and (sample) standard deviation 1.
# transpose turns the rows into one array per subject.
columns = rows.transpose
means = columns.map { |c| c.sum / n }
sds = columns.zip(means).map { |c, m| Math.sqrt(c.sum { |x| (x - m)**2 } / (n - 1)) }
z = rows.map { |row| row.zip(means, sds).map { |x, m, s| (x - m) / s } }
# The correlation matrix: the covariance of the standardized columns.
corr = Array.new(k) { |i| Array.new(k) { |j| z.sum { |r| r[i] * r[j] } / (n - 1) } }
value_1, pc_1 = top_eigenvector(corr)
# Deflation: subtract the first component's share, leaving the rest.
deflated = Array.new(k) { |i| Array.new(k) { |j| corr[i][j] - value_1 * pc_1[i] * pc_1[j] } }
value_2, pc_2 = top_eigenvector(deflated)
puts "count #{n}"
[[value_1, pc_1], [value_2, pc_2]].each.with_index(1) do |(value, pc), c|
puts format("eigenvalue_%d %.6f", c, value)
# Standardized columns each have variance 1, so the total is k.
puts format("explained_%d %.6f", c, value / k)
SUBJECTS.zip(pc).each { |subject, loading| puts format("pc%d_%s %.6f", c, subject, loading) }
end
puts format("student_1_pc1 %.6f", dot(z[0], pc_1))
puts format("student_1_pc2 %.6f", dot(z[0], pc_2))Vectors are plain arrays and the matrix is an array of rows, so each step is
one line: transpose turns the rows into columns for the means and standard
deviations, Array.new(k) { |i| ... } builds the correlation matrix, and
power iteration is STEPS.times { v = normalize(mat_vec(m, v)) }.
The built-in we skip is Matrix#eigensystem, from the bundled matrix gem.
On this correlation matrix it found 2.496154 and 1.110197 too, but listed the
eigenvalues smallest first. Its eigenvector for 2.496154 matched ours, and
its one for 1.110197 pointed the opposite way, starting at -0.497992. As in
the other Ruby implementations, we also skip the csv gem and split lines
by hand.
//! Principal component analysis by power iteration, using only the standard
//! library.
//!
//! Standardizes the four scores, builds their correlation matrix, and finds
//! the top two eigenvectors (the principal components): power iteration finds
//! the first, then deflation removes it from the matrix so the same steps find
//! the second.
use std::env;
use std::fs;
const SUBJECTS: [&str; K] = ["math", "physics", "reading", "writing"];
const K: usize = 4;
const STEPS: usize = 1000;
/// Fixed-size arrays: the length is part of the type, and they are `Copy`,
/// so passing one to a function hands over a copy.
type Vector = [f64; K];
type Matrix = [Vector; K];
/// There is no CSV parser in the standard library; the file has no quoted
/// fields, so commas suffice.
fn read(path: &str) -> Vec<Vector> {
let text = fs::read_to_string(path).expect("cannot read data file");
let mut lines = text.lines();
let header: Vec<&str> = lines.next().expect("empty file").split(',').collect();
let index = SUBJECTS.map(|name| {
header
.iter()
.position(|&h| h == name)
.expect("column not found")
});
lines
.map(|line| {
let fields: Vec<&str> = line.split(',').collect();
index.map(|i| fields[i].parse().unwrap())
})
.collect()
}
fn dot(u: &Vector, v: &Vector) -> f64 {
u.iter().zip(v).map(|(a, b)| a * b).sum()
}
fn mat_vec(m: &Matrix, v: &Vector) -> Vector {
m.map(|row| dot(&row, v))
}
fn normalize(v: Vector) -> Vector {
let length = dot(&v, &v).sqrt();
v.map(|x| x / length)
}
/// Power iteration: multiply by the matrix and rescale, STEPS times.
/// Returns the eigenvalue and the eigenvector.
fn top_eigenvector(m: &Matrix) -> (f64, Vector) {
let mut v = [0.0; K];
v[0] = 1.0;
for _ in 0..STEPS {
v = normalize(mat_vec(m, &v));
}
// An eigenvector can point either way; make its first entry positive.
if v[0] < 0.0 {
v = v.map(|x| -x);
}
(dot(&v, &mat_vec(m, &v)), v)
}
fn main() {
let path = env::args().nth(1).expect("usage: main <data.csv>");
let rows = read(&path);
let n = rows.len() as f64;
// Standardize each column to mean 0 and (sample) standard deviation 1.
// `columns` is [0, 1, 2, 3]; mapping over it builds a new fixed-size array.
let columns: [usize; K] = std::array::from_fn(|j| j);
let means = columns.map(|j| rows.iter().map(|r| r[j]).sum::<f64>() / n);
let sds = columns.map(|j| {
let squares: f64 = rows.iter().map(|r| (r[j] - means[j]).powi(2)).sum();
(squares / (n - 1.0)).sqrt()
});
let z: Vec<Vector> = rows
.iter()
.map(|r| columns.map(|j| (r[j] - means[j]) / sds[j]))
.collect();
// The correlation matrix: the covariance of the standardized columns.
let corr: Matrix =
columns.map(|i| columns.map(|j| z.iter().map(|r| r[i] * r[j]).sum::<f64>() / (n - 1.0)));
let (value_1, pc_1) = top_eigenvector(&corr);
// Deflation: subtract the first component's share, leaving the rest.
let deflated = columns.map(|i| columns.map(|j| corr[i][j] - value_1 * pc_1[i] * pc_1[j]));
let (value_2, pc_2) = top_eigenvector(&deflated);
println!("count {}", rows.len());
for (c, (value, pc)) in [(value_1, pc_1), (value_2, pc_2)].iter().enumerate() {
let c = c + 1;
println!("eigenvalue_{c} {value:.6}");
// Standardized columns each have variance 1, so the total is K.
println!("explained_{c} {:.6}", value / K as f64);
for (subject, loading) in SUBJECTS.iter().zip(pc) {
println!("pc{c}_{subject} {loading:.6}");
}
}
println!("student_1_pc1 {:.6}", dot(&z[0], &pc_1));
println!("student_1_pc2 {:.6}", dot(&z[0], &pc_2));
}The vector and matrix are fixed-size arrays, [f64; 4] and [[f64; 4]; 4],
named Vector and Matrix. The length is part of the type: passing a
three-element array to dot is a compile error, “expected an array with a
size of 4, found one with a size of 3”. Arrays of f64 are Copy, so
normalize takes its vector by value with no borrow, and array map returns
another array, so v.map(|x| -x) flips the sign without allocating a Vec.
The standard library has no linear algebra or eigenvalue routine, so there
was no built-in to skip.
// Principal component analysis by power iteration, using only Swift's
// standard library and Foundation, which ships with it.
//
// Standardizes the four scores, builds their correlation matrix, and finds
// the top two eigenvectors (the principal components): power iteration finds
// the first, then deflation removes it from the matrix so the same steps find
// the second.
import Foundation
let subjects = ["math", "physics", "reading", "writing"]
let k = subjects.count
let steps = 1000
// The named column of a CSV file, as numbers. A function that can fail is
// marked throws, and every call to it is marked try.
func readColumn(_ path: String, _ name: String) throws -> [Double] {
let text = try String(contentsOfFile: path, encoding: .utf8)
// split skips empty pieces, so a blank last line doesn't become a row.
var lines = text.split(separator: "\n")
let header = lines.removeFirst().split(separator: ",").map(String.init)
guard let index = header.firstIndex(of: name) else {
fatalError("no column \(name)")
}
// Double(_:) gives nil for text that isn't a number, and ! stops the
// program there instead of carrying on with a missing value.
return lines.map { line in
Double(line.split(separator: ",", omittingEmptySubsequences: false)[index])!
}
}
// Swift has no matrix type: a matrix is an array of rows, [[Double]].
// Arrays are values, so assigning one copies it.
func dot(_ u: [Double], _ v: [Double]) -> Double {
zip(u, v).map { $0 * $1 }.reduce(0, +)
}
func matVec(_ m: [[Double]], _ v: [Double]) -> [Double] {
m.map { row in dot(row, v) }
}
func normalize(_ v: [Double]) -> [Double] {
let length = dot(v, v).squareRoot()
return v.map { $0 / length }
}
// Power iteration: multiply by the matrix and rescale, steps times.
// Returns the eigenvalue and the eigenvector as a named pair (a tuple).
func topEigenvector(_ m: [[Double]]) -> (value: Double, vector: [Double]) {
var v = [Double](repeating: 0, count: k)
v[0] = 1
for _ in 0..<steps {
v = normalize(matVec(m, v))
}
// An eigenvector can point either way; make its first entry positive.
if v[0] < 0 {
v = v.map { -$0 }
}
return (dot(v, matVec(m, v)), v)
}
let path = CommandLine.arguments[1]
let columns = try subjects.map { try readColumn(path, $0) }
let n = columns[0].count
let count = Double(n)
// Standardize each column to mean 0 and (sample) standard deviation 1.
let z = columns.map { column in
let mean = column.reduce(0, +) / count
let squares = column.map { ($0 - mean) * ($0 - mean) }.reduce(0, +)
let sd = (squares / (count - 1)).squareRoot()
return column.map { ($0 - mean) / sd }
}
// The correlation matrix: the covariance of the standardized columns.
let corr = (0..<k).map { i in
(0..<k).map { j in dot(z[i], z[j]) / (count - 1) }
}
let first = topEigenvector(corr)
// Deflation: subtract the first component's share, leaving the rest.
var deflated = corr
for i in 0..<k {
for j in 0..<k {
deflated[i][j] -= first.value * first.vector[i] * first.vector[j]
}
}
let second = topEigenvector(deflated)
func show(_ name: String, _ value: Double) {
print(name, String(format: "%.6f", value))
}
print("count", n)
for (c, pc) in [first, second].enumerated() {
show("eigenvalue_\(c + 1)", pc.value)
// Standardized columns each have variance 1, so the total is k.
show("explained_\(c + 1)", pc.value / Double(k))
for (subject, loading) in zip(subjects, pc.vector) {
show("pc\(c + 1)_\(subject)", loading)
}
}
// The first student's standardized scores, one from each column.
let student = z.map { $0[0] }
show("student_1_pc1", dot(student, first.vector))
show("student_1_pc2", dot(student, second.vector))Swift has no matrix type, so the correlation matrix is a [[Double]], an
array of row arrays. matVec maps dot over the rows, and power
iteration is a loop that replaces v with normalize(matVec(m, v)).
Arrays are values, not references. We deflate by writing
var deflated = corr and then subtracting in place with
deflated[i][j] -= ..., and corr is left unchanged, rows and all. In
Java the same assignment would share the rows. The standard library has
no linear algebra or eigenvalue routine, so there was no built-in to skip,
and Apple’s Accelerate framework isn’t in this image: import Accelerate
fails with “no such module”.
# Principal component analysis by power iteration, using only Tcl's built-in
# commands.
#
# Standardizes the four scores, builds their correlation matrix, and finds
# the top two eigenvectors (the principal components): power iteration finds
# the first, then deflation removes it from the matrix so the same steps find
# the second. Tcl has no matrix type or eigen solver (Tcllib's
# math::linearalgebra has one, but it isn't installed), so a vector is a
# list and a matrix is a list of row lists. expr does the arithmetic; its
# {braces} keep Tcl from substituting twice.
set subjects {math physics reading writing}
set steps 1000
# Returns one list of numbers per row, in the order of the subjects.
proc read_rows {path subjects} {
set f [open $path r]
set header [split [gets $f] ,]
set cols {}
foreach name $subjects {
set col [lsearch -exact $header $name]
if {$col < 0} { error "column not found: $name" }
lappend cols $col
}
set rows {}
while {[gets $f line] >= 0} {
set fields [split $line ,]
lappend rows [lmap col $cols { lindex $fields $col }]
}
close $f
return $rows
}
# foreach can walk two lists side by side, one item from each per pass.
proc dot {u v} {
set total 0.0
foreach a $u b $v {
set total [expr {$total + $a * $b}]
}
return $total
}
# lmap builds a new list from what its body returns on each pass.
proc mat_vec {m v} {
return [lmap row $m { dot $row $v }]
}
proc normalize {v} {
set length [expr {sqrt([dot $v $v])}]
return [lmap x $v { expr {$x / $length} }]
}
# Power iteration: multiply by the matrix and rescale, again and again.
# Returns the eigenvalue and the eigenvector as a two-element list.
proc top_eigenvector {m steps} {
set v [lrepeat [llength $m] 0.0]
lset v 0 1.0
for {set i 0} {$i < $steps} {incr i} {
set v [normalize [mat_vec $m $v]]
}
# An eigenvector can point either way; make its first entry positive.
if {[lindex $v 0] < 0} {
set v [lmap x $v { expr {-$x} }]
}
return [list [dot $v [mat_vec $m $v]] $v]
}
proc say {name value} {
puts "$name [format %.6f $value]"
}
set rows [read_rows [lindex $argv 0] $subjects]
set n [llength $rows]
set k [llength $subjects]
# Standardize each column to mean 0 and (sample) standard deviation 1.
set means {}
set sds {}
for {set j 0} {$j < $k} {incr j} {
set total 0.0
foreach row $rows { set total [expr {$total + [lindex $row $j]}] }
set mean [expr {$total / $n}]
set sumsq 0.0
foreach row $rows {
set sumsq [expr {$sumsq + ([lindex $row $j] - $mean) ** 2}]
}
lappend means $mean
lappend sds [expr {sqrt($sumsq / ($n - 1))}]
}
set z [lmap row $rows {
lmap x $row m $means s $sds { expr {($x - $m) / $s} }
}]
# The correlation matrix: the covariance of the standardized columns.
set corr {}
for {set i 0} {$i < $k} {incr i} {
set matrix_row {}
for {set j 0} {$j < $k} {incr j} {
set total 0.0
foreach r $z {
set total [expr {$total + [lindex $r $i] * [lindex $r $j]}]
}
lappend matrix_row [expr {$total / ($n - 1)}]
}
lappend corr $matrix_row
}
lassign [top_eigenvector $corr $steps] value_1 pc_1
# Deflation: subtract the first component's share, leaving the rest.
# lindex takes several indexes, so [lindex $corr $i $j] is row i, column j.
set deflated {}
for {set i 0} {$i < $k} {incr i} {
set matrix_row {}
for {set j 0} {$j < $k} {incr j} {
lappend matrix_row [expr {[lindex $corr $i $j]
- $value_1 * [lindex $pc_1 $i] * [lindex $pc_1 $j]}]
}
lappend deflated $matrix_row
}
lassign [top_eigenvector $deflated $steps] value_2 pc_2
puts "count $n"
foreach c {1 2} value [list $value_1 $value_2] pc [list $pc_1 $pc_2] {
say eigenvalue_$c $value
# Standardized columns each have variance 1, so the total is k.
say explained_$c [expr {$value / $k}]
foreach subject $subjects loading $pc {
say pc${c}_$subject $loading
}
}
say student_1_pc1 [dot [lindex $z 0] $pc_1]
say student_1_pc2 [dot [lindex $z 0] $pc_2]A vector is a plain list and the matrix is a list of row lists, so
[lindex $corr $i $j] reads row i, column j, and mat_vec is one lmap
that takes the dot product of each row with the vector. Power iteration is a
for loop that sets v [normalize [mat_vec $m $v]] 1000 times, and
deflation builds a second matrix row by row with lappend.
The built-in we skip doesn’t exist here: Tcl’s own commands have no matrix
type or eigen solver. Tcllib’s math::linearalgebra package has one, but it
isn’t installed, so package require math::linearalgebra fails with “can’t
find package”.