1
$\begingroup$

I am currently trying to benchmark my code with a Haar circuit and I require to sample clifford gates in matrix form. I know a function "stim.Tableau.random(n)" which does that and gives me the tableau representation. Can I obtain the corresponding matrix form of the clifford unitary gate? Qiskit allows me to do it using the function

circuit = qiskit.quantum_info.random_clifford(n)

qiskit.circuit.Gate.to_matrix(circuit)

Is there any function available in stim which does that?

$\endgroup$

1 Answer 1

1
$\begingroup$

In Stim v1.9+ you can get the unitary matrix of a tableau by calling stim.Tableau.to_unitary_matrix.


(Old answer from before v1.9)

Stim v1.8 doesn't have a method to do this, but you can get around that by using the state channel duality and the tableau simulator's state_vector method. Note that, once you get past 8 qubit tableaus, this starts to take multiple seconds to finish.

import stim
import numpy as np


def tableau_to_unitary(tableau: stim.Tableau,
                       *,
                       canonical_global_phase: bool,
                       endian: str) -> np.ndarray:
    assert endian in ['little', 'big']
    sim = stim.TableauSimulator()
    n = len(tableau)

    # Create n Bell pairs.
    for q in range(n):
        sim.h(q)
        sim.cnot(q, q + n)

    # Take one qubit from each Bell pair, and apply the custom tableau to those qubits.
    # Operating on [0,n) instead of [n,2n) gets transposed unitaries (Y rotations are backwards).
    # Reversing the order switches the endian-ness.
    qubits = range(n, 2*n)
    if endian == 'big':
        qubits = qubits[::-1]
    state = sim.current_inverse_tableau()
    state.prepend(tableau**-1, qubits)
    sim.set_inverse_tableau(state)

    # Get the state and interpret it as a matrix via the state channel duality.
    result = sim.state_vector().reshape((2**n, 2**n)) * 2**(n/2)
    if canonical_global_phase:
        result = with_canonical_global_phase(result)
    return result


def with_canonical_global_phase(u: np.ndarray) -> np.ndarray:
    f = u.flat
    best_index = 0
    best_val = abs(f[0])
    for i in range(1, len(f)):
        v = abs(f[i])
        if v > best_val * 2:
            best_index = i
            best_val = v

    v = f[best_index]
    assert v != 0
    v /= abs(v)
    return u * (np.conj(v) / abs(v))

Testing that it works:

sqrt_y = tableau_to_unitary(stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("-Z"),
    ],
    zs=[
        stim.PauliString("X"),
    ],
), canonical_global_phase=True, endian='big')
np.testing.assert_allclose(
    sqrt_y,
    np.array([
        [1, -1],
        [1, 1],
    ]) / np.sqrt(2),
    atol=1e-6,
)


sim = stim.TableauSimulator()
sim.x(0)
sim.h(1)
x_tensor_h = tableau_to_unitary(sim.current_inverse_tableau()**-1,
                                canonical_global_phase=True,
                                endian='big')
np.testing.assert_allclose(
    x_tensor_h,
    np.array([
        [0, 0, 1, 1],
        [0, 0, 1, -1],
        [1, 1, 0, 0],
        [1, -1, 0, 0],
    ]) / np.sqrt(2),
    atol=1e-6,
)

cnot = tableau_to_unitary(stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("XX"),
        stim.PauliString("_X"),
    ],
    zs=[
        stim.PauliString("Z_"),
        stim.PauliString("ZZ"),
    ],
), canonical_global_phase=True, endian='big')

iswap = tableau_to_unitary(stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("ZY"),
        stim.PauliString("YZ"),
    ],
    zs=[
        stim.PauliString("_Z"),
        stim.PauliString("Z_"),
    ],
), canonical_global_phase=True, endian='big')

s = tableau_to_unitary(stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("Y"),
    ],
    zs=[
        stim.PauliString("Z"),
    ],
), canonical_global_phase=True, endian='big')
s_dag = tableau_to_unitary(stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("-Y"),
    ],
    zs=[
        stim.PauliString("Z"),
    ],
), canonical_global_phase=True, endian='big')

np.testing.assert_allclose(
    cnot,
    np.array([
        [1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 0, 1],
        [0, 0, 1, 0],
    ]),
    atol=1e-6,
)

np.testing.assert_allclose(
    iswap,
    np.array([
        [1, 0, 0, 0],
        [0, 0, 1j, 0],
        [0, 1j, 0, 0],
        [0, 0, 0, 1],
    ]),
    atol=1e-6,
)
np.testing.assert_allclose(s, np.diag([1, 1j]), atol=1e-6)
np.testing.assert_allclose(s_dag, np.diag([1, -1j]), atol=1e-6)
$\endgroup$
4
  • $\begingroup$ Thanks a lot for your detailed answer. The code is really helpful and works for me in sampling Clifford unitary in their expected matrix form. I only require it for small number of qubits in order to benchmark it with Haar code. I am also curious to know if their is any additional function which can perform single qubit measurement on a state and return a fixed state(say up spin) as an outcome? which is to similar to performing post selection. $\endgroup$
    – user21113
    Jun 30, 2022 at 3:18
  • $\begingroup$ @user21113 You can use stim.TableauSimulator.measure_kickback to get postselection-like effects, because it tells you how to get to the other post-measurement state. $\endgroup$ Jun 30, 2022 at 7:02
  • $\begingroup$ Thanks for suggesting the pseudo post measurement scheme and for all your answers regarding stim. If I may ask, Is their any way to reverse the output state vector of a tableau simulator from its default little endian to big endian order? $\endgroup$
    – user21113
    Jul 3, 2022 at 19:43
  • $\begingroup$ I found the answer in stim latest version after installing "pip install stim~=1.9.dev". $\endgroup$
    – user21113
    Jul 3, 2022 at 21:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.