cqlib_core/compile/commutation/
matrix.rs

1// This code is part of Cqlib.
2//
3// (C) Copyright China Telecom Quantum Group 2025-2026
4//
5// This code is licensed under the Apache License, Version 2.0. You may
6// obtain a copy of this license in the LICENSE.txt file in the root directory
7// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
8//
9// Any modifications or derivative works of this code must retain this
10// copyright notice, and modified files need to carry a notice indicating
11// that they have been altered from the originals.
12
13//! Small local-matrix fallback for commutation checks.
14//!
15//! This module proves commutation by embedding both operations into the matrix
16//! space of their combined qubit support, multiplying them in both orders, and
17//! comparing the products up to a global phase.  It is intended as a bounded
18//! fallback for concrete parameters after cheaper symbolic checks fail.
19//!
20//! The fallback does not attempt symbolic evaluation.  Any parameter that
21//! cannot be evaluated to a concrete `f64`, unsupported matrix construction, or
22//! support wider than the configured qubit bound produces `None`.
23
24use crate::circuit::{Instruction, Parameter, Qubit};
25use crate::compile::commutation::checker::{Commutation, CommutationResult};
26use crate::compile::{NUMERIC_ZERO_TOLERANCE, UNIT_PHASE_NORM_TOLERANCE};
27use ndarray::Array2;
28use num_complex::Complex64;
29use std::borrow::Cow;
30
31const MATRIX_TOLERANCE: f64 = 1e-10;
32
33/// Attempts matrix-based commutation on the union support of two operations.
34///
35/// The returned phase satisfies `lhs * rhs = exp(i * phase) * rhs * lhs` on the
36/// expanded local Hilbert space.  The caller controls the maximum support size
37/// to prevent exponential blow-up.
38pub fn matrix_commutation(
39    lhs_inst: &Instruction,
40    lhs_qubits: &[Qubit],
41    lhs_params: &[Parameter],
42    rhs_inst: &Instruction,
43    rhs_qubits: &[Qubit],
44    rhs_params: &[Parameter],
45    max_matrix_qubits: usize,
46) -> CommutationResult {
47    let mut combined_qubits = Vec::with_capacity(lhs_qubits.len() + rhs_qubits.len());
48    for &qubit in lhs_qubits.iter().chain(rhs_qubits) {
49        if !combined_qubits.contains(&qubit) {
50            combined_qubits.push(qubit);
51        }
52    }
53    combined_qubits.sort_unstable();
54
55    if combined_qubits.len() > max_matrix_qubits {
56        return None;
57    }
58
59    let lhs_values = lhs_params
60        .iter()
61        .map(|param| param.evaluate(&None).ok())
62        .collect::<Option<Vec<_>>>()?;
63    let rhs_values = rhs_params
64        .iter()
65        .map(|param| param.evaluate(&None).ok())
66        .collect::<Option<Vec<_>>>()?;
67    let lhs_matrix = lhs_inst.matrix(&lhs_values)?;
68    let rhs_matrix = rhs_inst.matrix(&rhs_values)?;
69
70    let lhs_positions = lhs_qubits
71        .iter()
72        .map(|qubit| combined_qubits.iter().position(|item| item == qubit))
73        .collect::<Option<Vec<_>>>()?;
74    let rhs_positions = rhs_qubits
75        .iter()
76        .map(|qubit| combined_qubits.iter().position(|item| item == qubit))
77        .collect::<Option<Vec<_>>>()?;
78
79    let expanded_lhs = expand_unitary(lhs_matrix, &lhs_positions, combined_qubits.len())?;
80    let expanded_rhs = expand_unitary(rhs_matrix, &rhs_positions, combined_qubits.len())?;
81
82    let lhs_rhs = expanded_lhs.dot(&expanded_rhs);
83    let rhs_lhs = expanded_rhs.dot(&expanded_lhs);
84
85    let phase = phase_between(&lhs_rhs, &rhs_lhs)?;
86
87    if phase.abs() <= MATRIX_TOLERANCE {
88        Some(Commutation::Exact)
89    } else {
90        let two_pi = 2.0 * std::f64::consts::PI;
91        let mut normalized =
92            (phase + std::f64::consts::PI).rem_euclid(two_pi) - std::f64::consts::PI;
93        if normalized <= -std::f64::consts::PI {
94            normalized += two_pi;
95        }
96        Some(Commutation::UpToGlobalPhase(Parameter::from(normalized)))
97    }
98}
99
100/// Expands a gate matrix to the full local space described by `total_qubits`.
101///
102/// `gate_positions` maps each gate-local qubit to its position in the sorted
103/// combined support used by [`matrix_commutation`].  Bit extraction follows the
104/// gate-local matrix convention: earlier qubits correspond to more
105/// significant bits in the local basis index.
106fn expand_unitary(
107    matrix: Cow<'_, Array2<Complex64>>,
108    gate_positions: &[usize],
109    total_qubits: usize,
110) -> Option<Array2<Complex64>> {
111    let dim = 1usize.checked_shl(total_qubits as u32)?;
112
113    if gate_positions.is_empty() {
114        if matrix.nrows() != 1 || matrix.ncols() != 1 {
115            return None;
116        }
117        let mut expanded = Array2::<Complex64>::eye(dim);
118        let scalar = matrix[[0, 0]];
119        expanded.iter_mut().for_each(|value| *value *= scalar);
120        return Some(expanded);
121    }
122
123    let expected = 1usize.checked_shl(gate_positions.len() as u32)?;
124    if matrix.nrows() != expected || matrix.ncols() != expected {
125        return None;
126    }
127    if gate_positions
128        .iter()
129        .any(|&position| position >= total_qubits)
130    {
131        return None;
132    }
133
134    let mut xor_mask = dim - 1;
135    for &position in gate_positions {
136        xor_mask ^= 1usize << (total_qubits - 1 - position);
137    }
138
139    let mut mapped_indices = vec![0usize; dim];
140    for (index, mapped) in mapped_indices.iter_mut().enumerate() {
141        for (gate_index, &position) in gate_positions.iter().enumerate() {
142            let bit = total_qubits - 1 - position;
143            if (index & (1usize << bit)) != 0 {
144                *mapped |= 1usize << (gate_positions.len() - 1 - gate_index);
145            }
146        }
147    }
148
149    let mut expanded = Array2::<Complex64>::zeros((dim, dim));
150    for row in 0..dim {
151        for col in 0..dim {
152            if (row & xor_mask) == (col & xor_mask) {
153                expanded[[row, col]] = matrix[[mapped_indices[row], mapped_indices[col]]];
154            }
155        }
156    }
157    Some(expanded)
158}
159
160/// Returns the global phase relating two matrices, if they differ only by one.
161///
162/// The first non-zero entry establishes the candidate unit complex ratio; every
163/// entry is then checked against that ratio with matrix tolerance.  Mismatched
164/// zero structure or a non-unit candidate ratio means the matrices are not equal
165/// up to global phase.
166fn phase_between(lhs: &Array2<Complex64>, rhs: &Array2<Complex64>) -> Option<f64> {
167    if lhs.dim() != rhs.dim() {
168        return None;
169    }
170
171    let mut ratio = None;
172    for (&lhs_value, &rhs_value) in lhs.iter().zip(rhs.iter()) {
173        let lhs_zero = lhs_value.norm() <= NUMERIC_ZERO_TOLERANCE;
174        let rhs_zero = rhs_value.norm() <= NUMERIC_ZERO_TOLERANCE;
175        if lhs_zero != rhs_zero {
176            return None;
177        }
178        if lhs_zero {
179            continue;
180        }
181
182        let candidate = lhs_value / rhs_value;
183        let norm = candidate.norm();
184        if !norm.is_finite() || (norm - 1.0).abs() > UNIT_PHASE_NORM_TOLERANCE {
185            return None;
186        }
187        ratio = Some(candidate);
188        break;
189    }
190
191    let ratio = ratio.unwrap_or_else(|| Complex64::new(1.0, 0.0));
192    for (&lhs_value, &rhs_value) in lhs.iter().zip(rhs.iter()) {
193        if (lhs_value - ratio * rhs_value).norm() > MATRIX_TOLERANCE {
194            return None;
195        }
196    }
197
198    Some(ratio.im.atan2(ratio.re))
199}
cqlib_core/compile/commutation/matrix.rs

cqlib_core/compile/commutation/
matrix.rs
