Tour

First steps

Let's start by constructing a Pauli polynomial.

>>> from fastfermion import poly
>>> A = poly("X0 X1 + (.3 X0 Z1 + .4 Y0 Y1) Z2")

The above expression creates a Pauli polynomial with three terms. One can get the number of terms in \(A\) using the len function, and individual monomial coefficients using the coefficient function.

>>> len(A)
3
>>> A.coefficient("X0 Z1 Z2")
(0.3+0j)

Polynomials can be added and multiplied using the usual operations. Several properties of a polynomial can be accessed via functions such as degree, norm, coefficient, etc.

>>> A += 2
>>> print(A)
X0 X1  +  0.3 X0 Z1 Z2  +  0.4 Y0 Y1 Z2  +  2
>>> A.norm('inf')
2.0
>>> A.degree()
3
>>> A.degree("Y")
2

To construct \(A\) we used the function poly which parses a string into a fastfermion polynomial. Alternatively, one can use paulis to get the generators of the Pauli algebra and use standard Python operations to build up a polynomial. For example, to construct the Heisenberg Hamiltonian on a line, one can proceed as follows.

>>> from fastfermion import paulis
>>> L = 4
>>> X,Y,Z = paulis(L)
>>> edges = [(i,i+1) for i in range(L-1)]
>>> H = sum(X[i]*X[j] + Y[i]*Y[j] + Z[i]*Z[j] for i,j in edges)
>>> print(H)
X0 X1  +  Y0 Y1  +  Z0 Z1  +  X1 X2  +  Y1 Y2  +  Z1 Z2  +  X2 X3  +  Y2 Y3  +  Z2 Z3

Importantly, fastfermion also supports fermionic polynomials in creation/annihilation operators. Fermionic creation/annihilation operators \(f_p^{\dagger}\), \(f_p\) obey the canonical anticommutation relations:

\[ \{f_p,f_q\} = 0, \quad \{f_p^{\dagger}, f_q\} = \delta_{pq}. \]

We can also use the poly function to input the polynomial as a string.

>>> from fastfermion import poly
>>> B = poly("f0 f0^ + 3.1 f2^ f0 + 1")
>>> print(B)
- f0^ f0  +  2  +  3.1 f2^ f0
>>> B.coefficient("f2^ f0")
(3.1+0j)

The notation f0^ indicates a creation operator and f0 (without a trailing ^) indicates an annihilation operator. Note that the expression \(f_0f_0^{\dagger}\) was automatically changed to \(1-f_0^{\dagger} f_0\) in the internal representation of \(B\). This is because fastfermion stores all fermionic polynomials in normal ordered form: all the creation operators appear to the left of annihilation operators, and creation and annihilation operators are ordered in decreasing order from left to right.

Fermionic polynomials can also be constructed by first getting generators of the algebra (annihilation operators) using the fermis command, and then building up a polynomial. The code below creates a quadratic polynomial \(\sum_{ij} f_i^{\dagger} f_{i+1}\):

>>> from fastfermion import fermis
>>> L = 5
>>> edges = [(i,i+1) for i in range(L-1)]
>>> f = fermis(L)
>>> H = sum(f[i].dagger() * f[j] for i,j in edges)
>>> print(H)
f0^ f1  +  f1^ f2  +  f2^ f3  +  f3^ f4

Finally, fastfermion also supports polynomials in Majorana operators \(m_p\), which are Hermitian operators that obey the following anticommutation relation:

\[ \{ m_p,m_q \} = 2 \delta_{pq}. \]

One can construct Majorana polynomials in the same way we constructed Pauli polynomials or Fermi polynomials.

>>> from fastfermion import poly, majoranas
>>> C = poly("(m0 m1 m4 + 0.2 m2 m1) m4")
>>> print(C)
m0 m1  -  0.2 m1 m2 m4
>>> m = majoranas(5)
>>> D = (m[0]*m[1]*m[4] + 0.2*m[2]*m[1])*m[4]
>>> C-D
0

Polynomials in Majorana operators are also automatically put in normal ordered form. In this case, normal form means that Majorana modes appear in increasing order from left to right. This is why the second term in \(C\), which was input \(m_2 m_1 m_4\) is represented as \(-m_1 m_2 m_4\) instead.

Transforms

fastfermion offers the possibility to transform polynomials from one type {Pauli, Fermi, Majorana} to any other type. The well-known Jordan-Wigner transform transforms polynomials in Fermionic creation/annihilation operators to a polynomial in Pauli operators using the rule:

