Circuit extraction

[matulni]

This PR introduces functionality to perform the circuit extraction algorithm presented in Ref. [1] which allows to extract a circuit without ancillas from a strongly deterministic pattern.

Context

The main result of the paper is a procedure to perform the transformation

$$U_{\mathcal{P}} \rightarrow C V(\mathbf{\theta}),$$

where $U_{\mathcal{P}}$ is the isometry implemented by the pattern $\mathcal{P}$ (a unitary when the pattern has the same number of inputs and outputs), $C$ is a Clifford map, and $V(\mathbf{\theta})$ is an ordered product of Pauli exponentials.

The structure of the extraction process reads as follows:

flowchart LR

n0["Pattern"]
n1["OpenGraph"]
n2("Focused PauliFlow")

subgraph s1["ExtractionResult"]
    s10("PauliExponentialDAG")
    s11("CliffordMap")
end

%% Junction node to merge arrows
j1(( ))

n3("Graphix circuit")

n0 --> n1
n1 --> n2
n2 --> s10
n2 --> s11

%% Merging into the junction
s10 --> j1
s11 --> j1
j1 --> n3

%% Styling the junction to be small
style j1 fill:#000,stroke-width:0px,width:10px

Loading

• The transformation Pattern -> OpenGraph -> Focused PauliFlow already exists in the current implementation of Graphix. The $O(N^3)$ Pauli-flow extraction algorithm in Ref. [2] always returns focused Pauli flows, contrary to the $O(N^5)$ version in Ref. [1] which requires an additional post-processing to focus the correction function.

• ExtractionResult is the output of the algorithm which is denoted "PdDAG" in [1]. In particular:

  • PauliExponentialDAG is a mapping between measured nodes in the open graph and Pauli exponentials, and a partial order between the measured node (the flow's partial order). Pauli exponentials are a pair of an Angle (the node's measurement angle) and a PauliString acting on the output nodes of the open graph. The corresponding Pauli string is obtained from the focused flow's correction function.

  • Clifford Map describes a linear transformation between the space of input qubits and the space of output qubits. It is encoded as a map from the Pauli-group generators ($X$ and $Z$) over the input nodes to Pauli strings over the output nodes.

    The $Z$-map is also obtained from the focused flow's correction function. The $X$-map is obtained from the focused flow's correction function of the extended open graph (see Def. C.3 in [1]).

Summary of changes

• Added the following to graphix.flow.core.PauliFlow:

  • is_focused : (method) Verifies if a Pauli flow is focused according to definition 4.3 in Ref. [1].
  • extract_circuit : (method) Returns a ExtractionResult object.
  • pauli_strings : (cached_property) Associates a Pauli string to every corrected node in the flow's correction function.

• Added new module graphix.circ_ext.extraction

  • This module contains the machinery to perfom the transformation Focused PauliFlow -> ExtractionResult.
  • We introduce the following (self-descriptive) classes:
    • ExtractionResult
    • PauliString
    • PauliExponential
    • PauliExponentialDAG
    • CliffordMap

• Added new module graphix.circ_ext.compilation

  • This module contains the machinery to perform the transformation ExtractionResult -> Graphix circuit.
  • It lays out a structure of abstract compilation passes which will allow in the future to experiment with different strategies. Currently:
    • The transformation PauliExponentialDAG -> Graphix circuit is done with the naive "star" construction (see [3]).
    • The transformation CliffordMap -> Graphix circuit relies on stim and currently only supports unitaries (i.e., patterns with the same number of inputa and outputs). To avoid introducing a dependency on stim, the StimCliffordPass is implemented in the external plugin graphix-stim-compiler. This compilation pass together with the tooling introduced in this PR allow to do a round-trip conversion circuit -> MBQC pattern -> circuit.

We note that all the transformations in the figure work with parametric objects as well.

References

[1] Simmons, 2021, arXiv:2109.05654.

[2] Mitosek and Backens, 2024, arXiv:2410.23439.

[3] https://quantumcomputing.stackexchange.com/questions/5567/circuit-construction-for-hamiltonian-simulation/11373#11373

[codecov]

Codecov Report

