Coverage for qml_essentials/coefficients.py: 99%
263 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-09-08 14:29 +0000
1from __future__ import annotations
2import numpy as np
3import math
4from collections import defaultdict
5from dataclasses import dataclass
6import pennylane as qml
7from pennylane.operation import Operator
8import pennylane.ops.op_math as qml_op
9from typing import List, Tuple, Optional, Any, Dict, Union
11from qml_essentials.model import Model
14class Coefficients:
15 @staticmethod
16 def get_spectrum(
17 model: Model,
18 mfs: int = 1,
19 mts: int = 1,
20 shift=False,
21 trim=False,
22 **kwargs,
23 ) -> np.ndarray:
24 """
25 Extracts the coefficients of a given model using a FFT (np-fft).
27 Note that the coefficients are complex numbers, but the imaginary part
28 of the coefficients should be very close to zero, since the expectation
29 values of the Pauli operators are real numbers.
31 It can perform oversampling in both the frequency and time domain
32 using the `mfs` and `mts` arguments.
34 Args:
35 model (Model): The model to sample.
36 mfs (int): Multiplicator for the highest frequency. Default is 2.
37 mts (int): Multiplicator for the number of time samples. Default is 1.
38 shift (bool): Whether to apply np-fftshift. Default is False.
39 trim (bool): Whether to remove the Nyquist frequency if spectrum is even.
40 Default is False.
41 kwargs (Any): Additional keyword arguments for the model function.
43 Returns:
44 np.ndarray: The sampled Fourier coefficients.
45 """
46 kwargs.setdefault("force_mean", True)
47 kwargs.setdefault("execution_type", "expval")
49 coeffs, freqs = Coefficients._fourier_transform(
50 model, mfs=mfs, mts=mts, **kwargs
51 )
53 if not np.isclose(np.sum(coeffs).imag, 0.0, rtol=1.0e-5):
54 raise ValueError(
55 f"Spectrum is not real. Imaginary part of coefficients is:\
56 {np.sum(coeffs).imag}"
57 )
59 if trim:
60 for ax in range(model.n_input_feat):
61 if coeffs.shape[ax] % 2 == 0:
62 coeffs = np.delete(coeffs, len(coeffs) // 2, axis=ax)
63 freqs = np.delete(freqs, len(freqs) // 2, axis=ax)
65 if shift:
66 return np.fft.fftshift(
67 coeffs, axes=list(range(model.n_input_feat))
68 ), np.fft.fftshift(freqs)
69 else:
70 return coeffs, freqs
72 @staticmethod
73 def _fourier_transform(
74 model: Model, mfs: int, mts: int, **kwargs: Any
75 ) -> np.ndarray:
76 # Create a frequency vector with as many frequencies as model degrees,
77 # oversampled by nfs
78 n_freqs: int = 2 * mfs * model.degree + 1
80 start, stop, step = 0, 2 * mts * np.pi, 2 * np.pi / n_freqs
81 # Stretch according to the number of frequencies
82 inputs: np.ndarray = np.arange(start, stop, step)
84 # permute with input dimensionality
85 nd_inputs = np.array(np.meshgrid(*[inputs] * model.n_input_feat)).T.reshape(
86 -1, model.n_input_feat
87 )
89 # Output vector is not necessarily the same length as input
90 outputs = model(inputs=nd_inputs, **kwargs)
91 outputs = outputs.reshape(*(inputs.shape * model.n_input_feat), -1).squeeze()
93 coeffs = np.fft.fftn(outputs, axes=list(range(model.n_input_feat)))
95 # assert (
96 # mts * n_freqs,
97 # ) * model.n_input_feat == coeffs.shape, f"Expected shape\
98 # {(mts * n_freqs,) * model.n_input_feat} but got {coeffs.shape}"
100 freqs = np.fft.fftfreq(mts * n_freqs, 1 / n_freqs)
102 # TODO: this could cause issues with multidim input
103 # FIXME: account for different frequencies in multidim input scenarios
104 # Run the fft and rearrange +
105 # normalize the output (using product if multidim)
106 return (
107 coeffs / np.prod(outputs.shape[0 : model.n_input_feat]),
108 freqs,
109 # np.repeat(freqs[:, np.newaxis], model.n_input_feat, axis=1).squeeze(),
110 )
112 @staticmethod
113 def get_psd(coeffs: np.ndarray) -> np.ndarray:
114 """
115 Calculates the power spectral density (PSD) from given Fourier coefficients.
117 Args:
118 coeffs (np.ndarray): The Fourier coefficients.
120 Returns:
121 np.ndarray: The power spectral density.
122 """
123 # TODO: if we apply trim=True in advance, this will be slightly wrong..
125 def abs2(x):
126 return x.real**2 + x.imag**2
128 scale = 2.0 / (len(coeffs) ** 2)
129 return scale * abs2(coeffs)
131 @staticmethod
132 def evaluate_Fourier_series(
133 coefficients: np.ndarray,
134 frequencies: np.ndarray,
135 inputs: Union[np.ndarray, list, float],
136 ) -> float:
137 """
138 Evaluate the function value of a Fourier series at one point.
140 Args:
141 coefficients (np.ndarray): Coefficients of the Fourier series.
142 frequencies (np.ndarray): Corresponding frequencies.
143 inputs (np.ndarray): Point at which to evaluate the function.
144 Returns:
145 float: The function value at the input point.
146 """
147 dims = len(coefficients.shape)
149 if not isinstance(inputs, (np.ndarray, list)):
150 inputs = [inputs]
152 frequencies = np.stack(np.meshgrid(*[frequencies] * dims)).T.reshape(-1, dims)
153 freq_inputs = np.einsum("...j,j->...", frequencies, inputs)
154 coeffs = coefficients.flatten()
155 freq_inputs = freq_inputs.flatten()
157 exp = 0.0
158 for omega_x, c in zip(freq_inputs, coeffs):
159 exp += c * np.exp(1j * omega_x)
161 return np.real_if_close(exp)
164class FourierTree:
165 """
166 Sine-cosine tree representation for the algorithm by Nemkov et al.
167 This tree can be used to obtain analytical Fourier coefficients for a given
168 Pauli-Clifford circuit.
169 """
171 class CoefficientsTreeNode:
172 """
173 Representation of a node in the coefficients tree for the algorithm by
174 Nemkov et al.
175 """
177 def __init__(
178 self,
179 parameter_idx: Optional[int],
180 observable: FourierTree.PauliOperator,
181 is_sine_factor: bool,
182 is_cosine_factor: bool,
183 left: Optional[FourierTree.CoefficientsTreeNode] = None,
184 right: Optional[FourierTree.CoefficientsTreeNode] = None,
185 ):
186 """
187 Coefficient tree node initialisation. Each node has information about
188 its creation context and it's children, i.e.:
190 Args:
191 parameter_idx (Optional[int]): Index of the corresp. param. index i.
192 observable (FourierTree.PauliOperator): The nodes observable to
193 obtain the expectation value that contributes to the constant
194 term.
195 is_sine_factor (bool): If this node belongs to a sine coefficient.
196 is_cosine_factor (bool): If this node belongs to a cosine coefficient.
197 left (Optional[CoefficientsTreeNode]): left child (if any).
198 right (Optional[CoefficientsTreeNode]): right child (if any).
199 """
200 self.parameter_idx = parameter_idx
202 assert not (
203 is_sine_factor and is_cosine_factor
204 ), "Cannot be sine and cosine at the same time"
205 self.is_sine_factor = is_sine_factor
206 self.is_cosine_factor = is_cosine_factor
208 # If the observable does not constist of only Z and I, the
209 # expectation (and therefore the constant node term) is zero
210 if np.logical_or(
211 observable.list_repr == 0, observable.list_repr == 1
212 ).any():
213 self.term = 0.0
214 else:
215 self.term = observable.phase
217 self.left = left
218 self.right = right
220 def evaluate(self, parameters: list[float]) -> float:
221 """
222 Recursive function to evaluate the expectation of the coefficient tree,
223 starting from the current node.
225 Args:
226 parameters (list[float]): The parameters, by which the circuit (and
227 therefore the tree) is parametrised.
229 Returns:
230 float: The expectation for the current node and it's children.
231 """
232 factor = (
233 parameters[self.parameter_idx]
234 if self.parameter_idx is not None
235 else 1.0
236 )
237 if self.is_sine_factor:
238 factor = 1j * np.sin(factor)
239 elif self.is_cosine_factor:
240 factor = np.cos(factor)
241 if not (self.left or self.right): # leaf
242 return factor * self.term
244 sum_children = 0.0
245 if self.left:
246 left = self.left.evaluate(parameters)
247 sum_children = sum_children + left
248 if self.right:
249 right = self.right.evaluate(parameters)
250 sum_children = sum_children + right
252 return factor * sum_children
254 def get_leafs(
255 self,
256 sin_list: List[int],
257 cos_list: List[int],
258 existing_leafs: List[FourierTree.TreeLeaf] = [],
259 ) -> List[FourierTree.TreeLeaf]:
260 """
261 Traverse the tree starting from the current node, to obtain the tree
262 leafs only.
263 The leafs correspond to the terms in the sine-cosine tree
264 representation that eventually are used to obtain coefficients and
265 frequencies.
266 Sine and cosine lists are recursively passed to the children until a
267 leaf is reached (top to bottom).
268 Leafs are then passed bottom to top to the caller.
270 Args:
271 sin_list (List[int]): Current number of sine contributions for each
272 parameter. Has the same length as the parameters, as each
273 position corresponds to one parameter.
274 cos_list (List[int]): Current number of cosine contributions for
275 each parameter. Has the same length as the parameters, as each
276 position corresponds to one parameter.
277 existing_leafs (List[TreeLeaf]): Current list of leaf nodes from
278 parents.
280 Returns:
281 List[TreeLeaf]: Updated list of leaf nodes.
282 """
284 if self.is_sine_factor:
285 sin_list[self.parameter_idx] += 1
286 if self.is_cosine_factor:
287 cos_list[self.parameter_idx] += 1
289 if not (self.left or self.right): # leaf
290 if self.term != 0.0:
291 return [FourierTree.TreeLeaf(sin_list, cos_list, self.term)]
292 else:
293 return []
295 if self.left:
296 leafs_left = self.left.get_leafs(
297 sin_list.copy(), cos_list.copy(), existing_leafs.copy()
298 )
299 else:
300 leafs_left = []
302 if self.right:
303 leafs_right = self.right.get_leafs(
304 sin_list.copy(), cos_list.copy(), existing_leafs.copy()
305 )
306 else:
307 leafs_right = []
309 existing_leafs.extend(leafs_left)
310 existing_leafs.extend(leafs_right)
311 return existing_leafs
313 @dataclass
314 class TreeLeaf:
315 """
316 Coefficient tree leafs according to the algorithm by Nemkov et al., which
317 correspond to the terms in the sine-cosine tree representation that
318 eventually are used to obtain coefficients and frequencies.
320 Args:
321 sin_indices (List[int]): Current number of sine contributions for each
322 parameter. Has the same length as the parameters, as each
323 position corresponds to one parameter.
324 cos_list (List[int]): Current number of cosine contributions for
325 each parameter. Has the same length as the parameters, as each
326 position corresponds to one parameter.
327 term (np.complex): Constant factor of the leaf, depending on the
328 expectation value of the observable, and a phase.
329 """
331 sin_indices: List[int]
332 cos_indices: List[int]
333 term: np.complex128
335 class PauliOperator:
336 """
337 Utility class for storing Pauli Rotations, the corresponding indices
338 in the XY-Space (whether there is a gate with X or Y generator at a
339 certain qubit) and the phase.
341 Args:
342 pauli (Union[Operator, np.ndarray[int]]): Pauli Rotation Operation
343 or list representation
344 n_qubits (int): Number of qubits in the circuit
345 prev_xy_indices (Optional[np.ndarray[bool]]): X/Y indices of the
346 previous Pauli sequence. Defaults to None.
347 is_observable (bool): If the operator is an observable. Defaults to
348 False.
349 is_init (bool): If this Pauli operator is initialised the first
350 time. Defaults to True.
351 phase (float): Phase of the operator. Defaults to 1.0
352 """
354 def __init__(
355 self,
356 pauli: Union[Operator, np.ndarray[int]],
357 n_qubits: int,
358 prev_xy_indices: Optional[np.ndarray[bool]] = None,
359 is_observable: bool = False,
360 is_init: bool = True,
361 phase: float = 1.0,
362 ):
363 self.is_observable = is_observable
364 self.phase = phase
366 if is_init:
367 if not is_observable:
368 pauli = pauli.generator()[0].base
369 self.list_repr = self._create_list_representation(pauli, n_qubits)
370 else:
371 assert isinstance(pauli, np.ndarray)
372 self.list_repr = pauli
374 if prev_xy_indices is None:
375 prev_xy_indices = np.zeros(n_qubits, dtype=bool)
376 self.xy_indices = np.logical_or(
377 prev_xy_indices,
378 self._compute_xy_indices(self.list_repr, rev=is_observable),
379 )
381 @staticmethod
382 def _compute_xy_indices(
383 op: np.ndarray[int], rev: bool = False
384 ) -> np.ndarray[bool]:
385 """
386 Computes the positions of X or Y gates to an one-hot encoded boolen
387 array.
389 Args:
390 op (np.ndarray[int]): Pauli-Operation list representation.
391 rev (bool): Whether to negate the array.
393 Returns:
394 np.ndarray[bool]: One hot encoded boolean array.
395 """
396 xy_indices = (op == 0) + (op == 1)
397 if rev:
398 xy_indices = ~xy_indices
399 return xy_indices
401 @staticmethod
402 def _create_list_representation(op: Operator, n_qubits: int) -> np.ndarray[int]:
403 """
404 Create list representation of a Pennylane Operator.
405 Generally, the list representation is a list of length n_qubits,
406 where at each position a Pauli Operator is encoded as such:
407 I: -1
408 X: 0
409 Y: 1
410 Z: 2
412 Args:
413 op (Operator): Pennylane Operator
414 n_qubits (int): number of qubits in the circuit
416 Returns:
417 np.ndarray[int]: List representation
418 """
419 pauli_repr = -np.ones(n_qubits, dtype=int)
420 op = op.terms()[1][0] if isinstance(op, qml_op.Prod) else op
421 op = op.base if isinstance(op, qml_op.SProd) else op
423 if isinstance(op, qml.PauliX):
424 pauli_repr[op.wires[0]] = 0
425 elif isinstance(op, qml.PauliY):
426 pauli_repr[op.wires[0]] = 1
427 elif isinstance(op, qml.PauliZ):
428 pauli_repr[op.wires[0]] = 2
429 else:
430 for o in op:
431 if isinstance(o, qml.PauliX):
432 pauli_repr[o.wires[0]] = 0
433 elif isinstance(o, qml.PauliY):
434 pauli_repr[o.wires[0]] = 1
435 elif isinstance(o, qml.PauliZ):
436 pauli_repr[o.wires[0]] = 2
437 return pauli_repr
439 def is_commuting(self, pauli: np.ndarray[int]) -> bool:
440 """
441 Computes if this Pauli commutes with another Pauli operator.
442 This computation is based on the fact that The commutator is zero
443 if and only if the number of anticommuting single-qubit Paulis is
444 even.
446 Args:
447 pauli (np.ndarray[int]): List representation of another Pauli
449 Returns:
450 bool: If the current and other Pauli are commuting.
451 """
452 anticommutator = np.where(
453 pauli < 0,
454 False,
455 np.where(
456 self.list_repr < 0,
457 False,
458 np.where(self.list_repr == pauli, False, True),
459 ),
460 )
461 return not (np.sum(anticommutator) % 2)
463 def tensor(self, pauli: np.ndarray[int]) -> FourierTree.PauliOperator:
464 """
465 Compute tensor product between the current Pauli and a given list
466 representation of another Pauli operator.
468 Args:
469 pauli (np.ndarray[int]): List representation of Pauli
471 Returns:
472 FourierTree.PauliOperator: New Pauli operator object, which
473 contains the tensor product
474 """
475 diff = (pauli - self.list_repr + 3) % 3
476 phase = np.where(
477 self.list_repr < 0,
478 1.0,
479 np.where(
480 pauli < 0,
481 1.0,
482 np.where(
483 diff == 2,
484 1.0j,
485 np.where(diff == 1, -1.0j, 1.0),
486 ),
487 ),
488 )
490 obs = np.where(
491 self.list_repr < 0,
492 pauli,
493 np.where(
494 pauli < 0,
495 self.list_repr,
496 np.where(
497 diff == 2,
498 (self.list_repr + 1) % 3,
499 np.where(diff == 1, (self.list_repr + 2) % 3, -1),
500 ),
501 ),
502 )
503 phase = self.phase * np.prod(phase)
504 return FourierTree.PauliOperator(
505 obs, phase=phase, n_qubits=obs.size, is_init=False, is_observable=True
506 )
508 def __init__(self, model: Model, inputs=1.0):
509 """
510 Tree initialisation, based on the Pauli-Clifford representation of a model.
511 Currently, only one input feature is supported.
513 **Usage**:
514 ```
515 # initialise a model
516 model = Model(...)
518 # initialise and build FourierTree
519 tree = FourierTree(model)
521 # get expectaion value
522 exp = tree()
524 # Get spectrum (for each observable, we have one list element)
525 coeff_list, freq_list = tree.spectrum()
526 ```
528 Args:
529 model (Model): The Model, for which to build the tree
530 inputs (bool, optional): Possible default inputs. Defaults to 1.0.
531 """
532 self.model = model
533 self.tree_roots = None
535 if not model.as_pauli_circuit:
536 model.as_pauli_circuit = True
538 inputs = (
539 self.model._inputs_validation(inputs)
540 if inputs is not None
541 else self.model._inputs_validation(1.0)
542 )
543 inputs.requires_grad = False
545 quantum_tape = qml.workflow.construct_tape(self.model.circuit)(
546 params=model.params, inputs=inputs
547 )
548 self.parameters = [np.squeeze(p) for p in quantum_tape.get_parameters()]
549 self.input_indices = [
550 i for (i, p) in enumerate(self.parameters) if not p.requires_grad
551 ]
553 self.observables = self._encode_observables(quantum_tape.observables)
555 pauli_rot = FourierTree.PauliOperator(
556 quantum_tape.operations[0],
557 self.model.n_qubits,
558 )
559 self.pauli_rotations = [pauli_rot]
560 for op in quantum_tape.operations[1:]:
561 pauli_rot = FourierTree.PauliOperator(
562 op, self.model.n_qubits, pauli_rot.xy_indices
563 )
564 self.pauli_rotations.append(pauli_rot)
566 self.tree_roots = self.build()
567 self.leafs: List[List[FourierTree.TreeLeaf]] = self._get_tree_leafs()
569 def __call__(
570 self,
571 params: Optional[np.ndarray] = None,
572 inputs: Optional[np.ndarray] = None,
573 **kwargs,
574 ) -> np.ndarray:
575 """
576 Evaluates the Fourier tree via sine-cosine terms sum. This is
577 equivalent to computing the expectation value of the observables with
578 respect to the corresponding circuit.
580 Args:
581 params (Optional[np.ndarray], optional): Parameters of the model.
582 Defaults to None.
583 inputs (Optional[np.ndarray], optional): Inputs to the circuit.
584 Defaults to None.
586 Returns:
587 np.ndarray: Expectation value of the tree.
589 Raises:
590 NotImplementedError: When using other "execution_type" as expval.
591 NotImplementedError: When using "noise_params"
594 """
595 params = (
596 self.model._params_validation(params)
597 if params is not None
598 else self.model.params
599 )
600 inputs = (
601 self.model._inputs_validation(inputs)
602 if inputs is not None
603 else self.model._inputs_validation(1.0)
604 )
605 inputs.requires_grad = False
607 if kwargs.get("execution_type", "expval") != "expval":
608 raise NotImplementedError(
609 f'Currently, only "expval" execution type is supported when '
610 f"building FourierTree. Got {kwargs.get('execution_type', 'expval')}."
611 )
612 if kwargs.get("noise_params", None) is not None:
613 raise NotImplementedError(
614 "Currently, noise is not supported when building FourierTree."
615 )
617 quantum_tape = qml.workflow.construct_tape(self.model.circuit)(
618 params=self.model.params, inputs=inputs
619 )
620 self.parameters = [np.squeeze(p) for p in quantum_tape.get_parameters()]
621 self.input_indices = [
622 i for (i, p) in enumerate(self.parameters) if not p.requires_grad
623 ]
625 results = np.zeros(len(self.tree_roots))
626 for i, root in enumerate(self.tree_roots):
627 results[i] = np.real_if_close(root.evaluate(self.parameters))
629 if kwargs.get("force_mean", False):
630 return np.mean(results)
631 else:
632 return results
634 def build(self) -> List[CoefficientsTreeNode]:
635 """
636 Creates the coefficient tree, i.e. it creates and initialises the tree
637 nodes.
638 Leafs can be obtained separately in _get_tree_leafs, once the tree is
639 set up.
641 Returns:
642 List[CoefficientsTreeNode]: The list of root nodes (one root for
643 each observable).
644 """
645 tree_roots = []
646 pauli_rotation_idx = len(self.pauli_rotations) - 1
647 for obs in self.observables:
648 root = self._create_tree_node(obs, pauli_rotation_idx)
649 tree_roots.append(root)
650 return tree_roots
652 def _encode_observables(
653 self, tape_obs: List[Operator]
654 ) -> List[FourierTree.PauliOperator]:
655 """
656 Encodes Pennylane observables from tape as FourierTree.PauliOperator
657 utility objects.
659 Args:
660 tape_obs (List[Operator]): Pennylane tape operations
662 Returns:
663 List[FourierTree.PauliOperator]: List of Pauli operators
664 """
665 observables = []
666 for obs in tape_obs:
667 observables.append(
668 FourierTree.PauliOperator(obs, self.model.n_qubits, is_observable=True)
669 )
670 return observables
672 def _get_tree_leafs(self) -> List[List[TreeLeaf]]:
673 """
674 Obtain all Leaf Nodes with its sine- and cosine terms.
676 Returns:
677 List[List[TreeLeaf]]: For each observable (root), the list of leaf
678 nodes.
679 """
680 leafs = []
681 for root in self.tree_roots:
682 sin_list = np.zeros(len(self.parameters), dtype=np.int32)
683 cos_list = np.zeros(len(self.parameters), dtype=np.int32)
684 leafs.append(root.get_leafs(sin_list, cos_list, []))
685 return leafs
687 def get_spectrum(
688 self, force_mean: bool = False
689 ) -> Tuple[List[np.ndarray], List[np.ndarray]]:
690 """
691 Computes the Fourier spectrum for the tree, consisting of the
692 frequencies and its corresponding coefficinets.
693 If the frag force_mean was set in the constructor, the mean coefficient
694 over all observables (roots) are computed.
696 Args:
697 force_mean (bool, optional): Whether to average over multiple
698 observables. Defaults to False.
700 Returns:
701 Tuple[List[np.ndarray], List[np.ndarray]]:
702 - List of frequencies, one list for each observable (root).
703 - List of corresponding coefficents, one list for each
704 observable (root).
705 """
706 parameter_indices = [
707 i for i in range(len(self.parameters)) if i not in self.input_indices
708 ]
710 coeffs = []
711 for leafs in self.leafs:
712 freq_terms = defaultdict(np.complex128)
713 for leaf in leafs:
714 leaf_factor, s, c = self._compute_leaf_factors(leaf, parameter_indices)
716 for a in range(s + 1):
717 for b in range(c + 1):
718 comb = math.comb(s, a) * math.comb(c, b) * (-1) ** (s - a)
719 freq_terms[2 * a + 2 * b - s - c] += comb * leaf_factor
721 coeffs.append(freq_terms)
723 frequencies, coefficients = self._freq_terms_to_coeffs(coeffs, force_mean)
724 return coefficients, frequencies
726 def _freq_terms_to_coeffs(
727 self, coeffs: List[Dict[int, np.ndarray]], force_mean: bool
728 ) -> Tuple[List[np.ndarray], List[np.ndarray]]:
729 """
730 Given a list of dictionaries of the form:
731 [
732 {
733 freq_obs1_1: coeff1,
734 freq_obs1_2: coeff2,
735 ...
736 },
737 {
738 freq_obs2_1: coeff3,
739 freq_obs2_2: coeff4,
740 ...
741 }
742 ...
743 ],
744 Compute two separate lists of frequencies and coefficients.
745 such that:
746 freqs: [
747 [freq_obs1_1, freq_obs1_1, ...],
748 [freq_obs2_1, freq_obs2_1, ...],
749 ...
750 ]
751 coeffs: [
752 [coeff1, coeff2, ...],
753 [coeff3, coeff4, ...],
754 ...
755 ]
757 If force_mean is set length of the resulting frequency and coefficent
758 list is 1.
760 Args:
761 coeffs (List[Dict[int, np.ndarray]]): Frequency->Coefficients
762 dictionary list, one dict for each observable (root).
763 force_mean (bool): Whether to average coefficients over multiple
764 observables.
766 Returns:
767 Tuple[List[np.ndarray], List[np.ndarray]]:
768 - List of frequencies, one list for each observable (root).
769 - List of corresponding coefficents, one list for each
770 observable (root).
771 """
772 frequencies = []
773 coefficients = []
774 if force_mean:
775 all_freqs = sorted(set([f for c in coeffs for f in c.keys()]))
776 coefficients.append(
777 np.array([np.mean([c.get(f, 0.0) for c in coeffs]) for f in all_freqs])
778 )
779 frequencies.append(np.array(all_freqs))
780 else:
781 for freq_terms in coeffs:
782 freq_terms = dict(sorted(freq_terms.items()))
783 frequencies.append(np.array(list(freq_terms.keys())))
784 coefficients.append(np.array(list(freq_terms.values())))
785 return frequencies, coefficients
787 def _compute_leaf_factors(
788 self, leaf: TreeLeaf, parameter_indices: List[int]
789 ) -> Tuple[float, int, int]:
790 """
791 Computes the constant coefficient factor for each leaf.
792 Additionally sine and cosine contributions of the input parameters for
793 this leaf are returned, which are required to obtain the corresponding
794 frequencies.
796 Args:
797 leaf (TreeLeaf): The leaf for which to compute the factor.
798 parameter_indices (List[int]): Variational parameter indices.
800 Returns:
801 Tuple[float, int, int]:
802 - float: the constant factor for the leaf
803 - int: number of sine contributions of the input
804 - int: number of cosine contributions of the input
805 """
806 leaf_factor = 1.0
807 for i in parameter_indices:
808 interm_factor = (
809 np.cos(self.parameters[i]) ** leaf.cos_indices[i]
810 * (1j * np.sin(self.parameters[i])) ** leaf.sin_indices[i]
811 )
812 leaf_factor = leaf_factor * interm_factor
814 # Get number of sine and cosine factors to which the input contributes
815 c = np.sum([leaf.cos_indices[k] for k in self.input_indices], dtype=np.int32)
816 s = np.sum([leaf.sin_indices[k] for k in self.input_indices], dtype=np.int32)
818 leaf_factor = leaf.term * leaf_factor * 0.5 ** (s + c)
820 return leaf_factor, s, c
822 def _early_stopping_possible(
823 self, pauli_rotation_idx: int, observable: FourierTree.PauliOperator
824 ):
825 """
826 Checks if a node for an observable can be discarded as all expecation
827 values that can result through further branching are zero.
828 The method is mentioned in the paper by Nemkov et al.: If the one-hot
829 encoded indices for X/Y operations in the Pauli-rotation generators are
830 a basis for that of the observable, the node must be processed further.
831 If not, it can be discarded.
833 Args:
834 pauli_rotation_idx (int): Index of remaining Pauli rotation gates.
835 Gates itself are attributes of the class.
836 observable (FourierTree.PauliOperator): Current observable
837 """
838 xy_indices_obs = np.logical_or(
839 observable.xy_indices, self.pauli_rotations[pauli_rotation_idx].xy_indices
840 ).all()
842 return not xy_indices_obs
844 def _create_tree_node(
845 self,
846 observable: FourierTree.PauliOperator,
847 pauli_rotation_idx: int,
848 parameter_idx: Optional[int] = None,
849 is_sine: bool = False,
850 is_cosine: bool = False,
851 ) -> Optional[CoefficientsTreeNode]:
852 """
853 Builds the Fourier-Tree according to the algorithm by Nemkov et al.
855 Args:
856 observable (FourierTree.PauliOperator): Current observable
857 pauli_rotation_idx (int): Index of remaining Pauli rotation gates.
858 Gates itself are attributes of the class.
859 parameter_idx (Optional[int]): Index of the current parameter.
860 Parameters itself are attributes of the class.
861 is_sine (bool): If the current node is a sine (left) node.
862 is_cosine (bool): If the current node is a cosine (right) node.
864 Returns:
865 Optional[CoefficientsTreeNode]: The resulting node. Children are set
866 recursively. The top level receives the tree root.
867 """
868 if self._early_stopping_possible(pauli_rotation_idx, observable):
869 return None
871 # remove commuting paulis
872 while pauli_rotation_idx >= 0:
873 last_pauli = self.pauli_rotations[pauli_rotation_idx]
874 if not observable.is_commuting(last_pauli.list_repr):
875 break
876 pauli_rotation_idx -= 1
877 else: # leaf
878 return FourierTree.CoefficientsTreeNode(
879 parameter_idx, observable, is_sine, is_cosine
880 )
882 last_pauli = self.pauli_rotations[pauli_rotation_idx]
884 left = self._create_tree_node(
885 observable,
886 pauli_rotation_idx - 1,
887 pauli_rotation_idx,
888 is_cosine=True,
889 )
891 next_observable = self._create_new_observable(last_pauli.list_repr, observable)
892 right = self._create_tree_node(
893 next_observable,
894 pauli_rotation_idx - 1,
895 pauli_rotation_idx,
896 is_sine=True,
897 )
899 return FourierTree.CoefficientsTreeNode(
900 parameter_idx,
901 observable,
902 is_sine,
903 is_cosine,
904 left,
905 right,
906 )
908 def _create_new_observable(
909 self, pauli: np.ndarray[int], observable: FourierTree.PauliOperator
910 ) -> FourierTree.PauliOperator:
911 """
912 Utility function to obtain the new observable for a tree node, if the
913 last Pauli and the observable do not commute.
915 Args:
916 pauli (np.ndarray[int]): The int array representation of the last
917 Pauli rotation in the operation sequence.
918 observable (FourierTree.PauliOperator): The current observable.
920 Returns:
921 FourierTree.PauliOperator: The new observable.
922 """
923 observable = observable.tensor(pauli)
924 return observable