Create tools to build operators

[a-corni]

I am taking build_operator out of Hamiltonian and QutipEmulator, such that it can be used to generate effective noise operators.
I associate the operators to the eigenstates of a SequenceSamples (introduced in #705). I provide:

• default_operators to generate the default operators associated to a SequenceSamples.
• build_projector to build a specific 1qudit projector.
• build_1qubit_operator to make linear combinations of 1qudit projectors (or other 1 qudit matrices provided in an optional dictionary).
• build_operator to make tensor product of 1 qudit operator (function that was already in place).
  What I can provide:
• a tutorial on how to use these operators
• a function that decomposes any operator into a tensor product of linear combinations of projectors.
  Other question:
• At the moment, doing sim.build_operator(("I", "global")) returns III + III + III if sim has 3 atoms in its register, is this really the wanted behaviour ? (sim.build_operator has always behaved this way)

[a-corni]

@a-quelle, it will be easier to discuss this PR :)

[HGSilveri]

• At the moment, doing sim.build_operator(("I", "global")) returns III + III + III if sim has 3 atoms in its register, is this really the wanted behaviour ? (sim.build_operator has always behaved this way)

I agree this looks weird (especially with the example you provided above) and we can probably change it, I doubt anyone is using it.

But just for context, it does do what we want when building the Hamiltonian with a global drive. Take the detuning term, for example:

-\delta(t) (\sum_i Z_i) = -delta(t) (ZII + IZI + IIZ)

This is why the "global" is there.

[a-corni]

• At the moment, doing sim.build_operator(("I", "global")) returns III + III + III if sim has 3 atoms in its register, is this really the wanted behaviour ? (sim.build_operator has always behaved this way)

I agree this looks weird (especially with the example you provided above) and we can probably change it, I doubt anyone is using it.

But just for context, it does do what we want when building the Hamiltonian with a global drive. Take the detuning term, for example:

-\delta(t) (\sum_i Z_i) = -delta(t) (ZII + IZI + IIZ)

This is why the "global" is there.

Sure, I have never doubted about the usefulness of the "global", I was just wondering about its behaviour with the identity - to me, it's natural that it behaves this way but user might forget about this.

[HGSilveri]

Sure, I have never doubted about the usefulness of the "global", I was just wondering about its behaviour with the identity - to me, it's natural that it behaves this way but user might forget about this.

Indeed, I don't think it's something we want to expose. We can just change the one place in the code base where it is being used and get rid of it.

[a-quelle]

could you change this example into, e.g.
[(sigma_gg, ["q0"]), (sigma_rr, ["q2"])] applies sigma_ggIsigma_rrI
to make explicit that the tensor product, rather than the sum is applied. I guess this could be inferred from

Returns the operator given by the tensor
product of {operator_i applied on qubits_i} and Id on the rest.

But that phrasing is just vague enough that I assumed it would sum over list elements, rather than do the product.

[a-quelle]

You feed the sequence samples to get the used basis, for which I see you've added a method to SequenceSamples. What is the reason you've added it here, and not to Sequence. Now you force the user to discretize the sequence whenever they want to define an operator, while in principle, that should be an operation the backend should do as part of its internal workings, no?

[a-quelle]

Alternatively, if all you care about from the SequenceSamples is to obtain a tuple like ('r','g') from it, why not ask the user to feed in the basis like that directly. In that case, the user can get it from either Sequence or SequenceSamples since I see methods in both of them for getting the basis. It has the added benefit of making it more explicit what information is actually used by these functions, and also that you can use it on an equal footing with Backends which take the Sequence and Simulators which take the SequenceSamples, since this API then no longer requires either of those.

[a-quelle]

The ArrayLike here makes this interface dependent on numpy, even though the abstract action of creating an operator should not need to know about numpy (specifically, I want to implement this function in my MPS backend that only uses torch). Is this type hint really necessary? The values of the dict are currently not used for example, do you envision a need for them in the future?

[a-corni]

I agree I should just ask for the qubit_ids, there is no need to have the qubit locations for this.

