microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
source/noisy_simulator/src/instrument.rs
162lines · modecode
| 1 | // Copyright (c) Microsoft Corporation. |
| 2 | // Licensed under the MIT License. |
| 3 | |
| 4 | //! This module contains the `Instrument` struct, used to make measurments |
| 5 | //! in a quantum system. |
| 6 | |
| 7 | #[cfg(test)] |
| 8 | mod tests; |
| 9 | |
| 10 | use crate::{operation::Operation, Error, SquareMatrix, TOLERANCE}; |
| 11 | use nalgebra::{DMatrix, DVector}; |
| 12 | |
| 13 | /// An instrument is the means by which we make measurements on a quantum system. |
| 14 | pub struct Instrument { |
| 15 | operations: Vec<Operation>, |
| 16 | summed_operation: SquareMatrix, |
| 17 | summed_effect: SquareMatrix, |
| 18 | summed_effect_transpose: SquareMatrix, |
| 19 | summed_kraus_operators: Vec<SquareMatrix>, |
| 20 | } |
| 21 | |
| 22 | impl Instrument { |
| 23 | /// Creates a new instrument. |
| 24 | pub fn new(operations: Vec<Operation>) -> Result<Self, Error> { |
| 25 | if operations.is_empty() { |
| 26 | return Err(Error::FailedToConstructInstrument( |
| 27 | "there should be at least one Operation".to_string(), |
| 28 | )); |
| 29 | } |
| 30 | |
| 31 | let number_of_qubits = operations |
| 32 | .first() |
| 33 | .expect("there should be at least an element") |
| 34 | .number_of_qubits(); |
| 35 | |
| 36 | if !operations |
| 37 | .iter() |
| 38 | .all(|op| number_of_qubits == op.number_of_qubits()) |
| 39 | { |
| 40 | return Err(Error::FailedToConstructInstrument( |
| 41 | "all Operations should target the same number of qubits".to_string(), |
| 42 | )); |
| 43 | } |
| 44 | |
| 45 | let summed_operation: SquareMatrix = operations.iter().map(Operation::matrix).sum(); |
| 46 | let summed_effect: SquareMatrix = operations.iter().map(Operation::effect_matrix).sum(); |
| 47 | let summed_effect_transpose = summed_effect.transpose(); |
| 48 | let summed_kraus_operators = summed_kraus_operators(&operations)?; |
| 49 | |
| 50 | Ok(Self { |
| 51 | operations, |
| 52 | summed_operation, |
| 53 | summed_effect, |
| 54 | summed_effect_transpose, |
| 55 | summed_kraus_operators, |
| 56 | }) |
| 57 | } |
| 58 | |
| 59 | /// Return operation corresponding to the i-th outcome. |
| 60 | #[must_use] |
| 61 | pub fn operation(&self, i: usize) -> &Operation { |
| 62 | &self.operations[i] |
| 63 | } |
| 64 | |
| 65 | /// Return the matrix corresponding to the sum over all |
| 66 | /// operations in this instrument: |
| 67 | /// Σᵢ Σₖ (Kᵢₖ ⊗ Kᵢₖ*) |
| 68 | /// where ⊗ denotes the Kronecker product |
| 69 | /// and * denotes the complex conjugate (entry-wise). |
| 70 | #[must_use] |
| 71 | pub fn non_selective_operation_matrix(&self) -> &SquareMatrix { |
| 72 | &self.summed_operation |
| 73 | } |
| 74 | |
| 75 | /// Return Kraus operators of the operation corresponding to non selective evolution. |
| 76 | #[must_use] |
| 77 | pub fn non_selective_kraus_operators(&self) -> &Vec<SquareMatrix> { |
| 78 | &self.summed_kraus_operators |
| 79 | } |
| 80 | |
| 81 | /// Return total effect Σᵢ Σₖ (Kᵢₖ† Kᵢₖ), where † denotes the adjoint. |
| 82 | #[must_use] |
| 83 | pub fn total_effect(&self) -> &SquareMatrix { |
| 84 | &self.summed_effect |
| 85 | } |
| 86 | |
| 87 | /// Return transposed total effect Σᵢ Σₖ (Kᵢₖ† Kᵢₖ)^T, where † denotes the adjoint. |
| 88 | #[must_use] |
| 89 | pub fn total_effect_transposed(&self) -> &SquareMatrix { |
| 90 | &self.summed_effect_transpose |
| 91 | } |
| 92 | |
| 93 | /// Return number of operations/outcomes in this instrument. |
| 94 | #[must_use] |
| 95 | pub fn num_operations(&self) -> usize { |
| 96 | self.operations.len() |
| 97 | } |
| 98 | } |
| 99 | |
| 100 | fn summed_kraus_operators(operations: &[Operation]) -> Result<Vec<SquareMatrix>, Error> { |
| 101 | let choi_matrix = compute_choi_matrix(operations); |
| 102 | let (choi_dim, _) = choi_matrix.shape(); |
| 103 | let eigen_decomposition = choi_matrix.symmetric_eigen(); |
| 104 | let eigenvectors = eigen_decomposition.eigenvectors; |
| 105 | let eigenvalues = eigen_decomposition.eigenvalues; |
| 106 | let mut summed_kraus_operators = Vec::with_capacity(eigenvalues.len()); |
| 107 | |
| 108 | let (krows, kcols) = operations[0].kraus_operators()[0].shape(); |
| 109 | for (col, eigenvalue) in eigenvalues.iter().enumerate() { |
| 110 | if *eigenvalue > 0. { |
| 111 | let mut kraus_operator = SquareMatrix::zeros(krows, kcols); |
| 112 | let sqrt_eigenvalue = eigenvalue.sqrt(); |
| 113 | |
| 114 | // Performance note: for performance reasons we transpose the |
| 115 | // Kraus operators before storing them. We want to store each column |
| 116 | // of the eigenvector matrix as a Kraus matrix in row major order, and |
| 117 | // then transpose it. But because nalgebra stores and iterates over |
| 118 | // its matrices in column major order, we can skip the transposition, |
| 119 | // since it is already implicit. |
| 120 | // |
| 121 | // See noisy_simulator/src/operation.rs/Operation::new for more details. |
| 122 | for row in 0..choi_dim { |
| 123 | let idx = kraus_operator.vector_to_matrix_index(row); |
| 124 | kraus_operator[idx] = sqrt_eigenvalue * eigenvectors[(row, col)]; |
| 125 | } |
| 126 | |
| 127 | summed_kraus_operators.push(kraus_operator); |
| 128 | } else if *eigenvalue < -TOLERANCE { |
| 129 | return Err(Error::FailedToConstructInstrument(format!( |
| 130 | "failed to decompose Choi matrix, lambda_i = {eigenvalue}" |
| 131 | ))); |
| 132 | } |
| 133 | } |
| 134 | |
| 135 | Ok(summed_kraus_operators) |
| 136 | } |
| 137 | |
| 138 | fn compute_choi_matrix(operations: &[Operation]) -> SquareMatrix { |
| 139 | operations |
| 140 | .iter() |
| 141 | .map(|op| { |
| 142 | op.kraus_operators() |
| 143 | .iter() |
| 144 | .map(|k| { |
| 145 | let vectorized_k = vectorize(k); |
| 146 | &vectorized_k * &vectorized_k.adjoint() |
| 147 | }) |
| 148 | .sum::<SquareMatrix>() |
| 149 | }) |
| 150 | .sum() |
| 151 | } |
| 152 | |
| 153 | /// Stacks the columns of `matrix` into a single column vector. |
| 154 | /// |
| 155 | /// Performance note: Typically vectorization stacks the |
| 156 | /// rows of a matrix into a single column vector. But since we |
| 157 | /// transposed all matrices until now, we stack the columns instead. |
| 158 | /// |
| 159 | /// See `noisy_simulator/src/operation.rs/Operation::new` for more details. |
| 160 | fn vectorize<T: nalgebra::Scalar + Copy>(matrix: &DMatrix<T>) -> DVector<T> { |
| 161 | DVector::<T>::from_iterator(matrix.len(), matrix.iter().copied()) |
| 162 | } |
| 163 | |