Implement Reverse Jordan-Wigner Transformation

[molar-volume]

What is the expected behavior?

While Jordan-Wigner transformation is implemented in the mappers framework (JordanWignerMapper), the reverse transformation is missing.

Inspiration can be drawn from https://github.com/quantumlib/OpenFermion/blob/master/src/openfermion/transforms/opconversions/reverse_jordan_wigner.py

PS: If maintainers agree with this issue, it will be solved by collaboration among Quantum Hackathon Korea participants.

[0sophy1]

Thank you for creating the issue! We are mentors of this event (Also Qiskit Advocates) - https://github.com/qiskit-community/quantum-hackathon-korea-21 and while handling the questions, we found there is no reverse Jordan Wigner function in Qiskit. We hope to contribute to Qiskit Nature by providing this function with our attendees.

[woodsp-ibm]

Do you have a use-case/example of situations where this reverse transformation would be useful for a particular algorithm or problem being solved?

[0sophy1]

it can be used for making anti-symmetric groups. also https://github.com/quantumlib/OpenFermion-FQE/blob/e4261fb7de19ae64a2a608e7c136503088976ca9/src/fqe/openfermion_utils.py this doc shows another use-case of this function in the openfermion side. In this time we've received a question because they want to reverse encode the computed qubit state into a fermionic state.

[BStar14]

@GyeonghumKim implemented mapper
https://github.com/molar-volume/qiskit-nature/blob/reverse_jw/qiskit_nature/mappers/second_quantization/reverse_jordan_wigner_mapper.py

and test file
https://github.com/molar-volume/qiskit-nature/blob/reverse_jw/test/mappers/second_quantization/test_reverse_jordan_wigner_mapper.py

We would improve these codes and consider making relevant examples.

[mrossinek]

#702 shows that we are in the process of migrating the FermionicOp into a new location.
As shown in #701 this will also relocate the mappers into a new location (namely qiskit_nature.second_quantization.operators.fermionic.mappers).

Maybe you can use this opportunity in a few weeks to open a PR which adds this reverse mapping? If you're still interested, I would suggest that you subscribe to #703 and once that is resolved, you add this mapper of yours to the aforementioned location. What do you think?

[0sophy1]

#702 shows that we are in the process of migrating the FermionicOp into a new location. As shown in #701 this will also relocate the mappers into a new location (namely qiskit_nature.second_quantization.operators.fermionic.mappers).

Maybe you can use this opportunity in a few weeks to open a PR which adds this reverse mapping? If you're still interested, I would suggest that you subscribe to #703 and once that is resolved, you add this mapper of yours to the aforementioned location. What do you think?

Hi @mrossinek why not! thank you so much! we will follow #703 and add our mapper!

[mrossinek]

Hi @0sophy1! It has been a while and Qiskit Nature just underwent yet another large refactoring (hopefully the last one as big as this for a long time).

Are you still interested in adding the reverse Jordan-Wigner mapping to Qiskit Nature?

I see that the implementation was basically already done (@BStar14). Instead of the version you proposed there I would suggest to add the reverse map to the existing JordanWignerMapper (because the QubitMapper class dictates a SparseLabelOp -> PauliSumOp direction which your mapper would not match). But that is a design detail which should be easy to figure out 👍

[MarcoBarroca]

Hey, I will likely need this feature in the near-ish future.

Would you guys mind if I tried adding it to the JordanWignerMappermyself?

[MarcoBarroca]

Was this implementation working previously? @BStar14 @GyeonghumKim @0sophy1

I was going through it and updating it to work with the recent updates. Mostly this means using PauliSparseOp and converting the result from dense labels to sparse labels so it's compatible with the new FermionicOp

It seems to map all the Pauli strings with Z correctly but it's missing terms on the way back. Only things I did were what I described above.

Here's a test script, there are some print statements commented for debugging:

from itertools import product
import math

import numpy as np

from qiskit.quantum_info.operators import Pauli
from qiskit.quantum_info import SparsePauliOp

from qiskit_nature.second_q.mappers.fermionic_mapper import FermionicMapper
from qiskit_nature.second_q.operators import FermionicOp

