unitaryHACK #1: Faster to_sparse_matrix method

[TimWeaving]

Goal: Write a faster to_sparse_matrix method for the PauliwordOp class in symmer.

Reward: $200

Details:

One of the core features of the symmer codebase is to provide a symbolic representation of the Pauli operators:

$$X =\begin{pmatrix} 0 & 1 \ 1 & 0 \end{pmatrix} $$

$$Y =\begin{pmatrix} 0 & -1i \ 1i & 0 \end{pmatrix} $$

$$Z =\begin{pmatrix} 1 & 0 \ 0 & -1 \end{pmatrix} $$

$$I =\begin{pmatrix} 1 & 0 \ 0 & 1 \end{pmatrix} $$

The PauliwordOp class gives different functions that can be applied on these operators and tensor products of them. The individual operators can be written as:

from symmer import PauliwordOp

# dicitonary: {string: coefficient}
X = PauliwordOp.from_dictionary({"X":1})
Y = PauliwordOp.from_dictionary({"Y":1})
Z = PauliwordOp.from_dictionary({"Z":1})
I = PauliwordOp.from_dictionary({"I":1})

print(X)
print(Y)
print(Z)
print(I)

We can also represent tensor products of these Pauli operators as:

from symmer import PauliwordOp
# dicitonary: {string: coefficient}
XYZ = PauliwordOp.from_dictionary({"XYZ":1})
print(XYZ)

This object represents:

$$XYZ = X \otimes Y \otimes Z =\begin{pmatrix} 0 & 1 \ 1 & 0 \end{pmatrix} \otimes \begin{pmatrix} 0 & -1i \ 1i & 0 \end{pmatrix} \otimes \begin{pmatrix} 1 & 0 \ 0 & -1 \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0& -1j& 0\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & +1j \ 0 & 0 & 0 & 0 & +1j & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & -1j & 0 & 0\ 0 & 0 & -1j & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & +1j & 0 & 0 & 0 & 0 \ +1j & 0 & 0 & 0 & 0 & 0 & 0 & 0\ 0 & -1j & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix} $$

and can be generated via:

from symmer import PauliwordOp
# dicitonary: {string: coefficient}
XYZ = PauliwordOp.from_dictionary({"XYZ":1})
print(XYZ.to_sparse_matrix.toarray())

Given the structure of Pauli operators $X, Y, Z, I$ rather than finding the kronecker product of each operator:

from scipy.linalg import kron
import numpy as np
from functools import reduce
X = np.array([[0, 1],
                       [1, 0]])
Y = np.array([[ 0., -1j],
                    [ 1j,  0.]])
Z = np.array([[1, 0],
                      [0, -1]])
I = np.eye(2)

XYZ = reduce(kron, [X,Y,Z])
print(XYZ)

We use a to_sparse_matrix method, that determines the indices of the output matrix from the structure of the operators. This increases the speed in finding the matrix respresentation of the operator over manually taking the individual kronecker products. However, there are still improvements that can be made hence this pull request. We want someone to speed up this functionality.

One approach is outlined in a qiskit issue: https://github.com/Qiskit/qiskit-terra/issues/8772 , that points to: https://github.com/chetmurthy/qrusty. This codebase is written in rust. It would be great to have someone convert this rust implementation into a python version for symmer.

Alternatively, if someone can find a different method that is faster than the current one this would also be a valid solution!

The winning PR will be awarded $200.

[metapunkkk]

Hi, could you please assign this task to me?

[AlexisRalli]

Hi @Nizdew! Of course you can work on this issue. Good luck and feel free to post any questions here.