How to write the $iSWAP$ unitary as a linear combination of tensor products between 1-qubit gates? [duplicate]

Question

As far as I understood, it should always be possible to decompose any $n$-qubits unitary $W$ into a linear combination of tensor products between $n$ single-qubit unitaries $U_i$: $$W = \sum_k \lambda_k \left( \bigotimes_{i=1}^{n} U_i \right)$$

For example, the unitary matrix of the $SWAP$ gate can be written as follows: $$SWAP = \frac{1}{2} \left( I \otimes I + X \otimes X + Y \otimes Y + Z \otimes Z\right)$$

How can I write the unitary of the $iSWAP$ gate in the same form? And, more in general, is there a mechanical procedure to find such a decomposition for any given $W$?

Answer 1

ISWAP =
    + II * (0.5+0j)
    + XX * 0.5j
    + YY * 0.5j
    + ZZ * (0.5+0j)

Output is from this code:

iswap = np.array([
    [1, 0, 0, 0],
    [0, 0, 1j, 0],
    [0, 1j, 0, 0],
    [0, 0, 0, 1],
])
iswap_terms = matrix_to_pauli_terms(iswap)
print("ISWAP =")
for iswap_pauli_string, iswap_coefficient in iswap_terms.items():
    print("    +", iswap_pauli_string, '*', repr(iswap_coefficient))
np.testing.assert_array_equal(pauli_terms_to_matrix(iswap_terms), iswap)

Using these utilities, which are not optimized but work on any matrix. The complexity here is $16^n$ where $n$ is the number of qubits. I think it should be possible to achieve $n 4^n$ instead:

import itertools
from typing import Dict, Iterable

import numpy as np

i = np.eye(2)
x = np.array([
    [0, 1],
    [1, 0],
])
z = np.array([
    [1, 0],
    [0, -1],
])
y = np.array([
    [0, -1j],
    [1j, 0],
])
paulis = {"I": i, "X": x, "Y": y, "Z": z}


def pauli_string_to_matrix(pauli_string: Iterable[str]) -> np.ndarray:
    t = np.eye(1)
    for p in pauli_string:
        t = np.kron(t, paulis[p])
    return t


def matrix_to_pauli_terms(matrix: np.ndarray) -> Dict[str, complex]:
    w, h = matrix.shape
    n = w.bit_length() - 1
    assert w == h and w == 2**n  # check matrix is square with power of 2 size

    terms: Dict[str, complex] = {}
    for pauli_tuple in itertools.product("IXYZ", repeat=n):
        pauli_string = ''.join(pauli_tuple)
        pauli_matrix = pauli_string_to_matrix(pauli_string)
        coefficient = complex(np.dot(matrix.flat, np.conj(pauli_matrix.flat))) / 2**n
        if coefficient:
            terms[pauli_string] = coefficient
    if not terms:
        terms["I" * n] = 0
    return terms


def pauli_terms_to_matrix(terms: Dict[str, complex]) -> np.ndarray:
    assert terms
    n = len(next(iter(terms.keys())))
    total = np.zeros(shape=(2**n, 2**n), dtype=np.complex64)
    for pauli_string, coefficient in terms.items():
        total += pauli_string_to_matrix(pauli_string) * coefficient
    return total

---

Update: here's the fast version that takes $O(n4^n)$ time to compute the terms

def matrix_to_pauli_terms_fast(matrix: np.ndarray) -> Dict[str, complex]:
    w, h = matrix.shape
    n = w.bit_length() - 1
    N = 2**n
    assert w == h == N  # check matrix is square with power of 2 size

    # Permute by xoring row coordinate into column coordinate
    term_matrix = np.empty(shape=matrix.shape, dtype=np.complex64)
    term_matrix[0, :] = matrix[0, :]
    indices = np.array(range(N))
    for k in range(1, N):
        indices ^= k ^ (k - 1)
        term_matrix[k, :] = matrix[k, indices]

    # Hadamard transform the columns and account for scalar phase from Y.
    term_matrix.shape = (2,) * (2 * n)
    for k in range(n):
        index: List[Union[slice, int]] = [slice(None)] * (2 * n)
        index[k] = 0
        a = tuple(index)
        index[k] = 1
        b = tuple(index)
        # Inplace hadamard transform.
        term_matrix[a] += term_matrix[b]
        term_matrix[b] *= -2
        term_matrix[b] += term_matrix[a]
        # Scalar phase from Y.
        index[k + n] = 1
        term_matrix[tuple(index)] *= 1j
    term_matrix /= 2**n
    term_matrix.shape = (N, N)

    # Convert from dense matrix representation to sparse dict representation.
    terms: Dict[str, complex] = {}
    for pauli_tuple in itertools.product("IXYZ", repeat=n):
        pauli_string = ''.join(pauli_tuple)
        xk = sum(2 **k * (pauli_string[n-k-1] in 'XY') for k in range(n))
        zk = sum(2 **k * (pauli_string[n-k-1] in 'YZ') for k in range(n))
        coefficient = term_matrix[zk, xk]
        if coefficient:
            terms[pauli_string] = coefficient
    if not terms:
        terms["I" * n] = 0

    return terms

The funny thing is that this fancy looking transformation of the matrix actually corresponds to a very familiar quantum circuit. Flattening the matrix into vector, what's happening is equivalent to applying the CNOT-and-H prefix of a bunch of bell basis measurements! So basically you convert from matrix to Pauli terms by Bell basis measurement.

Answer 2

You can check that $$\mathrm{iSWAP}=\frac12(I\otimes I+iX\otimes X+iY\otimes Y + Z\otimes Z)$$