Skip to content
Computational Lab.

Conjugate Gradient

Three solvers, one Poisson system.

Advanced / linear-solver / convergence / optimization

The same grounded-box Poisson system as the Electrostatic Potential experiment, handed to three iterative solvers at once. Jacobi and Gauss-Seidel creep; conjugate gradient builds search directions from the residuals and converges in far fewer iterations - plotted per iteration on a log axis, so the race is a measurement rather than a claim.

Tagged: CG, Jacobi, Gauss-Seidel, residual, sparse

Live simulation

Every tick advances all three solvers the same number of iterations.

Speed

Every tick advances all three solvers the same number of iterations.

Mathematics

For a symmetric positive-definite system Ax = b, conjugate gradient minimizes the quadratic f(x) = ½xᵀAx − bᵀx by building search directions from the residuals, each A-orthogonal to the last. In exact arithmetic it converges in n steps; with roundoff it converges whenever the residual says so.

α = (rᵀr) / (pᵀAp)

x ← x + α p; r ← r − α Ap

β = (r_newᵀr_new) / (rᵀr); p ← r + β p

The Poisson matrix 4I − adjacency is SPD, so CG applies. The plot shows relative residual ‖r‖/‖b‖ per iteration for Jacobi, Gauss-Seidel, and CG on identical systems: stationary methods decay geometrically at a rate set by the grid, while CG's polynomial optimizer crushes the low modes first and wins by orders of magnitude. That system is exactly the one Electrostatic Potential animates with SOR.

How it works

All three solvers advance the same number of iterations per tick on their own copies of the system. The operator is applied as a 5-point stencil - the matrix is never stored. Residuals are recorded each iteration into a history buffer and plotted on a log scale; a small inset shows the CG solution field with its equipotential banding.

The convergence test uses the relative residual against the selectable tolerance; a solver that meets it stops and its curve goes flat.

Implementation

