Skip to content

Commit 4bef428

Browse files
jakelishmanCryoris
jakelishman
and
Cryoris
authored
Add SparseObservable.compose (Qiskit#13766)
* Add `SparseObservable.compose` This adds matrix multiplication between `SparseObservable`s (or one `SparseObservable` and anything that can be coerced to an observable), using the same signature as similar methods in `quantum_info`. The current implementation attempts to be relatively memory efficient, though there are specialisations that may be implemented if we wanted to be faster: - if all the contained terms are Paulis, so each pairwise multiplication produces either zero or one alphabet letters, we currently have a redundant copy-in and copy-out of the Cartesian-product generator that effectively doubles the overhead. We could check whether each term contains only Paulis once, and then avoid that overhead in those cases. - alternatively, the Cartesian-product iterator could be taught to use the buffers of the output `SparseObservable` directly as its output space, rather than requiring copy-out of its buffers, which would remove the overhead from all calls, not just the single-Pauli ones, but requires a bit more bookkeeping, since the references would keep need to be updated. - the case of `qargs` being not `None` and `front=True` currently involves an extra copy step to simplify the logic. If this is an important case, we could add a right-matmul specialisation of `compose_map`. This commit currently does not contain tests, because we currently don't have suitable methods to test the validity of the output without coupling to the particular choice of factorisation of the matrix multiplication. We need either a way to convert to a Pauli-only representation (with the exponential explosion that entails) or to a matrix; either of these can be uniquely canonicalised. * Make `lookup::matmul` non-`const` Before Rust 1.83, `const` functions can't borrow `static`s, even immutably and without interior mutability.o There's no need for `matmul` to be `const`; it's always inlined anyway, and `matmul_generate` can be used in `const` contexts. * Add `BitTerm.label` property * Add tests of `SparseObservable.compose` * Fix typographical details Co-authored-by: Julien Gacon <[email protected]> --------- Co-authored-by: Julien Gacon <[email protected]>
1 parent 3ae0d14 commit 4bef428

File tree

4 files changed

+981
-0
lines changed

4 files changed

+981
-0
lines changed

‎crates/accelerate/src/sparse_observable/lookup.rs‎

Lines changed: 294 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,294 @@
1+
// This code is part of Qiskit.
2+
//
3+
// (C) Copyright IBM 2025
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+
use num_complex::Complex64;
14+
// The lookup tables are pretty illegible if we have all the syntactic noise of `BitTerm`.
15+
use super::BitTerm::{self, *};
16+
17+
/// Short-hand alias for [Complex64::new] that retains its ability to be used in `const` contexts.
18+
const fn c64(re: f64, im: f64) -> Complex64 {
19+
Complex64::new(re, im)
20+
}
21+
22+
/// The allowable items of `BitTerm`. This is used by the lookup expansion; we need a const-safe
23+
/// way of iterating through all of the variants.
24+
static BIT_TERMS: [BitTerm; 9] = [X, Plus, Minus, Y, Right, Left, Z, Zero, One];
25+
26+
/// Create a full binary lookup table for all valid combinations of `BitTerm`.
27+
macro_rules! expand_lookup_binary {
28+
(static $table:ident: [$ty:ty], $generate:ident, $dummy:expr) => {
29+
static $table: [$ty; 256] = {
30+
let mut out = [const { $dummy }; 256];
31+
let mut i = 0;
32+
while i < 9 {
33+
let left = BIT_TERMS[i];
34+
let mut j = 0;
35+
while j < 9 {
36+
let right = BIT_TERMS[j];
37+
out[((left as usize) << 4) | (right as usize)] = $generate(left, right);
38+
j += 1;
39+
}
40+
i += 1;
41+
}
42+
out
43+
};
44+
};
45+
}
46+
47+
expand_lookup_binary!(
48+
static MATMUL: [Option<&'static [(Complex64, BitTerm)]>],
49+
matmul_generate,
50+
None
51+
);
52+
/// Calculate the matrix multiplication of two [BitTerm]s.
53+
///
54+
/// This is done by a static lookup table.
55+
///
56+
/// In the output, `None` means that the result was zero. Beyond that, the slice should be
57+
/// interpreted as a sum of the coefficient multiplied by the [BitTerm] pairs. An empty slice means
58+
/// the identity.
59+
#[inline(always)]
60+
pub fn matmul(left: BitTerm, right: BitTerm) -> Option<&'static [(Complex64, BitTerm)]> {
61+
// This can be `const` from Rust 1.83 - before that, `const` functions can't refer to `static`s.
62+
MATMUL[((left as usize) << 4) | (right as usize)]
63+
}
64+
65+
/// The `const` version of [matmul].
66+
///
67+
/// This is less efficient at runtime than [matmul] (which inlines to a couple of bitwise operations
68+
/// and a single offset load), but can be used in `const` contexts.
69+
const fn matmul_generate(left: BitTerm, right: BitTerm) -> Option<&'static [(Complex64, BitTerm)]> {
70+
// This typically gets compiled to a non-inlinable (because the code size is too big) two
71+
// separate sets of operate-and-jump instructions, but we use it in the `const` context to build
72+
// a static lookup table which can easily be inlined and is a couple of bitwise operations
73+
// followed by an offset load of the result.
74+
//
75+
// These rules were generated by a brute-force search over increasing-n-tuples of (coeff, term)
76+
// pairs, where `coeff` was restricted to magnitudes `(0.25, 0.5, 1.0)` at complex phases
77+
// `(1, i, -1, -i)`. Under those conditions, the rules are minimal in number of sum terms, but
78+
// don't have any further logic to "curate" the choice of output factorisation.
79+
match (left, right) {
80+
(X, X) => Some(const { &[] }),
81+
(X, Plus) => Some(const { &[(c64(1., 0.), Plus)] }),
82+
(X, Minus) => Some(const { &[(c64(-1., 0.), Minus)] }),
83+
(X, Y) => Some(const { &[(c64(0., 1.), Z)] }),
84+
(X, Right) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Z)] }),
85+
(X, Left) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Z)] }),
86+
(X, Z) => Some(const { &[(c64(0., -1.), Y)] }),
87+
(X, Zero) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Y)] }),
88+
(X, One) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Y)] }),
89+
(Plus, X) => Some(const { &[(c64(1., 0.), Plus)] }),
90+
(Plus, Plus) => Some(const { &[(c64(1., 0.), Plus)] }),
91+
(Plus, Minus) => None,
92+
(Plus, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), Z)] }),
93+
(Plus, Right) => Some(
94+
const {
95+
&[
96+
(c64(0.25, 0.), X),
97+
(c64(0.5, 0.), Right),
98+
(c64(0., 0.25), Z),
99+
]
100+
},
101+
),
102+
(Plus, Left) => Some(
103+
const {
104+
&[
105+
(c64(0.25, 0.), X),
106+
(c64(0.5, 0.), Left),
107+
(c64(0., -0.25), Z),
108+
]
109+
},
110+
),
111+
(Plus, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), Y)] }),
112+
(Plus, Zero) => Some(
113+
const {
114+
&[
115+
(c64(0.25, 0.), X),
116+
(c64(0.5, 0.), Zero),
117+
(c64(0., -0.25), Y),
118+
]
119+
},
120+
),
121+
(Plus, One) => {
122+
Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., 0.25), Y)] })
123+
}
124+
(Minus, X) => Some(const { &[(c64(-1., 0.), Minus)] }),
125+
(Minus, Plus) => None,
126+
(Minus, Minus) => Some(const { &[(c64(1., 0.), Minus)] }),
127+
(Minus, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), Z)] }),
128+
(Minus, Right) => Some(
129+
const {
130+
&[
131+
(c64(0.25, 0.), Y),
132+
(c64(0.5, 0.), Minus),
133+
(c64(0., -0.25), Z),
134+
]
135+
},
136+
),
137+
(Minus, Left) => Some(
138+
const {
139+
&[
140+
(c64(-0.25, 0.), X),
141+
(c64(0.5, 0.), Left),
142+
(c64(0., 0.25), Z),
143+
]
144+
},
145+
),
146+
(Minus, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), Y)] }),
147+
(Minus, Zero) => Some(
148+
const {
149+
&[
150+
(c64(0.25, 0.), Z),
151+
(c64(0.5, 0.), Minus),
152+
(c64(0., 0.25), Y),
153+
]
154+
},
155+
),
156+
(Minus, One) => Some(
157+
const {
158+
&[
159+
(c64(-0.25, 0.), X),
160+
(c64(0.5, 0.), One),
161+
(c64(0., -0.25), Y),
162+
]
163+
},
164+
),
165+
(Y, X) => Some(const { &[(c64(0., -1.), Z)] }),
166+
(Y, Plus) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), Z)] }),
167+
(Y, Minus) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), Z)] }),
168+
(Y, Y) => Some(const { &[] }),
169+
(Y, Right) => Some(const { &[(c64(1., 0.), Right)] }),
170+
(Y, Left) => Some(const { &[(c64(-1., 0.), Left)] }),
171+
(Y, Z) => Some(const { &[(c64(0., 1.), X)] }),
172+
(Y, Zero) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), X)] }),
173+
(Y, One) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), X)] }),
174+
(Right, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Z)] }),
175+
(Right, Plus) => Some(
176+
const {
177+
&[
178+
(c64(0.25, 0.), X),
179+
(c64(0.5, 0.), Right),
180+
(c64(0., -0.25), Z),
181+
]
182+
},
183+
),
184+
(Right, Minus) => Some(
185+
const {
186+
&[
187+
(c64(0.25, 0.), Y),
188+
(c64(0.5, 0.), Minus),
189+
(c64(0., 0.25), Z),
190+
]
191+
},
192+
),
193+
(Right, Y) => Some(const { &[(c64(1., 0.), Right)] }),
194+
(Right, Right) => Some(const { &[(c64(1., 0.), Right)] }),
195+
(Right, Left) => None,
196+
(Right, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), X)] }),
197+
(Right, Zero) => {
198+
Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), Zero), (c64(0., 0.25), X)] })
199+
}
200+
(Right, One) => {
201+
Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., -0.25), X)] })
202+
}
203+
(Left, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Z)] }),
204+
(Left, Plus) => {
205+
Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), Left), (c64(0., 0.25), Z)] })
206+
}
207+
(Left, Minus) => Some(
208+
const {
209+
&[
210+
(c64(-0.25, 0.), X),
211+
(c64(0.5, 0.), Left),
212+
(c64(0., -0.25), Z),
213+
]
214+
},
215+
),
216+
(Left, Y) => Some(const { &[(c64(-1., 0.), Left)] }),
217+
(Left, Right) => None,
218+
(Left, Left) => Some(const { &[(c64(1., 0.), Left)] }),
219+
(Left, Z) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), X)] }),
220+
(Left, Zero) => Some(
221+
const {
222+
&[
223+
(c64(0.25, 0.), Z),
224+
(c64(0.5, 0.), Left),
225+
(c64(0., -0.25), X),
226+
]
227+
},
228+
),
229+
(Left, One) => {
230+
Some(const { &[(c64(-0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., 0.25), X)] })
231+
}
232+
(Z, X) => Some(const { &[(c64(0., 1.), Y)] }),
233+
(Z, Plus) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), Y)] }),
234+
(Z, Minus) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), Y)] }),
235+
(Z, Y) => Some(const { &[(c64(0., -1.), X)] }),
236+
(Z, Right) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., -0.5), X)] }),
237+
(Z, Left) => Some(const { &[(c64(0.5, 0.), Z), (c64(0., 0.5), X)] }),
238+
(Z, Z) => Some(const { &[] }),
239+
(Z, Zero) => Some(const { &[(c64(1., 0.), Zero)] }),
240+
(Z, One) => Some(const { &[(c64(-1., 0.), One)] }),
241+
(Zero, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., 0.5), Y)] }),
242+
(Zero, Plus) => {
243+
Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), Zero), (c64(0., 0.25), Y)] })
244+
}
245+
(Zero, Minus) => Some(
246+
const {
247+
&[
248+
(c64(0.25, 0.), Z),
249+
(c64(0.5, 0.), Minus),
250+
(c64(0., -0.25), Y),
251+
]
252+
},
253+
),
254+
(Zero, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., -0.5), X)] }),
255+
(Zero, Right) => Some(
256+
const {
257+
&[
258+
(c64(0.25, 0.), Y),
259+
(c64(0.5, 0.), Zero),
260+
(c64(0., -0.25), X),
261+
]
262+
},
263+
),
264+
(Zero, Left) => {
265+
Some(const { &[(c64(0.25, 0.), Z), (c64(0.5, 0.), Left), (c64(0., 0.25), X)] })
266+
}
267+
(Zero, Z) => Some(const { &[(c64(1., 0.), Zero)] }),
268+
(Zero, Zero) => Some(const { &[(c64(1., 0.), Zero)] }),
269+
(Zero, One) => None,
270+
(One, X) => Some(const { &[(c64(0.5, 0.), X), (c64(0., -0.5), Y)] }),
271+
(One, Plus) => {
272+
Some(const { &[(c64(0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., -0.25), Y)] })
273+
}
274+
(One, Minus) => {
275+
Some(const { &[(c64(-0.25, 0.), X), (c64(0.5, 0.), One), (c64(0., 0.25), Y)] })
276+
}
277+
(One, Y) => Some(const { &[(c64(0.5, 0.), Y), (c64(0., 0.5), X)] }),
278+
(One, Right) => {
279+
Some(const { &[(c64(0.25, 0.), Y), (c64(0.5, 0.), One), (c64(0., 0.25), X)] })
280+
}
281+
(One, Left) => Some(
282+
const {
283+
&[
284+
(c64(-0.25, 0.), Y),
285+
(c64(0.5, 0.), One),
286+
(c64(0., -0.25), X),
287+
]
288+
},
289+
),
290+
(One, Z) => Some(const { &[(c64(-1., 0.), One)] }),
291+
(One, Zero) => None,
292+
(One, One) => Some(const { &[(c64(1., 0.), One)] }),
293+
}
294+
}

0 commit comments

Comments
 (0)