¹

²

Unitary transformations lie at the foundation of quantum mechanics, governing the evolution of quantum states in closed systems [1], and are central to modern photonic computing and quantum information processing [2]–[4]. In linear optical systems, unitary matrices determine how multimode fields interfere and propagate through reconfigurable interferometric networks, allowing direct hardware implementations of matrix operations with ultrahigh speed, low energy consumption, and intrinsic parallelism [4], [5]. These capabilities make photonic platforms attractive for optical neural networks [6]–[11], where unitaries implement trainable linear layers, and for quantum computation [2], [12]–[15], where they realize universal linear-optical circuits. Rapid advances in integrated photonics have produced large-scale, low-loss, and programmable interferometer meshes [16]–[18], motivating scalable and systematic methods to realize arbitrary high-dimensional unitaries in hardware.

In practice, arbitrary \(N \times N\) unitary transformations are implemented by decomposing them into sequences of Mach–Zehnder interferometers (MZIs), which serve as universal two-dimensional beam-splitter gates [19], [20]. Mathematically, any matrix can be expressed via singular value decomposition (SVD) as \(W = U \Sigma V^\dagger\) [21], producing two unitary matrices \(U\) and \(V^\dagger\) and a diagonal matrix \(\Sigma\). Since diagonal matrices \(\Sigma\) can be directly realized with MZIs, the challenge reduces to efficiently decomposing the unitary components \(U\) and \(V^\dagger\). Classical schemes, such as the Reck triangular mesh [22] and the Clements rectangular mesh [23], provide hand-derived constructions using \(N(N-1)/2\) MZIs, forming the algorithmic backbone for universal linear optics and underpinning most present-day photonic processors [7], [12].

Despite their foundational role, these schemes reflect specific human-designed elimination strategies derived from analytic reasoning and occupy only a small region of the combinatorial space of possible elimination orderings, leaving open the question of whether alternative, equally fundamental strategies exist. To address this, we employ automated reasoning to systematically explore this landscape, uncovering strategies that extend beyond human intuition and yield generalizable constructions. In pursuit of this goal, we develop a synthesis engine—built on the DreamCoder framework [24] and extended to manipulate complex-valued linear algebra—that autonomously discovers strategies for decomposing unitary matrices. The synthesized programs are not merely confined to known triangular or rectangular meshes. Instead, the system identifies previously unreported elimination orderings that achieve the optimal \(N(N-1)/2\) MZIs. Analysis of the learned programs reveals that these constructions encode dimension-agnostic algorithmic regularities: strategies discovered on \(5 \times 5\) matrices generalize to arbitrary higher dimensions such as \(64 \times 64\) without retraining, indicating that the system has internalized domain-general patterns rather than dimension-specific heuristics. These results establish program synthesis as a means of autonomously identifying algorithmic structure beyond known human-designed schemes, enabling the discovery of interpretable and generalizable rules for universal unitary decompositions. Beyond universal constructions, the same synthesis engine adapts to structured inputs, automatically discovering structure-specific reductions independent of any prior assumptions about the matrix structure. For analytically structured matrices such as Householder reflectors, the system discovers a dimension-independent rule requiring only \(2N-3\) MZIs, reducing complexity from quadratic to linear scaling, and generalizes to arbitrary \(N\) without retraining. For matrices obtained from the SVD of sparse matrices, it uncovers structure-dependent reductions that tend to increase with sparsity, reaching up to \(38\%\) fewer MZIs than the universal theoretical bound \(N(N-1)/2\) at \(95\%\) sparsity. Together, these results position program synthesis as a unified framework for both discovering generalizable decomposition rules and exploiting matrix-specific structure, enabling scalable automated design of photonic unitary circuits.

As shown in Fig. 1(a), our learning objective is to discover programs that transform arbitrary unitary matrices into diagonal form. This formulation directly corresponds to photonic implementation: synthesizing an arbitrary unitary matrix \(U\) on a photonic chip decomposes into two steps—an MZI cascade network \(P\) and a diagonal matrix \(D\) such that \(U = P^\dagger D\). Since \(D\) requires no algorithmic search, our learning objective simplifies to finding the sequence \(P\) that performs the diagonalization. Training tasks consist of randomly generated unitary matrices \(U_k\) [25]. A program \(p\) transforms \(U_k\) into a diagonal matrix \(D_k = p(U_k)\), where \(D_k\) retains phase information \(e^{i\phi_j}\) on its diagonals. Success is determined solely by the diagonalization constraint: specifically, all off-diagonal elements (\(i \neq j\)) must satisfy \(|p(U_k)_{ij}| < \epsilon_{\text{elem}}\), where \(\epsilon_{\text{elem}} = 5 \times 10^{-4}\) in the simulation. The specific value of this threshold does not affect program acceptance, as successfully synthesized programs typically drive residuals down to machine precision (\(\sim 10^{-16}\)). Meanwhile, the diagonal phase values in \(D_k\) are not prescribed, but emerge naturally from the decomposition process.

