Skip to content
Computational Lab.

SVD & Image Compression

Keep the big singular values, discard the rest.

Advanced / linear-algebra / decomposition / approximation

Every matrix factors into rotation, scale, rotation. Keeping only the largest singular values gives the best possible rank-k approximation of the image, and the metrics show exactly what that trade costs: energy retained against compression ratio and measured reconstruction error. Computed with one-sided Jacobi orthogonalization, right on the page.

Tagged: SVD, rank-k approximation, compression, Jacobi

Live simulation

The rank slider is the compression dial; the view selector shows the error field.

Speed

The rank slider is the compression dial; the view selector shows the error field.

Mathematics

The singular value decomposition factors any matrix into orthogonal rotations and a diagonal scale: A = UΣVᵀ. Truncating to the k largest singular values gives the best rank-k approximation of A in the Frobenius norm - the Eckart-Young theorem, and the entire theory of the compression slider.

A = U Σ Vᵀ, σ₁ ≥ σ₂ ≥ … ≥ 0

A_k = Σ_{j≤k} σ_j u_j v_jᵀ minimizes ‖A − B‖_F over rank-k B

energy retained: Σ_{j≤k} σ_j² / Σ_j σ_j²

The decomposition here is computed by one-sided Jacobi orthogonalization: rotate column pairs until all columns of A are orthogonal; their norms are the singular values and the accumulated rotations form V. The reconstruction formula then collapses to a sum of column products. Zero iterations, no learning rate - the same problem least squares solves, seen from the spectral side.

How it works

The 64×64 image is a matrix. One-sided Jacobi sweeps rotate column pairs (p, q) by the angle that zeroes their inner product, accumulating the rotations into V; after a few sweeps the off-diagonal Gram values drop below tolerance and the columns are the scaled left singular vectors.

The rank-k reconstruction sums the top k columns against their V rows. The split view compares original and reconstruction; the error view maps the absolute difference through a token-built diverging colormap.

Implementation