src/engine/experiments/conjugate-gradient.ts · 300 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import { sampleInto } from "../render/colormaps";
3import type { ParameterDef, ThemeColors } from "../core/types";
4import { divergingFor } from "../render/diverging-cache";
5
6const N = 56;
7
8/**
9 * Three iterative solvers race on the SAME Poisson system as the
10 * Electrostatic Potential experiment (A = 4I - adjacency, Dirichlet zeros):
11 * Jacobi, Gauss-Seidel, and Conjugate Gradient. Residual norms are recorded
12 * per iteration and plotted on a log axis - CG's advantage is a measurement,
13 * not a slogan.
14 */
15export class ConjugateGradientExperiment extends BaseExperiment {
16 readonly id = "conjugate-gradient";
17
18 private b = new Float64Array(0);
19 private xJacobi = new Float64Array(0);
20 private xGauss = new Float64Array(0);
21 private xCG = new Float64Array(0);
22 private rCG = new Float64Array(0);
23 private pCG = new Float64Array(0);
24 private ap = new Float64Array(0);
25 private xTmp = new Float64Array(0);
26 private history: { jacobi: number; gauss: number; cg: number }[] = [];
27 private iterations = 0;
28 private residualJacobi = 1;
29 private residualGauss = 1;
30 private residualCG = 1;
31 private normB = 1;
32 private cgStarted = false;
33 private rho = 0;
34 private image: ImageData | null = null;
35 private scratch: HTMLCanvasElement | null = null;
36
37 protected params(): ParameterDef[] {
38 return [
39 { key: "iterations", label: "Iterations per tick", min: 1, max: 20, step: 1, defaultValue: 4 },
40 {
41 key: "tolerance",
42 label: "Target rel. residual",
43 options: [
44 { value: "1e-4", label: "1e-4" },
45 { value: "1e-8", label: "1e-8" },
46 { value: "1e-12", label: "1e-12" },
47 ],
48 defaultValue: "1e-8",
49 },
50 {
51 key: "preset",
52 label: "Charge system",
53 options: [
54 { value: "dipole", label: "Dipole (same as #45)" },
55 { value: "quad", label: "Quadrupole" },
56 { value: "single", label: "Single charge" },
57 ],
58 defaultValue: "dipole",
59 },
60 ];
61 }
62
63 protected onReset(): void {
64 const size = N * N;
65 this.b = new Float64Array(size);
66 this.xJacobi = new Float64Array(size);
67 this.xGauss = new Float64Array(size);
68 this.xCG = new Float64Array(size);
69 this.rCG = new Float64Array(size);
70 this.pCG = new Float64Array(size);
71 this.ap = new Float64Array(size);
72 this.xTmp = new Float64Array(size);
73 this.history = [];
74 this.iterations = 0;
75 this.stampCharges();
76 this.normB = this.norm(this.b);
77 this.residualJacobi = this.normB;
78 this.residualGauss = this.normB;
79 this.residualCG = this.normB;
80 this.cgStarted = false;
81 this.rho = 0;
82 }
83
84 private stampCharges(): void {
85 this.b.fill(0);
86 const place = (fx: number, fy: number, q: number) => {
87 const cx = Math.round(fx * N);
88 const cy = Math.round(fy * N);
89 const radius = 2.2;
90 for (let y = Math.max(1, Math.floor(cy - 4)); y < Math.min(N - 1, Math.ceil(cy + 4)); y++) {
91 for (let x = Math.max(1, Math.floor(cx - 4)); x < Math.min(N - 1, Math.ceil(cx + 4)); x++) {
92 const d = Math.hypot(x - cx, y - cy);
93 if (d <= 4) this.b[y * N + x] += q * 8 * Math.exp(-(d * d) / (radius * radius));
94 }
95 }
96 };
97 if (this.str("preset") === "dipole") {
98 place(0.3, 0.5, 1);
99 place(0.7, 0.5, -1);
100 } else if (this.str("preset") === "quad") {
101 place(0.3, 0.3, 1);
102 place(0.7, 0.3, -1);
103 place(0.7, 0.7, 1);
104 place(0.3, 0.7, -1);
105 } else {
106 place(0.5, 0.5, 1);
107 }
108 }
109
110 protected onParameterChange(key: string): void {
111 this.reset();
112 void key;
113 }
114
115 /** y = A x for the 5-point positive-definite Laplacian, Dirichlet zeros. */
116 private applyA(x: Float64Array, out: Float64Array): void {
117 out.fill(0);
118 for (let y = 1; y < N - 1; y++) {
119 for (let x2 = 1; x2 < N - 1; x2++) {
120 const i = y * N + x2;
121 out[i] = 4 * x[i] - x[i - 1] - x[i + 1] - x[i - N] - x[i + N];
122 }
123 }
124 }
125
126 private dot(a: Float64Array, b: Float64Array): number {
127 let sum = 0;
128 for (let i = 0; i < a.length; i++) sum += a[i] * b[i];
129 return sum;
130 }
131
132 private norm(v: Float64Array): number {
133 let sum = 0;
134 for (let i = 0; i < v.length; i++) sum += v[i] * v[i];
135 return Math.sqrt(sum);
136 }
137
138 protected onUpdate(): void {
139 const tolerance = Number(this.str("tolerance"));
140 const perTick = this.num("iterations");
141 for (let k = 0; k < perTick; k++) {
142 this.iterations += 1;
143 // Jacobi: simultaneous update through a scratch buffer.
144 if (this.residualJacobi > tolerance * this.normB) {
145 this.applyA(this.xJacobi, this.ap);
146 for (let i = 0; i < this.b.length; i++) {
147 this.xTmp[i] = this.xJacobi[i] + (this.b[i] - this.ap[i]) / 4;
148 }
149 this.xJacobi.set(this.xTmp);
150 for (let i = 0; i < this.b.length; i++) this.ap[i] = this.b[i];
151 this.applyA(this.xJacobi, this.xTmp);
152 for (let i = 0; i < this.b.length; i++) this.ap[i] -= this.xTmp[i];
153 this.residualJacobi = this.norm(this.ap);
154 }
155 // Gauss-Seidel: in-place sweep (red-black order irrelevant for the norm).
156 if (this.residualGauss > tolerance * this.normB) {
157 for (let y = 1; y < N - 1; y++) {
158 for (let x2 = 1; x2 < N - 1; x2++) {
159 const i = y * N + x2;
160 const sigma =
161 this.xGauss[i - 1] + this.xGauss[i + 1] + this.xGauss[i - N] + this.xGauss[i + N];
162 this.xGauss[i] = (this.b[i] + sigma) / 4;
163 }
164 }
165 this.applyA(this.xGauss, this.ap);
166 for (let i = 0; i < this.b.length; i++) this.ap[i] = this.b[i] - this.ap[i];
167 this.residualGauss = this.norm(this.ap);
168 }
169 // Conjugate gradient with the standard rho recurrence.
170 if (this.residualCG > tolerance * this.normB) {
171 if (!this.cgStarted) {
172 this.cgStarted = true;
173 this.applyA(this.xCG, this.rCG);
174 for (let i = 0; i < this.b.length; i++) this.rCG[i] = this.b[i] - this.rCG[i];
175 this.pCG.set(this.rCG);
176 this.rho = this.dot(this.rCG, this.rCG);
177 }
178 this.applyA(this.pCG, this.ap);
179 const pDotAp = this.dot(this.pCG, this.ap);
180 const alpha = this.rho / (pDotAp || 1e-300);
181 for (let i = 0; i < this.b.length; i++) {
182 this.xCG[i] += alpha * this.pCG[i];
183 this.rCG[i] -= alpha * this.ap[i];
184 }
185 const rhoNew = this.dot(this.rCG, this.rCG);
186 const beta = rhoNew / (this.rho || 1e-300);
187 for (let i = 0; i < this.b.length; i++) {
188 this.pCG[i] = this.rCG[i] + beta * this.pCG[i];
189 }
190 this.rho = rhoNew;
191 this.residualCG = Math.sqrt(rhoNew);
192 }
193 this.history.push({
194 jacobi: this.residualJacobi / this.normB,
195 gauss: this.residualGauss / this.normB,
196 cg: this.residualCG / this.normB,
197 });
198 if (this.history.length > 800) this.history.shift();
199 }
200 }
201
202 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
203 ctx.fillStyle = theme.bg;
204 ctx.fillRect(0, 0, this.width, this.height);
205 const pad = 26;
206 const plotH = this.height - 90;
207 // Log-residual grid: decades
208 ctx.strokeStyle = theme.line;
209 ctx.lineWidth = 1;
210 ctx.font = "10px ui-monospace, monospace";
211 ctx.fillStyle = theme.fgTertiary;
212 for (let decade = 0; decade >= -14; decade -= 2) {
213 const y = pad + (Math.log10(Math.max(1e-14, 1)) - decade) * 0; // gridlines by value below
214 void y;
215 }
216 const toY = (relResidual: number) => {
217 const clamped = Math.min(1, Math.max(1e-14, relResidual));
218 return pad + (-Math.log10(clamped) / 14) * plotH;
219 };
220 const toX = (iteration: number) => (iteration / 600) * (this.width - 2 * pad) + pad;
221 for (let decade = 0; decade >= -14; decade -= 2) {
222 const y = toY(Math.pow(10, decade));
223 ctx.beginPath();
224 ctx.moveTo(pad, y);
225 ctx.lineTo(this.width - pad, y);
226 ctx.stroke();
227 ctx.fillText(`1e${decade}`, 2, y + 3);
228 }
229 const plot = (pick: (h: { jacobi: number; gauss: number; cg: number }) => number, color: string) => {
230 if (this.history.length < 2) return;
231 ctx.strokeStyle = color;
232 ctx.lineWidth = 1.75;
233 ctx.beginPath();
234 for (let i = 0; i < this.history.length; i++) {
235 const x = toX(i);
236 const y = toY(pick(this.history[i]));
237 if (i === 0) ctx.moveTo(x, y);
238 else ctx.lineTo(x, y);
239 }
240 ctx.stroke();
241 };
242 plot((h) => h.jacobi, theme.viz[3]);
243 plot((h) => h.gauss, theme.viz[1]);
244 plot((h) => h.cg, theme.accent);
245 // Legend
246 ctx.font = "11px ui-monospace, monospace";
247 ctx.fillStyle = theme.viz[3];
248 ctx.fillText("JACOBI", pad, this.height - 52);
249 ctx.fillStyle = theme.viz[1];
250 ctx.fillText("GAUSS-SEIDEL", pad + 60, this.height - 52);
251 ctx.fillStyle = theme.accent;
252 ctx.fillText("CG", pad + 160, this.height - 52);
253 ctx.fillStyle = theme.fgTertiary;
254 ctx.fillText(`${this.iterations} iterations`, pad + 195, this.height - 52);
255 // Solution inset: CG field, banded
256 const inset = 92;
257 const insetX = this.width - inset - pad;
258 const insetY = this.height - inset - 16;
259 const diverging = divergingFor(theme);
260 if (!this.image || this.image.width !== N) {
261 this.image = new ImageData(N, N);
262 this.scratch = document.createElement("canvas");
263 this.scratch.width = N;
264 this.scratch.height = N;
265 }
266 let maxAbs = 1e-6;
267 for (let i = 0; i < this.xCG.length; i++) maxAbs = Math.max(maxAbs, Math.abs(this.xCG[i]));
268 const data = this.image.data;
269 for (let i = 0; i < this.xCG.length; i++) {
270 const banded = Math.round((this.xCG[i] / maxAbs) * 8) / 8;
271 sampleInto(diverging, Math.min(1, Math.max(0, 0.5 + banded * 0.5)), data, i * 4);
272 }
273 this.scratch!.getContext("2d")!.putImageData(this.image, 0, 0);
274 ctx.imageSmoothingEnabled = true;
275 ctx.drawImage(this.scratch!, insetX, insetY, inset, inset);
276 ctx.strokeStyle = theme.line;
277 ctx.strokeRect(insetX - 0.5, insetY - 0.5, inset + 1, inset + 1);
278 }
279
280 getMetrics() {
281 const fmt = (r: number) => (r / this.normB).toExponential(2);
282 const tolerance = Number(this.str("tolerance"));
283 return {
284 iterations: this.iterations,
285 jacobi: fmt(this.residualJacobi),
286 gaussSeidel: fmt(this.residualGauss),
287 cg: fmt(this.residualCG),
288 cgDone: this.residualCG <= tolerance * this.normB ? "yes" : "no",
289 };
290 }
291
292 describe(): string {
293 return `After ${this.iterations} iterations on the Poisson system, CG's relative residual is ${(this.residualCG / this.normB).toExponential(1)} versus ${(this.residualJacobi / this.normB).toExponential(1)} for Jacobi.`;
294 }
295
296 entities(): number {
297 return N * N * 3;
298 }
299}
300

Performance

Time per tick
O(cells) per iteration per solver (one stencil application for CG)
Space
Five grid-sized float arrays for CG, two per stationary solver
Notes
Watch the CG curve's staircase: superlinear convergence zones appear as low-frequency error modes are eliminated.

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.