To solve this diagonalization task, we define a set of \(N(N-1)\) MZI primitives \(\{R_{ij}, L_{ij}\}\) for \(i > j\), with subscripts denoting the target element position at row \(i\), column \(j\). Each primitive eliminates one matrix element. The decomposition strategy exploits the unitarity constraint \(U^\dagger U = I\), where systematic elimination of lower-triangular elements automatically enforces upper-triangular nullity, yielding a diagonal matrix. Each MZI primitive corresponds to a beam splitter transformation acting on adjacent modes, eliminating matrix element \(u_{ij}\) via adaptive phase shifts. For \(R_{ij}\)-type primitives, the transformation eliminates element \(u_{ij}\) by mixing columns \(j\) and \(j+1\) via right multiplication. From the nullification condition \((U R_{ij})_{ij} = 0\), the required phase \(\theta_{ij}\) and mixing angle \(\omega_{ij}\) are analytically derived: the ratio \(r = u_{ij}/u_{i,j+1}\) is computed, with \(\theta_{ij} = \arg(r)\) compensating the phase difference and \(\omega_{ij} = \arctan(|r|)\) determining the mixing ratio. This constructs an \(N \times N\) transformation matrix by embedding the \(2 \times 2\) MZI block \[M_{R,ij} = \begin{pmatrix} e^{-i\theta_{ij}}\cos\omega_{ij} & e^{-i\theta_{ij}}\sin\omega_{ij} \\ -\sin\omega_{ij} & \cos\omega_{ij} \end{pmatrix}\] at positions \((j, j+1)\) within an identity matrix, applying it via right multiplication \(U' = U R_{ij}\). \(L_{ij}\)-type primitives similarly mix rows \(i\) and \(i-1\) via left multiplication, deriving parameters from condition \((L_{ij} U)_{ij} = 0\) and embedding the \(2 \times 2\) MZI block \[M_{L,ij} = \begin{pmatrix} e^{i\theta_{ij}}\cos\omega_{ij} & -\sin\omega_{ij} \\ e^{i\theta_{ij}}\sin\omega_{ij} & \cos\omega_{ij} \end{pmatrix}\] at positions \((i-1, i)\) within an identity matrix, applying it via left multiplication \(U' = L_{ij} U\). The complete \(N \times N\) matrix forms of \(R_{ij}\) and \(L_{ij}\) primitives are provided in Supplemental Material (SM). Each primitive takes a complex matrix as input and returns a transformed matrix, so primitives can be freely chained into sequences: e.g., a \(3 \times 3\) decomposition might be (lambda (R21 (L10 (L20 $0)))) [applied in sequence to positions \((2,0)\), \((1,0)\), and \((2,1)\)], where $0 denotes the input matrix and operations are applied from innermost to outermost.

Given the task definition and primitive operations, we employ program synthesis to autonomously discover effective decomposition strategies. We adopt the DreamCoder framework [24], which operates through an iterative wake-sleep cycle, as illustrated in Fig. 1(b). During the wake phase, the system searches over candidate MZI operation sequences using the current primitive library and evaluates them against decomposition tasks. During the sleep phase, it identifies recurring patterns across successful solutions and abstracts them into new reusable primitives, expanding the library. The updated library is then carried into the next wake phase, where the richer primitive vocabulary enables more efficient search over deeper program structures [26]. We extend DreamCoder to the complex-valued matrix decomposition domain; implementation details are provided in SM. The MZI primitives \(\{R_{ij}, L_{ij}\}\) encode physically meaningful beam-splitter transformations with analytically derived phase and rotation angles. This design can be generalized to a broad class of quantum information processing tasks involving matrix operations, such as quantum circuit synthesis [27], [28].

Verification of the synthesized programs depends on their target application. Universal decompositions are evaluated on multiple randomly generated unitary matrices, whereas structure-specific programs (e.g., for sparse matrices) are tested directly on their original inputs. In all cases, successfully verified programs achieve off-diagonal residuals at machine precision (\(\sim 10^{-16}\)). Furthermore, we validate the cross-dimensional generalization of universal programs by applying the learned patterns to larger matrices, confirming their ability to scale to arbitrary \(N\).

Having established the framework, we now demonstrate its effectiveness on unitary decomposition tasks. To assess the system’s capacity for autonomous algorithmic discovery, we employ independent learning where the system constructs complete decompositions from scratch without prior knowledge of classical schemes. This approach ensures discovered strategies emerge from the search process itself rather than reflecting biases from known methods. We first investigate the \(5 \times 5\) case using independent learning within tractable search spaces. After nine wake-sleep iterations, the system autonomously discovered multiple distinct decomposition strategies, all achieving the theoretical minimum of 10 MZI operations. Figure 2 visualizes four random decomposition patterns, where the orange paths trace the sequence of operations through the matrix. Each pattern shows a different execution ordering that systematically eliminates lower-triangular elements, confirming the system learned flexible decomposition strategies rather than a single fixed pattern. These strategies exhibit novel elimination orderings distinct from both Reck and Clements schemes, as presented in Fig. 2. The learned decomposition strategies are not tied to specific unitary instances but apply universally to arbitrary unitary matrices, functioning as generalizable algorithmic rules. To the best of our knowledge, such strategies have not been reported previously. Additional representative examples are provided in Table S1 in SM.

