Operator apply permutation

[alexanderivrii]

Summary

Adds a method apply_permutation to Operator class that modifies operator's data by composing it with a permutation. The permutation can be composed either before or after the operator This is achieved by directly swapping the rows or the columns of the array, and is more efficient than the more general compose.

Closes #9402.

Details and comments

This still misses a few tests and release notes. I am also not sure what is the most general version of the Operator, i.e. maybe this code should only work for unitary matrices.

[qiskit-bot]

Thank you for opening a new pull request.

Before your PR can be merged it will first need to pass continuous integration tests and be reviewed. Sometimes the review process can be slow, so please be patient.

While you're waiting, please feel free to review other open PRs. While only a subset of people are authorized to approve pull requests for merging, everyone is encouraged to review open pull requests. Doing reviews helps reduce the burden on the core team and helps make the project's code better for everyone.

One or more of the the following people are requested to review this:

• @Qiskit/terra-core
• @ikkoham

[alexanderivrii]

Here are some fun numbers on the efficiency of the approach. Here for each number of qubits nq we generate a single random unitary operator, take 100 random permutations, and provide cumulative times for applying permutation to an operator using one of 4 different methods.

def method0(op, perm_pattern):
    nq = len(perm_pattern)
    qc = QuantumCircuit(nq)
    qc.append(UnitaryGate(op), range(nq))
    qc.append(Permutation(len(perm_pattern), perm_pattern), range(nq))
    return Operator(qc)

def method1(op, perm_pattern):
    perm_op = Operator(Permutation(len(perm_pattern), perm_pattern))
    op &= perm_op
    return op

def method2(op, perm_pattern):
    perm_op = PermutationGate(perm_pattern).__array__()
    op &= perm_op
    return op

def method3(op, perm_pattern):
    op = op.apply_permutation(perm_pattern, front=False)
    return op

permutation-back
nq = 2, time0 =   0.11,   time1 =   0.05,   time2 =   0.00,   time3 =   0.00
nq = 3, time0 =   0.12,   time1 =   0.06,   time2 =   0.00,   time3 =   0.00
nq = 4, time0 =   0.16,   time1 =   0.08,   time2 =   0.01,   time3 =   0.00
nq = 5, time0 =   0.19,   time1 =   0.10,   time2 =   0.01,   time3 =   0.01
nq = 6, time0 =   0.31,   time1 =   0.16,   time2 =   0.04,   time3 =   0.01
nq = 7, time0 =   0.88,   time1 =   0.25,   time2 =   0.07,   time3 =   0.03
nq = 8, time0 =   4.75,   time1 =   0.77,   time2 =   0.18,   time3 =   0.08
nq = 9, time0 =  35.53,   time1 =   3.88,   time2 =   1.00,   time3 =   0.22
nq = 10, time0 = 294.70,   time1 =  19.63,   time2 =   5.68,   time3 =   0.78

permutation-front
nq = 2, time0 =   0.10,   time1 =   0.05,   time2 =   0.00,   time3 =   0.00
nq = 3, time0 =   0.12,   time1 =   0.06,   time2 =   0.01,   time3 =   0.00
nq = 4, time0 =   0.15,   time1 =   0.09,   time2 =   0.01,   time3 =   0.00
nq = 5, time0 =   0.18,   time1 =   0.10,   time2 =   0.01,   time3 =   0.01
nq = 6, time0 =   0.31,   time1 =   0.16,   time2 =   0.04,   time3 =   0.02
nq = 7, time0 =   0.81,   time1 =   0.24,   time2 =   0.05,   time3 =   0.04
nq = 8, time0 =   4.89,   time1 =   0.85,   time2 =   0.18,   time3 =   0.15
nq = 9, time0 =  35.54,   time1 =   5.04,   time2 =   1.09,   time3 =   0.53
nq = 10, time0 = 288.06,   time1 =  19.10,   time2 =   5.59,   time3 =   2.70

method1 (for permutation-back) is what we have right now, method3 is what we want to apply instead.
This is not an invitation to create method4 in Rust which will be 100x faster still :).

