Synthesize diagonal unitaries with Walsh coefficients

Download Notebook - diagonal_walsh.ipynb

Turn a list of diagonal phases into commuting Pauli-Z rotations:

\[ D=\operatorname{diag}(e^{i\theta_0},\ldots,e^{i\theta_{2^n-1}}) =\prod_{k=0}^{2^n-1}e^{ia_kW_k}, \qquad W_k=\bigotimes_{r=0}^{n-1}Z_r^{k_r}. \]

The normalized Walsh coefficients are

\[ a_k=2^{-n}\sum_{j=0}^{2^n-1}(-1)^{\sum_r j_rk_r}\theta_j. \]

fast_walsh_hadamard_transform computes the unnormalized transform; divide by \(2^n\) to obtain \(a_k\). diagonal_unitary_walsh builds the rotations and omits the \(k=0\) global phase. Each remaining term uses a parity ladder, a Z rotation, and the inverse ladder.

Run this notebook from a source checkout with the development dependencies installed. The repository default is little endian.

Hide code cell source

import numpy as np
from guppyalgos.algorithms.state_preparation.diagonal import diagonal_unitary_walsh, fast_walsh_hadamard_transform
from guppyalgos.tests.helpers import Endianness, assert_allclose_ignorephase, get_unitary

1. One qubit: the S gate

\[\begin{split} S=\begin{pmatrix}1&0\\0&i\end{pmatrix} =e^{i\pi/4}e^{-i\pi Z/4}. \end{split}\]

The normalized coefficients are \((\pi/4,-\pi/4)\). After removing the global phase, one Z rotation remains.

diagonal_1q = np.array([1.0, 1j])
phases_1q = np.angle(diagonal_1q)
walsh_coeffs_1q = fast_walsh_hadamard_transform(phases_1q)
print("Phases:         ", phases_1q)
print("Walsh coeffs:   ", walsh_coeffs_1q)
print("Normalized a_k: ", walsh_coeffs_1q / len(diagonal_1q))
Phases:          [0.         1.57079633]
Walsh coeffs:    [ 1.57079633 -1.57079633]
Normalized a_k:  [ 0.78539816 -0.78539816]
circ_1q = diagonal_unitary_walsh(diagonal_1q)
unitary_1q = get_unitary(circ_1q, 1, endianness=Endianness.LITTLE)
print("Synthesized unitary:\n", unitary_1q)
assert_allclose_ignorephase(np.diag(diagonal_1q), unitary_1q)
print("Matches S gate up to global phase.")
Synthesized unitary:
 [[7.07106781e-01-7.07106781e-01j 4.32978028e-17+4.32978028e-17j]
 [0.00000000e+00+0.00000000e+00j 7.07106781e-01+7.07106781e-01j]]
Matches S gate up to global phase.

2. Two qubits: the CZ gate

\[ \mathrm{CZ}=\operatorname{diag}(1,1,1,-1) =e^{i\pi/4}e^{-i\pi Z_0/4}e^{-i\pi Z_1/4}e^{i\pi Z_0Z_1/4}. \]

The three non-global Walsh terms produce two single-qubit rotations and one two-qubit parity rotation. The check compares the synthesized matrix with CZ up to global phase.

diagonal_2q = np.array([1.0, 1.0, 1.0, -1.0])
phases_2q = np.angle(diagonal_2q)
walsh_coeffs_2q = fast_walsh_hadamard_transform(phases_2q)
n_terms = np.sum(np.abs(walsh_coeffs_2q / len(diagonal_2q)) > 1e-10)
print("Walsh coeffs: ", walsh_coeffs_2q)
print(f"Non-trivial terms (excluding global phase): {n_terms - 1}")
Walsh coeffs:  [ 3.14159265 -3.14159265 -3.14159265  3.14159265]
Non-trivial terms (excluding global phase): 3
circ_2q = diagonal_unitary_walsh(diagonal_2q)
unitary_2q = get_unitary(circ_2q, 2, endianness=Endianness.LITTLE)
print("Synthesized unitary:\n", np.round(unitary_2q, 6))
assert_allclose_ignorephase(np.diag(diagonal_2q), unitary_2q)
print("Matches CZ gate up to global phase.")
Synthesized unitary:
 [[-0.707107+0.707107j -0.      -0.j        0.      -0.j
  -0.      -0.j      ]
 [ 0.      +0.j       -0.707107+0.707107j  0.      +0.j
   0.      +0.j      ]
 [ 0.      +0.j       -0.      -0.j       -0.707107+0.707107j
  -0.      -0.j      ]
 [ 0.      +0.j        0.      -0.j        0.      +0.j
   0.707107-0.707107j]]
Matches CZ gate up to global phase.

3. Trade accuracy for fewer rotations

truncation_threshold drops small Walsh coefficients. Because the terms commute, the omitted coefficients bound the phase error:

\[ \|D-\widetilde D\|_2\leq\sum_{k\ \mathrm{omitted}}|a_k|, \]

using the same global-phase convention. The example below compares exact and truncated synthesis after aligning their global phases, and reports the Frobenius matrix error.

rng = np.random.default_rng(42)
diagonal_generic = np.exp(1j * rng.uniform(-np.pi, np.pi, 4))

circ_exact = diagonal_unitary_walsh(diagonal_generic, truncation_threshold=0.0)
circ_approx = diagonal_unitary_walsh(diagonal_generic, truncation_threshold=0.3)

u_exact = get_unitary(circ_exact, 2, endianness=Endianness.LITTLE)
u_approx = get_unitary(circ_approx, 2, endianness=Endianness.LITTLE)

phase = np.angle(np.vdot(u_exact, u_approx))
u_approx_aligned = u_approx * np.exp(-1j * phase)
error = np.linalg.norm(u_exact - u_approx_aligned)
print(f"Frobenius error from truncation: {error:.4f}")
Frobenius error from truncation: 0.5445