Among the discovered programs, the one ranked highest by the system, which favors solutions that are both correct and structurally simple (i.e., expressible more compactly in terms of the learned library operations) [29], [30] is (lambda (R10 (R21 (R20 (R32 (R31 (R30 (R43 (R42 (R41 (L40 $0))))))))))). This program corresponds to the matrix equation: \[D = L_{40} \, U \, R_{41} R_{42} R_{43} R_{30} R_{31} R_{32} R_{20} R_{21} R_{10}\, ,\] where \(D\) is the resulting diagonal matrix. Inverting this relation yields the unitary decomposition \[U = L_{40}^{-1} D \, R_{10}^{-1} R_{21}^{-1} R_{20}^{-1} R_{32}^{-1} R_{31}^{-1} R_{30}^{-1} R_{43}^{-1} R_{42}^{-1} R_{41}^{-1}\, .\] This program directly maps to a photonic circuit implementation, where the five optical waveguides correspond to matrix modes and each MZI block executes one primitive operation. The circuit processes light through a cascade of 10 MZI stages, achieving the theoretical minimum count. This demonstrates that the system has learned a flexible decomposition strategy, revealing that the systematic elimination of lower-triangular elements admits multiple valid operation sequences.

Critically, the learned patterns exhibit cross-dimensional reasoning. Analysis of discovered programs reveals dimension-agnostic principles: systematic elimination strategies that remain valid across matrix scales. For instance, algorithmic abstractions derived from \(5 \times 5\) matrices exhibit compositional generalization: they can be directly applied to higher-dimensional unitary matrices without retraining or modification. As an illustration, consider Pattern 4 shown in Fig. 2, which implements a row-pair interleaving strategy. Starting from the largest row index \(i = N-1\) and decrementing by two each step, the algorithm processes rows in pairs: it first emits primitive \(\mathtt{R}_{i,0}\), then interleaves primitives from two adjacent rows as \(\mathtt{R}_{i,j} \to \mathtt{R}_{i-1,\,j-1}\) for \(j = 1, \ldots, i-1\). When a single row remains (\(i=1\)), \(\mathtt{R}_{1,0}\) closes the sequence. This fixed rule, discovered on \(5 \times 5\) matrices, applies unchanged at any \(N\) and enumerates all \(N(N-1)/2\) operations exactly. Full pseudocode and generalization results are given in SM. We validated this cross-dimensional generalization by testing programs discovered on \(5 \times 5\) training instances against progressively larger matrices, including \(8 \times 8\), \(16 \times 16\), \(32 \times 32\), and \(64 \times 64\) dimensions. All test cases achieve successful diagonalization at machine precision (\(\sim 10^{-16}\)), confirming that the learned decomposition strategy scales correctly to arbitrary dimensions.

While our learned patterns exhibit dimension-agnostic principles that allow \(5 \times 5\) strategies to generalize, direct training on \(6 \times 6\) matrices remains essential to uncover richer decomposition orderings inaccessible from smaller examples. However, this objective confronts a fundamental bottleneck: the search space grows exponentially with matrix dimension, making the direct enumeration of deep 15-MZI sequences intractable. As matrix dimension grows, both the primitive vocabulary and program depth increase, causing the search space to explode combinatorially. To address this challenge, we employ curriculum learning [31], [32] by decomposing the problem into subtasks \(\{T_3, T_2, T_1, T_0\}\). This progression allows the model to first master short subsequences in early stages (\(T_3\)), rapidly abstracting reusable patterns into library primitives via the sleep phase. Consequently, these consolidated abstractions effectively reduce the complexity for later stages (\(T_0\)), enabling the synthesis of deeply nested programs. Specifically, our results demonstrate that after just seven wake-sleep iterations, the system not only converged to the theoretical minimum depth but also expanded its library from 30 to 65 hierarchical abstractions. The system successfully rediscovered the Reck decomposition and identified novel patterns distinct from \(5 \times 5\) independent learning. One such pattern can be expressed by the diagonalization relation \(D = L_{5,1}L_{4,0}L_{5,0}\,U\, R_{5,2}R_{5,3}R_{5,4}R_{4,1}R_{4,2}R_{3,0}R_{4,3} R_{3,1}R_{2,0}R_{3,2}R_{2,1}R_{1,0}\,\). This demonstrates that curriculum learning effectively navigates the expanded search space to uncover both classical and genuinely new decomposition strategies.

Beyond universal decompositions, the system also discovers structure-aware reductions for specific matrix classes, without requiring any specialized analytical rules. We identify two broad categories of such structure-aware reductions.

The first category comprises analytically structured matrices, with Householder matrices serving as a prime example. Defined as \(H = I - 2vv^\dagger\) (\(I\) is the identity matrix, \(v \in \mathbb{C}^N\), \(\|v\|=1\)), Householder matrices are unitary reflectors widely used in mathematics and physics for stable unitary transformations, including tridiagonalization of Hamiltonians and eigenvalue calculations in quantum and many-body systems [33]–[35]. By exploiting these inherent structural constraints, our system autonomously identifies a dimension-independent decomposition rule that strictly bounds the total MZI count to \(2N-3\). The recursive logic is as follows: extending to \(N+1\) dimensions requires prepending two operations, \(\mathtt{L}_{N,0}\) followed by \(\mathtt{R}_{N,N-1}\), to the \(N\)-dimensional sequence. At \(N=6\), the corresponding diagonalization relation is \(D = L_{4,0}L_{5,0}\,H\,R_{5,4}R_{4,3}R_{3,0}R_{3,1}R_{3,2}R_{2,1}R_{1,0}\) (9 operations \(= 2\times6-3\)). This strategy reduces the hardware scaling from the universal quadratic bound of \(N(N-1)/2\) to a linear complexity. When evaluated at \(N=64\), this generalized rule perfectly reconstructs the unitary using only 125 MZIs instead of 2016. This represents a substantial 93.8% reduction in hardware components, verified at machine precision without any retraining.