[a-corni]

[Draft]: I am not sure this is exactly what you had in mind @a-quelle, and I am not sure that's exactly what I want from this PR, but I have started introducing classes to ease the definition of 1-qudit and n-qudit operators. Do you think it helps the user ?

[a-quelle]

[Draft]: I am not sure this is exactly what you had in mind @a-quelle, and I am not sure that's exactly what I want from this PR, but I have started introducing classes to ease the definition of 1-qudit and n-qudit operators. Do you think it helps the user ?

I do think this workflow is nicer for the user, however for the implementation details, I think it might be a good idea if we sit down, so I can show you what I'm trying to do on my end, and why.

[a-corni]

I could use some guidance here.
I have just figured that with the new operators, we can build the Hamiltonian without Qutip. We could then take out the Hamiltonian class out of pulser-simulation such that it can be used by the other frameworks ?
In pulser-simulation, I would put a "QutipHamiltonian" that would convert the Hamiltonian into a QobjEvo.
However, this might go a bit beyond the scope of this PR and I think it should be left as an issue after this PR is merged
What do you think @HGSilveri @a-quelle ?
Could also mention @lvignoli whom I discussed this possibility with.

[a-quelle]

What is the need for this? I'm not familiar with pulser simulation, but I suspect it's a bit specific to there. I'm currently implementing something very similar to the Backend ABC, and I have a run method like

        coeff = 0.001  # Omega and delta are given in rad/ms, dt in ns
        dt = mps_config.dt
        omega_delta = _extract_omega_delta(sequence, dt, mps_config.with_modulation)
        emuct_register = get_qubit_positions(sequence.register)

        evolve_state = MPS(
            len(sequence.register.qubits),
            truncate=False,
            precision=mps_config.precision,
            max_bond_dim=mps_config.max_bond_dim,
            num_devices_to_use=mps_config.num_devices_to_use,
        ) 

        for step in range(omega_delta.shape[1]):
            t = step * dt
            mpo_t = make_H(
                emuct_register,
                omega_delta[0, step, :],
                omega_delta[1, step, :],
                sequence.device.interaction_coeff,
            )
            evolve_tdvp(-coeff * dt * 1j, evolve_state, mpo_t, mps_config.max_krylov_dim)

        return result

so there is no real reason to explicitly store a time series anywhere. If the user wants to measure an operator, the can create it from the Operator ABC that you created below.

[a-corni]

The idea with this Hamiltonian(TimeOperator) is that it would avoid you making _extract_omega_delta, get_qubit_positions, and make_H, you could just do Hamiltonian.from_sequence(sequence, sampling_times), that would make a list of mpo_t (as Operator) and time steps, and then you could convert these Operators into mpo_t (as matrix product states, with the converter you have defined).
That could enable us to share the same definition of the Pulser Hamiltonian ? I think this could be useful for other emulators as well.

[a-quelle]

For code efficiency reasons, I would avoid that code path even if it were available. Having an abstract representation for operators and states is unavoidable if you want a backend independent interface, but within the backend, everything goes. Constructing the Hamiltonians by going through the abstract format that we're now defining would be at least an order of magnitude slower than the bit fiddling I'm doing in make_H, as well as creating problems with differentiability due to svd operations in the general MPO addition algorithm that I took some pains to avoid in make_H.

[a-quelle]

I appreciate that you're trying to be helpful, but that is not your job. Rather it's my job to be helpful by removing the burden what logic is needed to actually efficiently simulate a pulse in a given backend from the pulser team.
For historical reasons, it doesn't quite work that way yet, but there is no time like the present to try and improve things.

[HGSilveri]

I also think this would severely complicate things for marginal benefit. If the Hamiltonian was a single operator I could see it, but the time-dependence seems like a deal breaker to me

[a-corni]

Thanks for your feedback @a-quelle and @HGSilveri, let's forget about this feature then.

[a-corni]

