Constant-depth preparation of matrix product states with dynamic circuits

Matrix product states (MPS)

Matrix product states (MPS) are an important class of quantum states. They can be described and manipulated conveniently on classical computers and are able to approximate relevant states in quantum many-body systems. In particular, MPS can efficiently describe ground states of (gapped local) one-dimensional Hamiltonians, which in addition can be found efficiently using density matrix renormalization group (DMRG) algorithms. Also check out our introduction to matrix product states, and see 3 for a review of MPS.

Following 1, we will look at translation-invariant MPS of a quantum \(N\)-body system where each body, or site, has local (physical) dimension \(d.\) We will further restrict to periodic boundary conditions, corresponding to translational invariance of the MPS on a closed loop. A general MPS with these properties is given by

where \(\vec{m}\) is a multi-index of \(N\) indices ranging over \(d,\) i.e., \(\vec{m}\in\{0, 1 \dots, d-1\}^N,\) and \(\{A^{m}\}_{m}\) are \(d\) square matrices of equal dimension \(D.\) (Note that we are using the same \(A^m\) for each physical site, due to the translational invariance.) This dimension \(D\) is called the bond dimension, which is a crucial quantity for the expressivity and complexity of the MPS. We see that \(|\Psi\rangle\) is fully specified by the \(d\times D\times D=dD^2\) numbers in the rank-3 tensor \(A.\) However, this specification is not unique. We remove some of the redundancies by assuming that \(A\) is in so-called left-canonical form, meaning that it satisfies

We will not discuss additional redundancies here but refer to 1 and the mentioned reviews for more details.

Example

As an example, consider a chain of \(N\) qubits \((d=2)\) and an MPS \(|\Psi(g)\rangle\) on this system with \(D=2,\) defined by the matrices

where \(\eta = \frac{1}{\sqrt{1-g}}\) and \(g\in[-1, 0]\) is a freely chosen parameter.

This MPS is a simple yet important example (also discussed in Sec. III C 1. of 1) because it can be tuned from the GHZ state \(|\Psi(0)\rangle=\frac{1}{\sqrt{2}}(|0\rangle^{\otimes N} + |1\rangle^{\otimes N})\) to the product state \(|\Psi(-1)\rangle=|+\rangle^{\otimes N},\) two states that differ greatly in their physical properties.

One such property is the correlation length. For the GHZ state, a correlator between sites \(0\) and \(j\) is

That is, the correlation does not decay with the distance \(j,\) but the correlation length is infinite. For the product state, on the other hand, we get

so that correlations decay “instantaneously.” The correlation length vanishes. Long-range correlated states require linear-depth circuits if we are restricted to unitary operations. Therefore, the constant-depth circuit for \(|\Psi(0)\rangle\) will demonstrate that dynamic quantum circuits, which are not purely unitary, are more powerful than unitary operations alone!

The MPS \(|\Psi(g)\rangle\) will be our running example throughout this demo. To warm up our coding, let’s define the tensor \(A\) and test that it is in left-canonical form.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch

import pennylane as qml


def A(g):
    """Compute the tensor A in form of d=2 matrices with shape (D, D) = (2, 2) each."""
    eta = 1 / np.sqrt(1 - g)
    A_0 = np.array([[1, 0], [np.sqrt(-g), 0]]) * eta
    A_1 = np.array([[0, -np.sqrt(-g)], [0, 1]]) * eta
    return (A_0, A_1)


g = -0.6
A_0, A_1 = A(g)
is_left_canonical = np.allclose(A_0.conj().T @ A_0 + A_1.conj().T @ A_1, np.eye(2))
print(f"The matrices A^m are in left-canonical form: {is_left_canonical}")

xi = 1 / abs(np.log((1 + g) / (1 - g)))
print(f"For {g=}, the theoretical correlation length is {xi=:.4f}")

The matrices A^m are in left-canonical form: True
For g=-0.6, the theoretical correlation length is xi=0.7213

Sequential preparation circuit

An MPS like the one above can be prepared in linear depth using an established technique by Schön et al. 2. We introduce this technique here because the new constant-depth construction uses it as a (parallelized) subroutine, and because we will later compare the two approaches.

