Decompose a unitary into a sequence of Givens rotation gates with phase shifts and a diagonal phase matrix.

This decomposition is based on the construction scheme given in Optica, 3, 1460 (2016), which allows one to write any unitary matrix \(U\) as:

\[U = D \left(\prod_{(m, n) \in G} T_{m, n}(\theta, \phi)\right),\]

where \(D\) is a diagonal phase matrix, \(T(\theta, \phi)\) is the Givens rotation gates with phase shifts and \(G\) defines the specific ordered sequence of the Givens rotation gates acting on wires \((m, n)\). The unitary for the \(T(\theta, \phi)\) gates appearing in the decomposition is of the following form:

\[\begin{split}T(\theta, \phi) = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & e^{i \phi} \cos(\theta) & -\sin(\theta) & 0 \\ 0 & e^{i \phi} \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix},\end{split}\]

where \(\theta \in [0, \pi/2]\) is the angle of rotation in the \(\{|01\rangle, |10\rangle \}\) subspace and \(\phi \in [0, 2 \pi]\) represents the phase shift at the first wire.


unitary = np.array([[ 0.73678+0.27511j, -0.5095 +0.10704j, -0.06847+0.32515j],
                    [-0.21271+0.34938j, -0.38853+0.36497j,  0.61467-0.41317j],
                    [ 0.41356-0.20765j, -0.00651-0.66689j,  0.32839-0.48293j]])

phase_mat, ordered_rotations = givens_decomposition(unitary)
>>> phase_mat
tensor([-0.20604358+0.9785369j , -0.82993272+0.55786114j,
        0.56230612-0.82692833j], requires_grad=True)
>>> ordered_rotations
[(tensor([[-0.65087861-0.63937521j, -0.40933651-0.j        ],
        [-0.29201359-0.28685265j,  0.91238348-0.j        ]], requires_grad=True),
(0, 1)),
(tensor([[ 0.47970366-0.33308926j, -0.8117487 -0.j        ],
        [ 0.66677093-0.46298215j,  0.5840069 -0.j        ]], requires_grad=True),
(1, 2)),
(tensor([[ 0.36147547+0.73779454j, -0.57008306-0.j        ],
        [ 0.2508207 +0.51194108j,  0.82158706-0.j        ]], requires_grad=True),
(0, 1))]

unitary (tensor) – unitary matrix on which decomposition will be performed


diagonal elements of the phase matrix \(D\) and Givens rotation matrix \(T\) with their indices.

Return type

(np.ndarray, list[(np.ndarray, tuple)])


ValueError – if the provided matrix is not square.

Givens Rotation

Applying the Givens rotation \(T(\theta, \phi)\) performs the following transformation of the basis states:

\[\begin{split}&|00\rangle \mapsto |00\rangle\\ &|01\rangle \mapsto e^{i \phi} \cos(\theta) |01\rangle - \sin(\theta) |10\rangle\\ &|10\rangle \mapsto e^{i \phi} \sin(\theta) |01\rangle + \cos(\theta) |10\rangle\\ &|11\rangle \mapsto |11\rangle.\end{split}\]


The algorithm that implements the decomposition is the following:

def givens_decomposition(U):
    for i in range(1, N):
        if i % 2:
            for j in range(0, i):
                # Find T^-1(i-j, i-j+1) matrix that nulls element (N-j, i-j) of U
                # Update U = U @ T^-1(i-j, i-j+1)
            for j in range(1, i):
                # Find T(N+j-i-1, N+j-i) matrix that nulls element (N+j-i, j) of U
                # Update U = T(N+j-i-1, N+j-i) @ U