❌ Patch coverage is 78.57143% with 42 lines in your changes missing coverage. Please review.
✅ Project coverage is 88.75%. Comparing base (06ba6a8) to head (6ec46ed).
⚠️ Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
graphix/circ_ext/compilation.py	62.50%	21 Missing ⚠️
graphix/circ_ext/extraction.py	86.79%	14 Missing ⚠️
graphix/flow/core.py	76.00%	6 Missing ⚠️
graphix/transpiler.py	88.88%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #445      +/-   ##
==========================================
- Coverage   89.06%   88.75%   -0.31%     
==========================================
  Files          44       46       +2     
  Lines        6547     6743     +196     
==========================================
+ Hits         5831     5985     +154     
- Misses        716      758      +42     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[thierry-martinez]

The following code is a bit simpler and faster if the graph is sparse (but slower if the graph is dense; however, it would be much faster in any case if we used rustworkx!):

Suggested change

	for edge in combinations(c_set, 2):
	negative_sign ^= edge in og.graph.edges()
	negative_sign ^= graph.subgraph(c_set).number_of_edges() % 2 == 1

[matulni]

Nice, thanks. Done in f99c285. Do you suggest moving to rustworks (adding this to the backlog)?

[matulni]

Hi @thierry-martinez, thanks for the review. I addressed all your comments, the corresponding modifications in the stim plugin are in commit matulni/graphix-stim-compiler@843601c. I think it's ready for another round!

[thierry-martinez]

Some minor comments.

[thierry-martinez]

Suggested change

	For every input node :math:`v`, a new node :math:`u` and edge :math:`(u, v)` are added to the open graph. Node :math:`u` is measured in plane :math:`XY` with angle :math:`\alpha = 0` and replaces :math:`v` in the open graph's sequence of input nodes.
	For every input node :math:`v`, a new node :math:`u` and edge :math:`(u, v)` are added to the open graph. Node :math:`u` is measured in Pauli axis `X` and replaces :math:`v` in the open graph's sequence of input nodes.

[matulni]

Fixed in 2b4a008

[thierry-martinez]

All tests pass with the following change:

Suggested change

	modified_qubits = [
	qubit
	for qubit in range(circuit.width)
	if qubit in pexp.pauli_string.x_nodes | pexp.pauli_string.y_nodes | pexp.pauli_string.z_nodes
	]
	modified_qubits = list(pexp.pauli_string.x_nodes | pexp.pauli_string.y_nodes | pexp.pauli_string.z_nodes)

I realize we discussed this earlier and the ladder transformation output does change with this modification. However, if the order matters, should we add a test to enforce it?

[matulni]

Hi @thierry-martinez, this is tricky.

From a code perspective, I think that in the current Python implementation, calling list on a set of (small?) integers always returns an ordered list:

assert [0, 1, 2] == list[{1, 0, 2}]

Ultimately, I think that this is because (small?) integers are hashed with their own value:

a, b = 0, 1
assert a.__hash__() == 0
assert b.__hash__() == 1

So in our case, since we are dealing with small integers (number of qubits), in both the PR's implementation and your suggestion, we get the same list of instructions (and hence all tests pass trivially).

Below, I paste a pexp_ladder_pass function where you can specify the qubit order by hand:

from itertools import pairwise, chain
from graphix.circ_ext.extraction import PauliExponential, PauliExponentialDAG
from graphix.transpiler import Circuit
from graphix.fundamentals import ANGLE_PI