I have investigated the implementation of multiplication by an operator or a qudit-state, and could generalize this to multi-qudit states, but I would like to know if this is necessary before moving forward @a-quelle @lvignoli

[a-quelle]

It's pretty cool that you're implementing an algebra structure on the abstract representation, but I do wonder if you're not making your life needlessly difficult. I have no immediate use for it, given that I can so far create the state/operator in the backend and add, multiply, apply afterwards. In addition I see the following dificulties:

I imagine that you want the contents of operators: list[str] to be essentially arbitrary, since then you can make use of the dict in build_operator below, but that means that you're going to get stuck with complicated things in your abstract representation for which implementing the algebra structure is hard, and also very slow compared to creating the state in the backend, and doing the operations when things are numbers, and not strings.

[a-corni]

So I have restricted the operators that can be defined in this QuditOperatorString to projectors, namely {sigma_ab with (a, b) in [u, d, r, g, h]²}, such that I can make this algebra easily. However, I do agree it might be quite slow, plus we don't see yet a clear use of this and it needs quite some work to develop, so I guess it's best to store this as a possible development (I would only need the addition and dagger operations if I were to develop the Hamiltonian(TimeOperator))

[a-corni]

This is surely not the most elegant way to store it, I can try to avoid storing one operator per qubit and try to store one (operator, [all the qubits it addresses]) I would just like to agree on the global structure before moving forward.

[a-quelle]

I don't know qutip quite well enough to understand the difference between a Qutip.QObj and a QutipOperator. Given that this part of the code is unrelated the the interface we're trying to define, I don't really feel qualified to comment.

[a-quelle]

Should this class not have a from_string abstract method or something, that forces backends to specify how to use the operator_string to actually create an instance? That way you can do all required operator releated things through this class.

[a-corni]

So in this implementation I have investigated a class that could be common to all backends, that all backends could inherit from and just implement a to_myobject() and a from_myobject() method - see QutipOperator in pulser-simulation/operators.
The advantage of this approach seemed at first that some checks could be shared with all the backends, but the issue is that to perform any operation you should:

• either implement your own algebra on Operator (solution I favored)
• convert the Operator into your objects [to_myobject()], make the operation with these objects, convert the resulting object into an Operator [from_myobject()]. This seems to be very costly, see my comment in pulser_simulation/operators.

[a-corni]

And so in the current implementation, the from_string is the init of the Operator class. I thought of this implementation when writing QutipOperator. I realized then that converting an object into a string is complex and takes time, whereas if you store the string you have it immediately when you want to send your operator to another backend. That's why I favoured this implementation rather than storing the object and converting the object into an OperatorString when needed (some operations might take more time (I don't think that addition is so long compared to addition in qutip) but the main operation being the conversion from object to OperatorString, this will take less time)

[a-corni]

If you still favor storing the object and having an abstract class method from_operator_string (and an abstract method to_operator_string), I am still open to this, just let me know if my arguments made sense :)

[a-quelle]

Yesterday I was very tired and in the train, and I did not look and think as closely as I should have. I think I understand what you're trying to do, and I don't think it's compatible with the intended use of an ABC.

ABCs are intended to enforce a software contract. In this case, we want all backends to enable the same workflow, and at this step we use an ABC to enforce that however an Operator works under the hood, the user can always use it the same way, with a few given methods. This is what @abstractmethod does. The moment an ABC has an abstractmethod you can no longer construct it, and any subclass must implement that method with identical signature.

This means you cannot make an Operator directly, and since __kron__ etc. are abstract, I must implement them in a subclass, and they will never be called unless I choose to do so via super(). This effectively makes the code in your implementations unreachable. Conversely, __add__ and __rmul__ are not abstract, so I'm not forced to implement them. Furthermore, the default implementations of __add__ and __rmul__ depend on operator_string etc. but in my own subclasses to Operator I have no incentive to actually populate that field.

It looks to me like you're trying to be helpful by already giving some default implementations for things. But since these implementations have to be general, they are unlikely to be appropriate to any specific backend, and you also have to sacrifice the ability of the ABC to enforce a software contract to do it.