from collections import defaultdict 

def transform_dict_keys(data):
    label_transformation = {
        "I": "",
        "N": "+_{} -_{}",
        "E": "-_{} +_{}",
        "+": "+_{}",
        "-": "-_{}",
    }

    new_data = {}
    for key, value in data.items():
        new_key_parts = []
        for i, char in enumerate(key):
            if char in label_transformation:
                # Apply transformation if the character is in the mapping
                transformed = label_transformation[char].format(i, i)
                if transformed:
                    new_key_parts.append(transformed)
        new_key = ' '.join(new_key_parts)  # Join the parts with spaces
        new_data[new_key] = value

    return new_data
    
def reverse_map(second_q_op: SparsePauliOp) -> FermionicOp:
        """Maps a class:`OperatorBase` to a `FermionicOp`.

        Args:
            second_q_op: the :class:`OperatorBase` to be mapped.

        Returns:
            The `FermionicOp` corresponding to the Hamiltonian in the Fermionic space.
        """

        num_qubits = second_q_op.num_qubits # get number of qubits from input second quantized operator
        fermionic_op = None
        for term in second_q_op:
            transform_list : List[Tuple(str, float)] = []  # list of tuple(pauli, coeff)
            coef_term = term.coeffs[0]  # Assuming coeffs is a list with a single element
            target_pauli_op = term.paulis[0]
            #print(target_pauli_op)

            for i in range(num_qubits):
                one_pauli = target_pauli_op[num_qubits - 1 - i]
                pauli_char=one_pauli.to_label()
                if pauli_char == 'Z': # dealing Pauli Z op
                    transform_list.append((('I',1),('N', -2))) # Zj -> I - 2*Nj => [('I', 1), ('N', -2)]
                elif pauli_char == 'X': # dealing Pauli X op
                    transform_list.append((("+", 1), ('-', 1))) # Xj -> aj_dag + aj => [('+', 1), ('-', 1)]
                    target_pauli_op &= Pauli("I" * (i+1) + "Z" * (num_qubits - i - 1)) # apply Z(j-1)Z(j-2) ... Z(0)
                elif one_pauli.to_label() == 'Y': # dealing Pauli Y op
                    transform_list.append((('+', -1j), ('-', 1j))) # Yj -> i(aj - aj_dag) => [('+', -1j), ('-', 1j)]
                    target_pauli_op &= Pauli("I" * (i+1) + "Z" * (num_qubits - i - 1)) # apply Z(j-1)Z(j-2) ... Z(0)
                else: 
                    # dealing Pauli I op
                    transform_list.append((('I', 0.5), ('I', 0.5))) # Ij -> Ij => [('I', 0.5), ('I', 0.5)]; split I into 0.5I + 0.5I for code consistency
            
            # dealing the phase
            if target_pauli_op.phase == 1:
                coef_term *= -1j
            elif target_pauli_op.phase == 2:
                coef_term *= -1
            elif target_pauli_op.phase == 3:
                coef_term *= 1j

            #print(transform_list)
            pauli_coefs = []
            pauli_strings = []
            # create fermionic operator for a term based on transform_list
            for idxes in product(*[[0, 1]]*num_qubits):
                pauli_coefs.append(math.prod([t[i][1] for t, i in zip(transform_list, idxes)]))
                pauli_strings.append("".join([t[i][0] for t, i in zip(transform_list, idxes)])[::-1])
            ferm_op=list(zip(pauli_strings, pauli_coefs))

            ferm_dict = defaultdict(float)

            # Iterate over each tuple in the list
            for key, value in ferm_op:
                # Sum the float values for tuples with the same string
                ferm_dict[key] += value

            ferm_dict_sparse=transform_dict_keys(ferm_dict)

            #print(ferm_dict)
            #print(ferm_dict_sparse)
            #print(coef_term)

            if not fermionic_op:
                fermionic_op = coef_term * FermionicOp(ferm_dict_sparse,num_spin_orbitals=num_qubits).simplify()
            else:
                fermionic_op += coef_term * FermionicOp(ferm_dict_sparse,num_spin_orbitals=num_qubits).simplify()

        return fermionic_op.simplify()

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver

driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.735",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

from IPython.display import display

problem = driver.run()

hamiltonian = problem.hamiltonian
second_q_op = hamiltonian.second_q_op()

mapper=JordanWignerMapper()
qubit_op=mapper.map(second_q_op)

reverse_op=reverse_map(qubit_op)

display("Fermioinic OP:")
display(print(second_q_op))
display("Qubit_OP:")
display(print(qubit_op))
display("Reversed_OP:")
display(print(reverse_op))

And the output:

'Fermioinic OP:'

Fermionic Operator
number spin orbitals=4, number terms=28
  (0.3322908651276482+0j) * ( +_0 +_1 -_1 -_0 )
+ (0.3378550774017583+0j) * ( +_0 +_2 -_2 -_0 )
+ (0.3322908651276482+0j) * ( +_0 +_3 -_3 -_0 )
+ (0.09046559989211572+0j) * ( +_0 +_1 -_0 -_1 )
+ (0.09046559989211572+0j) * ( +_0 +_2 -_3 -_1 )
+ (0.09046559989211572+0j) * ( +_0 +_3 -_2 -_1 )
+ (0.09046559989211572+0j) * ( +_1 +_0 -_1 -_0 )
+ (0.09046559989211572+0j) * ( +_1 +_2 -_3 -_0 )
+ (0.09046559989211572+0j) * ( +_1 +_3 -_2 -_0 )
+ (0.3322908651276482+0j) * ( +_1 +_0 -_0 -_1 )
+ (0.3322908651276482+0j) * ( +_1 +_2 -_2 -_1 )
+ (0.34928686136600906+0j) * ( +_1 +_3 -_3 -_1 )
+ (0.3378550774017583+0j) * ( +_2 +_0 -_0 -_2 )
+ (0.3322908651276482+0j) * ( +_2 +_1 -_1 -_2 )
+ (0.3322908651276482+0j) * ( +_2 +_3 -_3 -_2 )
+ (0.09046559989211572+0j) * ( +_2 +_0 -_1 -_3 )
+ (0.09046559989211572+0j) * ( +_2 +_1 -_0 -_3 )
+ (0.09046559989211572+0j) * ( +_2 +_3 -_2 -_3 )
+ (0.09046559989211572+0j) * ( +_3 +_0 -_1 -_2 )
+ (0.09046559989211572+0j) * ( +_3 +_1 -_0 -_2 )
+ (0.09046559989211572+0j) * ( +_3 +_2 -_3 -_2 )
+ (0.3322908651276482+0j) * ( +_3 +_0 -_0 -_3 )
+ (0.34928686136600906+0j) * ( +_3 +_1 -_1 -_3 )
+ (0.3322908651276482+0j) * ( +_3 +_2 -_2 -_3 )
+ (-1.25633907300325+0j) * ( +_0 -_0 )
+ (-0.47189600728114184+0j) * ( +_1 -_1 )
+ (-1.25633907300325+0j) * ( +_2 -_2 )
+ (-0.47189600728114184+0j) * ( +_3 -_3 )

None

'Qubit_OP:'

SparsePauliOp(['IIII', 'IIIZ', 'IIZI', 'IIZZ', 'IZII', 'IZIZ', 'ZIII', 'ZIIZ', 'YYYY', 'XXYY', 'YYXX', 'XXXX', 'IZZI', 'ZIZI', 'ZZII'],
              coeffs=[-0.81054798+0.j,  0.17218393+0.j, -0.22575349+0.j,  0.12091263+0.j,
  0.17218393+0.j,  0.16892754+0.j, -0.22575349+0.j,  0.16614543+0.j,
  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,
  0.16614543+0.j,  0.17464343+0.j,  0.12091263+0.j])

None

'Reversed_OP:'

Fermionic Operator
number spin orbitals=4, number terms=14
  (-1.25633907300325+0j) * ( +_0 -_0 )
