PhaseShift

Decomposition and approximation tools for linear optical unitaries.

Table of Contents

• About this project
• Package contents
• Installation
• Usage
• Documentation
• Citing this work
• License
• References

About This Project

PhaseShift is a Python package for performing various decompositions and approximations of unitary matrices into planar arrangements of simple optical components. These tools can be used to design and program universal multiport interferometers (UMIs), devices capable of implementing arbitrary linear transformations on multiple optical modes. Such devices have numerous applications in communication, imaging, and information processing.

The algorithms implemented in this package cover two main classes of planar UMI architectures:

• Networks of two-mode components
• Sequences of multichannel components

Two-mode Component Networks

This class of decompositions seeks to express an $N \times N$ unitary matrix as a planar mesh of configurable two-mode unit cells, typically realized using Mach–Zehnder interferometers (MZIs). The first design of this kind was proposed by Reck et al., 1994, who used a triangular mesh of asymmetric Mach–Zehnder interferometers to implement arbitrary unitary transformations. This design was later improved by Clements et al., 2016, who introduced a more compact rectangular mesh using the same unit cells. Bell et al., 2021 further compactified the Clements et al. design by using a symmetric Mach–Zehnder interferometer as unit cell, which helped reduce the optical depth of the interferometer.

<div align="center"> Rectangular mesh of Mach–Zehnder interferometers based on the Clements architecture for 6 modes </div>

Multichannel Component Sequences

This second class of decompositions aims to express an $N \times N$ unitary matrix as a sequence of configurable phase masks interleaved with a multichannel mixing layer, such as the discrete Fourier transform (DFT). Numerical evidence suggests that using $N+1$ layers of phase masks with any dense mixing layer is enough to result in a universal design (Saygin et al., 2020) (Zelaya et al., 2024). The first constructive design based on this approach was proposed by López Pastor et al., 2021 and generates a sequence of $6N + 1$ phase masks to implement a $N \times N$ unitary. We improved this design to reach $4N+1$ and $2N+5$ phase masks (Girouard et al., 2026).

<div align="center"> Sequence of phase masks interleaved with the discrete Fourier transform mixing layer for 6 modes </div>

Package Contents

PhaseShift provides tools to perform exact decompositions and numerical approximations of unitary matrices.

Exact Decompositions

PhaseShift includes four main modules to perform exact decompositions of unitary matrices:

• clements_interferometer: Implementation of the algorithm by Clements et al., 2016 to decompose $N \times N$ unitary matrices into a rectangular mesh of $N(N-1)/2$ asymmetric Mach–Zehnder interferometers.

• bell_interferometer: Implementation of the algorithm by Bell et al., 2021 to decompose $N \times N$ unitary matrices into a rectangular mesh of $N(N-1)/2$ symmetric Mach–Zehnder interferometers.

• lplm_interferometer: Implementation of the algorithm by López Pastor et al., 2021 to decompose $N \times N$ unitary matrices into a sequence of $6N+1$ phase masks interleaved with the DFT matrix.

• fourier_interferometer: Implementation of the Fourier decomposition and the compact Fourier decomposition (Girouard et al., 2026) to decompose $N \times N$ unitary matrices into sequences of $4N+1$ and $2N+5$ phase masks interleaved with DFT respectively.

Optimization Tools

In addition to exact decompositions, PhaseShift also has an optimization subpackage, which contains tools to approximate unitary matrices into a sequence of phase masks interleaved with a chosen mixing layer. The optimization subpackage has two modules:

• fourier_optimizer: Uses the basin-hopping algorithm from scipy.optimize to solve a global minimization problem, yielding the sequence of phase masks that minimizes the infidelity with respect to a target unitary.

• jax_optimizer: Uses Jax and Optax to perform gradient-based optimization of the phase masks with multiple restarts to minimize the infidelity or the geodesic distance (Álvarez-Vizoso et al.) with respect to a target unitary. This algorithm can run efficiently on CPU or GPU and is significantly faster than the SciPy-based implementation.

---

Note: For more detailed descriptions and usage examples, see the documentation of the individual modules.

Installation

Install from PyPI (recommended)

pip install phaseshift

Install from source

You can install PhaseShift from source as follows:

1. Clone the repository

git clone https://github.com/polyquantique/phaseshift.git
cd phaseshift

2. (Optional) Create and activate a virtual environment

• Linux / macOS:

  python3 -m venv venv
  source venv/bin/activate
  

• Windows (Command Prompt):

  python -m venv venv
  venv\Scripts\activate
  

3. Install the package and dependencies

• Standard installation:

  pip install .
  

• Editable (developer) installation:

  pip install -e .[dev]
  

4. (Optional) Run the tests

pip install pytest
pytest tests

Note

For GPU support with JAX on Linux, install the appropriate CUDA-enabled version:

pip install "jax[cuda]" -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html

Usage

1. clements_interferometer module

This first example shows how to use the clements_interferometer module to decompose a random unitary matrix.

>>> from phaseshift import clements_interferometer as ci
>>> from scipy.stats import unitary_group
>>> import numpy as np

The function clements_decomposition performs the decomposition described in Clements et al., 2016 on a random unitary matrix U.

>>> # Generate a random unitary matrix
>>> U = unitary_group(dim = 8, seed = 137).rvs() 

>>> # Compute the Clements decomposition
>>> decomposition = ci.clements_decomposition(U)

The output of this function is a Decomposition object, which has a circuit attribute. Decomposition.circuit is a list of MachZehnder objects that contain the parameters $\theta$ and $\phi$ of each unit cell in the mesh.

>>> # Extract the circuit from the decomposition
>>> circuit = decomposition.circuit

>>> # Print the parameters of the first unit cell in the circuit
>>> print(circuit[0])
MachZehnder(theta=np.float64(0.7697802543915319), phi=np.float64(3.8400842306814207), target=(5, 6))

The function circuit_reconstruction allows to compute the matrix that corresponds to a Decomposition object. This matrix can be compared to the original matrix.

>>> # Reconstruct the unitary matrix from the decomposition
>>> reconstructed_matrix = ci.circuit_reconstruction(decomposition)

>>> # Compare with the initial matrix
>>> print(np.allclose(U, reconstructed_matrix))
True

2. fourier_interferometer module

This example shows how to use the fourier_interferometer module to decompose a random unitary matrix.

>>> from phaseshift import fourier_interferometer as fi
>>> from scipy.stats import unitary_group
>>> import numpy as np

The compact_fourier_decomposition function decomposes a random $N \times N$ unitary matrix into a sequence of $2N+5$ phase masks and $2N+4$ DFT layers.

>>> # Generate a random unitary matrix
>>> U = unitary_group(dim = 8, seed = 137).rvs() 

>>> # Compute the Compact Fourier decomposition
>>> decomposition = fi.compact_fourier_decomposition(U)

The output is a FourierDecomp object, which contains the mask_sequence attribute. FourierDecomp.mask_sequence stores the sequence of phase masks to be interleaved with DFT matrices.

>>> # Extract the mask sequence from the decomposition
>>> mask_sequence = decomposition.mask_sequence

>>> # Print the first mask in the sequence
>>> print(mask_sequence[0].round(3))
[0.707-0.707j 0.707+0.707j 0.707-0.707j 0.707+0.707j 0.707-0.707j
 0.707+0.707j 0.707-0.707j 0.707+0.707j]

The function circuit_reconstruction computes the matrix given by a FourierDecomp object by inserting a DFT matrix between each phase mask. The matrix can then be compared to the initial matrix.

>>> # Reconstruct the unitary matrix from the decomposition
>>> reconstructed_ma