The MPS \(|\Psi\rangle\) above is given by the rank-3 tensor \(A.\) We can think of this as an operator that acts on a bond site and simultaneously creates a physical site. However, a unitary quantum circuit is not allowed to create qudits out of nowhere. To fix this, we construct a unitary operation from \(A\) that already takes a physical site as input. This means that we build a new rank-4 tensor \(U\) with an additional axis of dimension \(d.\) Then we demand that \(U\) acts like \(A\) if this new physical input axis is in the state \(|0\rangle.\) That is, after applying \(U\) to a physical site in \(|0\rangle\) and a bond site, we obtain the same result as applying \(A\) to the bond site alone.

We leave the action of \(U\) on other physical input states undefined, and only require it to make \(U\) unitary overall. In short, our construction can be written as

where \(C_\perp\) can be any operator making \(U\) unitary.

In the definition of the MPS, the bond site is passed on between copies of \(A,\) sequentially creating physical sites. Likewise, we will apply our unitary \(U\) to a sequence of physical sites (each in the initial state \(|0\rangle\)) while passing on a single bond site. This way, we chain up the unitaries on the bond site like beads on a string.

Alongside the correspondence between \(U\) acting on \(|0\rangle\) and \(A,\) these steps are visualized in the sketch above. Note that the sketch deviates from standard circuit diagrams: Usually each wire corresponds to one qudit in the algorithm and the position of the qudits is fixed through the wire position. Here, the two bond qudits start as the top two wires, and one of them is moved to the bottom throughout the circuit.

The sequential preparation circuit now consists of the following steps.

1. Start in the state \(|\psi_0\rangle = |0\rangle^{\otimes N}\otimes |00\rangle_D.\) The last two sites are non-physical bond sites, encoded by \(D\)-dimensional qudits.

2. Entangle the bond qudits into the state \(|I\rangle = \frac{1}{\sqrt{D}}\sum_j |jj\rangle_D.\)

3. Apply the unitary \(U\) to each of the \(N\) physical sites, with the \(D\)-dimensional tensor factor always acting on the first bond qudit. Denoting \(U\) acting on the \(i\)-th physical site and the first bond site by \(U^{(i)},\) this produces the following state.

   \[\begin{split}|\psi_1\rangle &= \prod_{i=N}^{1} U^{(i)} |0\rangle^{\otimes N} \frac{1}{\sqrt{D}} \sum_j|jj\rangle _D\\ &= \frac{1}{\sqrt{D}} \sum_j \prod_{i=N-1}^{1} U^{(i)} \sum_{m_N=0}^{d-1} |0\rangle^{\otimes N-1}|m_N\rangle A^{m_N}|jj\rangle _D\\ &= \frac{1}{\sqrt{D}} \sum_j \sum_{\vec{m}} |\vec{m}\rangle A^{m_1} A^{m_2}\cdots A^{m_N}|jj\rangle_D\\ &= \frac{1}{\sqrt{D}} \sum_{j,k} \sum_{\vec{m}} |\vec{m}\rangle \langle j|A^{m_1} A^{m_2}\cdots A^{m_N}|k\rangle |jk\rangle_D\end{split}\]

   We can see how each \(U^{(i)}\) contributes a factor of \(A\) and a corresponding state \(|m_i\rangle\) on the \(i\)-th physical site. The same would have come out when applying the 3-tensor \(A\) to the first bond site \(N\) times.

4. Measure the two bond qudits in the (generalized) Bell basis and postselect on the outcome being \(|I\rangle.\) Then discard the bond qudits, which collapses \(|\psi_1\rangle\) into the state

   \[|\psi_2\rangle = \sum_{\vec{m}} \operatorname{tr}[A^{m_1}A^{m_2}\cdots A^{m_N}]|\vec{m}\rangle = |\Psi\rangle.\]

   Note that this step is probabilistic and we only succeed to produce the state if we measure the state \(|I\rangle.\)