+ (-0.4718960072811419+0j) * ( +_1 -_1 )
+ (0.483650530471065+0j) * ( +_0 -_0 +_1 -_1 )
+ (-1.2563390730032502+0j) * ( +_2 -_2 )
+ (0.6757101548035166+0j) * ( +_0 -_0 +_2 -_2 )
+ (-0.4718960072811419+0j) * ( +_3 -_3 )
+ (0.6645817302552964+0j) * ( +_0 -_0 +_3 -_3 )
+ (0.18093119978423144+0j) * ( -_0 +_1 -_2 +_3 )
+ (-0.18093119978423144+0j) * ( +_0 -_1 -_2 +_3 )
+ (0.18093119978423144+0j) * ( +_0 -_1 +_2 -_3 )
+ (-0.18093119978423144+0j) * ( -_0 +_1 +_2 -_3 )
+ (0.6645817302552964+0j) * ( +_1 -_1 +_2 -_2 )
+ (0.6985737227320181+0j) * ( +_1 -_1 +_3 -_3 )
+ (0.483650530471065+0j) * ( +_2 -_2 +_3 -_3 )

None

[BStar14]

@MarcoBarroca, this issue was suspended because we couldn't find its usecase. I would be grateful if you could share your idea. It is also fine to update the code on your own. Maybe @GyeonghunKim can provide more help for it.

[MarcoBarroca]

I've modified the implementation to use the FermionicOp arithmetic instead of doing so only at the end. Seems simpler this way and I've recovered the same results from the original implementation. Still looks like it's missing something for Operators X and Y as the reverse still isn't the same as the original.

from itertools import product
import math

import numpy as np

from qiskit.quantum_info.operators import Pauli
from qiskit.quantum_info import SparsePauliOp

from qiskit_nature.second_q.mappers import JordanWignerMapper
from qiskit_nature.second_q.operators import FermionicOp

def invert_pauli_terms(sparse_pauli_op: SparsePauliOp) -> SparsePauliOp:
    """Inverts the order of Pauli operators in each term of a SparsePauliOp."""
    inverted_labels = [label[::-1] for label in sparse_pauli_op.paulis.to_labels()]
    # Create a new SparsePauliOp with the inverted labels but same coefficients
    inverted_sparse_pauli_op = SparsePauliOp(inverted_labels, sparse_pauli_op.coeffs)
    return inverted_sparse_pauli_op
    
def reverse_map(qubit_op: SparsePauliOp) -> FermionicOp:
        """Maps a class:`OperatorBase` to a `FermionicOp`.

        Args:
            second_q_op: the :class:`OperatorBase` to be mapped.

        Returns:
            The `FermionicOp` corresponding to the Hamiltonian in the Fermionic space.
        """

        num_qubits = qubit_op.num_qubits # get number of qubits from input second quantized operator
        qubit_op=invert_pauli_terms(qubit_op)
        total_fermionic_op = FermionicOp.zero()
        for term in qubit_op:
            coef_term = term.coeffs[0]  # Assuming coeffs is a list with a single element
            target_pauli_op = term.paulis[0]
            #print(target_pauli_op)
            #print(coef_term)
            ferm_term_ops=[]
            z_count=0
            for i in range(num_qubits):
                transform_dict={}
                parity_dict={}
                one_pauli = target_pauli_op[num_qubits - 1 - i]
                pauli_char=one_pauli.to_label()
                #print(pauli_char)
                if pauli_char == 'Z': # dealing Pauli Z op
                    transform_dict[f'']=1
                    transform_dict[f'+_{i} -_{i}']=-2
                    ferm_op_pauli=FermionicOp(transform_dict)
                elif pauli_char == 'X': # dealing Pauli X op
                    transform_dict[f'+_{i}']=1
                    transform_dict[f'-_{i}']=1
                    ferm_op_pauli=FermionicOp(transform_dict)
                    target_pauli_op &= Pauli("I" * (i+1) + "Z" * (num_qubits - i - 1)) # apply Z(j-1)Z(j-2) ... Z(0)
                elif one_pauli.to_label() == 'Y': # dealing Pauli Y op
                    phase_factor = (-1j) ** z_count
                    transform_dict[f'+_{i}']=1j
                    transform_dict[f'-_{i}']=-1j
                    ferm_op_pauli=FermionicOp(transform_dict)
                    target_pauli_op &= Pauli("I" * (i+1) + "Z" * (num_qubits - i - 1)) # apply Z(j-1)Z(j-2) ... Z(0)
                else:
                    ferm_op_pauli=FermionicOp.one()
                ferm_term_ops.append(ferm_op_pauli)
                #print(ferm_term_ops)
                
            term_fermionic_op = FermionicOp.one()
            for op in ferm_term_ops:
                term_fermionic_op = term_fermionic_op @ op
            #print(term_fermionic_op)

            #print(target_pauli_op)

            # dealing the phase
            if target_pauli_op.phase == 1:
                coef_term *= -1j
            elif target_pauli_op.phase == 2:
                coef_term *= -1
            elif target_pauli_op.phase == 3:
                coef_term *= 1j
            
            total_fermionic_op += coef_term * term_fermionic_op
        
            #print(coef_term * term_fermionic_op)
            #print(total_fermionic_op)
            

        return total_fermionic_op.simplify()

