SVD and Arbitrary Matrices
The mesh
theorem was a triumph with an asterisk: any unitary. But walk over to the
machine-learning building and ask what matrix they would like multiplied today, and the answer is a
trained weight layer — entries of every size and sign, rows nowhere near orthonormal, power
conservation violated in both directions at once. Unitary matrices are a soap-film-thin subset of
all matrices; if photonic computing stopped at the mesh theorem it would be a beautiful instrument
with no repertoire. This lesson removes the asterisk. The tool is the
singular value
decomposition, and the punchline deserves stating up front: any matrix is two
unitaries with a diagonal squeeze between them — photonically, two Clements meshes with a
row of attenuators in the corridor. Everything a neural network's linear layers will ever ask of a
photonic chip is buildable from parts this course already owns.
The factorisation, and what it looks like
The SVD writes any N \times N matrix (real or complex, invertible or
not) as
M \;=\; U\,\Sigma\,V^\dagger, \qquad
\Sigma = \mathrm{diag}(\sigma_1, \dots, \sigma_N),\quad
\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_N \ge 0,
with U and V unitary. Geometrically:
every linear map, however messy, is a rotation, then a stretch along perpendicular axes,
then another rotation. The stretching factors — the singular values — carry all the
non-unitariness; the two flanking matrices are exactly the kind of object a mesh can be programmed
to be.
A worked example that will follow us into the code: the shear
M = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix} — the very matrix the
beamsplitter
lesson proved no lossless device can implement. Form
M^\dagger M = \begin{pmatrix} 1 & 1 \\ 1 & 2 \end{pmatrix}; its
eigenvalues are \lambda_{\pm} = (3 \pm \sqrt 5)/2, so the singular
values are
\sigma_1 = \sqrt{\lambda_+} = \frac{1+\sqrt 5}{2} \approx 1.618, \qquad
\sigma_2 = \sqrt{\lambda_-} = \frac{\sqrt 5 - 1}{2} \approx 0.618 .
The golden ratio and its reciprocal, of all things — a pleasing coincidence of this particular
shear, but the structure is general: \sigma_1 \sigma_2 = |\det M| = 1,
so the shear amplifies along one axis by 1.618 while squeezing the perpendicular axis by
the same factor. That amplification is precisely why no passive chip can do it raw — and points
straight at the fix.
The sandwich architecture
Read M = U \Sigma V^\dagger right to left as light does, and the SVD is
not merely computable — it is a floorplan:
The division of labour is exquisite. The two meshes, being unitary, are implementable losslessly —
they are pure interference, all the machinery of the last three lessons. The entire departure from
unitarity is concentrated into N independent, diagonal, real,
non-negative numbers — no interference, no coupling, just per-channel dimming, which is
the single easiest operation in photonics: a
modulator
per waveguide (or an MZI used as a variable attenuator, dumping the surplus into an unused port).
Component bill for N \times N: 2 \times N(N-1)/2 = N(N-1)
MZIs plus N modulators.
- every matrix M factors as U \Sigma V^\dagger
with U, V unitary and \Sigma diagonal,
non-negative;
- a mesh–attenuators–mesh chip therefore implements M/\sigma_1
exactly: the largest singular value is normalised to transmission 1 and every other
channel is dimmed to \sigma_k/\sigma_1 (power transmission
(\sigma_k/\sigma_1)^2);
- the known scalar \sigma_1 is restored digitally at the readout —
a single multiplication;
- rank deficiency is free: a zero singular value is an attenuator turned fully off.
The 1/\sigma_1 normalisation is not a cheat but a physical necessity:
passive optics can redistribute and discard power, never mint it. A matrix that stretches must be
implemented as "shrink everything else", and the bookkeeping factor rides along as metadata.
Worked check: the golden shear, decomposed and rebuilt
The program computes the shear's SVD from scratch — eigen-decomposition of
M^\dagger M, singular values, both unitaries — then rebuilds
M from the factors and prints the attenuator settings a chip would use:
// SVD of the real 2×2 shear M = [[1,1],[0,1]] via eigenvectors of MᵀM.
const M = [ [1, 1], [0, 1] ];
// MᵀM is symmetric: [[a,b],[b,d]].
const a = M[0][0] * M[0][0] + M[1][0] * M[1][0];
const b = M[0][0] * M[0][1] + M[1][0] * M[1][1];
const d = M[0][1] * M[0][1] + M[1][1] * M[1][1];
const tr = a + d, det = a * d - b * b;
const disc = Math.sqrt(tr * tr - 4 * det);
const lam1 = (tr + disc) / 2, lam2 = (tr - disc) / 2;
const s1 = Math.sqrt(lam1), s2 = Math.sqrt(lam2);
console.log("singular values: " + s1.toFixed(4) + ", " + s2.toFixed(4));
console.log("their product (=|det M|): " + (s1 * s2).toFixed(4));
// Eigenvector of [[a,b],[b,d]] for eigenvalue lam: (b, lam − a), normalised.
function evec(lam: number): number[] {
const n = Math.hypot(b, lam - a);
return [b / n, (lam - a) / n];
}
const v1 = evec(lam1), v2 = evec(lam2); // columns of V
// Columns of U: u_k = M v_k / σ_k.
const u1 = [ (M[0][0] * v1[0] + M[0][1] * v1[1]) / s1, (M[1][0] * v1[0] + M[1][1] * v1[1]) / s1 ];
const u2 = [ (M[0][0] * v2[0] + M[0][1] * v2[1]) / s2, (M[1][0] * v2[0] + M[1][1] * v2[1]) / s2 ];
// Rebuild M = U Σ Vᵀ and measure the worst entry error.
let worst = 0;
for (let i = 0; i < 2; i++)
for (let j = 0; j < 2; j++) {
const U = [ [u1[0], u2[0]], [u1[1], u2[1]] ];
const V = [ [v1[0], v2[0]], [v1[1], v2[1]] ];
const rebuilt = s1 * U[i][0] * V[j][0] + s2 * U[i][1] * V[j][1];
worst = Math.max(worst, Math.abs(rebuilt - M[i][j]));
}
console.log("reconstruction error of U Σ Vᵀ: " + worst.toExponential(2));
// What the chip actually programs: Σ/σ₁ as amplitude (and power) transmissions.
console.log("attenuator amplitudes: 1.0000, " + (s2 / s1).toFixed(4));
console.log("attenuator powers: 1.0000, " + ((s2 / s1) ** 2).toFixed(4));
console.log("digital rescale factor σ₁: " + s1.toFixed(4));
Reconstruction error at machine precision; the "impossible" shear reduced to two rotations, one
channel dimmed to amplitude \sigma_2/\sigma_1 = 1/\varphi^2 \approx 0.382,
and a note to the DSP saying "multiply the answer by 1.618". Nothing was amplified; the matrix's
stretch was repackaged as everyone-else's shrink.
Why this is the neural-network theorem
Recall what a
forward
pass spends its time on: at every layer, y = Wx (then a
cheap nonlinearity). The weight matrix W is arbitrary — and now, up to
one scalar, photonically implementable in a single pass of light. This closes the loop
this module opened: interference computes unitaries; SVD stretches unitaries over all of linear
algebra; therefore interference computes linear algebra. Every photonic neural-network
accelerator in Module
5 is some engineering translation of the sandwich on this page (and its practical
variants — many chips implement W in other encodings, but the SVD mesh
is the conceptual anchor and the proof of possibility). Worth noting the price tags, though: two
meshes double the MZI count, depth, loss and calibration burden relative to a single unitary — one
reason some designs restrict layers to unitary weight matrices on purpose, trading a little
trainability for half the hardware.
If M = U \Sigma V^\dagger is a chip, is
M^{-1} = V \Sigma^{-1} U^\dagger just the chip read right-to-left?
Almost — and the "almost" is instructive. Reversing the meshes is easy: a unitary's inverse is its
conjugate transpose, another unitary, another programmable mesh. But inverting
\Sigma means replacing each transmission
\sigma_k/\sigma_1 by its reciprocal — the channels you dimmed
most now need the most gain, and passive optics has none to offer. The best a passive chip
can do is implement M^{-1} scaled by its own worst factor
\sigma_N, whereupon small singular values crush the whole chip's
throughput: the dynamic range you must squeeze through the attenuators is the condition number
\kappa = \sigma_1/\sigma_N. For genuinely ill-conditioned
M, the honest move is the
pseudoinverse:
clip the hopeless singular values to zero and invert only the healthy subspace — regularisation,
enforced by the laws of optics rather than by a numerical analyst's taste. The
calibration
lesson returns to condition numbers with error bars attached.
The diagram whispers a temptation: if the middle row can dim each channel, why not boost
a channel instead and skip the 1/\sigma_1 bookkeeping? Because gain is
not a passive operation. An amplifying element (a semiconductor optical amplifier, a pumped
section of gain medium) must be driven by an external power source — there goes the energy budget
— and, more fundamentally, amplification is never clean: every optical amplifier injects
amplified spontaneous emission, a quantum-mandated noise floor, into precisely the analog signal
you were trying to scale. The attenuation-only sandwich is not a limitation grudgingly accepted
but the correct engineering choice: implement M/\sigma_1
noiselessly and let a digital multiplier — which amplifies numbers, not photons — restore the
scale for free. When you meet a photonic-accelerator paper, check this detail: a quoted matrix
throughput that ignores the \sigma_1 normalisation, or silently
assumes noiseless gain, is quoting a machine that cannot be built.
Where this goes next
The hardware can now hold any matrix. The
next
lesson finally runs it: encode a vector on the inputs, let one transit of light
perform all N^2 multiply–accumulates at once, detect the result — and
do the honest accounting of latency, throughput and energy against the systolic arrays of
electronic accelerators, including the conversion tolls at the edges that photonic marketing
prefers not to mention. It is the lesson this module has been building toward, flagship simulation
included.