Paper deep dive
Program-Synthesis-Driven Autodesign of Universal Unitary Operators
Yifei Zhang, Dong Chen, Fan Wang, Wenrui Zhang, Yan Chen, Dingding Han, Jianmin Yuan, Xiangjin Kong, Yu-Gang Ma
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 95%
Last extracted: 7/18/2026, 11:57:34 AM
Summary
The paper demonstrates that AI-driven program synthesis, specifically using an extended DreamCoder framework, can autonomously discover strategies for decomposing unitary matrices into Mach-Zehnder interferometers (MZIs) for photonic networks. The system discovers dimension-agnostic universal decomposition rules that generalize across matrix sizes (e.g., from 5x5 to 64x64) and achieve the theoretical minimum of N(N-1)/2 MZIs, distinct from classical Reck and Clements architectures. Additionally, the system exploits matrix structure to reduce MZI counts below the universal bound, achieving linear scaling for Householder matrices (2N-3 MZIs) and significant reductions for sparse matrices, demonstrating a scalable paradigm for automated algorithm discovery and photonic circuit design.
Entities (8)
Relation Signals (8)
Program Synthesis → discovers → decomposition strategies
confidence 97% · We demonstrate that AI-driven program synthesis can autonomously discover fundamental strategies for decomposing unitary matrices...
Discovered Strategies → achieves → minimal N(N-1)/2 MZIs
confidence 96% · ...generates decomposition programs achieving the minimal N(N-1)/2 Mach-Zehnder interferometers...
Householder Matrices → requires → 2N-3 MZIs
confidence 96% · For Householder matrices, it discovers a dimension-independent rule that requires only 2N-3 MZIs.
Discovered Strategies → distinctfrom → Clements Architecture
confidence 95% · ...distinct from both Reck and Clements architectures.
Discovered Strategies → distinctfrom → Reck Architecture
confidence 95% · ...distinct from both Reck and Clements architectures.
DreamCoder → extends → complex-valued linear algebra
confidence 95% · By extending DreamCoder to complex-valued linear algebra, the system generates decomposition programs...
Discovered Strategies → generalizesto → higher dimensions
confidence 95% · strategies discovered for 5x5 matrices generalize to higher dimensions such as 64x64.
Sparse Matrices → exhibits →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We demonstrate that AI-driven program synthesis can autonomously discover fundamental strategies for decomposing unitary matrices in photonic networks. By extending DreamCoder to complex-valued linear algebra, the system generates decomposition programs achieving the minimal $N(N-1)/2$ Mach-Zehnder interferometers, distinct from both Reck and Clements architectures. Learned programs encode dimension-agnostic invariants: strategies discovered for $5 \times 5$ matrices generalize to higher dimensions such as $64 \times 64$. The discovered programs encode interpretable, dimension-agnostic construction rules. These rules generalize across matrix sizes without retraining, demonstrating that autonomous program synthesis can serve as a scalable paradigm for algorithm discovery and the automated design of universal unitary operators. Beyond universal decompositions, the system automatically exploits matrix structure to reduce the interferometer count below the universal theoretical bound. For instance, for Householder matrices, it discovers a dimension-independent rule that requires only $2N-3$ MZIs. This achieves linear, rather than quadratic, scaling and generalizes to arbitrary $N$ without retraining. For matrices obtained from the singular value decomposition of sparse matrices, reductions generally increase with sparsity, reaching up to 38% fewer MZIs than the universal theoretical bound $N(N-1)/2$ at 95% sparsity. These MZI reductions translate directly into practical hardware benefits for scalable photonic implementations. Taken together, the system functions as a single unified engine that discovers both universal decomposition rules and matrix-specific optimizations, without being provided with the structural or analytical properties of the input matrices.
Tags
Links
- Source: https://arxiv.org/abs/2607.10295v1
- Canonical: https://arxiv.org/abs/2607.10295v1
Trouble viewing inline? Open PDF directly →
Full Text
55,960 characters extracted from source content.
Expand or collapse full text
Program-Synthesis-Driven Autodesign of Universal Unitary Operators Yifei Zhang, 1, 2,∗ Dong Chen, 3,∗ Fan Wang, 4, 2 Wenrui Zhang, 3 Yan Chen, 5 Dingding Han, 6, 7 Jianmin Yuan, 8 Xiangjin Kong, 1, 2,† and Yu-Gang Ma 1, 2, 9,‡ 1 Key Laboratory of Nuclear Physics and Ion-beam Application (MOE), Institute of Modern Physics, Fudan University, Shanghai 200433, China 2 Research Center for Theoretical Nuclear Physics, NSFC and Fudan University, Shanghai 200438, China 3 Huawei Technologies Co., Ltd, Beijing 100095, China 4 Department of Industrial Engineering and Decision Analytics, Hong Kong University of Science and Technology, HongKong, China 5 Hunan Key Laboratory of Mechanism and Technology of Quantum Information, Changsha 410073, China 6 School of Information Science and Technology, Fudan University, Shanghai 200433, China 7 Research Institute of Intelligent Complex Systems, Fudan University, Shanghai 200433, China 8 Institute of Atomic and Molecular Physics, Jilin University, Changchun 130012, China 9 School of Physics, East China Normal University, Shanghai 200062, China We demonstrate that AI-driven program synthesis can autonomously discover fundamental strategies for de- composing unitary matrices in photonic networks. By extending DreamCoder to complex-valued linear algebra, the system generates decomposition programs achieving the minimal N(N−1)/2 Mach–Zehnder interferome- ters, distinct from both Reck and Clements architectures. Learned programs encode dimension-agnostic invari- ants: strategies discovered for 5× 5 matrices generalize to higher dimensions such as 64× 64. The discovered programs encode interpretable, dimension-agnostic construction rules. These rules generalize across matrix sizes without retraining, demonstrating that autonomous program synthesis can serve as a scalable paradigm for algorithm discovery and the automated design of universal unitary operators. Beyond universal decompositions, the system automatically exploits matrix structure to reduce the interferometer count below the universal theo- retical bound. For instance, for Householder matrices, it discovers a dimension-independent rule that requires only 2N− 3 MZIs. This achieves linear, rather than quadratic, scaling and generalizes to arbitrary N without retraining. For matrices obtained from the singular value decomposition of sparse matrices, reductions gener- ally increase with sparsity, reaching up to 38% fewer MZIs than the universal theoretical bound N(N− 1)/2 at 95% sparsity. These MZI reductions translate directly into practical hardware benefits for scalable photonic implementations. Taken together, the system functions as a single unified engine that discovers both universal decomposition rules and matrix-specific optimizations, without being provided with the structural or analytical properties of the input matrices. 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 com- puting and quantum information processing [2–4]. In linear optical systems, unitary matrices determine how multimode fields interfere and propagate through reconfigurable inter- ferometric networks, allowing direct hardware implementa- tions of matrix operations with ultrahigh speed, low energy consumption, and intrinsic parallelism [4, 5]. These capa- bilities make photonic platforms attractive for optical neu- ral networks [6–11], where unitaries implement trainable lin- ear 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], motivat- ing scalable and systematic methods to realize arbitrary high- dimensional unitaries in hardware. In practice, arbitrary N × 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]. Mathemati- cally, any matrix can be expressed via singular value decom- position (SVD) as W = U ΣV † [21], producing two unitary matrices U and V † and a diagonal matrix Σ. Since diagonal matrices Σ can be directly realized with MZIs, the challenge reduces to efficiently decomposing the unitary components U and V † . 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, form- ing the algorithmic backbone for universal linear optics and underpinning most present-day photonic processors [7, 12]. Despite their foundational role, these schemes reflect spe- cific human-designed elimination strategies derived from ana- lytic reasoning and occupy only a small region of the com- binatorial space of possible elimination orderings, leaving open the question of whether alternative, equally fundamen- tal strategies exist. To address this, we employ automated reasoning to systematically explore this landscape, uncover- ing strategies that extend beyond human intuition and yield generalizable constructions. In pursuit of this goal, we de- velop a synthesis engine—built on the DreamCoder frame- work [24] and extended to manipulate complex-valued linear algebra—that autonomously discovers strategies for decom- posing unitary matrices. The synthesized programs are not merely confined to known triangular or rectangular meshes. Instead, the system identifies previously unreported elimina- tion orderings that achieve the optimal N (N − 1)/2 MZIs. Analysis of the learned programs reveals that these con- structions encode dimension-agnostic algorithmic regulari- ties: strategies discovered on 5×5 matrices generalize to arbi- trary higher dimensions such as 64×64 without retraining, in- dicating that the system has internalized domain-general pat- arXiv:2607.10295v1 [physics.optics] 11 Jul 2026 2 (a) Task Series of MZI operations Unitary→Diagonal(Multiple Strategies) (b) programs decomposition patterns primitives Wake Decompose unitary matrices into MZIs Sleep Extract reusable decomposition patterns Library Growing primitive library Automated Discovery FIG. 1. Autonomous discovery of unitary decomposition strategies. (a) Task transformation: unitary matrices are decomposed into diagonal form through a series of MZI operations, with multiple decomposition strategies available. (b) DreamCoder’s abstraction loop: candidate MZI programs are enumerated and tested on decomposition tasks, recurring fragments from successful programs are compressed into reusable library primitives, and the enlarged library is reused to bias subsequent searches toward deeper decompositions. terns rather than dimension-specific heuristics. These results establish program synthesis as a means of autonomously iden- tifying algorithmic structure beyond known human-designed schemes, enabling the discovery of interpretable and gener- alizable 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 com- plexity 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 discover- ing generalizable decomposition rules and exploiting matrix- specific structure, enabling scalable automated design of pho- tonic unitary circuits. As shown in Fig. 1(a), our learning objective is to discover programs that transform arbitrary unitary matrices into diag- onal form. This formulation directly corresponds to photonic implementation: synthesizing an arbitrary unitary matrix U on a photonic chip decomposes into two steps—an MZI cas- cade network P and a diagonal matrix D such that U = P † D. Since D requires no algorithmic search, our learning objec- tive simplifies to finding the sequence P that performs the di- agonalization. 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 in- formation e iφ j on its diagonals. Success is determined solely by the diagonalization constraint: specifically, all off-diagonal elements (i ̸= j) must satisfy |p(U k ) ij | < ε elem , where ε elem = 5 × 10 −4 in the simulation. The specific value of this threshold does not affect program acceptance, as success- fully synthesized programs typically drive residuals down to machine precision (∼ 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 de- noting the target element position at row i, column j. Each primitive eliminates one matrix element. The decomposition strategy exploits the unitarity constraint U † U = I , where systematic elimination of lower-triangular elements automat- ically 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 (UR ij ) ij = 0, the required phase θ ij and mixing angle ω ij are analytically derived: the ratio r = u ij /u i,j+1 is computed, with θ ij = arg(r) compensating the phase dif- ference and ω ij = arctan(|r|) determining the mixing ratio. This constructs an N × N transformation matrix by embed- ding the 2× 2 MZI block M R,ij = e −iθ ij cosω ij e −iθ ij sinω ij − sinω ij cosω ij (1) at positions (j,j + 1) within an identity matrix, applying it via right multiplication U ′ = UR ij . L ij -type primitives similarly mix rows i and i− 1 via left multiplication, deriving param- eters from condition (L ij U ) ij = 0 and embedding the 2× 2 MZI block M L,ij = e iθ ij cosω ij − sinω ij e iθ ij sinω ij cosω ij (2) at positions (i− 1,i) within an identity matrix, applying it via left multiplication U ′ = L ij U . The complete N × N matrix forms of R ij and L ij primitives are provided in Supplemen- tal 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 × 3 decomposition might be (lambda (R21 (L10 (L20 $0)))) [applied 3 in sequence to positions (2, 0), (1, 0), and (2, 1)], where $0 denotes the input matrix and operations are applied from in- nermost to outermost. Given the task definition and primitive operations, we em- ploy program synthesis to autonomously discover effective decomposition strategies. We adopt the DreamCoder frame- work [24], which operates through an iterative wake-sleep cy- cle, as illustrated in Fig. 1(b). During the wake phase, the system searches over candidate MZI operation sequences us- ing the current primitive library and evaluates them against decomposition tasks. During the sleep phase, it identifies re- curring patterns across successful solutions and abstracts them into new reusable primitives, expanding the library. The up- dated 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; imple- mentation details are provided in SM. The MZI primitives R ij ,L ij encode physically meaningful beam-splitter trans- formations with analytically derived phase and rotation an- gles. This design can be generalized to a broad class of quan- tum 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, success- fully verified programs achieve off-diagonal residuals at ma- chine precision (∼ 10 −16 ). Furthermore, we validate the cross-dimensional generalization of universal programs by ap- plying 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 knowl- edge of classical schemes. This approach ensures discov- ered strategies emerge from the search process itself rather than reflecting biases from known methods.We first in- vestigate the 5 × 5 case using independent learning within tractable search spaces. After nine wake-sleep iterations, the system autonomously discovered multiple distinct decompo- sition strategies, all achieving the theoretical minimum of 10 MZI operations. Figure 2 visualizes four random decompo- sition patterns, where the orange paths trace the sequence of operations through the matrix. Each pattern shows a differ- ent 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 ma- trices, 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. × × × × × Pattern 1 × × × × × Pattern 2 × × × × × Reck × × × × × Pattern 3 × × × × × Pattern 4 × × × × × Clements FIG. 2. Decomposition patterns discovered for 5× 5 unitary ma- trices compared with classical reference schemes. Left four panels (orange paths, yellow background) show representative discovered patterns exhibiting diverse elimination orderings. Right two panels show classical reference schemes for comparison: Reck (top, blue path) uses bottom-up row-wise elimination, while Clements (bottom, green path) employs symmetric alternating diagonal pattern. All pat- terns achieve the theoretical minimum of 10 MZI operations. The discovered patterns demonstrate that the system learned flexible de- composition strategies different from classical approaches. Among the discovered programs, the one ranked high- est by the system, which favors solutions that are both cor- rect and structurally simple (i.e., expressible more com- pactly in terms of the learned library operations) [29, 30] is (lambda (R10 (R21 (R20 (R32 (R31 (R30 (R43 (R42 (R41 (L40 $0))))))))))). This pro- gram 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 ,(3) where D is the resulting diagonal matrix. Inverting this rela- tion yields the unitary decomposition U = L −1 40 DR −1 10 R −1 21 R −1 20 R −1 32 R −1 31 R −1 30 R −1 43 R −1 42 R −1 41 . (4) This program directly maps to a photonic circuit implemen- tation, where the five optical waveguides correspond to ma- trix modes and each MZI block executes one primitive op- eration. 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 decom- position strategy, revealing that the systematic elimination of lower-triangular elements admits multiple valid operation se- quences. Critically, the learned patterns exhibit cross-dimensional reasoning.Analysis of discovered programs reveals dimension-agnostic principles: systematic elimination strate- gies that remain valid across matrix scales. For instance, al- gorithmic abstractions derived from 5 × 5 matrices exhibit compositional generalization: they can be directly applied to higher-dimensional unitary matrices without retraining or 4 (a) 102030405060708090 NNZ-Controlled Sparsity (%) 70 75 80 85 90 95 100 MZI Count Ratio (%) MZI Ratio under NNZ-Controlled Sparsity (N=6) (b) 20406080 Sparsity Parameter p (%) 65 70 75 80 85 90 95 100 MZI Count Ratio (%) MZI Ratio under Bernoulli Sparsity Model (N=6) FIG. 3. MZI count ratio [normalized to the theoretical maximum N(N− 1)/2] as a function of sparsity level for N = 6 matrices under NNZ-controlled and Bernoulli sparsity models. Both models show a consistent overall decrease, confirming that the MZI reduction reflects an intrinsic property of sparse unitary matrices. Each data point is averaged over nearly 100 successfully decomposed matrices per sparsity level. 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 decrement- ing by two each step, the algorithm processes rows in pairs: it first emits primitiveR i,0 , then interleaves primitives from two adjacent rows asR i,j →R i−1,j−1 for j = 1,...,i− 1. When a single row remains (i = 1),R 1,0 closes the sequence. This fixed rule, discovered on 5× 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× 5 training instances against pro- gressively larger matrices, including 8× 8, 16× 16, 32× 32, and 64 × 64 dimensions. All test cases achieve successful diagonalization at machine precision (∼ 10 −16 ), confirming that the learned decomposition strategy scales correctly to ar- bitrary dimensions. While our learned patterns exhibit dimension-agnostic prin- ciples that allow 5× 5 strategies to generalize, direct train- ing on 6× 6 matrices remains essential to uncover richer de- composition 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 al- lows 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 con- solidated abstractions effectively reduce the complexity for later stages (T 0 ), enabling the synthesis of deeply nested pro- grams. 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 li- brary from 30 to 65 hierarchical abstractions. The system successfully rediscovered the Reck decomposition and iden- tified novel patterns distinct from 5× 5 independent learning. One such pattern can be expressed by the diagonalization re- lation 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 dis- covers 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 matri- ces, with Householder matrices serving as a prime exam- ple. Defined as H = I − 2v † (I is the identity matrix, v ∈C N , ∥v∥ = 1), Householder matrices are unitary reflec- tors widely used in mathematics and physics for stable uni- tary transformations, including tridiagonalization of Hamil- tonians 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 to- tal MZI count to 2N − 3. The recursive logic is as follows: extending to N + 1 dimensions requires prepending two op- erations,L N,0 followed byR N,N−1 , to the N -dimensional se- quence. 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 opera- tions = 2× 6− 3). This strategy reduces the hardware scaling from the universal quadratic bound of N (N − 1)/2 to a lin- ear 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% reduc- tion 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 5 (W = U ΣV † ), where the unitary factors U and V † in- herit 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 theo- retical bound N (N− 1)/2, and these reductions generally in- crease 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 arbi- trary 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 ben- efits: fewer components can help reduce optical insertion loss [23], calibration complexity [38], and power consump- tion [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 univer- sal unitaries—going beyond known schemes to reveal de- composition principles and resource optimizations that clas- sical Reck and Clements designs cannot capture. Built on a complex-valued matrix decomposition domain, our frame- work combines symbolic reasoning with physically grounded MZI primitives to systematically generate and verify decom- position programs, providing a scalable foundation that can be extended to broader quantum information tasks. These automatically generated constructions encode systematic con- struction methods, demonstrating the potential of AI-driven algorithm discovery to reveal fundamental principles and of- fering a powerful new approach for advancing photonic com- puting and quantum information processing. Acknowledgments—This work is supported by the Na- tional Key Research and Development Program of China un- der Contract No. 2024YFA1610900 and the National Nat- ural 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 Univer- sity. ∗ These authors contributed equally to this work. † kongxiangjin@fudan.edu.cn ‡ mayugang@fudan.edu.cn [1] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, Cam- bridge, England, 2010). [2] P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Rev. Mod. Phys. 79, 135 (2007). [3] B. J. Shastri, A. N. Tait, T. Ferreira de Lima, W. H. Pernice, H. Bhaskaran, C. D. Wright, and P. R. Prucnal, Photonics for artificial intelligence and neuromorphic computing, Nat. Pho- tonics 15, 102 (2021). [4] P. L. McMahon, The physics of optical computing, Nat. Rev. Phys. 5, 717 (2023). [5] E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature (London) 409, 46 (2001). [6] Y. Shen, N. C. Harris, S. Skirlo, M. Prabhu, T. Baehr-Jones, M. Hochberg, X. Sun, S. Zhao, H. Larochelle, D. Englund et al., Deep learning with coherent nanophotonic circuits, Nat. Pho- tonics 11, 441 (2017). [7] D. P ́ erez, I. Gasulla, L. Crudgington, D. J. Thomson, A. Z. Khokhar, K. Li, W. Cao, G. Z. Mashanovich, and J. Cap- many, Multipurpose silicon photonics signal processor core, Nat. Commun. 8, 636 (2017). [8] X. Lin, Y. Rivenson, N. T. Yardimci, M. Veli, Y. Luo, M. Jar- rahi, and A. Ozcan, All-optical machine learning using diffrac- tive deep neural networks, Science 361, 1004 (2018). [9] S. Pai, Z. Sun, T. W. Hughes, T. Park, B. Bartlett, I. A. Williamson, M. Minkov, M. Milanizadeh, N. Abebe, F. Morichetti et al., Experimentally realized in situ backprop- agation for deep learning in photonic neural networks, Science 380, 398 (2023). [10] Z. Xue, T. Zhou, Z. Xu, S. Yu, Q. Dai, and L. Fang, Fully for- ward mode training for optical neural networks, Nature (Lon- don) 632, 280 (2024). [11] T. Fu, J. Zhang, R. Sun, Y. Huang, W. Xu, S. Yang, Z. Zhu, and H. Chen, Optical neural networks: Progress and challenges, Light Sci. Appl. 13, 263 (2024). [12] W. Bogaerts, D. P ́ erez, J. Capmany, D. A. Miller, J. Poon, D. Englund, F. Morichetti, and A. Melloni, Programmable pho- tonic circuits, Nature (London) 586, 207 (2020). [13] J. L. O’brien, Optical quantum computing, Science 318, 1567 (2007). [14] S. Konno, W. Asavanant, F. Hanamura, H. Nagayoshi, K. Fukui, A. Sakaguchi, R. Ide, F. China, M. Yabuno, S. Miki et al., Log- ical states for fault-tolerant quantum computation with propa- gating light, Science 383, 289 (2024). [15] PsiQuantum team, A manufacturable platform for photonic quantum computing, Nature (London) 641, 876 (2025). [16] E. Pelucchi, G. Fagas, I. Aharonovich, D. Englund, E. Figueroa, Q. Gong, H. Hannes, J. Liu, C.-Y. Lu, N. Matsuda et al., The potential and global outlook of integrated photonics for quan- tum technologies, Nat. Rev. Phys. 4, 194 (2022). [17] I. Bente,S. Taheriniya,F. Lenzini,F. Br ̈ uckerhoff- Pl ̈ uckelmann, M. Kues, H. Bhaskaran, C. D. Wright, and W. Pernice, The potential of multidimensional photonic computing, Nat. Rev. Phys. 7, 439 (2025). [18] S. Shekhar, W. Bogaerts, L. Chrostowski, J. E. Bowers, M. Hochberg, R. Soref, and B. J. Shastri, Roadmapping the next generation of silicon photonics, Nat. Commun. 15, 751 (2024). 6 [19] D. A. Miller, Self-configuring universal linear optical compo- nent, Photonics Res. 1, 1 (2013). [20] J. Carolan, C. Harrold, C. Sparrow, E. Mart ́ ın-L ́ opez, N. J. Rus- sell, J. W. Silverstone, P. J. Shadbolt, N. Matsuda, M. Oguma, M. Itoh et al., Universal linear optics, Science 349, 711 (2015). [21] G. Golub and W. Kahan, Calculating the singular values and pseudo-inverse of a matrix, J. Soc. Ind. Appl. Math. 2, 205 (1965). [22] M. Reck, A. Zeilinger, H. J. Bernstein, and P. Bertani, Exper- imental realization of any discrete unitary operator, Phys. Rev. Lett. 73, 58 (1994). [23] W. R. Clements, P. C. Humphreys, B. J. Metcalf, W. S. Kolthammer, and I. A. Walmsley, Optimal design for universal multiport interferometers, Optica 3, 1460 (2016). [24] K. Ellis, C. Wong, M. Nye, M. Sabl ́ e-Meyer, L. Morales, L. He- witt, L. Cary, A. Solar-Lezama, and J. B. Tenenbaum, Dream- coder: Bootstrapping inductive program synthesis with wake- sleep library learning, in Proceedings of the 42nd ACM SIG- PLAN International Conference on Programming Language Design and Implementation (Association for Computing Ma- chinery (ACM), New York, 2021), p. 835–850. [25] F. Mezzadri, How to generate random matrices from the classi- cal compact groups, Not. Am. Math. Soc. 54, 592 (2007). [26] G. E. Hinton, P. Dayan, B. J. Frey, and R. M. Neal, The “wake- sleep” algorithm for unsupervised neural networks, Science 268, 1158 (1995). [27] L. Sarra, K. Ellis, and F. Marquardt, Discovering quantum cir- cuit components with program synthesis, Mach. Learn. Sci. Technol. 5, 025029 (2024). [28] F. J. Ruiz, T. Laakkonen, J. Bausch, M. Balog, M. Barekatain, F. J. Heras, A. Novikov, N. Fitzpatrick, B. Romera-Paredes, J. Van De Wetering et al., Quantum circuit optimization with alphatensor, Nat. Mach. Intell. 7, 374 (2025). [29] B. M. Lake, R. Salakhutdinov, and J. B. Tenenbaum, Human- level concept learning through probabilistic program induction, Science 350, 1332 (2015). [30] P. Liang, M. I. Jordan, and D. Klein, Learning programs: A hi- erarchical Bayesian approach, in ICML (Omnipress, Madison, 2010), Vol. 10, p. 639–646. [31] Y. Bengio, J. Louradour, R. Collobert, and J. Weston, Curricu- lum learning, in Proceedings of the 26th annual international conference on machine learning (Omnipress, Madison, 2009), p. 41–48. [32] X. Wang, Y. Chen, and W. Zhu, A survey on curriculum learn- ing, IEEE Trans. Pattern Anal. Mach. Intell. 44, 4555 (2021). [33] G. H. Golub and C. F. Van Loan, Matrix Computations (JHU Press, Baltimore, Maryland, 2013). [34] P. A. Ivanov, E. Kyoseva, and N. V. Vitanov, Engineering of ar- bitrary U (N) transformations by quantum householder reflec- tions, Phys. Rev. A 74, 022323 (2006). [35] C. Lu, X. Wang, and G. Ma, Experimental realization of special-unitary operations in classical mechanics by nonadia- batic evolutions, Phys. Rev. Lett. 135, 027201 (2025). [36] H. M. Markowitz, The elimination form of the inverse and its application to linear programming, Manage. Sci. 3, 255 (1957). [37] J. R. Gilbert, C. Moler, and R. Schreiber, Sparse matrices in matlab: Design and implementation, SIAM J. Matrix Anal. Appl. 13, 333 (1992). [38] S. Lin, J. Zeng, S. Lin, S. Yu, and Y. Zhang, High-fidelity and compact topology architecture for large-scale reconfig- urable linear optical networks, Adv. Photonics Nexus 4, 066012 (2025). [39] K. H. R. Mojaver, B. Zhao, E. Leung, S. M. R. Safaee, and O. Liboiron-Ladouceur, Addressing the programming chal- lenges of practical interferometric mesh based optical proces- sors, Opt. Express 31, 23851 (2023). [40] H. Zhou, H. Yang, N. Gangi, Z. R. Huang, H. Ren, and J. Gu, Apollo: Automated routing-informed placement for large-scale photonic integrated circuits, in 2025 IEEE/ACM International Conference On Computer Aided Design (ICCAD) (IEEE, New York, 2025), p. 1–9. [41] H. Zhou, K. Zhu, and J. Gu, Automated curvy waveguide routing for large-scale photonic integrated circuits, arXiv: 2410.01260. 1 SUPPLEMENTAL MATERIAL FOR “PROGRAM-SYNTHESIS–DRIVEN AUTO-DESIGN OF UNIVERSAL UNITARY OPERATORS” 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 decomposi- tion domain implementation. The DreamCoder framework operates through alternating wake and sleep phases. During the wake phase, the system per- forms 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 probabil- ities 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 ) ∝ 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 programsp 1 ,...,p n , the algorithm identifies a new abstraction f and rewritten programsp ′ 1 ,...,p ′ n that minimize the total cost cost(f ) + P i 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 matricesU 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 | < ε elem = 5× 10 −4 for all k and all i̸= j. In practice, successfully synthesized programs drive all off-diagonal elements to near machine precision (∼ 10 −16 ), so the specific value of ε elem does not affect which programs are accepted. Programs passing this check are then ranked by posterior probability P (p|T ) ∝ P (T|p)P (p), where the likelihood P (T|p) is measured by the average off- diagonal Frobenius norm 1 K P k ∥p(U k )− 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 logP (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 × 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: r = u ij /u i,j+1 , θ ij = arg(r), ω ij = arctan(|r|). The 2× 2 MZI transformation matrix is: M R,ij = e −iθ ij cosω ij e −iθ ij sinω ij − sinω ij cosω ij .(S.1) 2 Matrix type Define unitary matrix type for program synthesis Matrix representation Encode complex matrix as nested list of (real, imag) pairs MZI primitives GenerateN(N−1) beam-splitter operationsR ij ,L ij forN×N matrices Task definition Define training tasks as random unitary matrices paired with target diagonal matrices Parse input matrix Receive complex matrix data from Python frontend Execute MZI operations Apply each beam-splitter transformation to the matrix Evaluate program Accept if all off-diagonal ele- ments satisfy|D ij |< elem ; score by decomposition accuracy Search strategy Enumerate MZI programs ranked by accuracy and simplicity; re- turn top-Kvalid decompositions Python Frontend OCaml Backend JSON FIG. S1. Program synthesis framework for unitary matrix decomposition. The system learns from unitary matrices U and seeks to transform them into diagonal form, discovering MZI operation sequences via wake-sleep cycles. During waking, candidate MZI sequences are enumer- ated and evaluated against the diagonalization criterion. During sleep, recurring sub-sequences are abstracted into reusable building blocks and added to a growing library, enabling discovery of increasingly complex decompositions in later cycles. This submatrix is embedded in an N × N identity matrix, replacing columns (j,j + 1): G R ij = 1 · 000 · 0 . . . . . . . . . . . . . . . . . . . . . 0 · 100 · 0 0 · 0 e −iθ ij cosω ij e −iθ ij sinω ij · 0 0 · 0 − sinω ij cosω ij · 0 0 · 000 · 0 . . . . . . . . . . . . . . . . . . . . . 0 · 000 · 1 (S.2) Right multiplication U ′ = U · G R ij mixes columns j and j + 1, zeroing u ′ ij while preserving unitarity U ′† U ′ = I . The L ij primitive operates analogously via left multiplication, using element u ij and its upper neighbor u i−1,j : r = u ij /u i−1,j , θ ij = arg(r), ω ij = arctan(−|r|). The 2× 2 MZI transformation matrix is: M L,ij = e iθ ij cosω ij − sinω ij e iθ ij sinω ij cosω ij .(S.3) 3 This submatrix is embedded in an N × N identity matrix, replacing rows (i− 1,i): G L ij = 1 · 000 · 0 . . . . . . . . . . . . . . . . . . . . . 0 · 100 · 0 0 · 0 e iθ ij cosω ij − sinω ij · 0 0 · 0 e iθ ij sinω ij cosω ij · 0 0 · 000 · 0 . . . . . . . . . . . . . . . . . . . . . 0 · 000 · 1 .(S.4) Left multiplication U ′ = G L ij · U mixes rows i− 1 and i, zeroing u ′ ij while preserving unitarity U ′† U ′ = I . All MZI primitives share the uniform type signature tmatrix → 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 | < ε elem for i ̸= 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 ) ∝ 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× 5 matrices, independent learning successfully discovers decompositions with optimal MZI counts by constructing complete solutions from scratch. However, scaling to 6× 6 and larger dimensions benefits from curriculum learning, a training strategy that progressively increases task difficulty. The curriculum consists of subtasksT k ,T k−1 ,...,T 0 of increasing com- plexity, 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× 6 and larger matrices, where the lengthy search space proves intractable without hierarchical abstraction. Tasks are gener- ated 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× 6 Matrices Scaling to 6× 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 S2 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. S2 shows rapid growth 4 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. FIG. S2. New primitives added per iteration during curriculum-assisted learning for 6× 6 matrices. 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 vocab- ulary. REPRESENTATIVE SYNTHESIZED PROGRAMS Tables S1 and S2 present representative decomposition programs discovered by the DreamCoder framework under indepen- dent 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 | < ε elem for i̸= j. TABLE S1. 5× 5 programs (independent learning). All programs are validated on random unitary matrices. #Log-PosteriorExecution Order 1 −0.69L 40 ,R 41 ,R 42 ,R 43 ,R 30 ,R 31 ,R 32 ,R 20 ,R 21 ,R 10 2 −11.6L 40 ,L 30 ,R 41 ,R 42 ,R 43 ,R 31 ,R 32 ,R 20 ,R 21 ,R 10 3 −11.6L 40 ,R 41 ,R 42 ,L 30 ,R 31 ,R 43 ,R 32 ,R 20 ,R 21 ,R 10 4 −11.6L 40 ,R 41 ,R 30 ,R 42 ,R 43 ,R 31 ,R 32 ,R 20 ,R 21 ,R 10 5 −12.3L 40 ,L 30 ,L 41 ,R 42 ,R 43 ,R 31 ,R 32 ,R 20 ,R 21 ,R 10 TABLE S2. 6× 6 programs (curriculum learning). All programs are validated on random unitary matrices. #Log-PosteriorExecution Order 1 −1.39R 50 ,R 51 ,R 52 ,R 53 ,R 54 ,R 40 ,R 41 ,R 42 ,R 43 ,R 30 ,R 31 ,R 20 ,R 32 ,R 21 ,R 10 2 −1.39L 50 ,R 51 ,R 52 ,R 53 ,R 54 ,R 40 ,R 41 ,R 42 ,R 43 ,R 30 ,R 31 ,R 20 ,R 32 ,R 21 ,R 10 3 −15.5R 50 ,L 40 ,R 51 ,R 52 ,R 53 ,R 54 ,R 41 ,R 42 ,R 30 ,R 43 ,R 31 ,R 20 ,R 32 ,R 21 ,R 10 4 −20.1R 50 ,R 51 ,R 52 ,R 53 ,R 54 ,L 40 ,R 41 ,R 42 ,R 43 ,R 30 ,R 31 ,R 32 ,R 20 ,R 21 ,R 10 5 −23.1L 50 ,L 40 ,L 51 ,R 52 ,R 53 ,R 54 ,R 41 ,R 42 ,R 30 ,R 43 ,R 31 ,R 20 ,R 32 ,R 21 ,R 10 5 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 − 2v † (wherev ∈C N ,∥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 × 5 − 3), with the representative program L 4,0 ,R 4,3 ,R 3,0 ,R 3,1 ,R 3,2 ,R 2,1 ,R 1,0 . For N = 6, every task is solved in exactly 9 operations (= 2× 6− 3), with the representative programL 5,0 ,R 5,4 ,L 4,0 ,R 4,3 ,R 3,0 ,R 3,1 ,R 3,2 ,R 2,1 ,R 1,0 . The N → N + 1 extension prependsL N,0 andR 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 ΣV † yields unitary factors U and V † 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 = −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.000i0.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.000i0.000+0.000i 0.000+0.000i 0.000+0.000i 0.896−0.443i 0.000+0.000i 0.000+0.000i0.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 . 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 (∼ 10 −16 ). Bernoulli sparsity. We consider a representative example with source W : 86.1% sparse, 5/36 nonzero; each entry indepen- dently zeroed with probability 80%. The input U matrix is U = 0.000+0.000i0.000+0.000i 1.000+0.000i 0.000+0.000i0.000+0.000i 0.000+0.000i 0.336+0.000i0.942+0.000i 0.000+0.000i 0.000+0.000i0.000+0.000i 0.000+0.000i −0.532−0.643i 0.190+0.229i 0.000+0.000i 0.173−0.430i0.000+0.000i 0.000+0.000i 0.000+0.000i0.000+0.000i 0.000+0.000i 0.000+0.000i −1.000+0.015i 0.000+0.000i 0.000+0.000i0.000+0.000i 0.000+0.000i 0.000+0.000i0.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 . 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 (∼ 10 −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. 6 How low-dimensional patterns are expanded to high dimensions. Example 1 — Sequential row elimination (Table S1, Program 1):the N=5sequenceis L 4,0 ,R 4,1 ,R 4,2 ,R 4,3 ,R 3,0 ,R 3,1 ,R 3,2 ,R 2,0 ,R 2,1 ,R 1,0 .The rule is: for the top row (i = N − 1), applyL i,0 followed by R i,1 ,...,R i,i−1 . For all lower rows (i = N − 2 down to 1), applyR i,0 ,R i,1 ,...,R i,i−1 in order. Each row is fully eliminated before moving to the next. The operation count follows directly: P N−1 i=1 i = N (N− 1)/2, yielding 2016 operations at N = 64. Example 2 — Two-row leading elimination (Table S1, Program 2):the N=5sequenceis L 4,0 ,L 3,0 ,R 4,1 ,R 4,2 ,R 4,3 ,R 3,1 ,R 3,2 ,R 2,0 ,R 2,1 ,R 1,0 . The rule is: first applyL N−1,0 andL N−2,0 ; then complete row N − 1 with R N−1,1 ,...,R N−1,N−2 ; then row N− 2 withR N−2,1 ,...,R N−2,N−3 (starting from j = 1 since j = 0 was handled byL N−2,0 ); then all remaining rows i = N − 3 down to 1 withR i,0 ,R i,1 ,...,R i,i−1 . The same operation count N (N − 1)/2 is preserved. We verified Programs 1, 2, and 3 from Table S1 across N = 5, 6, 7, 8, 16, 32, 64; all achieve machine-precision accuracy (∼ 10 −16 ) with 100% success rate, as shown in Table S3. This confirms that the decompositions are mathematically exact: the residual error reflects only floating-point rounding and does not grow systematically with N . TABLE S3. Numerical precision verification for three universal programs across different dimensions. Off-diagonal error averaged over 5 random unitary matrices. NProgram 1Program 2Program 3 5 2.61× 10 −16 2.92× 10 −16 2.73× 10 −16 6 3.10× 10 −16 3.32× 10 −16 3.12× 10 −16 7 4.38× 10 −16 4.65× 10 −16 4.23× 10 −16 8 4.20× 10 −16 3.92× 10 −16 4.66× 10 −16 16 4.98× 10 −16 4.73× 10 −16 4.67× 10 −16 32 4.71× 10 −16 4.84× 10 −16 5.38× 10 −16 64 5.69× 10 −16 6.04× 10 −16 6.24× 10 −16 We next present the detailed algorithmic procedure for Pattern 4 (row-pair interleaving, Fig. 2; distinct from the programs listed in Table S1) 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 primitiveR i,0 , then interleaving primitives from two adjacent rows asR i,j →R i−1,j−1 for j = 1,...,i− 1. When a single row remains (i = 1), 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 ×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 7 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× 5) to deployment dimensions (e.g., 64× 64 or larger). The generalization validation protocol tests learned programs on matrices of progressively increasing dimensions: 8 × 8, 16 × 16, 32 × 32, and 64 × 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 | < ε elem for all k and all i ̸= j, with residuals at machine precision (∼ 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.