[coveralls]

Pull Request Test Coverage Report for Build 4747761829

• 29 of 29 (100.0%) changed or added relevant lines in 1 file are covered.
• 79 unchanged lines in 5 files lost coverage.
• Overall coverage decreased (-0.009%) to 85.923%

---

Files with Coverage Reduction	New Missed Lines	%
qiskit/qasm3/init.py	2	92.0%
crates/qasm2/src/lex.rs	5	90.63%
qiskit/transpiler/passes/routing/stochastic_swap.py	7	96.53%
crates/qasm2/src/parse.rs	12	97.11%
qiskit/dagcircuit/dagcircuit.py	53	89.55%

Totals
Change from base Build 4737986268:	-0.009%
Covered Lines:	71102
Relevant Lines:	82751

---

💛 - Coveralls

[alexanderivrii]

Here are a few questions. The function apply_permutation that I wrote accepts some permutation of qubits as input, lifts this permutation to a permutation of Operator's data, and then suitably permutes rows or columns of this data. However, this only makes sense when the operator is "composed of" qubits and not say of qubits and qutrits. Would it make more sense to have this function in a different file and not directly in operator.py? What is the best way to check that the Operator is only "composed of" qubits?

Update: here is a tutorial on Operators that explains about subsystems and gives an example of an Operator representing a system with one qubit and one qutrit: https://qiskit.org/documentation/tutorials/circuits_advanced/02_operators_overview.html

[ajavadia]

let's aim for consistency with this PR which implemented reverse_qargs for Operators but also works for States. And also works for qudits (because Qiskit operators don't assume qubits).
#5970

[ikkoham]

This PR proposes to calculate PO or OP where O is an operator and P is a permutation matrix. It is not a good idea to calculate permutations at the operator level, and in many cases they should be calculated on the statevector. But is there any specific use case?

I sometimes would like to calculate the conjugate action POP^\dagger or P^\dagger OP.
In this case, I think it could be naturally extended to qudit.

[jakelishman]

Ikko: why do you say it's not a good idea? Sure, it's more efficient to do on a statevector, but in the case of permutations it's not a terrible difference, since we can just apply a row/column-reordering rather than doing matrix-matrix multiplication.

The motivation here for allowing one-sided permutations on Operator is so we can more easily create operators from transpiled circuits, whose virtual-physical qubit mappings change between the start and the end of the circuit because of virtual swaps. We get several bugs from people complaining that they can't get a matrix representation of their transpiled circuit (or more precisely, that they say it's wrong), and in general to embed a transpiled circuit in a unitary operator, this is something we want too.

[alexanderivrii]

So, the current code only works for 2^n x 2^n operators, I don't think it immediately works for qutrits, etc.. I tried to think if it could be implemented along the lines of #5970 (cleverly using numpy.reshape, numpy.transpose and maybe numpy.moveaxis, etc.), and maybe this is possible, but I don't have enough experience reasoning about multi-dimensional numpy arrays.

The current code is more naive. For each row/column index from 0 to 2^n, we compute where it would move after applying the permutation, by first converting it to a bitstream: bit = bin(r)[2:].zfill(nq)[::-1]. However, if the operator's dimensions are (2x3) x (2x3), then we should work with row/column indices from 0 to 5 and convert these to mixed-based representation (*). Is there a simpler way than this? Is there a standard python code for converting integers to mixed base representations and back?

Note (*): this is a perfect opportunity to cite the paper by Georg Cantor Über einfache zahlensysteme. Zeitschrift für Mathematik und Physik, 14:121–128, 1869 (note the year!), and yes it's in German.

[ikkoham]

I see. I didn't know that use case. Then I think it is a good idea.
I tend to think of operators as observables or density matrices, so forgot about unitary operators. So I said conjugate action is better, but for the transpiler this proposal is nice.

Thank you. Almost LGTM! Could you add the release note?

[ikkoham]

How about inlace option? Or if this is destructive method, it would be nice to return nothing. If this is non-destructive, copy is needed in the first.