from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.drivers import PySCFDriver

driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.735",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
)

from IPython.display import display

problem = driver.run()

hamiltonian = problem.hamiltonian
second_q_op = hamiltonian.second_q_op()

mapper=JordanWignerMapper()
qubit_op=mapper.map(second_q_op)

reverse_op=reverse_map(qubit_op)

display("Fermioinic OP:")
display(print(second_q_op.simplify().simplify()))
display("Qubit_OP:")
display(print(qubit_op))
display("Reversed_OP:")
display(print(reverse_op))

Output:

'Fermioinic OP:'

Fermionic Operator
number spin orbitals=4, number terms=28
  (0.3322908651276482+0j) * ( +_0 +_1 -_1 -_0 )
+ (0.3378550774017583+0j) * ( +_0 +_2 -_2 -_0 )
+ (0.3322908651276482+0j) * ( +_0 +_3 -_3 -_0 )
+ (0.09046559989211572+0j) * ( +_0 +_1 -_0 -_1 )
+ (0.09046559989211572+0j) * ( +_0 +_2 -_3 -_1 )
+ (0.09046559989211572+0j) * ( +_0 +_3 -_2 -_1 )
+ (0.09046559989211572+0j) * ( +_1 +_0 -_1 -_0 )
+ (0.09046559989211572+0j) * ( +_1 +_2 -_3 -_0 )
+ (0.09046559989211572+0j) * ( +_1 +_3 -_2 -_0 )
+ (0.3322908651276482+0j) * ( +_1 +_0 -_0 -_1 )
+ (0.3322908651276482+0j) * ( +_1 +_2 -_2 -_1 )
+ (0.34928686136600906+0j) * ( +_1 +_3 -_3 -_1 )
+ (0.3378550774017583+0j) * ( +_2 +_0 -_0 -_2 )
+ (0.3322908651276482+0j) * ( +_2 +_1 -_1 -_2 )
+ (0.3322908651276482+0j) * ( +_2 +_3 -_3 -_2 )
+ (0.09046559989211572+0j) * ( +_2 +_0 -_1 -_3 )
+ (0.09046559989211572+0j) * ( +_2 +_1 -_0 -_3 )
+ (0.09046559989211572+0j) * ( +_2 +_3 -_2 -_3 )
+ (0.09046559989211572+0j) * ( +_3 +_0 -_1 -_2 )
+ (0.09046559989211572+0j) * ( +_3 +_1 -_0 -_2 )
+ (0.09046559989211572+0j) * ( +_3 +_2 -_3 -_2 )
+ (0.3322908651276482+0j) * ( +_3 +_0 -_0 -_3 )
+ (0.34928686136600906+0j) * ( +_3 +_1 -_1 -_3 )
+ (0.3322908651276482+0j) * ( +_3 +_2 -_2 -_3 )
+ (-1.25633907300325+0j) * ( +_0 -_0 )
+ (-0.47189600728114184+0j) * ( +_1 -_1 )
+ (-1.25633907300325+0j) * ( +_2 -_2 )
+ (-0.47189600728114184+0j) * ( +_3 -_3 )

