Coverage for qml_essentials/coefficients.py: 96%
389 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-10-02 13:10 +0000
1from __future__ import annotations
2import math
3from collections import defaultdict
4from dataclasses import dataclass
5import pennylane as qml
6import pennylane.numpy as np
7import numpy as nnp
8from pennylane.operation import Operator
9import pennylane.ops.op_math as qml_op
10from scipy.stats import rankdata
11from functools import reduce
12from typing import List, Tuple, Optional, Any, Dict, Union
14from qml_essentials.model import Model
16import logging
18log = logging.getLogger(__name__)
21class Coefficients:
22 @staticmethod
23 def get_spectrum(
24 model: Model,
25 mfs: int = 1,
26 mts: int = 1,
27 shift=False,
28 trim=False,
29 **kwargs,
30 ) -> Tuple[np.ndarray, np.ndarray]:
31 """
32 Extracts the coefficients of a given model using a FFT (np-fft).
34 Note that the coefficients are complex numbers, but the imaginary part
35 of the coefficients should be very close to zero, since the expectation
36 values of the Pauli operators are real numbers.
38 It can perform oversampling in both the frequency and time domain
39 using the `mfs` and `mts` arguments.
41 Args:
42 model (Model): The model to sample.
43 mfs (int): Multiplicator for the highest frequency. Default is 2.
44 mts (int): Multiplicator for the number of time samples. Default is 1.
45 shift (bool): Whether to apply np-fftshift. Default is False.
46 trim (bool): Whether to remove the Nyquist frequency if spectrum is even.
47 Default is False.
48 kwargs (Any): Additional keyword arguments for the model function.
50 Returns:
51 Tuple[np.ndarray, np.ndarray]: Tuple containing the coefficients
52 and frequencies.
53 """
54 kwargs.setdefault("force_mean", True)
55 kwargs.setdefault("execution_type", "expval")
57 coeffs, freqs = Coefficients._fourier_transform(
58 model, mfs=mfs, mts=mts, **kwargs
59 )
61 if not np.isclose(np.sum(coeffs).imag, 0.0, rtol=1.0e-5):
62 raise ValueError(
63 f"Spectrum is not real. Imaginary part of coefficients is:\
64 {np.sum(coeffs).imag}"
65 )
67 if trim:
68 for ax in range(model.n_input_feat):
69 if coeffs.shape[ax] % 2 == 0:
70 coeffs = np.delete(coeffs, len(coeffs) // 2, axis=ax)
71 freqs = np.delete(freqs, len(freqs) // 2, axis=ax)
73 if shift:
74 coeffs = np.fft.fftshift(coeffs, axes=list(range(model.n_input_feat)))
75 freqs = np.fft.fftshift(freqs)
77 return coeffs, freqs
79 @staticmethod
80 def _fourier_transform(
81 model: Model, mfs: int, mts: int, **kwargs: Any
82 ) -> np.ndarray:
83 # Create a frequency vector with as many frequencies as model degrees,
84 # oversampled by nfs
85 n_freqs: int = 2 * mfs * model.degree + 1
87 start, stop, step = 0, 2 * mts * np.pi, 2 * np.pi / n_freqs
88 # Stretch according to the number of frequencies
89 inputs: np.ndarray = np.arange(start, stop, step)
91 # permute with input dimensionality
92 nd_inputs = np.array(np.meshgrid(*[inputs] * model.n_input_feat)).T.reshape(
93 -1, model.n_input_feat
94 )
96 # Output vector is not necessarily the same length as input
97 outputs = model(inputs=nd_inputs, **kwargs)
98 outputs = outputs.reshape(*(inputs.shape * model.n_input_feat), -1).squeeze()
100 coeffs = np.fft.fftn(outputs, axes=list(range(model.n_input_feat)))
102 # TODO: in the future, this should take into account that there can be a
103 # different number of frequencies per dimension
104 freqs = [
105 np.fft.fftfreq(mts * n_freqs, 1 / n_freqs)
106 for _ in range(model.n_input_feat)
107 ]
108 # freqs = np.fft.fftfreq(mts * n_freqs, 1 / n_freqs)
110 # TODO: this could cause issues with multidim input
111 # FIXME: account for different frequencies in multidim input scenarios
112 # Run the fft and rearrange +
113 # normalize the output (using product if multidim)
114 return (
115 coeffs / np.prod(outputs.shape[0 : model.n_input_feat]),
116 np.array(freqs).squeeze(),
117 )
119 @staticmethod
120 def get_psd(coeffs: np.ndarray) -> np.ndarray:
121 """
122 Calculates the power spectral density (PSD) from given Fourier coefficients.
124 Args:
125 coeffs (np.ndarray): The Fourier coefficients.
127 Returns:
128 np.ndarray: The power spectral density.
129 """
130 # TODO: if we apply trim=True in advance, this will be slightly wrong..
132 def abs2(x):
133 return x.real**2 + x.imag**2
135 scale = 2.0 / (len(coeffs) ** 2)
136 return scale * abs2(coeffs)
138 @staticmethod
139 def evaluate_Fourier_series(
140 coefficients: np.ndarray,
141 frequencies: np.ndarray,
142 inputs: Union[np.ndarray, list, float],
143 ) -> float:
144 """
145 Evaluate the function value of a Fourier series at one point.
147 Args:
148 coefficients (np.ndarray): Coefficients of the Fourier series.
149 frequencies (np.ndarray): Corresponding frequencies.
150 inputs (np.ndarray): Point at which to evaluate the function.
151 Returns:
152 float: The function value at the input point.
153 """
154 dims = len(coefficients.shape)
156 if not isinstance(inputs, (np.ndarray, list)):
157 inputs = [inputs]
159 frequencies = np.stack(np.meshgrid(*frequencies)).T.reshape(-1, dims)
160 freq_inputs = np.einsum("...j,j->...", frequencies, inputs)
161 coeffs = coefficients.flatten()
162 freq_inputs = freq_inputs.flatten()
164 exp = 0.0
165 for omega_x, c in zip(freq_inputs, coeffs):
166 exp += c * np.exp(1j * omega_x)
168 return np.real_if_close(exp)
171class FourierTree:
172 """
173 Sine-cosine tree representation for the algorithm by Nemkov et al.
174 This tree can be used to obtain analytical Fourier coefficients for a given
175 Pauli-Clifford circuit.
176 """
178 class CoefficientsTreeNode:
179 """
180 Representation of a node in the coefficients tree for the algorithm by
181 Nemkov et al.
182 """
184 def __init__(
185 self,
186 parameter_idx: Optional[int],
187 observable: FourierTree.PauliOperator,
188 is_sine_factor: bool,
189 is_cosine_factor: bool,
190 left: Optional[FourierTree.CoefficientsTreeNode] = None,
191 right: Optional[FourierTree.CoefficientsTreeNode] = None,
192 ):
193 """
194 Coefficient tree node initialisation. Each node has information about
195 its creation context and it's children, i.e.:
197 Args:
198 parameter_idx (Optional[int]): Index of the corresp. param. index i.
199 observable (FourierTree.PauliOperator): The nodes observable to
200 obtain the expectation value that contributes to the constant
201 term.
202 is_sine_factor (bool): If this node belongs to a sine coefficient.
203 is_cosine_factor (bool): If this node belongs to a cosine coefficient.
204 left (Optional[CoefficientsTreeNode]): left child (if any).
205 right (Optional[CoefficientsTreeNode]): right child (if any).
206 """
207 self.parameter_idx = parameter_idx
209 assert not (
210 is_sine_factor and is_cosine_factor
211 ), "Cannot be sine and cosine at the same time"
212 self.is_sine_factor = is_sine_factor
213 self.is_cosine_factor = is_cosine_factor
215 # If the observable does not constist of only Z and I, the
216 # expectation (and therefore the constant node term) is zero
217 if np.logical_or(
218 observable.list_repr == 0, observable.list_repr == 1
219 ).any():
220 self.term = 0.0
221 else:
222 self.term = observable.phase
224 self.left = left
225 self.right = right
227 def evaluate(self, parameters: list[float]) -> float:
228 """
229 Recursive function to evaluate the expectation of the coefficient tree,
230 starting from the current node.
232 Args:
233 parameters (list[float]): The parameters, by which the circuit (and
234 therefore the tree) is parametrised.
236 Returns:
237 float: The expectation for the current node and it's children.
238 """
239 factor = (
240 parameters[self.parameter_idx]
241 if self.parameter_idx is not None
242 else 1.0
243 )
244 if self.is_sine_factor:
245 factor = 1j * np.sin(factor)
246 elif self.is_cosine_factor:
247 factor = np.cos(factor)
248 if not (self.left or self.right): # leaf
249 return factor * self.term
251 sum_children = 0.0
252 if self.left:
253 left = self.left.evaluate(parameters)
254 sum_children = sum_children + left
255 if self.right:
256 right = self.right.evaluate(parameters)
257 sum_children = sum_children + right
259 return factor * sum_children
261 def get_leafs(
262 self,
263 sin_list: List[int],
264 cos_list: List[int],
265 existing_leafs: List[FourierTree.TreeLeaf] = [],
266 ) -> List[FourierTree.TreeLeaf]:
267 """
268 Traverse the tree starting from the current node, to obtain the tree
269 leafs only.
270 The leafs correspond to the terms in the sine-cosine tree
271 representation that eventually are used to obtain coefficients and
272 frequencies.
273 Sine and cosine lists are recursively passed to the children until a
274 leaf is reached (top to bottom).
275 Leafs are then passed bottom to top to the caller.
277 Args:
278 sin_list (List[int]): Current number of sine contributions for each
279 parameter. Has the same length as the parameters, as each
280 position corresponds to one parameter.
281 cos_list (List[int]): Current number of cosine contributions for
282 each parameter. Has the same length as the parameters, as each
283 position corresponds to one parameter.
284 existing_leafs (List[TreeLeaf]): Current list of leaf nodes from
285 parents.
287 Returns:
288 List[TreeLeaf]: Updated list of leaf nodes.
289 """
291 if self.is_sine_factor:
292 sin_list[self.parameter_idx] += 1
293 if self.is_cosine_factor:
294 cos_list[self.parameter_idx] += 1
296 if not (self.left or self.right): # leaf
297 if self.term != 0.0:
298 return [FourierTree.TreeLeaf(sin_list, cos_list, self.term)]
299 else:
300 return []
302 if self.left:
303 leafs_left = self.left.get_leafs(
304 sin_list.copy(), cos_list.copy(), existing_leafs.copy()
305 )
306 else:
307 leafs_left = []
309 if self.right:
310 leafs_right = self.right.get_leafs(
311 sin_list.copy(), cos_list.copy(), existing_leafs.copy()
312 )
313 else:
314 leafs_right = []
316 existing_leafs.extend(leafs_left)
317 existing_leafs.extend(leafs_right)
318 return existing_leafs
320 @dataclass
321 class TreeLeaf:
322 """
323 Coefficient tree leafs according to the algorithm by Nemkov et al., which
324 correspond to the terms in the sine-cosine tree representation that
325 eventually are used to obtain coefficients and frequencies.
327 Args:
328 sin_indices (List[int]): Current number of sine contributions for each
329 parameter. Has the same length as the parameters, as each
330 position corresponds to one parameter.
331 cos_list (List[int]): Current number of cosine contributions for
332 each parameter. Has the same length as the parameters, as each
333 position corresponds to one parameter.
334 term (np.complex): Constant factor of the leaf, depending on the
335 expectation value of the observable, and a phase.
336 """
338 sin_indices: List[int]
339 cos_indices: List[int]
340 term: np.complex128
342 class PauliOperator:
343 """
344 Utility class for storing Pauli Rotations, the corresponding indices
345 in the XY-Space (whether there is a gate with X or Y generator at a
346 certain qubit) and the phase.
348 Args:
349 pauli (Union[Operator, np.ndarray[int]]): Pauli Rotation Operation
350 or list representation
351 n_qubits (int): Number of qubits in the circuit
352 prev_xy_indices (Optional[np.ndarray[bool]]): X/Y indices of the
353 previous Pauli sequence. Defaults to None.
354 is_observable (bool): If the operator is an observable. Defaults to
355 False.
356 is_init (bool): If this Pauli operator is initialised the first
357 time. Defaults to True.
358 phase (float): Phase of the operator. Defaults to 1.0
359 """
361 def __init__(
362 self,
363 pauli: Union[Operator, np.ndarray[int]],
364 n_qubits: int,
365 prev_xy_indices: Optional[np.ndarray[bool]] = None,
366 is_observable: bool = False,
367 is_init: bool = True,
368 phase: float = 1.0,
369 ):
370 self.is_observable = is_observable
371 self.phase = phase
373 if is_init:
374 if not is_observable:
375 pauli = pauli.generator()[0].base
376 self.list_repr = self._create_list_representation(pauli, n_qubits)
377 else:
378 assert isinstance(pauli, np.ndarray)
379 self.list_repr = pauli
381 if prev_xy_indices is None:
382 prev_xy_indices = np.zeros(n_qubits, dtype=bool)
383 self.xy_indices = np.logical_or(
384 prev_xy_indices,
385 self._compute_xy_indices(self.list_repr, rev=is_observable),
386 )
388 @staticmethod
389 def _compute_xy_indices(
390 op: np.ndarray[int], rev: bool = False
391 ) -> np.ndarray[bool]:
392 """
393 Computes the positions of X or Y gates to an one-hot encoded boolen
394 array.
396 Args:
397 op (np.ndarray[int]): Pauli-Operation list representation.
398 rev (bool): Whether to negate the array.
400 Returns:
401 np.ndarray[bool]: One hot encoded boolean array.
402 """
403 xy_indices = (op == 0) + (op == 1)
404 if rev:
405 xy_indices = ~xy_indices
406 return xy_indices
408 @staticmethod
409 def _create_list_representation(op: Operator, n_qubits: int) -> np.ndarray[int]:
410 """
411 Create list representation of a Pennylane Operator.
412 Generally, the list representation is a list of length n_qubits,
413 where at each position a Pauli Operator is encoded as such:
414 I: -1
415 X: 0
416 Y: 1
417 Z: 2
419 Args:
420 op (Operator): Pennylane Operator
421 n_qubits (int): number of qubits in the circuit
423 Returns:
424 np.ndarray[int]: List representation
425 """
426 pauli_repr = -np.ones(n_qubits, dtype=int)
427 op = op.terms()[1][0] if isinstance(op, qml_op.Prod) else op
428 op = op.base if isinstance(op, qml_op.SProd) else op
430 if isinstance(op, qml.PauliX):
431 pauli_repr[op.wires[0]] = 0
432 elif isinstance(op, qml.PauliY):
433 pauli_repr[op.wires[0]] = 1
434 elif isinstance(op, qml.PauliZ):
435 pauli_repr[op.wires[0]] = 2
436 else:
437 for o in op:
438 if isinstance(o, qml.PauliX):
439 pauli_repr[o.wires[0]] = 0
440 elif isinstance(o, qml.PauliY):
441 pauli_repr[o.wires[0]] = 1
442 elif isinstance(o, qml.PauliZ):
443 pauli_repr[o.wires[0]] = 2
444 return pauli_repr
446 def is_commuting(self, pauli: np.ndarray[int]) -> bool:
447 """
448 Computes if this Pauli commutes with another Pauli operator.
449 This computation is based on the fact that The commutator is zero
450 if and only if the number of anticommuting single-qubit Paulis is
451 even.
453 Args:
454 pauli (np.ndarray[int]): List representation of another Pauli
456 Returns:
457 bool: If the current and other Pauli are commuting.
458 """
459 anticommutator = np.where(
460 pauli < 0,
461 False,
462 np.where(
463 self.list_repr < 0,
464 False,
465 np.where(self.list_repr == pauli, False, True),
466 ),
467 )
468 return not (np.sum(anticommutator) % 2)
470 def tensor(self, pauli: np.ndarray[int]) -> FourierTree.PauliOperator:
471 """
472 Compute tensor product between the current Pauli and a given list
473 representation of another Pauli operator.
475 Args:
476 pauli (np.ndarray[int]): List representation of Pauli
478 Returns:
479 FourierTree.PauliOperator: New Pauli operator object, which
480 contains the tensor product
481 """
482 diff = (pauli - self.list_repr + 3) % 3
483 phase = np.where(
484 self.list_repr < 0,
485 1.0,
486 np.where(
487 pauli < 0,
488 1.0,
489 np.where(
490 diff == 2,
491 1.0j,
492 np.where(diff == 1, -1.0j, 1.0),
493 ),
494 ),
495 )
497 obs = np.where(
498 self.list_repr < 0,
499 pauli,
500 np.where(
501 pauli < 0,
502 self.list_repr,
503 np.where(
504 diff == 2,
505 (self.list_repr + 1) % 3,
506 np.where(diff == 1, (self.list_repr + 2) % 3, -1),
507 ),
508 ),
509 )
510 phase = self.phase * np.prod(phase)
511 return FourierTree.PauliOperator(
512 obs, phase=phase, n_qubits=obs.size, is_init=False, is_observable=True
513 )
515 def __init__(self, model: Model, inputs=1.0):
516 """
517 Tree initialisation, based on the Pauli-Clifford representation of a model.
518 Currently, only one input feature is supported.
520 **Usage**:
521 ```
522 # initialise a model
523 model = Model(...)
525 # initialise and build FourierTree
526 tree = FourierTree(model)
528 # get expectaion value
529 exp = tree()
531 # Get spectrum (for each observable, we have one list element)
532 coeff_list, freq_list = tree.spectrum()
533 ```
535 Args:
536 model (Model): The Model, for which to build the tree
537 inputs (bool, optional): Possible default inputs. Defaults to 1.0.
538 """
539 self.model = model
540 self.tree_roots = None
542 if not model.as_pauli_circuit:
543 model.as_pauli_circuit = True
545 inputs = (
546 self.model._inputs_validation(inputs)
547 if inputs is not None
548 else self.model._inputs_validation(1.0)
549 )
550 inputs.requires_grad = False
552 quantum_tape = qml.workflow.construct_tape(self.model.circuit)(
553 params=model.params, inputs=inputs
554 )
555 self.parameters = [np.squeeze(p) for p in quantum_tape.get_parameters()]
556 self.input_indices = [
557 i for (i, p) in enumerate(self.parameters) if not p.requires_grad
558 ]
560 self.observables = self._encode_observables(quantum_tape.observables)
562 pauli_rot = FourierTree.PauliOperator(
563 quantum_tape.operations[0],
564 self.model.n_qubits,
565 )
566 self.pauli_rotations = [pauli_rot]
567 for op in quantum_tape.operations[1:]:
568 pauli_rot = FourierTree.PauliOperator(
569 op, self.model.n_qubits, pauli_rot.xy_indices
570 )
571 self.pauli_rotations.append(pauli_rot)
573 self.tree_roots = self.build()
574 self.leafs: List[List[FourierTree.TreeLeaf]] = self._get_tree_leafs()
576 def __call__(
577 self,
578 params: Optional[np.ndarray] = None,
579 inputs: Optional[np.ndarray] = None,
580 **kwargs,
581 ) -> np.ndarray:
582 """
583 Evaluates the Fourier tree via sine-cosine terms sum. This is
584 equivalent to computing the expectation value of the observables with
585 respect to the corresponding circuit.
587 Args:
588 params (Optional[np.ndarray], optional): Parameters of the model.
589 Defaults to None.
590 inputs (Optional[np.ndarray], optional): Inputs to the circuit.
591 Defaults to None.
593 Returns:
594 np.ndarray: Expectation value of the tree.
596 Raises:
597 NotImplementedError: When using other "execution_type" as expval.
598 NotImplementedError: When using "noise_params"
601 """
602 params = (
603 self.model._params_validation(params)
604 if params is not None
605 else self.model.params
606 )
607 inputs = (
608 self.model._inputs_validation(inputs)
609 if inputs is not None
610 else self.model._inputs_validation(1.0)
611 )
612 inputs.requires_grad = False
614 if kwargs.get("execution_type", "expval") != "expval":
615 raise NotImplementedError(
616 f'Currently, only "expval" execution type is supported when '
617 f"building FourierTree. Got {kwargs.get('execution_type', 'expval')}."
618 )
619 if kwargs.get("noise_params", None) is not None:
620 raise NotImplementedError(
621 "Currently, noise is not supported when building FourierTree."
622 )
624 quantum_tape = qml.workflow.construct_tape(self.model.circuit)(
625 params=self.model.params, inputs=inputs
626 )
627 self.parameters = [np.squeeze(p) for p in quantum_tape.get_parameters()]
628 self.input_indices = [
629 i for (i, p) in enumerate(self.parameters) if not p.requires_grad
630 ]
632 results = np.zeros(len(self.tree_roots))
633 for i, root in enumerate(self.tree_roots):
634 results[i] = np.real_if_close(root.evaluate(self.parameters))
636 if kwargs.get("force_mean", False):
637 return np.mean(results)
638 else:
639 return results
641 def build(self) -> List[CoefficientsTreeNode]:
642 """
643 Creates the coefficient tree, i.e. it creates and initialises the tree
644 nodes.
645 Leafs can be obtained separately in _get_tree_leafs, once the tree is
646 set up.
648 Returns:
649 List[CoefficientsTreeNode]: The list of root nodes (one root for
650 each observable).
651 """
652 tree_roots = []
653 pauli_rotation_idx = len(self.pauli_rotations) - 1
654 for obs in self.observables:
655 root = self._create_tree_node(obs, pauli_rotation_idx)
656 tree_roots.append(root)
657 return tree_roots
659 def _encode_observables(
660 self, tape_obs: List[Operator]
661 ) -> List[FourierTree.PauliOperator]:
662 """
663 Encodes Pennylane observables from tape as FourierTree.PauliOperator
664 utility objects.
666 Args:
667 tape_obs (List[Operator]): Pennylane tape operations
669 Returns:
670 List[FourierTree.PauliOperator]: List of Pauli operators
671 """
672 observables = []
673 for obs in tape_obs:
674 observables.append(
675 FourierTree.PauliOperator(obs, self.model.n_qubits, is_observable=True)
676 )
677 return observables
679 def _get_tree_leafs(self) -> List[List[TreeLeaf]]:
680 """
681 Obtain all Leaf Nodes with its sine- and cosine terms.
683 Returns:
684 List[List[TreeLeaf]]: For each observable (root), the list of leaf
685 nodes.
686 """
687 leafs = []
688 for root in self.tree_roots:
689 sin_list = np.zeros(len(self.parameters), dtype=np.int32)
690 cos_list = np.zeros(len(self.parameters), dtype=np.int32)
691 leafs.append(root.get_leafs(sin_list, cos_list, []))
692 return leafs
694 def get_spectrum(
695 self, force_mean: bool = False
696 ) -> Tuple[List[np.ndarray], List[np.ndarray]]:
697 """
698 Computes the Fourier spectrum for the tree, consisting of the
699 frequencies and its corresponding coefficinets.
700 If the frag force_mean was set in the constructor, the mean coefficient
701 over all observables (roots) are computed.
703 Args:
704 force_mean (bool, optional): Whether to average over multiple
705 observables. Defaults to False.
707 Returns:
708 Tuple[List[np.ndarray], List[np.ndarray]]:
709 - List of frequencies, one list for each observable (root).
710 - List of corresponding coefficents, one list for each
711 observable (root).
712 """
713 parameter_indices = [
714 i for i in range(len(self.parameters)) if i not in self.input_indices
715 ]
717 coeffs = []
718 for leafs in self.leafs:
719 freq_terms = defaultdict(np.complex128)
720 for leaf in leafs:
721 leaf_factor, s, c = self._compute_leaf_factors(leaf, parameter_indices)
723 for a in range(s + 1):
724 for b in range(c + 1):
725 comb = math.comb(s, a) * math.comb(c, b) * (-1) ** (s - a)
726 freq_terms[2 * a + 2 * b - s - c] += comb * leaf_factor
728 coeffs.append(freq_terms)
730 frequencies, coefficients = self._freq_terms_to_coeffs(coeffs, force_mean)
731 return coefficients, frequencies
733 def _freq_terms_to_coeffs(
734 self, coeffs: List[Dict[int, np.ndarray]], force_mean: bool
735 ) -> Tuple[List[np.ndarray], List[np.ndarray]]:
736 """
737 Given a list of dictionaries of the form:
738 [
739 {
740 freq_obs1_1: coeff1,
741 freq_obs1_2: coeff2,
742 ...
743 },
744 {
745 freq_obs2_1: coeff3,
746 freq_obs2_2: coeff4,
747 ...
748 }
749 ...
750 ],
751 Compute two separate lists of frequencies and coefficients.
752 such that:
753 freqs: [
754 [freq_obs1_1, freq_obs1_1, ...],
755 [freq_obs2_1, freq_obs2_1, ...],
756 ...
757 ]
758 coeffs: [
759 [coeff1, coeff2, ...],
760 [coeff3, coeff4, ...],
761 ...
762 ]
764 If force_mean is set length of the resulting frequency and coefficent
765 list is 1.
767 Args:
768 coeffs (List[Dict[int, np.ndarray]]): Frequency->Coefficients
769 dictionary list, one dict for each observable (root).
770 force_mean (bool): Whether to average coefficients over multiple
771 observables.
773 Returns:
774 Tuple[List[np.ndarray], List[np.ndarray]]:
775 - List of frequencies, one list for each observable (root).
776 - List of corresponding coefficents, one list for each
777 observable (root).
778 """
779 frequencies = []
780 coefficients = []
781 if force_mean:
782 all_freqs = sorted(set([f for c in coeffs for f in c.keys()]))
783 coefficients.append(
784 np.array([np.mean([c.get(f, 0.0) for c in coeffs]) for f in all_freqs])
785 )
786 frequencies.append(np.array(all_freqs))
787 else:
788 for freq_terms in coeffs:
789 freq_terms = dict(sorted(freq_terms.items()))
790 frequencies.append(np.array(list(freq_terms.keys())))
791 coefficients.append(np.array(list(freq_terms.values())))
792 return frequencies, coefficients
794 def _compute_leaf_factors(
795 self, leaf: TreeLeaf, parameter_indices: List[int]
796 ) -> Tuple[float, int, int]:
797 """
798 Computes the constant coefficient factor for each leaf.
799 Additionally sine and cosine contributions of the input parameters for
800 this leaf are returned, which are required to obtain the corresponding
801 frequencies.
803 Args:
804 leaf (TreeLeaf): The leaf for which to compute the factor.
805 parameter_indices (List[int]): Variational parameter indices.
807 Returns:
808 Tuple[float, int, int]:
809 - float: the constant factor for the leaf
810 - int: number of sine contributions of the input
811 - int: number of cosine contributions of the input
812 """
813 leaf_factor = 1.0
814 for i in parameter_indices:
815 interm_factor = (
816 np.cos(self.parameters[i]) ** leaf.cos_indices[i]
817 * (1j * np.sin(self.parameters[i])) ** leaf.sin_indices[i]
818 )
819 leaf_factor = leaf_factor * interm_factor
821 # Get number of sine and cosine factors to which the input contributes
822 c = np.sum([leaf.cos_indices[k] for k in self.input_indices], dtype=np.int32)
823 s = np.sum([leaf.sin_indices[k] for k in self.input_indices], dtype=np.int32)
825 leaf_factor = leaf.term * leaf_factor * 0.5 ** (s + c)
827 return leaf_factor, s, c
829 def _early_stopping_possible(
830 self, pauli_rotation_idx: int, observable: FourierTree.PauliOperator
831 ):
832 """
833 Checks if a node for an observable can be discarded as all expecation
834 values that can result through further branching are zero.
835 The method is mentioned in the paper by Nemkov et al.: If the one-hot
836 encoded indices for X/Y operations in the Pauli-rotation generators are
837 a basis for that of the observable, the node must be processed further.
838 If not, it can be discarded.
840 Args:
841 pauli_rotation_idx (int): Index of remaining Pauli rotation gates.
842 Gates itself are attributes of the class.
843 observable (FourierTree.PauliOperator): Current observable
844 """
845 xy_indices_obs = np.logical_or(
846 observable.xy_indices, self.pauli_rotations[pauli_rotation_idx].xy_indices
847 ).all()
849 return not xy_indices_obs
851 def _create_tree_node(
852 self,
853 observable: FourierTree.PauliOperator,
854 pauli_rotation_idx: int,
855 parameter_idx: Optional[int] = None,
856 is_sine: bool = False,
857 is_cosine: bool = False,
858 ) -> Optional[CoefficientsTreeNode]:
859 """
860 Builds the Fourier-Tree according to the algorithm by Nemkov et al.
862 Args:
863 observable (FourierTree.PauliOperator): Current observable
864 pauli_rotation_idx (int): Index of remaining Pauli rotation gates.
865 Gates itself are attributes of the class.
866 parameter_idx (Optional[int]): Index of the current parameter.
867 Parameters itself are attributes of the class.
868 is_sine (bool): If the current node is a sine (left) node.
869 is_cosine (bool): If the current node is a cosine (right) node.
871 Returns:
872 Optional[CoefficientsTreeNode]: The resulting node. Children are set
873 recursively. The top level receives the tree root.
874 """
875 if self._early_stopping_possible(pauli_rotation_idx, observable):
876 return None
878 # remove commuting paulis
879 while pauli_rotation_idx >= 0:
880 last_pauli = self.pauli_rotations[pauli_rotation_idx]
881 if not observable.is_commuting(last_pauli.list_repr):
882 break
883 pauli_rotation_idx -= 1
884 else: # leaf
885 return FourierTree.CoefficientsTreeNode(
886 parameter_idx, observable, is_sine, is_cosine
887 )
889 last_pauli = self.pauli_rotations[pauli_rotation_idx]
891 left = self._create_tree_node(
892 observable,
893 pauli_rotation_idx - 1,
894 pauli_rotation_idx,
895 is_cosine=True,
896 )
898 next_observable = self._create_new_observable(last_pauli.list_repr, observable)
899 right = self._create_tree_node(
900 next_observable,
901 pauli_rotation_idx - 1,
902 pauli_rotation_idx,
903 is_sine=True,
904 )
906 return FourierTree.CoefficientsTreeNode(
907 parameter_idx,
908 observable,
909 is_sine,
910 is_cosine,
911 left,
912 right,
913 )
915 def _create_new_observable(
916 self, pauli: np.ndarray[int], observable: FourierTree.PauliOperator
917 ) -> FourierTree.PauliOperator:
918 """
919 Utility function to obtain the new observable for a tree node, if the
920 last Pauli and the observable do not commute.
922 Args:
923 pauli (np.ndarray[int]): The int array representation of the last
924 Pauli rotation in the operation sequence.
925 observable (FourierTree.PauliOperator): The current observable.
927 Returns:
928 FourierTree.PauliOperator: The new observable.
929 """
930 observable = observable.tensor(pauli)
931 return observable
934class FCC:
935 @staticmethod
936 def get_fcc(
937 model: Model,
938 n_samples: int,
939 seed: int,
940 method: Optional[str] = "pearson",
941 scale: Optional[bool] = False,
942 weight: Optional[bool] = False,
943 trim_redundant: Optional[bool] = True,
944 **kwargs,
945 ) -> float:
946 """
947 Shortcut method to get just the FCC.
948 This includes
949 1. What is done in `get_fourier_fingerprint`:
950 1. Calculating the coefficients (using `n_samples` and `seed`)
951 2. Correlating the result from 1) using `method`
952 3. Weighting the correlation matrix (if `weight` is True)
953 4. Remove redundancies
954 2. What is done in `calculate_fcc`:
955 1. Absolute of the fingerprint
956 2. Average
958 Args:
959 model (Model): The QFM model
960 n_samples (int): Number of samples to calculate average of coefficients
961 seed (int): Seed to initialize random parameters
962 method (Optional[str], optional): Correlation method. Defaults to "pearson".
963 scale (Optional[bool], optional): Whether to scale the number of samples.
964 Defaults to False.
965 weight (Optional[bool], optional): Whether to weight the correlation matrix.
966 Defaults to False.
967 trim_redundant (Optional[bool], optional): Whether to remove redundant
968 correlations. Defaults to False.
969 **kwargs: Additional keyword arguments for the model function.
971 Returns:
972 float: The FCC
973 """
974 fourier_fingerprint, _ = FCC.get_fourier_fingerprint(
975 model,
976 n_samples,
977 seed,
978 method,
979 scale,
980 weight,
981 trim_redundant=trim_redundant,
982 **kwargs,
983 )
985 return FCC.calculate_fcc(fourier_fingerprint)
987 def get_fourier_fingerprint(
988 model: Model,
989 n_samples: int,
990 seed: int,
991 method: Optional[str] = "pearson",
992 scale: Optional[bool] = False,
993 weight: Optional[bool] = False,
994 trim_redundant: Optional[bool] = True,
995 **kwargs,
996 ) -> Tuple[np.ndarray, np.ndarray]:
997 """
998 Shortcut method to get just the fourier fingerprint.
999 This includes
1000 1. Calculating the coefficients (using `n_samples` and `seed`)
1001 2. Correlating the result from 1) using `method`
1002 3. Weighting the correlation matrix (if `weight` is True)
1003 4. Remove redundancies (if `trim_redundant` is True)
1005 Args:
1006 model (Model): The QFM model
1007 n_samples (int): Number of samples to calculate average of coefficients
1008 seed (int): Seed to initialize random parameters
1009 method (Optional[str], optional): Correlation method. Defaults to "pearson".
1010 scale (Optional[bool], optional): Whether to scale the number of samples.
1011 Defaults to False.
1012 weight (Optional[bool], optional): Whether to weight the correlation matrix.
1013 Defaults to False.
1014 trim_redundant (Optional[bool], optional): Whether to remove redundant
1015 correlations. Defaults to True.
1016 **kwargs: Additional keyword arguments for the model function.
1018 Returns:
1019 Tuple[np.ndarray, np.ndarray]: The fourier fingerprint
1020 and the frequency indices
1021 """
1022 _, coeffs, freqs = FCC._calculate_coefficients(
1023 model, n_samples, seed, scale, **kwargs
1024 )
1025 fourier_fingerprint = FCC._correlate(coeffs.transpose(), method=method)
1027 # perform weighting if requested
1028 fourier_fingerprint = (
1029 FCC._weighting(fourier_fingerprint) if weight else fourier_fingerprint
1030 )
1032 if trim_redundant:
1033 mask = FCC._calculate_mask(freqs)
1035 # apply the mask on the fingerprint
1036 fourier_fingerprint = mask * fourier_fingerprint
1038 row_mask = np.any(np.isfinite(fourier_fingerprint), axis=1)
1039 col_mask = np.any(np.isfinite(fourier_fingerprint), axis=0)
1041 fourier_fingerprint = fourier_fingerprint[row_mask][:, col_mask]
1043 return fourier_fingerprint, freqs
1045 @staticmethod
1046 def calculate_fcc(
1047 fourier_fingerprint: np.ndarray,
1048 ) -> float:
1049 """
1050 Method to calculate the FCC based on an existing correlation matrix.
1051 Calculate absolute and then the average over this matrix.
1052 The Fingerprint can be obtained via `get_fourier_fingerprint`
1054 Args:
1055 coeff_coeff_correlation (np.ndarray): Correlation matrix of coefficients
1056 Returns:
1057 float: The FCC
1058 """
1059 # apply the mask on the fingerprint
1060 return np.nanmean(np.abs(fourier_fingerprint))
1062 def _calculate_mask(freqs: np.ndarray) -> np.ndarray:
1063 """
1064 Method to calculate a mask filtering out redundant elements
1065 of the Fourier correlation matrix, based on the provided frequency vector.
1066 It does so by 'simulating' the operations that would be performed
1067 by `_correlate`.
1069 Args:
1070 freqs (np.ndarray): Array of frequencies
1072 Returns:
1073 np.ndarray: The mask
1074 """
1075 # TODO: this part can be heavily optimized, by e.g. using a "positive_only"
1076 # flag when calculating the coefficients.
1077 # However this would change the numerical values
1078 # (while the order should be still the same).
1080 # disregard all the negativ frequencies
1081 freqs[freqs < 0] = np.nan
1082 # compute the outer product of the frequency vectors for arbitrary dimensions
1083 # or just use the existing frequency vector if it is 1D
1084 nd_freqs = (
1085 reduce(np.multiply, np.ix_(*freqs)) if len(freqs.shape) > 1 else freqs
1086 )
1087 # TODO: could prevent this if we're not using .squeeze()..
1089 # "simulate" what would happen on correlating the coefficients
1090 corr_freqs = np.outer(nd_freqs, nd_freqs)
1091 # mask all frequencies that are nan now
1092 # (i.e. all correlations with a negative frequency component)
1093 corr_mask = np.where(np.isnan(corr_freqs), corr_freqs, 1)
1094 # from this, disregard all the other redundant correlations (i.e. c_0_1 = c_1_0)
1095 corr_mask[np.triu_indices(corr_mask.shape[0], 0)] = np.nan
1097 return corr_mask
1099 @staticmethod
1100 def _calculate_coefficients(
1101 model: Model,
1102 n_samples: int,
1103 seed: int,
1104 scale: bool = False,
1105 **kwargs,
1106 ) -> Tuple[np.ndarray, np.ndarray]:
1107 """
1108 Calculates the Fourier coefficients of a given model
1109 using `n_samples` and `seed`.
1110 Optionally, `noise_params` can be passed to perform noisy simulation.
1112 Args:
1113 model (Model): The QFM model
1114 n_samples (int): Number of samples to calculate average of coefficients
1115 seed (int): Seed to initialize random parameters
1116 scale (bool, optional): Whether to scale the number of samples.
1117 Defaults to False.
1118 **kwargs: Additional keyword arguments for the model function.
1120 Returns:
1121 Tuple[np.ndarray, np.ndarray]: Parameters and Coefficients of size NxK
1122 """
1123 if n_samples > 0:
1124 if scale:
1125 total_samples = int(
1126 np.power(2, model.n_qubits) * n_samples * model.n_input_feat
1127 )
1128 log.info(f"Using {total_samples} samples.")
1129 else:
1130 total_samples = n_samples
1131 rng = np.random.default_rng(seed)
1132 model.initialize_params(rng=rng, repeat=total_samples)
1133 else:
1134 total_samples = 1
1136 coeffs, freqs = Coefficients.get_spectrum(
1137 model, shift=True, trim=True, **kwargs
1138 )
1140 return model.params, coeffs, freqs
1142 @staticmethod
1143 def _correlate(mat: np.ndarray, method: str = "pearson") -> np.ndarray:
1144 """
1145 Correlates two arrays using `method`.
1146 Currently, `pearson` and `spearman` are supported.
1148 Args:
1149 mat (np.ndarray): Array of shape (N, K)
1150 method (str, optional): Correlation method. Defaults to "pearson".
1152 Raises:
1153 ValueError: If the method is not supported.
1155 Returns:
1156 np.ndarray: Correlation matrix of `a` and `b`.
1157 """
1158 assert len(mat.shape) >= 2, "Input matrix must have at least 2 dimensions"
1160 # Note that for the general n-D case, we have to flatten along
1161 # the first axis (last one is batch).
1162 # Note that the order here is important so we can easily filter out
1163 # negative coefficients later.
1164 # Consider the following example: [[1,2,3],[4,5,6],[7,8,9]]
1165 # we want to get [1, 4, 7, 2, 5, 8, 3, 6, 9]
1166 # such that after correlation, all positive indexed coefficients
1167 # will be in the bottom right quadrant
1168 if method == "pearson":
1169 result = FCC._pearson(mat.reshape(mat.shape[0], -1))
1170 # result = FCC._pearson(mat.reshape(mat.shape[-1], -1, order="F"))
1171 elif method == "spearman":
1172 result = FCC._spearman(mat.reshape(mat.shape[0], -1))
1173 # result = FCC._spearman(mat.reshape(mat.shape[-1], -1, order="F"))
1174 else:
1175 raise ValueError(
1176 f"Unknown correlation method: {method}. \
1177 Must be 'pearson' or 'spearman'."
1178 )
1180 return result
1182 @staticmethod
1183 def _pearson(
1184 mat: np.ndarray, cov: Optional[bool] = False, minp: Optional[int] = 1
1185 ) -> np.ndarray:
1186 """
1187 Based on Pandas correlation method as implemented here:
1188 https://github.com/pandas-dev/pandas/blob/main/pandas/_libs/algos.pyx
1190 Compute Pearson correlation between columns of `mat`,
1191 permitting missing values (NaN or ±Inf).
1193 Args:
1194 mat : array_like, shape (N, K)
1195 Input data.
1196 minp : int, optional
1197 Minimum number of paired observations required to form a correlation.
1198 If the number of valid pairs for (i, j) is < minp, the result is NaN.
1200 Returns:
1201 corr : ndarray, shape (K, K)
1202 Pearson correlation matrix.
1203 """
1205 mat = np.asarray(mat, dtype=np.float64)
1206 N, K = mat.shape
1208 # pre‐compute finite‐mask
1209 mask = np.isfinite(mat)
1211 # output
1212 result = np.empty((K, K), dtype=np.float64)
1214 # TODO: optimize in future iterations
1215 # loop over column‐pairs
1216 for i in range(K):
1217 for j in range(i + 1):
1218 # find rows where both columns are finite
1219 m = mask[:, i] & mask[:, j]
1220 n = np.count_nonzero(m)
1221 if n < minp:
1222 # too few pairs
1223 value = np.nan
1224 else:
1225 x = mat[m, i]
1226 y = mat[m, j]
1228 # compute means
1229 mean_x = x.mean()
1230 mean_y = y.mean()
1232 # demeaned data
1233 dx = x - mean_x
1234 dy = y - mean_y
1236 # sum of squares and cross‐prod
1237 ssx = np.dot(dx, dx)
1238 ssy = np.dot(dy, dy)
1239 cxy = np.dot(dx, dy)
1241 if cov:
1242 # sample covariance (denominator n−1)
1243 value = cxy / (n - 1) if n > 1 else np.nan
1244 else:
1245 # Pearson r = cov / (σx σy)
1246 denom = np.sqrt(ssx * ssy)
1247 if denom == 0.0:
1248 value = np.nan
1249 else:
1250 value = cxy / denom
1251 # clip numerical drift
1252 if value > 1.0:
1253 value = 1.0
1254 elif value < -1.0:
1255 value = -1.0
1257 result[i, j] = result[j, i] = value
1259 return result
1261 def _spearman(mat: np.ndarray, minp: Optional[int] = 1) -> np.ndarray:
1262 """
1263 Based on Pandas correlation method as implemented here:
1264 https://github.com/pandas-dev/pandas/blob/main/pandas/_libs/algos.pyx
1266 Compute Spearman correlation between columns of `mat`,
1267 permitting missing values (NaN or ±Inf).
1269 Args:
1270 mat : array_like, shape (N, K)
1271 Input data.
1272 minp : int, optional
1273 Minimum number of paired observations required to form a correlation.
1274 If the number of valid pairs for (i, j) is < minp, the result is NaN.
1276 Returns:
1277 corr : ndarray, shape (K, K)
1278 Spearman correlation matrix.
1279 """
1280 N, K = mat.shape
1281 # trivial all-NaN answer if too few rows
1282 if N < minp:
1283 return np.full((K, K), np.nan, dtype=float)
1285 # mask of finite entries
1286 mask = np.isfinite(mat) # shape (N, K), dtype=bool
1288 # precompute ranks column-wise ignoring NaNs
1289 ranks = np.full((N, K), np.nan, dtype=float)
1290 for j in range(K):
1291 valid = mask[:, j]
1292 if valid.any():
1293 # rankdata by default gives average ranks for ties
1294 ranks[valid, j] = rankdata(mat[valid, j], method="average")
1296 # allocate result
1297 result = np.empty((K, K), dtype=float)
1299 # TODO: optimize in future iterations
1300 # loop lower triangle (including diagonal)
1301 for i in range(K):
1302 for j in range(i + 1):
1303 # find rows where both columns are finite
1304 valid = mask[:, i] & mask[:, j]
1305 nobs = valid.sum()
1307 if nobs < minp:
1308 rho = np.nan
1309 else:
1310 xi = ranks[valid, i]
1311 yj = ranks[valid, j]
1312 # subtract means
1313 xi = xi - xi.mean()
1314 yj = yj - yj.mean()
1315 num = np.dot(xi, yj)
1316 den = np.sqrt(np.dot(xi, xi) * np.dot(yj, yj))
1317 rho = num / den if den > 0 else np.nan
1319 result[i, j] = rho
1320 result[j, i] = rho
1322 return result
1324 @staticmethod
1325 def _weighting(fourier_fingerprint: np.ndarray) -> np.ndarray:
1326 """
1327 Performs weighting on the given correlation matrix.
1328 Here, low-frequent coefficients are weighted more heavily.
1330 Args:
1331 correlation (np.ndarray): Correlation matrix
1332 """
1333 # TODO: in Future iterations, this can be optimized by computing
1334 # on the trimmed matrix instead.
1336 assert (
1337 fourier_fingerprint.shape[0] % 2 != 0
1338 and fourier_fingerprint.shape[1] % 2 != 0
1339 ), "Correlation matrix must have odd dimensions. \
1340 Hint: use `trim` argument when calling `get_spectrum`."
1341 assert (
1342 fourier_fingerprint.shape[0] == fourier_fingerprint.shape[1]
1343 ), "Correlation matrix must be square."
1345 def quadrant_to_matrix(a: np.ndarray) -> np.ndarray:
1346 """
1347 Transforms [[1,2],[3,4]] to
1348 [[1,2,1],[3,4,3],[1,2,1]]
1350 Args:
1351 a (np.ndarray): _description_
1353 Returns:
1354 np.ndarray: _description_
1355 """
1356 # rotates a from [[1,2],[3,4]] to [[3,4],[1,2]]
1357 a_rot = np.rot90(a)
1358 # merge the two matrices
1359 left = np.concat([a, a_rot])
1360 # merges left and right (left flipped)
1361 b = np.concat(
1362 [left, np.flip(left)],
1363 axis=1,
1364 )
1365 # remove the middle column and row
1366 return np.delete(
1367 np.delete(b, (b.shape[0] // 2), axis=0), (b.shape[1] // 2), axis=1
1368 )
1370 nc = fourier_fingerprint.shape[0] // 2 + 1
1371 weights = nnp.mgrid[0:nc:1, 0:nc:1].sum(axis=0) / ((nc - 1) * 2)
1372 weights_matrix = quadrant_to_matrix(weights)
1374 return fourier_fingerprint * weights_matrix
1375 raise NotImplementedError("Weighting method is not implemented")