Too many subscripts in einsum

Hi,

I have a technical question. We want to compute the purity of a quantum state after adding depolarization noise. However, we get the error:

`ValueError: too many subscripts in einsum`

Is there a way to handle the subscripts or circumvent the problem?

To reproduce the error, I provide a minimum working example as well as the full error message below:

``````import pennylane as qml
from pennylane import numpy as npq
import numpy as np

# number of qubits
n_qubits = 4**2

# coefficients of cost hamiltonian
coeffs = [211.5, -24.25, -23.5, -23.5, -24.25, -24.5, -28.0, -28.0, -24.5, -23.25, -25.5, -25.5, -23.25, -28.0, -26.0, -26.0, -28.0, 5.0, 5.0, 5.0, 5.0, 0.25, 5.0, 1.25, 5.0, 0.25, 5.0, 5.0, 0.25, 5.0, 0.25, 1.25, 5.0, 1.25]

# Pauli operators of cost hamiltonian
obs = [qml.Identity(wires=[0]), qml.PauliZ(wires=[0]), qml.PauliZ(wires=[1]), qml.PauliZ(wires=[2]), qml.PauliZ(wires=[3]), qml.PauliZ(wires=[4]), qml.PauliZ(wires=[5]), qml.PauliZ(wires=[6]), qml.PauliZ(wires=[7]), qml.PauliZ(wires=[8]), qml.PauliZ(wires=[9]), qml.PauliZ(wires=[10]), qml.PauliZ(wires=[11]), qml.PauliZ(wires=[12]), qml.PauliZ(wires=[13]), qml.PauliZ(wires=[14]), qml.PauliZ(wires=[15]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[1]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[2]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[3]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[4]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[5]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[8]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[9]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[12]), qml.PauliZ(wires=[0]) @ qml.PauliZ(wires=[13]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[2]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[3]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[4]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[5]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[6]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[8]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[9]), qml.PauliZ(wires=[1]) @ qml.PauliZ(wires=[10])]

# cost hamiltonian
cost_h = qml.Hamiltonian(coeffs, obs)

single_qubit_coeffs = {0: -15.0, 1: -20, 2: -20, 3: -15.0, 4: -19.0, 5: -20, 6: -20, 7: -19.0, 8: -19.0, 9: -20, 10: -20, 11: -19.0, 12: -10, 13: -20, 14: -20, 15: -10}

two_qubit_coeffs = {(0, 1): 20, (0, 2): 20, (0, 3): 20, (0, 4): 20, (0, 5): 1.0, (0, 6): 0, (0, 7): 0, (0, 8): 20, (0, 9): 5.0, (0, 10): 0, (0, 11): 0, (0, 12): 20, (0, 13): 1.0, (0, 14): 0, (0, 15): 0, (1, 2): 20, (1, 3): 20, (1, 4): 1.0, (1, 5): 20, (1, 6): 1.0, (1, 7): 0, (1, 8): 5.0, (1, 9): 20, (1, 10): 5.0, (1, 11): 0, (1, 12): 1.0, (1, 13): 20, (1, 14): 1.0, (1, 15): 0, (2, 3): 20, (2, 4): 0, (2, 5): 1.0, (2, 6): 20, (2, 7): 1.0, (2, 8): 0, (2, 9): 5.0, (2, 10): 20, (2, 11): 5.0, (2, 12): 0, (2, 13): 1.0, (2, 14): 20, (2, 15): 1.0, (3, 4): 0, (3, 5): 0, (3, 6): 1.0, (3, 7): 20, (3, 8): 0, (3, 9): 0, (3, 10): 5.0, (3, 11): 20, (3, 12): 0, (3, 13): 0, (3, 14): 1.0, (3, 15): 20, (4, 5): 20, (4, 6): 20, (4, 7): 20, (4, 8): 20, (4, 9): 5.0, (4, 10): 0, (4, 11): 0, (4, 12): 20, (4, 13): 10.0, (4, 14): 0, (4, 15): 0, (5, 6): 20, (5, 7): 20, (5, 8): 5.0, (5, 9): 20, (5, 10): 5.0, (5, 11): 0, (5, 12): 10.0, (5, 13): 20, (5, 14): 10.0, (5, 15): 0, (6, 7): 20, (6, 8): 0, (6, 9): 5.0, (6, 10): 20, (6, 11): 5, (6, 12): 0, (6, 13): 10, (6, 14): 20, (6, 15): 10, (7, 8): 0, (7, 9): 0, (7, 10): 5, (7, 11): 20, (7, 12): 0, (7, 13): 0, (7, 14): 10, (7, 15): 20, (8, 9): 20, (8, 10): 20, (8, 11): 20, (8, 12): 20, (8, 13): 1, (8, 14): 0, (8, 15): 0, (9, 10): 20, (9, 11): 20, (9, 12): 1, (9, 13): 20, (9, 14): 1, (9, 15): 0, (10, 11): 20, (10, 12): 0, (10, 13): 1, (10, 14): 20, (10, 15): 1, (11, 12): 0, (11, 13): 0, (11, 14): 1.0, (11, 15): 20, (12, 13): 20, (12, 14): 20, (12, 15): 20, (13, 14): 20, (13, 15): 20, (14, 15): 20}