src/engine/experiments/svd-compression.ts · 246 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import { MAGMA, sampleInto } from "../render/colormaps";
3import type { ParameterDef, ThemeColors } from "../core/types";
4import { divergingFor } from "../render/diverging-cache";
5
6const SIZE = 64;
7
8function makeImage(pattern: number): Float64Array {
9 const a = new Float64Array(SIZE * SIZE);
10 for (let y = 0; y < SIZE; y++) {
11 for (let x = 0; x < SIZE; x++) {
12 const u = x / SIZE;
13 const v = y / SIZE;
14 let value: number;
15 if (pattern === 0) {
16 // Structured shapes: gradient + disk + band, low effective rank.
17 value =
18 0.35 * u +
19 (Math.hypot(u - 0.35, v - 0.4) < 0.18 ? 0.5 : 0) +
20 (Math.abs(v - 0.5 + 0.5 * u) < 0.05 ? 0.3 : 0);
21 } else {
22 // Smooth peaks: sum of three Gaussians, medium rank.
23 const g = (cx: number, cy: number, sx: number, sy: number, amp: number) =>
24 amp * Math.exp(-(((u - cx) / sx) ** 2 + ((v - cy) / sy) ** 2));
25 value =
26 g(0.3, 0.35, 0.18, 0.22, 0.9) +
27 g(0.7, 0.6, 0.2, 0.16, 0.7) +
28 g(0.5, 0.15, 0.3, 0.1, 0.4);
29 }
30 a[y * SIZE + x] = Math.min(1, Math.max(0, value));
31 }
32 }
33 return a;
34}
35
36/**
37 * Truncated SVD image compression, computed for real: one-sided Jacobi
38 * orthogonalization (Hestenes) yields singular values, U, and V for the
39 * 64x64 matrix; the rank-k reconstruction is compared against the original
40 * with measured error and the true storage ratio.
41 */
42export class SvdCompressionExperiment extends BaseExperiment {
43 readonly id = "svd-compression";
44
45 private original: Float64Array = new Float64Array(SIZE * SIZE);
46 private reconstruction: Float64Array = new Float64Array(SIZE * SIZE);
47 /** Descending singular values. */
48 private singularValues: number[] = [];
49 /** Column indices of `uColumns`/`vColumns` in descending-sigma order. */
50 private order: number[] = [];
51 private uColumns = new Float64Array(0) as Float64Array<ArrayBuffer>;
52 private vColumns = new Float64Array(0) as Float64Array<ArrayBuffer>;
53 private rank = 0;
54 private image: ImageData | null = null;
55 private scratch: HTMLCanvasElement | null = null;
56
57 protected params(): ParameterDef[] {
58 return [
59 { key: "rank", label: "Rank k", min: 1, max: 32, step: 1, defaultValue: 6 },
60 {
61 key: "image",
62 label: "Image",
63 options: [
64 { value: "shapes", label: "Shapes (low rank)" },
65 { value: "peaks", label: "Peaks (medium rank)" },
66 ],
67 defaultValue: "shapes",
68 },
69 {
70 key: "view",
71 label: "View",
72 options: [
73 { value: "split", label: "Original | rank-k" },
74 { value: "error", label: "Absolute error" },
75 ],
76 defaultValue: "split",
77 },
78 ];
79 }
80
81 protected onReset(): void {
82 this.original = makeImage(this.str("image") === "peaks" ? 1 : 0);
83 this.solve();
84 }
85
86 /** Static image: recomputes only when rank or image changes. */
87 protected onUpdate(): void {}
88
89 protected onParameterChange(key: string): void {
90 if (key === "rank") this.reconstruct();
91 else if (key === "image") this.reset();
92 }
93
94 private solve(): void {
95 const n = SIZE;
96 // One-sided Jacobi: rotating A's column pairs to orthogonality leaves
97 // A = U * Sigma with V as the accumulated rotation product.
98 const a = new Float64Array(this.original);
99 const v = new Float64Array(n * n);
100 for (let i = 0; i < n; i++) v[i * n + i] = 1;
101 for (let sweep = 0; sweep < 60; sweep++) {
102 let off = 0;
103 for (let p = 0; p < n - 1; p++) {
104 for (let q = p + 1; q < n; q++) {
105 let alpha = 0;
106 let beta = 0;
107 let gamma = 0;
108 for (let k = 0; k < n; k++) {
109 const ap = a[k * n + p];
110 const aq = a[k * n + q];
111 alpha += ap * ap;
112 beta += aq * aq;
113 gamma += ap * aq;
114 }
115 if (Math.abs(gamma) < 1e-14) continue;
116 off = Math.max(off, Math.abs(gamma) / Math.sqrt(alpha * beta));
117 const zeta = (beta - alpha) / (2 * gamma);
118 const t = Math.sign(zeta || 1) / (Math.abs(zeta) + Math.sqrt(1 + zeta * zeta));
119 const c = 1 / Math.sqrt(1 + t * t);
120 const s = c * t;
121 for (let k = 0; k < n; k++) {
122 const ap = a[k * n + p];
123 const aq = a[k * n + q];
124 a[k * n + p] = c * ap - s * aq;
125 a[k * n + q] = s * ap + c * aq;
126 const vp = v[k * n + p];
127 const vq = v[k * n + q];
128 v[k * n + p] = c * vp - s * vq;
129 v[k * n + q] = s * vp + c * vq;
130 }
131 }
132 }
133 if (off < 1e-10) break;
134 }
135 const columns: { sigma: number; index: number }[] = [];
136 for (let j = 0; j < n; j++) {
137 let sigmaSq = 0;
138 for (let k = 0; k < n; k++) sigmaSq += a[k * n + j] ** 2;
139 columns.push({ sigma: Math.sqrt(sigmaSq), index: j });
140 }
141 columns.sort((x, y) => y.sigma - x.sigma);
142 this.singularValues = columns.map((c) => c.sigma);
143 this.order = columns.map((c) => c.index);
144 this.uColumns = a;
145 this.vColumns = v;
146 this.reconstruct();
147 }
148
149 /** A_k = sum over the top k of (sigma_j * u_j) v_j^T = sum of a_col_j * v_col_j^T.
150 * Note V is stored row-major: v_j[col] lives at vColumns[col * n + j]. */
151 private reconstruct(): void {
152 const n = SIZE;
153 const k = Math.min(this.num("rank"), this.order.length);
154 this.reconstruction.fill(0);
155 for (let ki = 0; ki < k; ki++) {
156 const j = this.order[ki];
157 if (this.singularValues[ki] < 1e-12) break;
158 for (let row = 0; row < n; row++) {
159 const uSigma = this.uColumns[row * n + j];
160 if (uSigma === 0) continue;
161 for (let col = 0; col < n; col++) {
162 this.reconstruction[row * n + col] += uSigma * this.vColumns[col * n + j];
163 }
164 }
165 }
166 this.rank = k;
167 }
168
169 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
170 ctx.fillStyle = theme.bg;
171 ctx.fillRect(0, 0, this.width, this.height);
172 if (!this.image || this.image.width !== SIZE || this.image.height !== SIZE) {
173 this.image = new ImageData(SIZE, SIZE);
174 this.scratch = null;
175 }
176 const diverging = divergingFor(theme);
177 const data = this.image.data;
178 const view = this.str("view");
179 const half = SIZE / 2;
180 for (let y = 0; y < SIZE; y++) {
181 for (let x = 0; x < SIZE; x++) {
182 const i = y * SIZE + x;
183 if (view === "error") {
184 sampleInto(
185 diverging,
186 Math.min(1, Math.max(0, 0.5 + (this.reconstruction[i] - this.original[i]) * 2.5)),
187 data,
188 i * 4,
189 );
190 } else {
191 const value = x < half ? this.original[i] : this.reconstruction[i];
192 sampleInto(MAGMA, value, data, i * 4);
193 }
194 }
195 }
196 if (!this.scratch) {
197 this.scratch = document.createElement("canvas");
198 this.scratch.width = SIZE;
199 this.scratch.height = SIZE;
200 }
201 const offCtx = this.scratch.getContext("2d")!;
202 offCtx.putImageData(this.image, 0, 0);
203 ctx.imageSmoothingEnabled = true;
204 ctx.drawImage(this.scratch, 0, 0, this.width, this.height);
205 if (view === "split") {
206 ctx.strokeStyle = theme.fgSecondary;
207 ctx.beginPath();
208 ctx.moveTo(this.width / 2, 0);
209 ctx.lineTo(this.width / 2, this.height);
210 ctx.stroke();
211 }
212 }
213
214 getMetrics() {
215 let retained = 0;
216 let totalEnergy = 0;
217 for (let i = 0; i < this.singularValues.length; i++) {
218 totalEnergy += this.singularValues[i] ** 2;
219 if (i < this.rank) retained += this.singularValues[i] ** 2;
220 }
221 let errSq = 0;
222 let totalSq = 0;
223 for (let i = 0; i < this.original.length; i++) {
224 const d = this.reconstruction[i] - this.original[i];
225 errSq += d * d;
226 totalSq += this.original[i] * this.original[i];
227 }
228 const stored = this.rank * (2 * SIZE + 1);
229 return {
230 rank: this.rank,
231 compression: `${(SIZE * SIZE / stored).toFixed(1)} : 1`,
232 energyRetained: `${((retained / totalEnergy) * 100).toFixed(1)} %`,
233 relError: Math.sqrt(errSq / Math.max(1e-12, totalSq)).toExponential(2),
234 };
235 }
236
237 describe(): string {
238 const metrics = this.getMetrics();
239 return `Rank-${this.rank} reconstruction keeps ${metrics.energyRetained} of the image energy at a ${metrics.compression} compression ratio.`;
240 }
241
242 entities(): number {
243 return SIZE * SIZE;
244 }
245}
246

Performance

Time per tick
O(sweeps · n³) for the solve, once per image; O(n²k) per rank change
Space
Three 64×64 float matrices
Notes
The compression ratio is the true storage count: k(2n + 1) numbers against n² pixels.

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