Operator.compose should allow indices

[molar-volume]

resolves #1088
resolves #1143

This PR is created to verify the solution for expanding operator in #1088 (permute part of issue is not addressed yet)

The method expand for expanding the operator with identity is defined in OperatorBase,
because it enhances compose (defined for OperatorBase),
and uses tensor (defined for OperatorBase), and new method identity (abstract in OperatorBase).

identity is implemented for each subclass of 'OperatorBase',

but I am unsure of implementation for DictStateFn and VectorStateFn (instead of identity, zeros are used).

The new features are demonstrated intest_expand_on_state_fn and test_compose_op_of_different_dim

Summary

Progress for #1088 and solution for #1143

Details and comments

[Cryoris]

Thanks for the commitment @molar-volume! I think we need to clarify what #1088 should do and how we could do it. The idea was that you could pass concrete indices (not only ops of different sizes) and e.g. do

op = (X ^ 5)
op.compose( Z ^ Y, qubits=[0, 3], inplace=True)
print(op)  # ComposedOp( [X ^ X ^ X ^ X ^ X], [Z ^ I ^ I ^ Y ^ I] )

So it is closer to the same method of the quantum circuit, see QuantumCircuit.compose. (You could also add the front feature if you want! 🙂 ).

The approach in here does only allow for consecutive qubit indices (e.g. [1,2,3]) not for others (such as [0, 3]) if I understand it correctly. I have two concerns about the methods you added:

• expand already exists in the quantum_info operators but has a different signature. There you can call Operator.expand(other_operator_of_same_dim), I think we should keep the same syntax if possible and therefore not call it expand here. If this turns out to be necessary we could possibly call it expand_to_dim or pad_with_identities or so.
• Adding an identity method is maybe not the best approach since it is not clear to me what you would expect upon calling CircuitOp.identity? On a circuit this adds an identity gate, here it does something different. You could hide this inside pad_with_identities.

A possible approach would be to add a pad_with_identities (or whatever you want to call it) method looking like

def pad_with_identities(self, new_dim, qargs=None):
    """Extend this operator to ``new_dim`` where this operator acts on ``qargs``.

    Args:
        new_dim: The new dimension, must be larger equal than ``self.num_qubits``.
        qargs: The qubit indices this operator acts on, must be a list of length ``self.num_qubits``.
    """

And then in compose you can use the above to expand other to the correct size and compose it. What do you think?

[molar-volume]

Hi @Cryoris thank you for the feedback !
Yes, let’s rename expand to expand_to_dim.

Now identity (or we should rather name it identity_operator) is OperatorBase with primitive corresponding to identity matrix, so for CircuitOp, it returns (I ^ num_qubits).to_circuit_op() .

I think it has some advantages, for example you can define expand_to_dim for all subclasses in OperatorBase as

def expand_to_dim(self, num_qubits: int) -> OperatorBase:
	return self.tensor(self.identity_operator(num_qubits))

I understand that we need a more general method with indices, just like you defined

def pad_with_identities(self, new_dim, qargs=None)

For the CircuitOp, that functionality can be already provided using permute (as suggested by @dongreenberg )

circ_op = (X ^ Y).to_circuit_op().permute([0,2]) @ (X ^ Y ^ Z ^ X).to_circuit_op()

print(circ_op.to_circuit())
     ┌───┐     
q_0: ┤ X ├─────
     ├───┤┌───┐
q_1: ┤ Z ├┤ Y ├
     ├───┤└───┘
q_2: ┤ Y ├─────
     ├───┤┌───┐
q_3: ┤ X ├┤ X ├
     └───┘└───┘

I think I know how to implement permute for PauliOp, but do you have any suggestion how to do it for MatrixOp?

[molar-volume]

Permute implemented for PauliOp and MatrixOp as well.

import numpy as np
from qiskit.aqua.operators import MatrixOp, X, Y, Z, I, CircuitOp

pauli_op = (X ^ Z).permute([0, 2]) @ (X ^ Y ^ Z ^ X)

matrix_op = (X ^ Z).to_matrix_op().permute([0, 2]) @ (X ^ Y ^ Z ^ X).to_matrix_op()

circuit_op = (X ^ Z).to_circuit_op().permute([0, 2]) @ (X ^ Y ^ Z ^ X).to_circuit_op()
print(circuit_op.to_circuit())

circ_eq_pauli = np.allclose(circuit_op.to_matrix(), pauli_op.to_matrix())  # circ_eq_pauli: True
circ_eq_matrix = np.allclose(circuit_op.to_matrix(), matrix_op.to_matrix())  # circ_eq_matrix: True

     ┌───┐┌───┐
q_0: ┤ X ├┤ I ├
     ├───┤├───┤
q_1: ┤ Z ├┤ Z ├
     ├───┤└───┘
q_2: ┤ Y ├─────
     ├───┤┌───┐
q_3: ┤ X ├┤ X ├
     └───┘└───┘

Do you want to provide compose with indices (or permute) for ListOps as well @Cryoris ?
The implementation of permute for ListOps might go in line with MatrixOp.permute

[woodsp-ibm]

Is there are good reason for the rename of identity to identity_operator? Since this is a breaking change on any existing code, if we needed to do this, then we ought to be keep identity still (which would call over to the new function) but issue a deprecation warning with an appropriate message to use the method by its new name.

[molar-volume]

@woodsp-ibm identity in OperatorBase was defined by me as part of this PR, I changed it to identity_operator, because @Cryoris pointed out that it is unclear what to expect from calling CircuitOp.identity

[molar-volume]

after reading #1086 I note that expand_to_dim will be overridden for SummedOp and ComposedOp, because e.g SummedOp.num_qubits is supposed to return the List[int]

[woodsp-ibm]

Ah ok, identity is a new method that was added by this PR. so you just renamed your new method. From my perspective I am not sure that the name change does so much for me since they are all operators according to OperatorBase. Maybe the method anyway should be some class or even module level function to return a specific operator of said type i.e. an identity in this case. I do not know why it would be an instance method since it has no real need of self (or am I mistaken again) - unless it were an identity spanning the same number of qubits as the operator in self in which case I could somehow see it.

[molar-volume]

@woodsp-ibm self.identity_operator is called in method expand_to_dim, which is defined in OperatorBase.
So based on the self the specific implementation of identity_operator is called, see how it is used in PrimitiveOp._check_zero_for_composition_and_expand

[woodsp-ibm]

Yes, I see some implementations of identity_operator do use self (when they need to copy across a property like is_measurement) while some do not. Do you see this as function that would be called in general by a user of the API - I guess you can if you want an instance that is an identity across that number of qubits - but of course with it being tied to self one needs to do that off an existing object which just seems weird. Otherwise its just a really some helper method to facilitate the one expand_to_dim implementation in the base class rather than have it implemented in each of the classes.

[Cryoris]

LGTM thanks for all the work!

[woodsp-ibm]

I should add thanks for all your hard work and persistence to see this through. It has been a while, with a lot of discussion, so thanks again for your patience and perseverance in making this happen!

[molar-volume]

My pleasure @woodsp-ibm ! I know the discussion was necessary and helped me a lot, so thanks for all your comments. Thanks @Cryoris !