# unitary operator u_drive with parameter beta
def u_drive(
beta: float,
n_qubits: int
):
for qubit in range(n_qubits):
qml.RX(2*beta, wires=qubit)

# unitary operator u_problem with parameter gamma
def u_problem(
gamma: float,
n_qubits: float,
single_qubit_coeffs: dict,
two_qubit_coeffs: dict,
p: float
):
for qubit in range(n_qubits):
sq_coeff = single_qubit_coeffs[qubit]
qml.RZ(2*gamma*sq_coeff/2, wires=qubit)

for qubit1 in range(n_qubits):
for qubit2 in range(qubit1 + 1, n_qubits):
tq_coeff = two_qubit_coeffs[(qubit1, qubit2)]
qml.RZ(2*gamma*tq_coeff/4, wires=qubit1)
qml.RZ(2*gamma*tq_coeff/4, wires=qubit2)

if p == 0:
qml.CNOT(wires=[qubit1, qubit2])
qml.RZ(2*gamma*tq_coeff/4, wires=qubit2)
qml.CNOT(wires=[qubit1, qubit2])

elif p != 0:
qml.CNOT(wires=[qubit1, qubit2])
qml.DepolarizingChannel(p, wires=qubit1)
qml.DepolarizingChannel(p, wires=qubit2)
qml.RZ(2*gamma*tq_coeff/4, wires=qubit2)
qml.CNOT(wires=[qubit1, qubit2])
qml.DepolarizingChannel(p, wires=qubit1)
qml.DepolarizingChannel(p, wires=qubit2)

# device initialization
dev_state = qml.device("default.mixed", wires=n_qubits)

# circuit
@qml.qnode(dev_state)
def circuit_state(gammas, betas, depth=1, **kwargs):
n_qubits = dev_state.num_wires
single_qubit_coeffs = kwargs["single_qubit_coeffs"]
two_qubit_coeffs = kwargs["two_qubit_coeffs"]
p = kwargs["p"]
for wire in range(n_qubits):
for i in range(depth):
u_problem(gammas[i], n_qubits, single_qubit_coeffs, two_qubit_coeffs, p)
u_drive(betas[i], n_qubits)
return qml.probs(wires=range(n_qubits))

# calculate final purity
def purity_func(
p: float,
depth: int,
single_qubit_coeffs: dict,
two_qubit_coeffs: dict
) -> float:
output = circuit_state(params_init[0], params_init[1], depth=depth,
single_qubit_coeffs=single_qubit_coeffs, two_qubit_coeffs=two_qubit_coeffs, p=p)
purity = float(np.real(np.trace(output.dot(output))))
return purity