[jakelishman]

Yeah, I think this probably should act on a copy - that seems to be in line with what the rest of the Operator methods do, and they don't have inplace options. If we need inplace=True versions of things in the future, we might want to revisit on a larger scale.

[alexanderivrii]

Thank you! I have reimplemented the functionality based on the suggestions; there is no longer a need for the inplace option.

[jakelishman]

If the method I've sketched out does turn out to be the correct way to handle arbitrary-dimensioned tensors (I think it does work, at least up to the limit that Numpy only allows 31 dimensions in its arrays, so we can only have 15 spaces...), then we'd also need to add a couple of tests of that.

[jakelishman]

I think this can be done with numpy.transpose in order to support the more general n-dimensional case. It looks something vaguely like this:

import functools
import numpy as np

a, b, c = np.random.rand(4, 4), np.random.rand(3, 3), np.random.rand(2, 2)
_012 = functools.reduce(np.kron, (a, b, c))
_120 = functools.reduce(np.kron, (b, c, a))
_210 = functools.reduce(np.kron, (c, b, a))

def permute(matrix, shape_r, order):
    orig_shape = matrix.shape
    split_shape = tuple(shape_r) * 2
    new_order = tuple(order) + tuple(len(order) + x for x in order)
    return matrix.reshape(split_shape).transpose(*new_order).reshape(orig_shape)

[
    np.allclose(permute(_012, (4, 3, 2), (2, 1, 0)), _210),
    np.allclose(permute(_012, (4, 3, 2), (1, 2, 0)), _120),
]
# [True, True]

