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):
qml.Hadamard(wires=wire)
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:
params_init = npq.zeros([2, depth], requires_grad=True)
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,
668 gradient_fn=self.gradient_fn,
669 interface=self.interface,
670 gradient_kwargs=self.gradient_kwargs,
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