The second category involves sparse matrices. In our framework, an arbitrary input \(W\) is decomposed via SVD (\(W = U\Sigma V^\dagger\)), where the unitary factors \(U\) and \(V^\dagger\) inherit the sparsity of \(W\) and are further decomposed into MZI sequences. Our system exploits this inherited structure to discover highly efficient MZI sequences. We validate this on \(N=6\) matrices using two sparsity models: an NNZ-controlled setting (NNZ: number of nonzero elements) with a fixed number of nonzero elements in \(W\), and a Bernoulli model where each element of \(W\) is independently zeroed with probability \(p\) [36], [37]. As shown in Fig. 3, the system achieves significant reductions relative to the universal theoretical bound \(N(N-1)/2\), and these reductions generally increase with matrix sparsity. Under the NNZ-controlled model, reductions begin around \(39\%\) sparsity and reach up to \(31\%\) at \(89\%\) sparsity. Under the Bernoulli model, reductions emerge around \(40\%\) sparsity and grow more steeply at high sparsity, reaching up to \(38\%\) at \(95\%\) sparsity. This overall trend is consistent across both models, confirming that the efficiency gains reflect intrinsic structural properties of sparse unitaries rather than artifacts of specific matrix generation. For arbitrary unstructured sparse matrices, the discovered patterns are instance-specific, as there is no well-defined correspondence between sparsity patterns at different dimensions.

These MZI reductions can provide practical hardware benefits: fewer components can help reduce optical insertion loss [23], calibration complexity [38], and power consumption [39]. Although synthesized topologies may be irregular, modern photonic design flows increasingly support automated placement and routing under physical constraints such as bend radius and waveguide crossings [40], [41].

In conclusion, our results show that program synthesis can uncover generalizable algorithmic structures for universal unitaries—going beyond known schemes to reveal decomposition principles and resource optimizations that classical Reck and Clements designs cannot capture. Built on a complex-valued matrix decomposition domain, our framework combines symbolic reasoning with physically grounded MZI primitives to systematically generate and verify decomposition programs, providing a scalable foundation that can be extended to broader quantum information tasks. These automatically generated constructions encode systematic construction methods, demonstrating the potential of AI-driven algorithm discovery to reveal fundamental principles and offering a powerful new approach for advancing photonic computing and quantum information processing.

Acknowledgments—This work is supported by the National Key Research and Development Program of China under Contract No. 2024YFA1610900 and the National Natural Science Foundation of China (NSFC) under Contract No. 12447106, No. 12541501, No. 12547102 and No. 12450404. D. D. H. also acknowledges the support of NSFC under No. 11875133 and No. 11075057, the National Key Research and Development Program of China under Contract No. 2018YFB2101302. The computations in this research were performed using the CFFF platform of Fudan University.

In this Supplemental Material, we present more details on the framework implementation, curriculum learning strategy, and synthesized programs.

FRAMEWORK AND IMPLEMENTATION DETAILS↩︎

This section provides a detailed exposition of the DreamCoder program synthesis framework [24] and our matrix decomposition domain implementation.

The DreamCoder framework operates through alternating wake and sleep phases. During the wake phase, the system performs enumerative search over the space of programs expressible in the current grammar. Programs are enumerated in order of decreasing prior probability under a learned probabilistic context-free grammar (PCFG). This grammar assigns higher probabilities to operation sequences that have proven useful in previous iterations, guiding the search toward efficient decompositions. Programs are then evaluated against training tasks using a posterior-based scoring function. For each task, the search identifies programs that maximize the posterior probability \(P(p|T) \propto P(T|p)P(p)\), where \(P(T|p)\) measures how well program \(p\) solves task \(T\) and \(P(p)\) is the prior probability assigned by the grammar (favoring programs with shorter description length under the current library). The wake phase produces a corpus of successful programs paired with the tasks they solve.