(the arrays aren't bit-for-bit identical because I think the internal transpose is done by abstract tensor-contraction, so involves a few FLOPs per matrix entry as opposed to doing the Kronecker product directly)

It'd want to test that for sure - I just did that super roughly. Mine is a double-sided permutation - you'd just tweak the left/right side of new_order to be range(len(shape_r))/range(len(shape_r), 2*len(shape_r)) to do a back/front-only permutation, I think, but I didn't totally test.

Mine currently assumes that the operator is square, but I think if you're only accepting one-sided permutations, you don't even strictly need that assumption. The shape variables I've used here can probably be read out of Operator._op_shape in some form - I don't remember the details off the top of my head.

[jakelishman]

Also fwiw I didn't bother to check the endianness conventions when I wrote that - it was just a proof of principle. I might have done it backwards.

[alexanderivrii]

Thanks @ikkoham and @jakelishman! I am trying to make @jakelishman's solution work (also many thanks to @gadial, who has a very similar conceptual solution), but all this business with endianness, interpreting permutation order, etc. is really confusing.

[alexanderivrii]

I have reimplemented the code based on @jakelishman's and @gadial's code snippets. Took me a while to get all the indexing right. The new function can handle general qudits and one-sided permutations.

[jakelishman]

In tests, we should be able to just import everything at the top of the file.

[alexanderivrii]

Done.

[jakelishman]

Also, in terms of efficiencies - if my method is right, it likely wouldn't benefit much (if at all) from being moved to Rust, since all the work would already be happening in compiled Numpy code, and their einsum tensor-contraction stuff is super optimised, I believe. Though I suppose in the case where you know the contraction you're doing is just row/col re-ordering, maybe you can do it faster with vectorised strided memcpy stuff.

[alexanderivrii]

Also, in terms of efficiencies - if my method is right, it likely wouldn't benefit much (if at all) from being moved to Rust, since all the work would already be happening in compiled Numpy code, and their einsum tensor-contraction stuff is super optimised, I believe. Though I suppose in the case where you know the contraction you're doing is just row/col re-ordering, maybe you can do it faster with vectorised strided memcpy stuff.

I have done a very quick experiment, and it seems that the implementation based on numpy's reshape and transpose is even faster than the previous version.

[jakelishman]

Thanks for sticking with this Sasha!

[jakelishman]

Pretty minor, but since this is an if/else block, it's marginally clearer to reason about if the condition is "positive" (i.e. if front / else). I get that not front is the default, though.

[alexanderivrii]

changed that, thanks

[jakelishman]

I've got to say that this feels wrong to me - we seem to be splitting the matrix into a tensor with dimensions that don't match the actual input. Could you explain why this should be the case? Intuitively to me, it feels like the only thing that should vary based on the permutation is the rearrangement of the axes in the middle (the argument to transpose).

Assuming you're right, I think a code comment explaining the algorithm would be very helpful here.

Just as a minor legibility note - lines 230 to 233 might be clearer as:

raw_shape_l = self._op_shape.dims_l()
n_dims_l = len(raw_shape_l)
shape_l = tuple(raw_shape_l[n_dims_l - n - 1] for n in reversed(inv_perm))

rather than reversing the tuple twice. The whole thing still feels mathematically odd to me, as described above, though.

[alexanderivrii]

I have applied these and the following suggestions, making the code simpler and more readable, thanks!

[jakelishman]

It's hard to tell which parts of this code are logically shared between the two branches and which parts are different, so it's a bit tricky to read. I think you could re-organise this to

inv_perm = np.argsort(perm)
if front:
    if len(inv_perm) != len(...):
        ...
    shape_l = ...
    shape_r = ...
    axes_l = ...
    axes_r = ...
    new_shape_l = ...
    new_shape_r = ...
else:
    if len(inv_perm) != len(...):
        ...
    shape_l = ...
    ...
split_shape = shape_l + shape_r
axes_order = axes_l + axes_r
new_matrix = self._data.reshape(split_shape).transpile(axes_order).reshape(self._op_shape.shape)
return Operator(new_matrix, input_dims=new_shape_r, output_dims=new_shape_l)

To me at least, this makes the logic clearer (assuming I'm correct) - the difference between front and not front is more obviously localised just to the construction of the shape / permutation axes, and the actual re-arrangement and construction logic is the same from then on.

[alexanderivrii]

Done.

[jakelishman]

It'd be good to do this loop via ddt parametrisation instead, so we can see individual test runs.

[alexanderivrii]

Done.

[jakelishman]

We probably don't need to test these here, since the test is about the exceptions. We could test their full validity in a separate test, though.

[alexanderivrii]

Done.

[jakelishman]

I think we also need an explicit test that applying only a one-sided permutation to an Operator with heterogeneous tensored qudit spaces applies correctly (although I've got to say, I don't have a huge amount of intuition for what a clear test of that looks like).

[alexanderivrii]

Done. Also added tests for applying the permutation [1, 0] to operator with dimensions (2, 3) x (2, 3), where things can be easily worked out explicitly.

[alexanderivrii]

For completeness, here is Gadi's construction + proof of correctness (though please note that his notations are slightly different from Qiskit's).

# Explicit Example:
P = (1,2,0)
P_inv = (2,0,1)
extended_P = (4,5,3)
extended_P_inv = (5,3,4)
id = (0,1,2)
extended_id = (3,4,5)
shape = (2,2,3)
shape_P = (2,3,2)
shape_P_inv = (3,2,2)

# For OP:
O_matrix.reshape(shape + shape_P_inv).transpose(id + extended_P_inv).reshape((12,12))

# For PO:
O_matrix.reshape(shape_P + shape).transpose(P + extended_id).reshape((12,12))

What's important is that for math to work out, the shape has to be also permuted (something that I have also independently found out after countless experiments).

[jakelishman]

Thanks for all this work, Sasha and Gadi. I'm going to approve this now - I'm happy that the maths we're using works out, though I don't think we're coding it in an intuitive way right now. The proof in this PR doesn't account for the Kronecker-product structure of the permutation that we're actually applying, while the code form we're doing is inherently related to that via transpose. I'm content that the maths is correct for a pure row/col permutation, but it's currently unsatisfying to me about why we have to permute the tensor shape before we split it right now.

That said, I didn't have time to run my own maths on this to find a form I found more pleasant, and the tests seem to indicate that this is correct, so I don't want to hold up a correct implementation over me not fully clicking with the maths.