\[ \begin{aligned} f_j &\leftarrow Z_0 \dots Z_{j-1} (X_j + \text{i} Y_j) / 2\\ f_j^{\dagger} &\leftarrow Z_0 \dots Z_{j-1} (X_j - \text{i} Y_j) / 2. \end{aligned} \]

The Jordan-Wigner transform can be computed via the function jw. Furthermore the reverse Jordan-Wigner transform can be computed using the function rjw.

>>> from fastfermion import poly, jw, rjw
>>> B = poly("f0 f0^ + 3.1 f2^ f0 + 1")
>>> jw(B)
1.5  +  0.5 Z0  +  0.775j Y0 Z1 X2  +  0.775 X0 Z1 X2  +  0.775 Y0 Z1 Y2  -  0.775j X0 Z1 Y2
>>> rjw(jw(B)) - B
0

One can also convert any Fermi polynomial to a Majorana polynomial via the following transform

\[ \begin{aligned} f_j &\leftarrow (m_{2j} + \text{i} m_{2j+1})/2\\ f_j^{\dagger} &\leftarrow (m_{2j} - \text{i} m_{2j+1}) / 2. \end{aligned} \]

>>> B.tomajorana()
1.5  -  0.5j m0 m1  -  0.775 m0 m4  -  0.775j m1 m4  +  0.775j m0 m5  -  0.775 m1 m5

The functions topauli(), tofermi(), tomajorana() can be used to go from one type of polynomial to another.

Sparse matrix representations

Given any polynomial in fastfermion, one can use the sparse function to get a SciPy sparse matrix representation of the corresponding polynomial.

>>> from fastfermion import poly
>>> A = poly("X0*X1 + Z0 + 3")
>>> A.sparse()
<Compressed Sparse Column sparse matrix of dtype 'complex128'
        with 8 stored elements and shape (4, 4)>
>>> A.sparse().toarray()
array([[4.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],
       [0.+0.j, 4.+0.j, 1.+0.j, 0.+0.j],
       [0.+0.j, 1.+0.j, 2.+0.j, 0.+0.j],
       [1.+0.j, 0.+0.j, 0.+0.j, 2.+0.j]])

One can check that the above matrix is indeed equal to \(\sigma_x \otimes \sigma_x + \sigma_z \otimes I_2 + 3 I_2 \otimes I_2\), where \(\sigma_x\), \(\sigma_y\) and \(\sigma_z\) are the standard \(2\times 2\) Pauli matrices.

The sparse() function can also be used for Fermi polynomials and Majorana polynomials.

>>> from fastfermion import poly
>>> poly("f1^ f1 f0^ - 2 f1 f1^ + 3").sparse()
<Compressed Sparse Column sparse matrix of dtype 'complex128'
        with 5 stored elements and shape (4, 4)>
>>> poly("m0 m1 - m2 m0").sparse()
<Compressed Sparse Column sparse matrix of dtype 'complex128'
        with 8 stored elements and shape (4, 4)>

Note that a Fermi polynomial with \(n\) modes is represented (in the occupation number basis) as a \(2^n \times 2^n\) matrix, and a Majorana polynomial on \(n\) variables is represented as a \(2^{\lceil n/2 \rceil} \times 2^{\lceil n/2 \rceil}\) matrix.

Subspaces

Often, we may be interested in the sparse matrix representation of a polynomial, restricted to a certain subspace left invariant by the operator. A typical example for Pauli polynomials is the subspace corresponding to a fixed spin value. Getting the corresponding block is supported by fastfermion.

>>> from fastfermion import paulis
>>> n = 10
>>> sx, sy, sz = paulis(n)
>>> edge_list = [(i,(i+1)%n) for i in range(n)]
>>> H = sum(sx[i]*sx[j] + sy[i]*sy[j] + sz[i]*sz[j] for i,j in edge_list) # Heisenberg hamiltonian on a line
>>> H.sparse()
<Compressed Sparse Column sparse matrix of ... shape (1024, 1024)>
>>> H.sparse(nup=5)
<Compressed Sparse Column sparse matrix of ... shape (252, 252)>

With the parameter nup, the function sparse will only output the block corresponding to the subspace of states having nup spins equal to +1. The dimension of this subspace is \(\binom{n}{n_{up}}\).

Note: If the supplied polynomial does not have the required symmetry, an error is raised. For example, poly("X0+X1").sparse(nup=1) raises an error since the operator \((X_0+X_1)\) does not leave the appropriate subspace invariant.