During the sleep phase, the system analyzes the corpus of programs discovered during wake to identify frequently recurring substructures. It employs a compression-based abstraction algorithm that searches for program fragments whose factorization into reusable library functions would maximize the overall description length reduction across the entire corpus. Given a set of programs \(\{p_1, ..., p_n\}\), the algorithm identifies a new abstraction \(f\) and rewritten programs \(\{p'_1, ..., p'_n\}\) that minimize the total cost \(\text{cost}(f) + \sum_i \text{cost}(p'_i)\), where cost is measured by program size under the grammar. Successfully abstracted fragments are added to the library as new primitives, enabling the next wake phase to build upon these higher-level building blocks. The updated library is then used in the next wake phase, closing the learning loop and enabling progressively more efficient search over successive iterations.

Program evaluation with hard constraints. For matrix decomposition tasks, we extend DreamCoder’s evaluation mechanism from exact matching to constraint-based validation. Each task consists of unitary input matrices \(\{U_k\}\), with success defined by transforming each input into diagonal form. A program \(p\) succeeds when all off-diagonal elements satisfy \(|p(U_k)_{ij}| < \epsilon_{\text{elem}} = 5 \times 10^{-4}\) for all \(k\) and all \(i \neq j\). In practice, successfully synthesized programs drive all off-diagonal elements to near machine precision (\(\sim 10^{-16}\)), so the specific value of \(\epsilon_{\text{elem}}\) does not affect which programs are accepted. Programs passing this check are then ranked by posterior probability \(P(p|T) \propto P(T|p)P(p)\), where the likelihood \(P(T|p)\) is measured by the average off-diagonal Frobenius norm \(\frac{1}{K}\sum_k \|p(U_k) - \mathrm{diag}(p(U_k))\|_F\), serving as a continuous score reflecting decomposition quality, and \(P(p)\) is the prior from the learned grammar favoring shorter programs. In the tables of discovered programs, we report the log-posterior \(\log P(p|T)\): a less negative value indicates a program that is both more accurate and has a shorter description length under the learned library grammar. The program with the shortest description length under the learned library that correctly solves all tasks receives the highest (least negative) log-posterior and is selected as the representative solution.

The matrix decomposition domain extends DreamCoder through a two-component software interface. The Python frontend defines the MZI primitives, generates training tasks, and coordinates the learning loop. The OCaml backend performs the combinatorial search over program sequences, leveraging its efficient multicore enumeration capability. The two components communicate by passing matrix data and candidate programs as structured messages in JSON format (a lightweight, human-readable data interchange format). The Python component defines matrix decomposition tasks, constructs a parameterizable primitive library for \(N \times N\) matrices, and supports curriculum-style task generation, while the OCaml backend implements the corresponding primitive semantics for matrix evaluation and verification and carries out program search based on these task and primitive specifications.

For matrix element elimination, two primitive families are defined: \(R_{ij}\) and \(L_{ij}\) for \(i > j\). The \(R_{ij}\) primitive targets element \(u_{ij}\) via right multiplication (column operations), while \(L_{ij}\) uses left multiplication (row operations). Both are parameterized by row and column indices through partial application.

The \(R_{ij}\) primitive constructs its transformation from element \(u_{ij}\) and its right neighbor \(u_{i,j+1}\), computing the complex ratio, phase angle, and rotation angle: \[\begin{align} r &= u_{ij}/u_{i,j+1}, \quad \theta_{ij} = \arg(r), \quad \omega_{ij} = \arctan(|r|). \end{align}\] The \(2 \times 2\) MZI transformation matrix is: \[M_{R,ij} = \begin{bmatrix} e^{-i\theta_{ij}}\cos\omega_{ij} & e^{-i\theta_{ij}}\sin\omega_{ij} \\ -\sin\omega_{ij} & \cos\omega_{ij} \end{bmatrix}.\] This submatrix is embedded in an \(N \times N\) identity matrix, replacing columns \((j, j+1)\): \[G_{R_{ij}} = \begin{bmatrix} 1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\ 0 & \cdots & 0 & e^{-i\theta_{ij}}\cos\omega_{ij} & e^{-i\theta_{ij}}\sin\omega_{ij} & \cdots & 0 \\ 0 & \cdots & 0 & -\sin\omega_{ij} & \cos\omega_{ij} & \cdots & 0 \\ 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & \cdots & 1 \end{bmatrix}\] Right multiplication \(U' = U \cdot G_{R_{ij}}\) mixes columns \(j\) and \(j+1\), zeroing \(u'_{ij}\) while preserving unitarity \(U'^{\dagger}U' = I\).

The \(L_{ij}\) primitive operates analogously via left multiplication, using element \(u_{ij}\) and its upper neighbor \(u_{i-1,j}\): \[\begin{align} r &= u_{ij}/u_{i-1,j}, \quad \theta_{ij} = \arg(r), \quad \omega_{ij} = \arctan(-|r|). \end{align}\] The \(2 \times 2\) MZI transformation matrix is: \[M_{L,ij} = \begin{bmatrix} e^{i\theta_{ij}}\cos\omega_{ij} & -\sin\omega_{ij} \\ e^{i\theta_{ij}}\sin\omega_{ij} & \cos\omega_{ij} \end{bmatrix}.\] This submatrix is embedded in an \(N \times N\) identity matrix, replacing rows \((i-1, i)\): \[G_{L_{ij}} = \begin{bmatrix} 1 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 1 & 0 & 0 & \cdots & 0 \\ 0 & \cdots & 0 & e^{i\theta_{ij}}\cos\omega_{ij} & -\sin\omega_{ij} & \cdots & 0 \\ 0 & \cdots & 0 & e^{i\theta_{ij}}\sin\omega_{ij} & \cos\omega_{ij} & \cdots & 0 \\ 0 & \cdots & 0 & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & \cdots & 1 \end{bmatrix}.\] Left multiplication \(U' = G_{L_{ij}} \cdot U\) mixes rows \(i-1\) and \(i\), zeroing \(u'_{ij}\) while preserving unitarity \(U'^{\dagger}U' = I\).

All MZI primitives share the uniform type signature \(\text{tmatrix} \to \text{tmatrix}\), so they can be freely chained into sequences. We present complete decompositions as ordered products of MZI primitives: each \(L_{ij}\) left-multiplies the current matrix and each \(R_{ij}\) right-multiplies it, so an ordered sequence such as \(L_{4,1}\) followed by \(R_{2,1}\) is written as \(D=L_{4,1}\,U\,R_{2,1}\). The enumerative search explores different orderings and combinations of these primitives, constrained by type checking and guided by the learned grammar probabilities.

Task generation creates unitary input matrices \(\{U_k\}\) via QR decomposition of complex Gaussian matrices, with each task requiring the input to be transformed into diagonal form. Tasks support curriculum learning through a factory mechanism that applies a specified sequence of input transformation steps before presenting the matrix to the solver, enabling progressive learning of decomposition sub-problems.

The search engine evaluates each candidate MZI sequence by checking two criteria: (1) whether all off-diagonal elements satisfy \(|u_{ij}| < \epsilon_{\text{elem}}\) for \(i \neq j\) (the hard diagonalization constraint), and (2) its posterior score \(P(p|T)\) as defined above. Programs failing the hard constraint are immediately discarded. For programs passing it, an accuracy score is computed and combined with a simplicity prior (favoring shorter sequences) to form the posterior probability \(P(p|T) \propto P(T|p)P(p)\). The search enumerates candidate sequences in order of decreasing prior probability under the learned grammar, so that concise and accurate decompositions are found first. Search depth limits are set according to matrix size to control the enumeration space.

CURRICULUM LEARNING STRATEGY↩︎

For \(5 \times 5\) matrices, independent learning successfully discovers decompositions with optimal MZI counts by constructing complete solutions from scratch. However, scaling to \(6 \times 6\) and larger dimensions benefits from curriculum learning, a training strategy that progressively increases task difficulty. The curriculum consists of subtasks \(\{T_k, T_{k-1}, ..., T_0\}\) of increasing complexity, where task \(T_i\) provides partially decomposed matrices with the first \(i\) steps of a reference elimination sequence applied purely to constrain the initial search space. For instance, \(T_3\) presents matrices where the first 3 steps have already been executed, requiring the system to discover the remaining sequence. Early curriculum stages (large \(i\)) present tractable search spaces with more pre-applied steps and fewer remaining steps, enabling rapid discovery of short decomposition subsequences. The sleep phase abstracts these solutions into reusable library functions, which later stages leverage to construct longer programs through composition, avoiding exhaustive enumeration of deeply nested primitive sequences. This bootstrapping process is essential for \(6 \times 6\) and larger matrices, where the lengthy search space proves intractable without hierarchical abstraction. Tasks are generated by applying a specified sequence of MZI transformations to randomly generated unitary matrices before presenting them to the solver. In our experiments, we use a curriculum depth of 3, training on tasks \(\{T_3, T_2, T_1, T_0\}\) sequentially, progressing from simpler subtasks (fewer remaining steps) to more complex ones (more remaining steps). Deeper curricula can be used to discover decompositions for larger matrices directly, though this is not necessary since patterns learned from smaller matrices generalize compositionally. The wake-sleep cycle iterates within each curriculum stage until convergence criteria are met (e.g., solving all training instances or reaching iteration limits or timeout), then advances to the next stage.

Results for \(6 \times 6\) Matrices↩︎

Scaling to \(6 \times 6\) matrices required curriculum learning, progressively training on subtasks of increasing difficulty to bootstrap library construction. After 7 iterations, the system discovered programs achieving the theoretical 15 MZI minimum. Table 2 lists five representative optimal decompositions discovered by the system.

The evolution of the learned library reveals how the system builds hierarchical abstractions. Starting from 30 hand-defined MZI primitives, the grammar expanded to 65 primitives over 10 iterations via wake-sleep cycles. Fig. 5 shows rapid growth in early iterations with 8 and 7 new primitives in iterations 1-2, followed by stabilization with 1-2 new primitives per iteration after iteration 7, reflecting the transition from exploratory pattern discovery to consolidation of a compact, reusable vocabulary for expressing decompositions.

This growth pattern explains curriculum learning’s effectiveness: richer libraries enable reuse of extracted transformations, constructing deeper decomposition programs with less redundancy and revealing hierarchical organization in the learned vocabulary.

REPRESENTATIVE SYNTHESIZED PROGRAMS↩︎

Tables 1 and 2 present representative decomposition programs discovered by the DreamCoder framework under independent learning and curriculum learning configurations, respectively. All programs are validated to successfully reduce randomly generated unitary matrices to diagonal form with all off-diagonal elements satisfying \(|D_{ij}| < \epsilon_{\text{elem}}\) for \(i \neq j\).

10pt

10pt

STRUCTURE-AWARE DECOMPOSITIONS↩︎

When the input matrix carries structural constraints—such as Householder form or sparsity-induced structural constraints—the synthesized programs automatically exploit them, discovering far fewer MZI operations than the general-case bound. We present two families of examples.

Householder reflectors. A Householder matrix \(H = I - 2\mathbf{v}\mathbf{v}^\dagger\) (where \(\mathbf{v} \in \mathbb{C}^N\), \(\|\mathbf{v}\|=1\)) can be decomposed using only \(2N-3\) MZI operations, compared to the \(N(N-1)/2\) required for a general unitary. The synthesis system discovers programs achieving this reduction without being provided with the analytical formula of Householder reflectors.

For \(N=5\), every Householder task is solved in exactly 7 operations (\(=2\times5-3\)), with the representative program \(\mathtt{L}_{4,0},\mathtt{R}_{4,3},\mathtt{R}_{3,0},\mathtt{R}_{3,1},\mathtt{R}_{3,2},\mathtt{R}_{2,1},\mathtt{R}_{1,0}\). For \(N=6\), every task is solved in exactly 9 operations (\(=2\times6-3\)), with the representative program \(\mathtt{L}_{5,0},\mathtt{R}_{5,4},\mathtt{L}_{4,0},\mathtt{R}_{4,3},\mathtt{R}_{3,0},\mathtt{R}_{3,1},\mathtt{R}_{3,2},\mathtt{R}_{2,1},\mathtt{R}_{1,0}\). The \(N\to N+1\) extension prepends \(\mathtt{L}_{N,0}\) and \(\mathtt{R}_{N,N-1}\) as the innermost pair, consistent with \(2N-3\) scaling.

Sparse SVD decompositions. When the source matrix \(W\) is sparse, its SVD \(W=U\Sigma V^\dagger\) yields unitary factors \(U\) and \(V^\dagger\) with internal structure that reduces the number of MZI operations needed for decomposition. The system discovers and exploits this structure automatically from examples, without being provided with information about the sparsity level or structure of the input.

NNZ-controlled sparsity. We consider a representative example with source \(W\): \(81\%\) sparse, \(7/36\) nonzero entries controlled by NNZ. The input \(U\) matrix is shown below to three decimal places for readability, with all computations performed at full floating-point precision. \[U = \begin{pmatrix} {-}0.686{-}0.039i & {-}0.696{-}0.039i & 0.000{+}0.000i & 0.204{+}0.011i & 0.000{+}0.000i & 0.000{+}0.000i \\ 0.000{+}0.000i & 0.000{+}0.000i & 0.444{-}0.896i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i \\ {-}0.184{-}0.078i & 0.416{+}0.177i & 0.000{+}0.000i & 0.800{+}0.340i & 0.000{+}0.000i & 0.000{+}0.000i \\ 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.896{-}0.443i & 0.000{+}0.000i \\ 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.251{+}0.968i \\ {-}0.400{-}0.572i & 0.319{+}0.456i & 0.000{+}0.000i & {-}0.258{-}0.369i & 0.000{+}0.000i & 0.000{+}0.000i \end{pmatrix}.\] The corresponding diagonalization relation is \[D = L_{3,0}L_{4,0}L_{5,0}\,U\,R_{3,1}R_{3,2}R_{2,0}R_{2,1}R_{1,0},\] achieving a \(47\%\) reduction from the universal theoretical bound of 15 MZIs, verified at machine precision (\(\sim10^{-16}\)).

Bernoulli sparsity. We consider a representative example with source \(W\): \(86.1\%\) sparse, \(5/36\) nonzero; each entry independently zeroed with probability \(80\%\). The input \(U\) matrix is \[U = \begin{pmatrix} 0.000{+}0.000i & 0.000{+}0.000i & 1.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i \\ 0.336{+}0.000i & 0.942{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i \\ {-}0.532{-}0.643i & 0.190{+}0.229i & 0.000{+}0.000i & 0.173{-}0.430i & 0.000{+}0.000i & 0.000{+}0.000i \\ 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & {-}1.000{+}0.015i & 0.000{+}0.000i \\ 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.000{+}0.000i & 0.970{+}0.243i \\ 0.111{-}0.422i & {-}0.040{+}0.151i & 0.000{+}0.000i & {-}0.861{+}0.212i & 0.000{+}0.000i & 0.000{+}0.000i \end{pmatrix}.\] The corresponding diagonalization relation is \[D = L_{4,0}L_{5,0}\,U\,R_{3,1}R_{1,0}R_{2,1}R_{3,2},\] achieving a 60% reduction from the universal theoretical bound of 15 MZIs, verified at machine precision (\(\sim10^{-16}\)).

Across both structure classes, the synthesis system discovers these reductions automatically, without being provided with the analytical forms or sparsity properties of the input matrices, demonstrating that structure-aware efficiency emerges naturally from the program synthesis process itself.

CROSS-DIMENSIONAL GENERALIZATION: EXTRACTING DIMENSION-INDEPENDENT CONSTRUCTION RULES FROM SYNTHESIZED PROGRAMS↩︎

A critical capability of the discovered decomposition programs is their ability to generalize across matrix dimensions without retraining. The key mechanism is that the synthesis system does not memorize a fixed operation sequence at low \(N\), but extracts a dimension-independent construction rule. Once such a rule is identified, substituting any target \(N\) directly generates the corresponding execution sequence without retraining. We illustrate this with two concrete examples before presenting the detailed pseudocode for Pattern 4.

How low-dimensional patterns are expanded to high dimensions.

Example 1 — Sequential row elimination (Table 1, Program 1): the \(N=5\) sequence is \(\mathtt{L}_{4,0}, \mathtt{R}_{4,1}, \mathtt{R}_{4,2}, \mathtt{R}_{4,3}, \mathtt{R}_{3,0}, \mathtt{R}_{3,1}, \mathtt{R}_{3,2}, \mathtt{R}_{2,0}, \mathtt{R}_{2,1}, \mathtt{R}_{1,0}\). The rule is: for the top row (\(i=N-1\)), apply \(\mathtt{L}_{i,0}\) followed by \(\mathtt{R}_{i,1}, \ldots, \mathtt{R}_{i,i-1}\). For all lower rows (\(i = N-2\) down to \(1\)), apply \(\mathtt{R}_{i,0}, \mathtt{R}_{i,1}, \ldots, \mathtt{R}_{i,i-1}\) in order. Each row is fully eliminated before moving to the next. The operation count follows directly: \(\sum_{i=1}^{N-1} i = N(N-1)/2\), yielding 2016 operations at \(N=64\).

Example 2 — Two-row leading elimination (Table 1, Program 2): the \(N=5\) sequence is \(\mathtt{L}_{4,0}, \mathtt{L}_{3,0}, \mathtt{R}_{4,1}, \mathtt{R}_{4,2}, \mathtt{R}_{4,3}, \mathtt{R}_{3,1}, \mathtt{R}_{3,2}, \mathtt{R}_{2,0}, \mathtt{R}_{2,1}, \mathtt{R}_{1,0}\). The rule is: first apply \(\mathtt{L}_{N-1,0}\) and \(\mathtt{L}_{N-2,0}\); then complete row \(N-1\) with \(\mathtt{R}_{N-1,1}, \ldots, \mathtt{R}_{N-1,N-2}\); then row \(N-2\) with \(\mathtt{R}_{N-2,1}, \ldots, \mathtt{R}_{N-2,N-3}\) (starting from \(j=1\) since \(j=0\) was handled by \(\mathtt{L}_{N-2,0}\)); then all remaining rows \(i = N-3\) down to \(1\) with \(\mathtt{R}_{i,0}, \mathtt{R}_{i,1}, \ldots, \mathtt{R}_{i,i-1}\). The same operation count \(N(N-1)/2\) is preserved.

We verified Programs 1, 2, and 3 from Table 1 across \(N = 5, 6, 7, 8, 16, 32, 64\); all achieve machine-precision accuracy (\(\sim 10^{-16}\)) with 100% success rate, as shown in Table 3. This confirms that the decompositions are mathematically exact: the residual error reflects only floating-point rounding and does not grow systematically with \(N\).

We next present the detailed algorithmic procedure for Pattern 4 (row-pair interleaving, Fig. 2; distinct from the programs listed in Table 1) as a representative example, with explicit pseudocode for constructing the decomposition sequence at arbitrary \(N\).

The learned decomposition strategies encode systematic elimination orderings that can be expressed through regular patterns in the primitive sequence. Specifically, Pattern 4 implements a row-pair interleaving strategy: starting from the largest row index \(i = N-1\) and decrementing by two each step, the algorithm processes rows in pairs by first emitting primitive \(\mathtt{R}_{i,0}\), then interleaving primitives from two adjacent rows as \(\mathtt{R}_{i,j} \to \mathtt{R}_{i-1,\,j-1}\) for \(j = 1, \ldots, i-1\). When a single row remains (\(i=1\)), \(\mathtt{R}_{1,0}\) closes the sequence. This fixed rule enumerates all \(N(N-1)/2\) operations exactly and generalizes unchanged to arbitrary dimensions.

The following pseudocode formalizes this generalization procedure for constructing lambda expressions representing \(N \times N\) unitary decompositions:

    Algorithm: GenerateDimensionAgnosticDecomposition(N)
    Input: Matrix dimension N
    Output: Lambda expression for NxN unitary decomposition
    
    1: Initialize operation_sequence := []
    2: i := N - 1  // Start from the largest row index
    3:
    4: while i >= 1 do
    5:     if i = 1 then
    6:         // Base case: single remaining row
    7:         operation_sequence.append((1, 0))
    8:         break
    9:     end if
    10:
    11:    // Process row pair (i, i-1) with interleaving
    12:    // First emit the leading element of row i
    13:    operation_sequence.append((i, 0))
    14:
    15:    // Then interleave elements from rows i and i-1
    16:    for j = 1 to i-1 do
    17:        operation_sequence.append((i, j))      // Row i element
    18:        operation_sequence.append((i-1, j-1))  // Row i-1 element
    19:    end for
    20:
    21:    i := i - 2  // Move to next row pair
    22: end while
    23:
    24: // Construct lambda expression from innermost to outermost
    25: lambda_expr := "$0"  // Input matrix placeholder
    26:
    27: for each (row, col) in operation_sequence do
    28:     // Wrap with primitive R_{row,col}
    29:     lambda_expr := "(R_{" + row + col + "} " + lambda_expr + ")"
    30: end for
    31:
    32: lambda_expr := "(lambda " + lambda_expr + ")"
    33: return lambda_expr

This construction ensures that the total number of primitives equals \(N(N-1)/2\), the theoretical minimum for universal unitary decomposition. The dimension-agnostic nature of this pattern arises from its reliance on structural invariants—specifically, the row-pair interleaving rule and the fixed progression from larger to smaller row indices—rather than dimension-specific heuristics. Consequently, the same algorithmic template applies to matrices of any size, enabling automatic scaling from training dimensions (e.g., \(5 \times 5\)) to deployment dimensions (e.g., \(64 \times 64\) or larger).

The generalization validation protocol tests learned programs on matrices of progressively increasing dimensions: \(8 \times 8\), \(16 \times 16\), \(32 \times 32\), and \(64 \times 64\). For each test dimension, we verify that the compositionally extended program achieves successful diagonalization on random unitary matrices with all off-diagonal elements satisfying \(|p(U_k)_{ij}| < \epsilon_{\text{elem}}\) for all \(k\) and all \(i \neq j\), with residuals at machine precision (\(\sim 10^{-16}\)). Successful verification across all test dimensions confirms that the discovered decomposition programs encode dimension-agnostic decomposition strategies rather than memorized dimension-specific patterns, establishing the foundation for scalable photonic circuit synthesis.

References↩︎