_sanitize_monomials does not bring to canonical form

[apozas]

When inputting a monomial not in canonical form, the result of the function is not in canonical form either. For instance, inputting A_1_1_0*A_1_0_0 gives CompoundMonomial(A_1_1_0 A_1_0_0) instead of CompoundMonomial(A_1_0_0 A_1_1_0) (see the commit message). This happens both when the input is a CompoundMonomial, an array, or a string.

[eliewolfe]

Great call @apozas. Equality comparisons of compound monomial should work if all the atoms in MON1 are the conjugates of the atoms in MON2. I need to added this to the CompoundMonomial class.

• [x] Fix comparison between CompoundMonomial instances to account for conjugation
• [x] Adjust hash & signature & name of CompoundMonomial to minimize under conjugation.
• [x] Add interpret_compound_string and interpret_atomic_string functions to supplement _sanitize_monomial.

[eliewolfe]

So, I got most of it working.

Unfortunately, we are forced to take a major performance hit, as we can no longer use the orbits to assist in mapping atoms to their canonical form under inflation symmetries. This is because the orbits (and the whole initialization of the unsymmetrized momentmatrix) is manually set to enforce conjugation symmetry. I see two options for restoring the performance benefit:

1. Check to see if all the operators within the monomial commute with each other. (Easily done via the new commuting_operator_sequence_test test introduced in https://github.com/ecboghiu/inflation/commit/64e4aaaeb481866897c8a8aca2d463c63acf273e). If yes, then we can use the orbits as before to find the representative. This works for many monomials, but a diminishing fraction as the NPA level increases.
2. Overhaul both calculate_momentmatrix and apply_inflation_symmetries to NOT impose conjugation symmetry. This will allow us to construct an all-encompassing cache of where things go under inflation symmetry but not with conjugation symmetry. We would then have to call a new function (apply_conjugation_symmetry) which would then yield the final symmetrized-and-symmetric moment matrix. More work, but maybe worth it?

[eliewolfe]

I just implemented the first (easy) option per my previous comment.

[apozas]

I risk looking like a fool because I am not completely aware of what the current processing is, but couldn't it be a matter of adding the hash of mon_v2 to canonical_mon_to_idx_dict to look it up wherever appropriate?

[eliewolfe]

@apozas Let me explain why that won't work with an explicit example. Suppose we have the following four noncommuting operator pairings:

1. A^11_1 A^12_0 [by which I mean ]
2. A^12_0 A^11_1 [Note that #2 is the conjugate of #1]
3. A^11_0 A^12_1
4. A^12_1 A^11_0 [Note that #4 is the conjugate of #3] Note that #1 is TRULY equivalent to #4 under source relabeling (i.e., inflation symmetry), as #2 is truly equivalent to #3 under inflation symmetry.

During calculate_momentmatrix from fast_npa.py we currently create a dictionary (called unsymidx_to_unsym_monarray_dict) that says (for example) {...., index 25: op_seq #1, index 57: op_seq #3, ....} where the hashes of #2 and #4 simply do not appear as keys, because they are "less minimal" than #1 and #3 respectively. The inversion of this dictionary is called _unsymidx_from_hash_dict and it looks like: {...., hash(op_seq #1): index 25, hash(op_seq #3): index 57, ....}

Now, after applying inflation symmetry, we would recognize that index #25 and index #57 are in the same orbit.

We would then create a dictionary (canonsym_ndarray_from_hash_cache) which says that any operator sequence hash formerly associated with EITHER index #25 OR index #57 should be associated to the operators sequence of what was FORMERLY index #25. This dictionary was utilized as follows:

• If we encountered a numpy array representation that matched operator sequence #1 above, instead of calling to_representative from general_tools.py we would "cheat" by recognizing that hash(op_seq #1) is a key in the (inversion of) the dictionary constructed after calculate_momentmatrix.
• We then use canonsym_ndarray_from_hash_cache to return its representation under source relabelling.

Can you see the problem? This mean that by trying to cleverly bypass to_representative we would actually be returning the conjugate representative of op_seq #1, unintentionally.

[eliewolfe]

All tests now pass, issue can be closed, or we can talk about refining my implementation.

[apozas]

But isn't it the point that the expectation values of the four monomials that you described are, precisely, the same (under the assumption that we can take our moment matrix to be symmetric), and thus should all be associated to the same variable in the moment matrix? A different story is, for instance, monomials such as A^11_1 A^11_0 C^11_0 C^11_1, which should be identified with A^11_0 A^11_1 C^11_1 C^11_0 but neither with A^11_0 A^11_1 C^11_0 C^11_1 nor with A^11_1 A^11_0 C^11_1 C^11_0.

I wouldn't feel safe closing the issue, at least, until we see that we recover 3.085 for Mermin's inequality. Tests were passing before and yet the code had a bug. Which indicates that we need at least one more test! ;)

[eliewolfe]

Actually... I'm a tiny bit concerned that even my second try to sometimes bypass calling to_representative by using the orbits (namely, the implementation of https://github.com/ecboghiu/inflation/commit/073ebf090daae2f95fe61b9aae8233e9cc55137a, which limits the use of orbits to bypass to_representative only when the monomial's canonical representation must be invariant under conjugation because all the operators in the monomial commute) might still be slightly problematic. When choosing a representative via orbits, we select the monomial to represent the orbit based on it having the smallest index in the unsymmetrized momentmatrix among all monomials in the orbit. By contrast, when choosing a representative via to_representative we select based on the smallest lexicographic ordering. The inconsistency is the definition of "minimal" could potentially lead to the following problem:

Example: Consider the following physical three operator sequences (everything commutes)

• mon_1 = A^11 B^11
• mon_2 = A^12 B^11
• mon_3 = A^13 B^11

Suppose mon_1 and mon_2 are inside the moment matrix. It is plausible (given a manually specification of the generating monomials) that under apply_inflation_symmetries we find that mon_2 is selected to represent the orbit of {mon_1, mon_2}. However, if we call to_representative on mon_3, it will yield mon_1, as mon_1 is the true minimum under lexocgraphic ordering.

This would lead to AtomicMonomials being initialized for both mon_1 and mon_2 which would have different hashes (and would notably also be unrelated by conjugation). So, bottom line our SDP code would fail to treat these two things as the same variable.

As bugs go, this is low priority. 1) This seems like an edge case, not easy to come up with an example 2) It doesn't affect the validity of the SDP, it just means we are working with an outer approximation of the quantum set that could be made tighter.