depth = 1 # always 1 in this example
p = 0.03 # probability of depolarizing error

purity = purity_func(p, depth, single_qubit_coeffs, two_qubit_coeffs)

print("Purity of output state = {}".format(purity))
print("Renormalized purity of output state = {}".format(purity**(1/n_qubits)))
``````

Complete error message:

``````ValueError                                Traceback (most recent call last)
Cell In [8], line 94
91 depth = 1 # always 1 in this example
92 p = 0.03 # probability of depolarizing error
---> 94 purity = purity_func(p, depth, single_qubit_coeffs, two_qubit_coeffs)
96 print("Purity of output state = {}".format(purity))
97 print("Renormalized purity of output state = {}".format(purity**(1/n_qubits)))

Cell In [8], line 86, in purity_func(p, depth, single_qubit_coeffs, two_qubit_coeffs)
79 def purity_func(
80     p: float,
81     depth: int,
82     single_qubit_coeffs: dict,
83     two_qubit_coeffs: dict
84 ) -> float:
85     params_init = npq.zeros([2, depth], requires_grad=True)
---> 86     output = circuit_state(params_init[0], params_init[1], depth=depth,
87                            single_qubit_coeffs=single_qubit_coeffs, two_qubit_coeffs=two_qubit_coeffs, p=p)
88     purity = float(np.real(np.trace(output.dot(output))))
89     return purity

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/qnode.py:665, in QNode.__call__(self, *args, **kwargs)
661     self._update_original_device()
663     return res
--> 665 res = qml.execute(
666     [self.tape],
667     device=self.device,
669     interface=self.interface,
671     override_shots=override_shots,
672     **self.execute_kwargs,
673 )
675 if old_interface == "auto":
676     self.interface = "auto"

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/interfaces/execution.py:645, in execute(tapes, device, gradient_fn, interface, mode, gradient_kwargs, cache, cachesize, max_diff, override_shots, expand_fn, max_expansion, device_batch_transform)
641     return batch_fn(res)
643 if gradient_fn == "backprop" or interface is None:
644     return batch_fn(
--> 645         qml.interfaces.cache_execute(
646             batch_execute, cache, return_tuple=False, expand_fn=expand_fn
647         )(tapes)
648     )
650 # the default execution function is batch_execute
651 execute_fn = qml.interfaces.cache_execute(batch_execute, cache, expand_fn=expand_fn)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/interfaces/execution.py:206, in cache_execute.<locals>.wrapper(tapes, **kwargs)
202         return (res, []) if return_tuple else res
204 else:
205     # execute all unique tapes that do not exist in the cache
--> 206     res = fn(execution_tapes.values(), **kwargs)
208 final_res = []
210 for i, tape in enumerate(tapes):

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/interfaces/execution.py:131, in cache_execute.<locals>.fn(tapes, **kwargs)
129 def fn(tapes: Sequence[QuantumTape], **kwargs):  # pylint: disable=function-redefined
130     tapes = [expand_fn(tape) for tape in tapes]
--> 131     return original_fn(tapes, **kwargs)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/contextlib.py:79, in ContextDecorator.__call__.<locals>.inner(*args, **kwds)
76 @wraps(func)
77 def inner(*args, **kwds):
78     with self._recreate_cm():
---> 79         return func(*args, **kwds)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/_qubit_device.py:624, in QubitDevice.batch_execute(self, circuits)
621     self.reset()
623     # TODO: Insert control on value here
--> 624     res = self.execute(circuit)
625     results.append(res)
627 if self.tracker.active:

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/devices/default_mixed.py:594, in DefaultMixed.execute(self, circuit, **kwargs)
592         wires_list.append(obs.wires)
593     self.measured_wires = qml.wires.Wires.all_wires(wires_list)
--> 594 return super().execute(circuit, **kwargs)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/_qubit_device.py:375, in QubitDevice.execute(self, circuit, **kwargs)
372 self.check_validity(circuit.operations, circuit.observables)
374 # apply all circuit operations
--> 375 self.apply(circuit.operations, rotations=circuit.diagonalizing_gates, **kwargs)
377 # generate computational basis samples
378 if self.shots is not None or circuit.is_sampled:

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/devices/default_mixed.py:610, in DefaultMixed.apply(self, operations, rotations, **kwargs)
604         raise DeviceError(
605             f"Operation {operation.name} cannot be used after other Operations have already been applied "
606             f"on a {self.short_name} device."
607         )
609 for operation in operations:
--> 610     self._apply_operation(operation)
612 # store the pre-rotated state
613 self._pre_rotated_state = self._state

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/devices/default_mixed.py:539, in DefaultMixed._apply_operation(self, operation)
537     self._apply_diagonal_unitary(matrices, wires)
538 else:
--> 539     self._apply_channel(matrices, wires)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/devices/default_mixed.py:328, in DefaultMixed._apply_channel(self, kraus, wires)
322 # index mapping for einsum, e.g., 'iga,abcdef,idh->gbchef'
323 einsum_indices = (
324     f"{kraus_index}{new_row_indices}{row_indices}, {state_indices},"
325     f"{kraus_index}{col_indices}{new_col_indices}->{new_state_indices}"
326 )
--> 328 self._state = qnp.einsum(einsum_indices, kraus, self._state, kraus_dagger)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/pennylane/math/multi_dispatch.py:543, in einsum(indices, like, *operands)
541     like = get_interface(*operands)
542 operands = np.coerce(operands, like=like)
--> 543 return np.einsum(indices, *operands, like=like)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/autoray/autoray.py:80, in do(fn, like, *args, **kwargs)
31 """Do function named ``fn`` on ``(*args, **kwargs)``, peforming single
32 dispatch to retrieve ``fn`` based on whichever library defines the class of
33 the ``args[0]``, or the ``like`` keyword argument if specified.
(...)
77     <tf.Tensor: id=91, shape=(3, 3), dtype=float32>
78 """
79 backend = choose_backend(fn, *args, like=like, **kwargs)
---> 80 return get_lib_fn(backend, fn)(*args, **kwargs)