None

'Qubit_OP:'

SparsePauliOp(['IIII', 'IIIZ', 'IIZI', 'IIZZ', 'IZII', 'IZIZ', 'ZIII', 'ZIIZ', 'YYYY', 'XXYY', 'YYXX', 'XXXX', 'IZZI', 'ZIZI', 'ZZII'],
              coeffs=[-0.81054798+0.j,  0.17218393+0.j, -0.22575349+0.j,  0.12091263+0.j,
  0.17218393+0.j,  0.16892754+0.j, -0.22575349+0.j,  0.16614543+0.j,
  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,
  0.16614543+0.j,  0.17464343+0.j,  0.12091263+0.j])

None

'Reversed_OP:'

Fermionic Operator
number spin orbitals=4, number terms=14
  (-1.25633907300325+0j) * ( +_0 -_0 )
+ (-0.4718960072811419+0j) * ( +_1 -_1 )
+ (0.483650530471065+0j) * ( +_0 -_0 +_1 -_1 )
+ (-1.2563390730032502+0j) * ( +_2 -_2 )
+ (0.6757101548035166+0j) * ( +_0 -_0 +_2 -_2 )
+ (-0.4718960072811419+0j) * ( +_3 -_3 )
+ (0.6645817302552964+0j) * ( +_0 -_0 +_3 -_3 )
+ (0.18093119978423144+0j) * ( -_0 +_1 -_2 +_3 )
+ (-0.18093119978423144+0j) * ( -_0 +_1 +_2 -_3 )
+ (0.18093119978423144+0j) * ( +_0 -_1 +_2 -_3 )
+ (-0.18093119978423144+0j) * ( +_0 -_1 -_2 +_3 )
+ (0.6645817302552964+0j) * ( +_1 -_1 +_2 -_2 )
+ (0.6985737227320181+0j) * ( +_1 -_1 +_3 -_3 )
+ (0.483650530471065+0j) * ( +_2 -_2 +_3 -_3 )

None

[MarcoBarroca]

Found the issue! Both implementations work fine! Just do a print(second_q_op.normal_order()-reverse_op.normal_order()) and you'll see they match up to e-17!

Fermionic Operator
number spin orbitals=4, number terms=14
  0j * ( +_0 +_1 -_0 -_1 )
+ 0j * ( +_0 +_2 -_0 -_2 )
+ 0j * ( +_0 +_3 -_0 -_3 )
+ 0j * ( +_0 +_2 -_1 -_3 )
+ 0j * ( +_0 +_3 -_1 -_2 )
+ 0j * ( +_1 +_2 -_0 -_3 )
+ 0j * ( +_1 +_3 -_0 -_2 )
+ 0j * ( +_1 +_2 -_1 -_2 )
+ 0j * ( +_1 +_3 -_1 -_3 )
+ 0j * ( +_2 +_3 -_2 -_3 )
+ 0j * ( +_0 -_0 )
+ (5.551115123125783e-17+0j) * ( +_1 -_1 )
+ (2.220446049250313e-16+0j) * ( +_2 -_2 )
+ (5.551115123125783e-17+0j) * ( +_3 -_3 )

second_q_op.normal_order().equiv(reverse_op.normal_order()) also returns True.

I'll work to make a reverse_map() method in the JordanWignerMapper() that implements this and submit a PR. I'll probably use the second implementation I posted just because it seems easier to read. I'll try to include a unit test as well in the test_jordan_wigner_mapper.

[MarcoBarroca]

Done!

I think reverse JordanWigner will be enough for what I’m working on but I wonder if there shouldn’t be a ReverseQubitMapper() or modify/rename FermionicMapper() class in the future to handle reverse transformations for all mappers.

One other question. Why doesn’t FermionicOp.equiv() put all the operators in .normal_mode() before checking equivalence? I took a while to figure out I had recovered the same operator already.