Now, I can quickly write a patch that will take care of the aforementioned example. Namely, after to_representative does its thing, we can first check to see if NOW the ndarray can be processed by the dictionary of hashes. If so, we apply the dictionary instead of just returning it. (I'll write this patch now!) However, it still leaves open the extreme edge case where mon_2 above appears in the moment matrix before symmetrization while mon_1 does not. orbits would lead mon_2 to itself, but mon_3 to mon_1 .

Now, one way we could fix this is by changing apply_inflation_symmetries to use lexicographic minimality instead of index minimality. But I think that is a terrible idea, and we should not do it!! The problem is that apply_inflation_symmetries is already a computational bottleneck; it is slow on large problems. Making it compute lexocographic minimality would be a major performance hit. Totally not worth it, IMHO. To be discussed. Opinion @ecboghiu and @apozas ?

[eliewolfe]

But isn't it the point that the expectation values of the four monomials that you described are, precisely, the same (under the assumption that we can take our moment matrix to be symmetric), and thus should all be associated to the same variable in the moment matrix? A different story is, for instance, monomials such as A^11_1 A^11_0 C^11_0 C^11_1, which should be identified with A^11_0 A^11_1 C^11_1 C^11_0 but neither with A^11_0 A^11_1 C^11_0 C^11_1 nor with A^11_1 A^11_0 C^11_1 C^11_0.

I wouldn't feel safe closing the issue, at least, until we see that we recover 3.085 for Mermin's inequality. Tests were passing before and yet the code had a bug. Which indicates that we need at least one more test! ;)

That is NOT the point. The point is, as you say, that A^11_0 A^11_1 C^11_1 C^11_0 should not be identified with A^11_0 A^11_1 C^11_0 C^11_1. See, to_representative is exclusively used (since the introduction of the CompoundMonomial class) to obtain the representative of atomic factors. So, my goal is create a memoized bypass of to_represenative which does not consider conjugation symmetry.