File <__array_function__ internals>:180, in einsum(*args, **kwargs)

File /usr/local/anaconda3/envs/myenv/lib/python3.9/site-packages/numpy/core/einsumfunc.py:1371, in einsum(out, optimize, *operands, **kwargs)
1369     if specified_out:
1370         kwargs['out'] = out
-> 1371     return c_einsum(*operands, **kwargs)
1373 # Check the kwargs to avoid a more cryptic error later, without having to
1374 # repeat default values here
1375 valid_einsum_kwargs = ['dtype', 'order', 'casting']

ValueError: too many subscripts in einsum``````
1 Like

Hi @leolettuce thank you for your question!
Iâ€™m looking into it.

Hi @leolettuce !

Here are some hints on what may be happening:

1. If you donâ€™t specify the interface to use, the default interface is NumPy. Our version of NumPy comes from Autograd, and it has some restrictions. For instance, you should use np.dot(output,output) instead of output.dot(output). You can learn more about this here.
2. `np.dot(output,output)` is a single number, and getting the trace of a single number isnâ€™t possible so `purity = float(np.real(np.trace(output.dot(output))))` will throw an error. To bypass this Iâ€™ve fixed `purity = 1` just to be able to test other problems.
3. I tried this with `default.qubit` (commenting the lines which add the depolarizing channel) and I no longer get the error. However when I change `default.qubit` to `default.mixed` (keeping everything exactly the same) I get the error you get. By looking into the details of the error trace it looks like â€śbatch_executeâ€ť doesnâ€™t work well with `default.mixed` and it would seem like you are requiring the use of this function for some reason.

Next steps:
I would encourage you to try to simplify your code as much as possible so that youâ€™re not accidentally running things in batches. If you can reduce your code to a minimal example and the problem persists please post your code again here and I will take another look to see whatâ€™s happening!