def pexp_ladder_pass(pexp_dag: PauliExponentialDAG, circuit: Circuit, modified_qubits: list[int]) -> None:
    r"""Add a Pauli exponential DAG to a circuit by using a ladder decomposition.

    The input circuit is modified in-place. This function assumes that the Pauli exponential DAG has been remap, i.e., its Pauli strings are defined on qubit indices instead of output nodes. See :meth:`PauliString.remap` for additional information.

    Parameters
    ----------
    pexp_dag: PauliExponentialDAG
        The Pauli exponential rotation to be added to the circuit. Its Pauli strings are assumed to be defined on qubit indices.
    circuit : Circuit
        The circuit to which the operation is added. The input circuit is assumed to be compatible with ``pexp_dag.output_nodes``.

    Notes
    -----
    Pauli exponentials in the DAG are compiled sequentially following an arbitrary total order compatible with the DAG. Each Pauli exponential is decomposed into a sequence of basis changes, CNOT gates, and a single :math:`R_Z` rotation:

    .. math::

        R_Z(\phi) = \exp \left(-i \frac{\phi}{2} Z \right),

    with effective angle :math:`\phi = -2\alpha`, where :math:`\alpha` is the angle encoded in `self.angle`. Basis changes map :math:`X` and :math:`Y` operators to the :math:`Z` basis before entangling the qubits in a CNOT ladder.

    Gate set: H, CNOT, RZ, RY

    See https://quantumcomputing.stackexchange.com/questions/5567/circuit-construction-for-hamiltonian-simulation/11373#11373 for additional information.
    """

    def add_pexp(pexp: PauliExponential, circuit: Circuit) -> None:
        r"""Add the Pauli exponential unitary to a quantum circuit.

        This function modifies the input circuit in-place.

        Parameters
        ----------
        circuit : Circuit
            The quantum circuit to which the Pauli exponential is added.

        Notes
        -----
        It is assumed that the ``x``, ``y``, and ``z`` node sets of the Pauli string in the exponential are well-formed, i.e., contain valid qubit indices and are pairwise disjoint.
        """
        if pexp.angle == 0:  # No rotation
            return

        angle = -2 * pexp.angle * pexp.pauli_string.sign

        if len(modified_qubits) == 0:  # Identity
            return

        q0 = modified_qubits[0]

        if len(modified_qubits) == 1:
            if q0 in pexp.pauli_string.x_nodes:
                circuit.rx(q0, angle)
            elif q0 in pexp.pauli_string.y_nodes:
                circuit.ry(q0, angle)
            else:
                circuit.rz(q0, angle)
            return

        add_basis_change(pexp, q0, circuit)

        for q1, q2 in pairwise(modified_qubits):
            add_basis_change(pexp, q2, circuit)
            circuit.cnot(control=q1, target=q2)

        circuit.rz(modified_qubits[-1], angle)

        for q2, q1 in pairwise(modified_qubits[::-1]):
            circuit.cnot(control=q1, target=q2)
            add_basis_change(pexp, q2, circuit)

        add_basis_change(pexp, modified_qubits[0], circuit)

    def add_basis_change(pexp: PauliExponential, qubit: int, circuit: Circuit) -> None:
        """Apply an X or a Y basis change to a given qubit if required by the Pauli string.

        This function modifies the input circuit in-place.

        Parameters
        ----------
        pexp : PauliExponential
            The Pauli exponential under consideration.
        qubit : int
            The qubit on which the basis-change operation is performed.
        circuit : Circuit
            The quantum circuit to which the basis change is added.
        """
        if qubit in pexp.pauli_string.x_nodes:
            circuit.h(qubit)
        elif qubit in pexp.pauli_string.y_nodes:
            add_hy(qubit, circuit)

    def add_hy(qubit: int, circuit: Circuit) -> None:
        """Add a pi rotation around the z + y axis.

        This function modifies the input circuit in-place.

        Parameters
        ----------
        qubit : int
            The qubit on which the basis-change operation is performed.
        circuit : Circuit
            The quantum circuit to which the basis change is added.
        """
        circuit.rz(qubit, ANGLE_PI / 2)
        circuit.ry(qubit, ANGLE_PI / 2)
        circuit.rz(qubit, ANGLE_PI / 2)

    for node in chain(*reversed(pexp_dag.partial_order_layers[1:])):
        pexp = pexp_dag.pauli_exponentials[node]
        add_pexp(pexp, circuit)