[eliewolfe]

Update: https://github.com/ecboghiu/inflation/commit/f7faa8f304f7e74afd6628e52a33476aaf7ac6e6 provides protection against the first contrived example in https://github.com/ecboghiu/inflation/issues/75#issuecomment-1264002622.

[apozas]

So, my goal is create a memoized bypass of to_represenative which does not consider conjugation symmetry.

Wait, one step back. How is all this helping in solving the issue described in the first post?

In fact, I believe that the issue has a much simpler solution: check if the input (after transformation to whichever suitable format) is in the list of indices, and if not, compute the conjugate and check if this one is! This is particularly easy with the array formulation: you just reverse the list of components of the monomial and you sort by parties. So why do we have to "take a major performance hit" or "overhaul calculate_momentmatrix"?

[eliewolfe]

Because the "list of indices" you talk about (which maps unsymmetrized monomials to their orbit representative) already performs conjugation (sometimes) unintentionally! And we are blind as to whether or not any given key is mapped to its representative under JUST inflation symmetries or if its mapped to its representative under inflation symmetries AND conjugation symmetry. Re-read my explicit example in https://github.com/ecboghiu/inflation/issues/75#issuecomment-1263953725.

[eliewolfe]

Notes to self: the function to_representative_pairs in general_tools.py has a lot of room for serious performance improvements.

• [x] Employ caching: record all the permutation variants of a monomial while it is being walked to its lexicographic minimum under source relabelling. Save hashes of these in a dictionary with the final result, so that any time a variant is encountered we don't need to revisit it.
• [x] Don't use the inflation levels of the scenario, use only the inflation levels within that monomial (in to_repr_swap_plus_commutation)
• [x] Pass to_canonical_memoized into the function, so it doesn't have to rely on the slower to_canonical from fast_npa.py. (This is probably best accompished by converting it to a method of InflationSDP, instead of a function in general_tools.py.)
• [x] Improve to_canonical from fast_npa.py so that if all the operators in the monomial are mutually noncommuting (as is not super uncommon) we return the monomial immediately without further processing. (This also means no need to call to_repr_swap_plus_commutation after to_repr_lower_copy_indices_with_swaps.)

If we implement these tweaks, I think we can achieve high performance without trying to use the orbits to bypass to_representative, which we are seeing is problematic.

[apozas]

(Off-topic): How do you expect to find these notes to self once we close the issue? Please, open a new issue, a project, or anything else. I will proceed to do so on Tuesday morning in Europe if it has not been done by then. In any case, I'd say that right now the 100% of our effort must be put in making a ready-to-release code.

[apozas]

Because the "list of indices" you talk about (which maps unsymmetrized monomials to their orbit representative) already performs conjugation (sometimes) unintentionally!

This is not a problem, and in fact it is desired, if we are talking about complete monomials instead of atomic factors. The elements of the moment matrix are expectation values, and when using real moment matrices the expectation value of a cell and its conjugate coincide. Indeed this is not the case when we talk about atomic monomials. When we speak of the moment matrix, ALL the four monomials in https://github.com/ecboghiu/inflation/issues/75#issuecomment-1263953725 MUST be mapped to the same variable.

This was easily done in pre-Monomial code, namely:

1. Generation of the unsymmetrized moment matrix: Generate the monomial in a cell, check if it is already in the moment matrix and put the corresponding index if we have. Otherwise generate the conjugate of the cell, and check if this one appears in the moment matrix already. If yes, add a new entry mapping the non-conjugated form to the index that is already in the moment matrix (which corresponds to the conjugated monomial). Otherwise, add both the original monomial and its conjugate to the dictionary of indices, with the same index as target.
2. Apply inflation symmetries: In the unsymmetrized moment matrix, apply the symmetries by comparing $\Gamma{i,j}$ with $\Gamma{\pi(i),\pi(j)}$, and store which index went to which in the orbits dictionary.
3. Write a dictionary that maps every unsymmetrized monomial (at the time, in string form, which is what we had) to its final index in the symmetrized moment matrix.