We see that we prepared \(|\Psi\rangle\) with an entangling operation on the bond qudits, one unitary per physical qubit, and a final basis change to measure the Bell basis. Overall, this amounts to a linear operation count and circuit depth.

Example

For our example MPS \(|\Psi(g)\rangle,\) let us first find the unitary \(U.\) For this, consider the fixed part of \(U.\)

E_00 = np.array([[1, 0], [0, 0]])  # |0><0|
E_10 = np.array([[0, 0], [1, 0]])  # |1><0|
U_first_term = np.kron(A_0, E_00) + np.kron(A_1, E_10)

print(np.round(U_first_term, 5))
print(np.linalg.norm(U_first_term, axis=0))

[[ 0.79057  0.       0.       0.     ]
 [ 0.       0.      -0.61237  0.     ]
 [ 0.61237  0.       0.       0.     ]
 [ 0.       0.       0.79057  0.     ]]
[1. 0. 1. 0.]

We see that this fixed part has two norm-1 columns already, i.e., it maps the input states \(|00\rangle\) and \(|10\rangle\) to normalized, orthogonal quantum states. We can complement this by mapping \(|01\rangle\) and \(|11\rangle\) to two other normalized states that are mutually orthogonal again. This turns out to be easy for the small example here and leads to the following:

E_01 = np.array([[0, 1], [0, 0]])  # |0><1|
E_11 = np.array([[0, 0], [0, 1]])  # |1><1|
C_perp = np.kron(A_1, E_01) + np.kron(A_0, E_11)

U = U_first_term + C_perp

print(np.round(U, 5))
print(f"\nU is unitary: {np.allclose(U.conj().T @ U, np.eye(4))}")

[[ 0.79057  0.       0.      -0.61237]
 [ 0.       0.79057 -0.61237  0.     ]
 [ 0.61237  0.       0.       0.79057]
 [ 0.       0.61237  0.79057  0.     ]]

U is unitary: True

Great! This means that we already found the unitary \(U\) for this MPS. We will implement it as the matrix that it is via QubitUnitary. If you’re curious, you can try to decompose it into elementary gates (hint: you only need two!).

Now let’s code up the entire sequential preparation circuit (we will be applying the unitary \(U\) to the second instead of the first bond qubit). We implement the measurement and postselection step as a separate function project_measure for easier reusability later on.

def label_fn(self, *_, **__):
    """Report the physical wire this operation acts on in its label."""
    return f"$U^{{({self.wires[1]})}}$"


qml.QubitUnitary.label = label_fn


def sequential_preparation(g, wires):
    """Prepare the example MPS Ψ(g) on N qubits where N is the length
    of the passed wires minus 2. The bond qubits are still entangled."""
    eta = 1 / np.sqrt(1 - g)
    # Define bond qubits [0, N+1] and physical qubits [1, 2, ... N]
    bond0, bond1 = wires[0], wires[-1]
    phys_wires = wires[1:-1]

    # Prepare bond qubits in the Bell state (|00>+|11>)/sqrt(2)
    qml.Hadamard(bond0)
    qml.CNOT([bond0, bond1])

    # Apply unitary U to second bond qubit and each physical qubit
    for phys_wire in phys_wires:
        u = qml.QubitUnitary(U, wires=[bond1, phys_wire])


def project_measure(wire_0, wire_1):
    """Measure two qubits in the Bell basis and postselect on (|00>+|11>)/sqrt(2)."""
    # Move bond qubits from Bell basis to computational basis
    qml.CNOT([wire_0, wire_1])
    qml.Hadamard(wire_0)
    # Measure in computational basis and postselect |00>
    qml.measure(wire_0, postselect=0)
    qml.measure(wire_1, postselect=0)


dev = qml.device("default.qubit")


@qml.qnode(dev)
def sequential_circuit(N, g):
    """Run the preparation circuit and projectively measure the bond qubits."""
    wires = list(range(N + 2))
    sequential_preparation(g, wires)
    project_measure(0, N + 1)
    return qml.probs(wires=wires[1 : N + 1])

We will use this circuit later to compare the constant-depth circuit against it. For now, let’s draw it.