If you compile a Pauli exponential with different modified_qubit order, you get a different circuit (with the PR's implementation or your commit you always get the same set of instructions):

from graphix.circ_ext.compilation import pexp_ladder_pass
from graphix.circ_ext.extraction import PauliString, PauliExponential, PauliExponentialDAG
from graphix.transpiler import Circuit
from graphix.sim.base_backend import NodeIndex
from graphix.fundamentals import ANGLE_PI

alpha = 0.1

pexp = PauliExponentialDAG(
                    pauli_exponentials={
                        0: PauliExponential(alpha / 2, PauliString(x_nodes={4,3}, y_nodes={1, 2})),
                    },
                    partial_order_layers=[{1, 2, 3, 4}, {0}],
                    output_nodes=[1, 2, 3, 4],
                )


outputs_mapping = NodeIndex()
outputs_mapping.extend(pexp.output_nodes)
pexp_rmp = pexp.remap(outputs_mapping.index)

qc1 = Circuit(4)
pexp_ladder_pass_2(pexp_rmp, qc1, [0, 1, 2, 3])

qc2 = Circuit(4)
pexp_ladder_pass_2(pexp_rmp, qc2, [0, 2, 1, 3])
assert qc1.instruction != qc2.instruction

But to my surprise, they seem to represent the same unitary, up to a phase:

s1 = qc1.simulate_statevector().statevec
s2 = qc2.simulate_statevector().statevec
assert s1.isclose(s2)

In summary:

• I think that Python's implementation makes it difficult to check the different between ordering the qubits because of how integers are hashed.
• I haven't managed to find an example where the order matters when we specify it by hand.
• I'm not convinced that the order does not matter in general.

[matulni]

Edit: not even up to a phase:

import numpy as np
assert(np.all(np.isclose(s1.flatten(),s2.flatten())))

[matulni]

Ok, I convinced myself that the order on which we iterate over the qubits when synthesising the Pauli exponential does not matter.

Consider the ZZ-Pauli exponential:

$$U(\alpha) = \exp \left(i \frac{\alpha}{2} Z_0 Z_1 \right) =\text{diag}\left[e^{i \alpha/2}, e^{-i \alpha/2}, e^{-i \alpha/2}, e^{i \alpha/2}\right].$$

The star-pass synthesis in the order $[q_0, q_1]$ yields

q_0: ──■─────────────■──
     ┌─┴─┐┌───────┐┌─┴─┐
q_1: ┤ X ├┤ Rz(α) ├┤ X ├
     └───┘└───────┘└───┘

whereas the synthesis in the order $[q_1, q_0]$ yields

     ┌───┐┌───────┐┌───┐
q_0: ┤ X ├┤ Rz(α) ├┤ X ├
     └─┬─┘└───────┘└─┬─┘
q_1: ──■─────────────■──

and these circuits represent the same unitary. To see that, first, note that they are equivalent up to a relabelling $q_0 \leftrightarrow q_1$. Then, notice that according to Eq. (1) the computational-basis states that are sensitive to a relabelling, i.e., $|01\rangle$ and $|10\rangle$, pick up the same phase, $e^{-i \alpha/2}$. In other words,

$$ U(\alpha) = \text{SWAP} \cdot U(\alpha) \cdot \text{SWAP}^{\dagger}.$$

This is true for $N$-qubit case too $U_N(\alpha) = \exp \left(i \frac{\alpha}{2} Z_0 \dots Z_N \right)$, since all the computational basis states that are connected by a permutation have the same number of non-zero bits, and therefore pick up the same phase, $U(\alpha) |n\rangle = e^{(-1)^m i \alpha/2} |n\rangle$, with $m$ the number of 1s in bin(n). If there are $M$-identities in the Pauli string, we go back to the case $U_{N-M}(\alpha)$, so we can prove it by induction.

The argument does not change if we have $X$ or $Y$ gates in the Pauli string, since we simply proceed by doing a change of basis (i.e., applying $H$ or $H_Y$ on the corresponding qubit), which is independent of the synthesis order.

[matulni]

Unnecessary sorting removed in 2b4a008

[thierry-martinez]

Thank you very much for the explanation! Regarding the points you mentioned in the comment above, hash is indeed the identity over $[0, 2^{60}]$. However, the iteration order of a set is unspecified, and there is no guarantee that elements are enumerated in order of their hash values, nor that this order is stable from one Python version to another: for instance, list(set([1,2])) == [1,2] whereas list(set([1,8])) == [8,1].

[thierry-martinez]

Since x_nodes, y_nodes, and z_nodes are disjoint, I propose representing PauliString as Mapping[int, Axis], which is more canonical. See commit thierry-martinez@828f96d (and in thierry-martinez/graphix-stim-compiler@a743335 for graphix-stim-compiler).

[thierry-martinez]

(By "more canonical", I mean that disjointness is structurally enforced.)

[matulni]

Hi @thierry-martinez, thanks for this comment, indeed it bothered me that disjointness was not structurally enforced. I like your mapping idea, but I fear that using Axis may be a bit confusing: in Graphix this object is associated with measurements (or at least, with a notion of space), and I'm not sure if it makes much sense in this context.
What do you think about using graphix.Pauli objects ? Another advantage would be that we can explicitly map nodes to the identity (now, they are implicitly mapped by not being in the sets (or in the dictionary keys in your commit)). Including identity operators will certainly be less efficient, but all the Pauli-string processing routines will remain with the same complexity, so I think it's not a big problem (and it may improve readability). If you agree, I'll adapt your commit accordingly.

[thierry-martinez]

I think that Axis expresses precisely what we want to enforce structurally (at least for the value part; keys are not required to be actual nodes). Axis is defined in fundamentals and is a subtype of IXYZ, which is a component of Pauli. As you mentioned, using IXYZ or Pauli would allow nodes to be mapped to I; therefore, we lose canonicality because mapping to I is equivalent to not being mapped. This should not be a serious problem, but the invariant that the set of keys equals precisely the nodes not mapped to I is nice to have, and defining a rotation over I by a given angle is a bit strange (it would be the identity, of course). Pauli is a pair of IXYZ and Sign, and we will not use the Sign either, which is also a loss of canonicality.

[thierry-martinez]

(In other words, Pauli is a Pauli string of length 1.)

[matulni]

Fair enough, I'm convinced, thanks for the comment :)

[pranav97nair]

Small nitpick - I would rename this variable pexp_dag to match the convention you use everywhere else in the codebase and to minimize the chances of it being confused with the PauliExponential object.

[matulni]

Done in 9137a7d

[pranav97nair]

Great job on this PR, Mateo! I have one general comment regarding the assumption in several class functions in graphix.circ_ext.extraction that their flow argument is a focused PauliFlow. I agree broadly with the approach of leaving the responsibility of ensuring that the Pauli flow is focused to the user and allowing them to obtain an incorrect extraction result if it is not. But I wonder if it would be a good feature to warn the user if this is the case, similarly to how running Pauli pre-simulation on a pattern issues a warning if the pattern has non-inferred Pauli measurements. Perhaps it is worth discussing whether such a warning would be worth the cost of calling PauliFlow.isfocused() in some strategic places.

[thierry-martinez]

LGTM! Some minor remarks.

[thierry-martinez]

Suggested change

	This routine associates a two Pauli strings (one per generator of the Pauli group, X and Z) to each input node in ``flow.og``.
	This routine associates two Pauli strings (one per generator of the Pauli group, X and Z) to each input node in ``flow.og``.

[thierry-martinez]

Suggested change

	Compilation pass to synthetize a Clifford map. If ``None`` (default), a `ValueError` is raised since there is still no default pass for Clifford map integrated in Graphix.
	Compilation pass to synthetize a Clifford map. If ``None`` (default), a ``ValueError`` is raised since there is still no default pass for Clifford map integrated in Graphix.

[thierry-martinez]

Suggested change

	The identity operator is ommitted in this representation, which means that in general it is not possible to infer the size of the Hilbert space from an instance of `PauliString` alone.
	The identity operator is ommitted in this representation, which means that in general it is not possible to infer the size of the Hilbert space from an instance of ``PauliString`` alone.

[thierry-martinez]

Suggested change

	Primary extraction string associated to the input measured nodes. The sets in the returned `PauliString` instance are disjoint.
	Primary extraction string associated to the input measured nodes. The sets in the returned ``PauliString`` instance are disjoint.

[thierry-martinez]

Suggested change

	Partial order between the Pauli exponentials in a layer form. The set `layers[i]` comprises the nodes in layer `i`. Nodes in layer `i` are "larger" in the partial order than nodes in layer `i+1`. The pattern's output nodes are always in layer 0.
	Partial order between the Pauli exponentials in a layer form. The set ``layers[i]`` comprises the nodes in layer ``i``. Nodes in layer ``i`` are "larger" in the partial order than nodes in layer ``i+1``. The pattern's output nodes are always in layer 0.

[thierry-martinez]

@pranav97nair You make a good point: PauliExponentialDAG.from_focused_flow requires the flow to be focused but does not verify this. When the flow is obtained from OpenGraph.extract_pauli_flow, we don't need such a check since Pauli-flow extraction always returns a focused Pauli flow by construction, but this may not hold for arbitrary user-constructed flows. Note that we don't even check that the flow is well-formed. The current approach merely relies on the function name from_focused_flow to indicate that the flow should be focused (and well-formed, of course). Another option would be to run is_well_formed and is_focused systematically, but that would lead to unnecessary checks in the common cases where flows come from extract_pauli_flow. I would tend to propose introducing FocusedPauliFlow and FocusedGFlow classes to indicate that a flow is focused at the type level, but I think we can offload this discussion from the current PR and open a new issue on that topic once this PR is merged.