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=(2n, 2n), 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.