[a-quelle]

For this reason, i would recommend a minimal implementation like so:

class Operator(ABC):
    @abstractmethod
    def __add__(self, operator: Operator) -> Operator:
        pass
     @abstractmethod
    def __rmul__(self, scalar: Number) -> Operator:
        pass
    @abstractmethod
    def __mul__(self, state: State) -> State:
        pass
    @abstractmethod
    def __matmul__(self, operator: Operator) -> Operator:
        pass
    @abstractmethod
    def kron(self, operator: Operator) -> Operator:
        pass
    @abstractmethod
    def from_string(self, string:OperatorString) -> Operator:
        pass
    @abstractmethod
    def to_string(self) -> OperatorString:
        pass

this implies we need some State(ABC) similar to how we have the OperatorString and the Operator. If we want to implement the algebra on OperatorString we could have it inherit from Operator.

[a-quelle]

I guess we could just have OperatorString inherit from Operator anyway, since it implements some of the methods, and for the others we could throw an exception

[a-corni]

Thanks for your feedback @a-quelle. I realized yesterday that I was indeed not using ABC correctly, and I now fully understand the way to go.

[a-quelle]

tuple(zip(...)) or change the type hint since zip returns an iterator

[a-quelle]

This implementation suffers from the same issue we discussed with Operator, namely, any state can be written as a sum of kron, where the kron is over basis elements, but not every state can be written as a kron of sum where the sum is over basis elements (these last are the product states, which have no entanglement entropy for any bipartition of the qubits). For this reason, I would do something like

QuditString = str

class MultiQuditString:
    """weighted sum of QuditString."""
    coefficients: list[Number]
    states: list[QuditString]
    def as_tuple(self) -> tuple[tuple[Number, QuditString]]:
        pass #too lazy to write the implementation

I'm not entirely sure how STATES_RANK works, but I assume you can properly check for consistency even in this situation? Maybe QuditString should not simply be a type alias but an actual class with a post_init?

[a-corni]

Are QuditString and MultiQuditString really necessary ? I implemented this to implement the mul operation in the current Operator, but in the new implementation I am not sure it is...

Otherwise I completely agree with your comment.

[HGSilveri]

So, if I understand your motivation correctly, you were defining these so that the user could have a way to define operators without a backend specific support, just from strings.
I understand the need for not being dependent on qutip to define your operators (let's say you only want to interact with Emu-T, it doesn't make sense to have a qutip dependency), but perhaps it would be simpler to define operators with a numpy support instead for that case, rather than working out a string based algebra.

[a-corni]

If I were to implement states - I want to be 100% sure this is needed beforehand. I could do (to keep the spirit of Operator): QuditString (1-qudit state)- TargetedQuditString (multi-qudit state as kron of qudit states) - MultiQuditString (sum of kron) - MultiQudit(ABC) (abstract class like Operator) - NumpyMultiQudit(MultiQudit) (stores a numpy array, makes operation with this array, convert this array into string and from string into this array)

[HGSilveri]

I feel like this is getting out a hand, I would suggest we meet and come up with a plan that keeps things as simple as possible while delivering the functionality we are looking for.

[a-corni]

I think QuditString is too far fetched and should not be implemented, nor the operations on OperatorString. I definitely think we should go with something minimal. I will resume my work on this PR tomorrow.

[a-quelle]

I would like to have an abstract representation for strings, so that users can measure coefficients of the wavefunction via the inner product on states. I don't care about the vector space structure and interaction with operators very much, if you can create a state from a string like "rgrgrgrgrgr" then it's already good enough for me, since I could do something like MPS.from_string("rgrg) and then use the vector space structure implemented for MPS to make any state I'd need. So at the very least, this means having a State ABC similar to Operator with a from_string abstractmethod, but in the interest of a minimal implementation I don't really think we'd need the QuditString, and as far as I care, we don't need to implement State in this PR either.

What I would like to know already though, is if a from_string method to make a state would require qubit_ids and eigenbasis just like in the Operator case.