This procedure leads to a moment matrix where all four cases in your example are mapped to the same index (as they should, because the four expectation values are the same in real moment matrices), but where $A^{11}_0 A^{11}_1 C^{11}_1 C^{11}_0$ and $A^{11}_0 A^{11}_1 C^{11}_0 C^{11}_1$ are NOT mapped to the same index. Note that, in this whole process, to_representative is ONLY used to find a representative of the monomials AFTER the whole process above has been completed, and it is applied AT THE LEVEL OF WHOLE MONOMIALS and not at the level of atoms. This has two consequences, namely that using to_representative IS SUPERFLUOUS and it is just for ease-of-use/sanitization reasons (but without it the code would work equally well), and that monomials corresponding to different indices in the symmetrized moment matrix are NEVER mapped to the same representative, $A^{11}_0 A^{11}_1 C^{11}_1 C^{11}_0$ and $A^{11}_0 A^{11}_1 C^{11}_0 C^{11}_1$.

I really think that replicating this behavior does not need so many new lines of code that have been added referencing this issue, nor change core functions (calculate_momentmatrix and apply_inflation_symmetries) that are working the way as intended. It should "just" be a matter of thinking everything in terms of expectation values instead of operators.

[apozas]

EDIT: I have been able to go through the code in a more detailed way, and I think that the problem is a design choice. With the assumption that our monomials represent expectation values, a CompoundMonomial is not just a product of AtomicMonomials, because we need extra information about the global ordering of the operators. I thus think that the core object should be a Monomial (which, can and must have properties such as .factors). Keeping AtomicMonomial as the basic building block, and patching the code to make CompoundMonomial work as intended will make our code much harder to maintain and extend. I am working on the code refactoring, which I will put in a different branch.

[eliewolfe]

@apozas you are absolutely right that this "problem" is stemming from a design choice, and it is a design choice I will defend. Although our monomials represent expectation values, it is extremely valuable to first and foremost consider them to JUST be products of AtomicMonomials. This is useful for many reasons. 1) set_values works based on factors. This is a BIG DEAL. A basic task that comes up over and over again is trying to recognize when a subset of factors inside a CompoundMonomial is equivalent (after inflation symmetry relabelling) to some other CompoundMonomial. This is super straightforward because of this design choice. 2) there is no actual problem. Now that each AtomicMonomials carries with it the information about its conjugate (at essentially zero added overhead) everything works just as we want it to. 3) I worry that I may have given the wrong impression when commenting earlier in this thread about a performance hit which occurs when we do not aggressively make use of the orbits to identify monomials. The first use of orbits to improve the performance of to_representative only showed up in our codebase (added by your truly) only about a month ago. Clearly the code was working just fine without this tweak. It was merely a performance tweak I thought was particularly clever. Even if we had NO benefits from this tweak (because of design choice) it would probably not be noticeable to you. AND: 3a) we DO still get benefits from this tweak, as per https://github.com/ecboghiu/inflation/commit/073ebf090daae2f95fe61b9aae8233e9cc55137a, and additionally we now also have very high performance WITHOUT the tweak per https://github.com/ecboghiu/inflation/commit/e9e8f3e94124a41990747bf704d099a0e7471e15.

[apozas]

I am not saying that we have to refuse to use something like AtomicMonomial because it is indeed very useful. To me, the problem is that CompoundMonomials are built from AtomicMonomials. Look at how much code had to be added/fixed in order to build information about the conjugate. And moreover, since we use real moment matrices, having conjugates this way is only confusing.

I have been thinking through this problem today, and I think I have found a solution that should address all of our problems. Essentially, we would be building objects close to CompoundMonomials directly from the symmetrized moment matrix and the list of numpy arrays inside it, and each of them will have a number of Factors (similar to AtomicMonomials, but among others clearly not self-adjoint). I will explain tomorrow.

[eliewolfe]

Well, for now I'm still on a mission to convince you @apozas otherwise. I have tried to significantly further simplify the codebase for easier maintenance in https://github.com/ecboghiu/inflation/commit/f2478281ff469119195d9e36d89f710c14761588. Looking forward to talking about it over Zoom.

[apozas]

I think that after our discussion today, this issue can be finally closed. Thank you for the work on this.