N = 7
_ = qml.draw_mpl(sequential_circuit)(N, g)

The sequential preparation circuit already uses mid-circuit measurements, as is visible from the measurement steps on the first and last qubit, which are set to postselect the outcome \(|0\rangle.\) However, there is no feedforward control that modifies the circuit dynamically based on measured values. This will change in the following sections.

Fusion of MPS states with mid-circuit measurements

The next ingredient for the constant-depth MPS preparation circuit is to fuse the product state of two MPS states together into one (entangled) MPS. Recall that a single sequential creation of this final MPS would pass a bond site between the first and second group of physical sites, chaining them like beads on the same string. The fusion step allows us to replace this connection with dynamic circuit components after creating two independent MPS. This is like knotting together two strings on which we chained beads independently before. As a matter of fact, chaining beads on smaller strings and the fact that we can knot them together simultaneously is at the heart of going from linear to constant circuit depth.

Before we dive into the mathematical details, we visualize the idea behind the fusion step schematically:

Consider a state \(|\Phi\rangle=|\Phi_1\rangle\otimes|\Phi_2\rangle,\) where

for \(r=1,2\) are MPS on \(N_r\) qudits, respectively. We assume that they have been prepared with the sequential technique from above, but that the bond qudits have not been measured yet. In particular, the postselection step has not been performed yet. Using the form of \(|\psi_1\rangle\) from the recipe above, we can write this state as

Now we want to measure two bond qudits in an entangled basis, one of \(|\Phi_1\rangle\) and \(|\Phi_2\rangle\) each (e.g., the ones indexed by \(k\) and \(\ell\) above). Assume this basis to be given in the form of (orthogonal) states

which in turn are characterized by \(D\times D\) matrices \(B^k.\) We can change the basis on the bond qudits to be measured by applying the unitary \(V=\sum_k |k\rangle\!\langle B^k|,\) which leads to the state

where we used that \(\sum_{j,\ell} |j\rangle B^k_{j\ell}\langle\ell|=B^k.\) We now can measure the two bond qudits and will know that for a given measurement outcome \(k,\) we obtained the MPS state on \(K+L\) qubits, up to a defect, namely the matrix \(B^k\) that depends to the outcome.

In full generality, the fusion strategy may seem complicated, but it can take a very simple form, as we will now see in our example.

Example

We will use the Bell basis to perform the mid-circuit measurement, i.e.,

This choice must seem arbitrary at this point, but will become clear later on. There, we will see that in general the matrices \(B^k\) and the MPS tensor \(A\) have to “play well” together. As an example, this choice leads to \(B^1_{j\ell} = 1\) for \((j,\ell)=(0, 1)\) and \((j, \ell)=(1,0),\) and \(B^1_{j\ell}=0\) otherwise (note that the factor \(\frac{1}{\sqrt{2}}\) is not part of \(B^k\) in the general equation above). This means that \(B^1=X,\) the Pauli matrix! Similarly, we find

