stabilizer_ch_form_rust/form/
inner_product.rs

1use crate::StabilizerCHForm;
2use crate::error::{Error, Result};
3use crate::form::types::{InternalGate, PhaseFactor};
4use num_complex::Complex64;
5
6impl StabilizerCHForm {
7    /// Computes the inner product 〈self|other〉.
8    ///
9    /// ## Arguments
10    /// * `other` - The other StabilizerCHForm to compute the inner product with.
11    ///
12    /// ## Returns
13    /// A [`Result`] containing the complex inner product value.
14    pub fn inner_product(&self, other: &StabilizerCHForm) -> Result<Complex64> {
15        // TODO: Implement batch inner product calculation since the result of
16        // `_self._get_normalize_to_zero_ops()` can be reused.
17        if self.n != other.n {
18            return Err(Error::QubitCountMismatch {
19                operation: "calculating inner product",
20                left: self.n,
21                right: other.n,
22            });
23        }
24
25        // Get operations to transform `self` to |0...0>
26        // i.e. U_{ops} |self> = global_phase * phase * |0...0>
27        let (ops, phase) = self.get_normalize_to_zero_ops()?;
28
29        // Apply the same operations to `other`
30        // i.e. U_{ops} |other> = |transformed_other>
31        let transformed_other = other.get_ops_applied_state(&ops)?;
32
33        // Get the amplitude of |0...0> in `transformed_other`
34        // i.e. res = <0...0| U_{ops} |other>
35        let res = transformed_other.amplitude_at_zero()?;
36
37        // Combine the results
38        // The inner product is <self|other> = <self|U_dag U|other>
39        // Since U|self> = omega * phase * |0...0>, then <self|U_dag = omega.conj() * phase.conj() * <0...0|
40        // So, <self|other> = omega.conj() * phase.conj() * <0...0|U|other>
41        // <0...0|U|other> is `res.to_complex() / other.omega`, because _get_ops_applied_state carries over the omega.
42        // The total omega of transformed_other is other.omega, but _amplitude_at_zero doesn't include it.
43        // So res.to_complex() is the part without the original omega.
44        let inner_product_val = (res * phase.conjugated()).to_complex();
45
46        // We need to account for the global phases
47        Ok(self.global_phase() * other.global_phase() * inner_product_val)
48    }
49
50    /// Returns the sequence of operations needed to transform the current state to |0...0>
51    /// along with the phase factor of the resulting state.
52    fn get_normalize_to_zero_ops(&self) -> Result<(Vec<InternalGate>, PhaseFactor)> {
53        let mut ops = Vec::new();
54        let mut self_clone = self.clone();
55        let n = self_clone.n;
56
57        // Step 1: Convert G to identity matrix using left CNOTs
58        // NOTE: When G is converted to identity, F also becomes identity
59        for j in 0..n {
60            let mut pivot_row = j;
61            if !self_clone.mat_g[[j, j]] {
62                if let Some(k) = (j + 1..n).find(|&k| self_clone.mat_g[[k, j]]) {
63                    pivot_row = k;
64                } else {
65                    // Unreachable if the state is valid
66                    unreachable!("G matrix is not full rank, invalid stabilizer state.");
67                }
68            }
69
70            if pivot_row != j {
71                // Swap rows j and pivot_row using CNOTs: (k,j), (j,k), (k,j)
72                ops.push(InternalGate::CX(pivot_row, j));
73                self_clone.left_multiply_cx(pivot_row, j)?;
74                ops.push(InternalGate::CX(j, pivot_row));
75                self_clone.left_multiply_cx(j, pivot_row)?;
76                ops.push(InternalGate::CX(pivot_row, j));
77                self_clone.left_multiply_cx(pivot_row, j)?;
78            }
79
80            for i in 0..n {
81                if i != j && self_clone.mat_g[[i, j]] {
82                    ops.push(InternalGate::CX(j, i));
83                    self_clone.left_multiply_cx(j, i)?;
84                }
85            }
86        }
87
88        // Step 2-1: Convert off-diagonal M to zero using left CZs
89        for r in 0..n {
90            for c in (r + 1)..n {
91                if self_clone.mat_m[[r, c]] {
92                    ops.push(InternalGate::CZ(r, c));
93                    self_clone.left_multiply_cz(r, c)?;
94                }
95            }
96        }
97
98        // Step 2-2: Convert diagonal M to zero using left Sdg gates
99        for q in 0..n {
100            if self_clone.mat_m[[q, q]] {
101                ops.push(InternalGate::Sdg(q));
102                self_clone.left_multiply_sdg(q)?;
103            }
104        }
105
106        // Step 3: Convert vec_v to zero using Hs
107        for i in 0..n {
108            if self_clone.vec_v[i] {
109                ops.push(InternalGate::H(i));
110                self_clone.left_multiply_h(i)?;
111            }
112        }
113
114        // Step 4: Convert vec_s to zero using Xs
115        for i in 0..n {
116            if self_clone.vec_s[i] {
117                ops.push(InternalGate::X(i));
118                self_clone.left_multiply_x(i)?;
119            }
120        }
121
122        Ok((ops, self_clone.phase_factor))
123    }
124
125    fn get_ops_applied_state(&self, ops: &[InternalGate]) -> Result<StabilizerCHForm> {
126        let mut new_state = self.clone();
127        for op in ops {
128            match op {
129                InternalGate::H(q) => new_state.left_multiply_h(*q)?,
130                // InternalGate::S(q) => new_state._left_multiply_s(*q),
131                InternalGate::Sdg(q) => new_state.left_multiply_sdg(*q)?,
132                InternalGate::X(q) => new_state.left_multiply_x(*q)?,
133                InternalGate::CX(c, t) => new_state.left_multiply_cx(*c, *t)?,
134                InternalGate::CZ(c, t) => new_state.left_multiply_cz(*c, *t)?,
135            }
136        }
137        Ok(new_state)
138    }
139}
140
141#[cfg(test)]
142mod tests {
143    use crate::circuit::CliffordCircuit;
144
145    use super::*;
146
147    #[test]
148    fn test_random_state_inner_product() {
149        let num_qubits = 4;
150        for i in 0..10 {
151            let state1 = StabilizerCHForm::from_clifford_circuit(
152                &CliffordCircuit::random_clifford(num_qubits, Some([i + 12; 32])),
153            )
154            .unwrap();
155            let state2 = StabilizerCHForm::from_clifford_circuit(
156                &CliffordCircuit::random_clifford(num_qubits, Some([i + 34; 32])),
157            )
158            .unwrap();
159
160            let inner_product = state1.inner_product(&state2).unwrap();
161
162            let statevector1 = state1.to_statevector().unwrap();
163            let statevector2 = state2.to_statevector().unwrap();
164            let expected_inner_product = statevector1
165                .iter()
166                .zip(statevector2.iter())
167                .map(|(a, b)| a.conj() * b)
168                .sum::<Complex64>();
169
170            assert!((inner_product - expected_inner_product).norm() < 1e-8);
171        }
172    }
173}
stabilizer_ch_form_rust/form/inner_product.rs

stabilizer_ch_form_rust/form/
inner_product.rs