For Fermi polynomials, one can supply the total number of particles to the sparse function to return only the block corresponding to the supplied total occupation number. Of course, if the Fermi polynomial does not preserve total particle number, an error is raised.

>>> H = poly("f1^ f0 + .2 * f1^ f0^ f2 f1")
>>> H.sparse()
<Compressed Sparse Column sparse matrix of ... shape (8, 8)>
>>> H.sparse(nocc=1)
<Compressed Sparse Column sparse matrix of ... shape (3, 3)>

If the Fermi polynomial is defined on \(n\) modes, then H.sparse(nocc=k) returns a sparse matrix of size \(\binom{n}{k} \times \binom{n}{k}\).

Unitaries and evolution

A unitary is any operator \(U\) such that \(U^{\dagger} U = 1\). Such a unitary acts on polynomials by conjugation: \(P \to U^{\dagger} P U\). fastfermion currently supports a limited number of built-in unitaries for Pauli polynomials and Majorana polynomials. For example:

>>> from fastfermion import poly, CNOT
>>> A = poly("X0 X1 + (.3 X0 Z1 + .4 Y0 Y1) Z2")
>>> CNOT(0,1)(A) # Apply the CNOT(0,1) unitary to A
X0  -  0.3 Y0 Y1 Z2  -  0.4 X0 Z1 Z2

A circuit is an ordered collection of unitaries \(C = (U_1,\ldots,U_d)\) where \(d\) is the depth of the circuit. Starting from a polynomial \(P\), the evolution of \(P\) through the circuit is given by

\[ U_1^{\dagger} \dots U_d^{\dagger} P U_d \dots U_1. \]

Pauli propagation

The propagate function is a direct method to propagate a Pauli polynomial through a circuit.

>>> from fastfermion import poly, CNOT, H, S, ROT, propagate
>>> P = poly('X0')
>>> circuit = [CNOT(0,1), H(1), ROT("ZZ",[0,1],0.785)]
>>> Q = propagate(circuit, P)
>>> Q
0.7073882691671998 X0 X1 - 0.706825181105366 Y0

The propagate function is functionally equivalent to the following simple for loop:

def propagate(circuit, poly):
    for gate in reversed(circuit):
        poly = gate(poly)
    return poly

For large circuits, one can specify an additional parameter to propagate that truncates the degree of the polynomial after the action of non-Clifford gates (namely ROT gates).

>>> S = propagate(circuit, P, maxdegree=3)

Majorana propagation

The propagate function can also be used to propagate an arbitrary Majorana polynomial through a sequence of Majorana rotations. A Majorana rotation is a unitary of the form

\[ U = e^{-i \theta/2 M} \]

where \(M\) is a Hermitian Majorana monomial of the form \(M = c m_{i_1} \dots m_{i_k}\) where \(i_1 < \dots < i_k\) and \(c\) is real if \(m_{i_1} \dots m_{i_k}\) is Hermitian, otherwise it is pure imaginary. Majorana rotations can be constructed by the MROT function.

>>> from fastfermion import MROT, MajoranaString
>>> U = MROT(MajoranaString([0,1]),1j,1.0)
>>> print(U)
MROT(1j m0 m1, theta=1.000000) = e^{0.500000 m0 m1}

A Majorana circuit is a sequence of Majorana rotations.

>>> circuit = [MROT(MajoranaString([0,1]),1j,0.1) , MROT(MajoranaString([0,1,2,3]),1,0.1)]
>>> obs = poly("m0 m2")
>>> Q = propagate(circuit, obs, maxdegree=3)

Interface with OpenFermion and Cirq

The functions from_openfermion and to_openfermion can be used to convert from/to OpenFermion operators. Also the function from_cirq can be used to convert a Cirq circuit to a list of FastFermion gates.

>>> import cirq
>>> from fastfermion import from_cirq
>>> q0, q1 = cirq.LineQubit.range(2)
>>> cirq_circuit = cirq.Circuit(cirq.H(q0),cirq.CNOT(q0, q1), cirq.Rx(rads=0.23)(q0))
>>> ff_circuit = from_cirq(cirq_circuit)
>>> ff_circuit
[H(0), CNOT(0,1), ROT(X0,0.230000)]

Currently, only the Cirq gates cirq.H, cirq.S, cirq.CZ, cirq.SWAP, cirq.CNOT, cirq.Rx, cirq.Ry, cirq.Ry are supported.