Implement flexible operator substitutions

[apozas]

If I am not mistaken, we still depend on ncpol2sdpa to generate problems with arbitrary substitution rules. I guess that the final goal will be to remove this dependence completely, mainly for speed reasons. We already have implementations of the generation of the moment matrix for non-commuting parties and for all commuting operators using numba, and the aim would be to substitute it with a framework as versatile as ncpol2sdpa, since we want to leave flexibility to the users in choosing their substitution rules. What is what is preventing allowing arbitrary substitutions with numba?

[ecboghiu]

Some quick thoughts. I think this depends on what kind of substitutions.

If it is simply arbitrary commutation rules, then such commutation rules could be encoded as a network whose (i,j) entry says if operator i commutes with operator j. This can be used to determine if two arbitrary operators commute, and it could replace functions such as this one. The user can input a dictionary, and we can translate it into this type of matrix. This is one idea. Since everything is just looking up elements of a matrix, it should be very fast.

If it is more than commutations, then this can still be done if they are simple enough, for example, AB<=C would simply consist on looking for the product AB in a monomial and substituting it with C, and then simplifying/canonicalizing again.

If it is polynomial substitutions, if they are linear, it is still doable. A-> B+C is doable, but we need to distinguish: 1) doing this substitution at the level of operators or 2) doing this at the level of moments.

First, for 2), this can be done at the last step. If we substitute $\langle A \rangle$ with $\langle B \rangle$ + $\langle C \rangle$ then this simply requires updating the sparse matrices $F_0, \ldots,F_i$ that define the moment matrix through

$F_0 + \sum_i F_i x_i$

For a), we would need to encode the moment matrix FROM THE BEGINNING in the notation $F_0 + \sum_i F_i x_i$. If we have this description, then we can write entries (i,j) of the moment as linear combinations of $x_i$ easily. Do you agree?

So as a summary, simple substitutions are easy, polynomials substitutions require a bit of restructuring.

[apozas]

Polynomial substitutions are going to be quite hard. I was thinking of those that you define as "more than commutations, but still simple enough". In the end this is what we will use for using sources of different type or separable measurements. The procedure that you describe is precisely the one that ncpol2sdpa uses, and thus I was wondering what was preventing us of doing the same.

[ecboghiu]

ncpol2sdpa was removed in atomic-complete branch, the lastest version https://github.com/ecboghiu/inflation/commit/8c3f814efee735081ca54280096f19dfa99a8d06. Given that this issue is more general than that, I'm updating the name.

[ecboghiu]

One important step next is to implement flexible substitution rules. Currently we hardcode the following two types of substitutions:

• Projective rules, i.e., $A^2=A$. This is done with remove_projector_squares.
• Orthogonality rules, i.e., $A{xa} A{xa} = \delta_{aa'}$. This is done with mon_is_zero. After lexicographic ordering, it checks whether there is such a product present. Note: it is crucial to have the outputs last in our ordering convention, else when bringing to canonical form, products of orthogonal projectors might not end up together.
• Particular commutation rules. These are done through the function is_commuting.

What I propose for the generalisation of the 3rd rule above is to process commutation-type-rule in the substitutions dictionary into a matrix of commutations. This would change is_commuting in that instead of computing whether to operators commute, it would hash these operators and look up in the commutations matrix whether they commute or not:

def to_canonical(op1, op2, commutation_graph):
    i, j = to_mat_index(op1), to_mat_index(op2)
    return commutation_graph[i, j]

One can also think of other ways to implement this. Whatever way we choose it must be compilable with numba. This is because this is at the core of to_canonical, and this function will be called an extremely high number of times. This is why I suggest storing the boolean value of whether two operators commute in a matrix, and is_commuting to just access the value in this matrix.

One idea for to_idx_mat(o) would be to do a linear lookup, but this is O(n), simply have a matrix

np.array([[0, o1], [1, o2], [2, o3], [4, 04]...])

and loop through until o==o1 and return the linear index. Another idea would be to convert the tuple of numbers to a single numeric digit:

@njit()
def nb_ravel_multi_index(op, dims):
    summ = op[0]
    prod = 1
    for i in range(1, op.shape[0]):
        prod *= dims[i-1]
        summ += op[i]*prod
    return summ

This function is very fast, around 300ns on my pc. Once we have this hash, we could use a linear search, but now instead of comparing arrays of numbers, we just compare integers, which probably leads to a speedup.

In any case the tasks I propose would be:

• [ ] Convert the commutation-type-substitution rules into a boolean matrix indexed by operators, and change is_commuting to simply look up values in this matrix in a Numba compatible way.
• [ ] Then implementing the commutation rules for hybrid networks and separable measurements amounts to simply writing the substiutions dictionary correctly.

[eliewolfe]

Can we close this now, @apozas ?

[apozas]

Well, I guess that we can understand that flexible substitutions are there, and implementing custom commutation rules for hybrid scenarios is something different, so sure.