i.e., all standard Pauli matrices. The measurement procedure can then be coded up in the following short function, which is similar to the project_measure function from the sequential preparation routine. Note, however, that instead of postselecting, we record the measurement outcome in the form of two bits representing the index \(k\) of the operator :math:`B^k, which we will use further below.

def fuse(wire_0, wire_1):
    """Measure two qubits in the Bell basis and return the outcome
    encoded in two bits."""
    qml.CNOT([wire_0, wire_1])
    qml.Hadamard(wire_0)
    return np.array([qml.measure(w) for w in [wire_0, wire_1]])

Fusing two MPS together that have been prepared by the sequential routine now is really simple:

def two_qubit_mps_by_fusion(g):
    sequential_preparation(g, [0, 1, 2])
    sequential_preparation(g, [3, 4, 5])
    mcms = fuse(2, 3)

However, as mentioned above, the produced state is not quite the MPS state on two qubits, because of the impurities caused by the fusion measurement. We can check this by undoing the fusion-based preparation with the two-qubit sequential circuit and measuring an exemplary expectation value. If the two separately prepared MPS and the fusion step did prepare the correct MPS already, the test measurement would just be \(\langle 00 | Z_0| 00\rangle=1.\)

@qml.qnode(qml.device("default.qubit", wires=6))
def prepare_and_unprepare(g):
    two_qubit_mps_by_fusion(g)
    # The bond qubits for a sequential preparation are just 0 and 5
    # The bond qubits 2 and 3 have been measured out in the fusion protocol
    qml.adjoint(sequential_preparation)(g, [0, 1, 4, 5])
    return qml.expval(qml.PauliZ(1))


test = prepare_and_unprepare(g)
print(f"Test measurement of fusion preparation + sequential unpreparation is {test:.2f}")

Test measurement of fusion preparation + sequential unpreparation is 0.00

This means we still need a way to remove the defect matrices \(B^k\) from the state that ended up in the prepared state when we performed the fusion measurement. Operator pushing will allow us to do exactly that!

Operator pushing

The last building block we need is so-called operator pushing. In fact, whether or not an MPS can be prepared with the constant-depth circuit depends on whether certain operator pushing relations are satisfied between the matrices \(B^k\) from the fusion step and the MPS tensor \(A.\) We will not go into detail about these conditions but assume here that we work with compatible MPS and measurements.

Operator pushing then allows us to push a defect operator \(B^k\) from one bond axis through the tensor \(A\) of the MPS to the other bond axis and to the physical axis. In general, doing so will change the operator, i.e., we end up with different operators \(C^k\) and \(D^k\) on the bond and physical axes, respectively. Visually, we find:

For simplicity, we denote a push by \(B^k\mapsto (C^k, D^k).\)

Example

In the fusion step above we picked the Bell basis as measurement basis, without further specifying why. As we saw, this basis leads to Pauli matrices as impurities \(B^k.\) It turns out that those matrices satisfy simple pushing relations together with the tensor \(A\) of our example MPS \(|\Psi(g)\rangle,\) explaining this choice of measurement basis. In particular, the relations are

That is, in the notation above we have \(X\mapsto (Y, Y),\) \(Y\mapsto (Y, X),\) and \(Z\mapsto (I, Z).\) Note that the bond axis operator either is trivial (the identity), or a Pauli \(Y\) operator. These pushing relations allow us to push the bond axis operator through to an open boundary bond axis of the fused MPS. In addition, we can remove the operators on the physical axes by applying the correct Pauli operation to the physical site, conditioned on the fusion measurement outcomes. If we have an MPS on \(L\) sites, pushing through a Pauli will have the following effect:

• If it is a Pauli \(Y,\) all physical sites need to be corrected by a Pauli \(X,\) and a Pauli \(Y\) is pushed to the opposite bond site.

• If it is a Pauli \(X,\) a Pauli \(Y\) is applied to the first physical qubit, and afterwards the case above is recovered, leading to applying Pauli \(X\) and pushing through a Pauli \(Y.\)

• If it is a Pauli \(Z,\) a Pauli \(Z\) is applied to the first physical qubit and we are done. No operator is pushed through on the bond axis.

We can recast this condition into the following dynamic circuit instructions:

• If the first/second bit of the measurement outcome bitstring is active, apply a Pauli \(Z/Y\) to the first physical qubit.

• If the second bit is active, apply a Pauli \(X\) to all other physical qubits.

This pushing and correction step is captured by the following function.

def push_and_correct(op_id, phys_wires):
    """Push an operator from left to right along the bond sites of an MPS and
    conditionally apply correcting operations on the corresponding physical sites.
    - The operator is given by two bits indicating the Pauli operator type.
    - The physical MPS sites are given by phys_wires
    """
    w = phys_wires[0]

    # Apply Z if input is Z or Y
    qml.cond(op_id[0], qml.Z)(w)
    # Apply Y if input is X or Y
    qml.cond(op_id[1], qml.Y)(w)
    # Apply X on other physical sites if input is X or Y
    for i in phys_wires[1:]:
        qml.cond(op_id[1], qml.X)(i)
    # Push through Y if input is X or Y
    return np.array([op_id[1], op_id[1]])