From ee497574fc037bc5c7234e67bb038f827404f9e4 Mon Sep 17 00:00:00 2001 From: Craig Gidney Date: Thu, 14 Mar 2024 01:32:30 -0700 Subject: [PATCH 1/5] baseline --- src/stim/stabilizers/conversions.perf.cc | 48 ++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/src/stim/stabilizers/conversions.perf.cc b/src/stim/stabilizers/conversions.perf.cc index 0461355a7..abc8bc6ef 100644 --- a/src/stim/stabilizers/conversions.perf.cc +++ b/src/stim/stabilizers/conversions.perf.cc @@ -67,3 +67,51 @@ BENCHMARK(independent_to_disjoint_xyz_errors) { std::cout << "data dependence"; } } + +BENCHMARK(stabilizers_to_tableau) { + std::vector> offsets{ + {1, 0}, + {-1, 0}, + {0, 1}, + {0, -1}, + {3, 6}, + {-6, 3}, + }; + size_t w = 24; + size_t h = 12; + + auto normalize = [&](std::complex c) -> std::complex { + return {fmodf(c.real() + w*10, w), fmodf(c.imag() + h*10, h)}; + }; + auto q2i = [&](std::complex c) -> size_t { + c = normalize(c); + return (int)c.real() / 2 + c.imag() * (w / 2); + }; + + std::vector> stabilizers; + for (size_t x = 0; x < w; x++) { + for (size_t y = x % 2; y < h; y += 2) { + std::complex s{x % 2 ? -1.0f : +1.0f, 0.0f}; + std::complex c{(float)x, (float)y}; + stim::PauliString<128> ps(w * h / 2); + for (const auto &offset : offsets) { + size_t i = q2i(c + offset * s); + if (x % 2 == 0) { + ps.xs[i] = 1; + } else { + ps.zs[i] = 1; + } + } + stabilizers.push_back(ps); + } + } + + size_t dep = 0; + benchmark_go([&]() { + Tableau<128> t = stabilizers_to_tableau(stabilizers, true, true, false); + dep += t.xs[0].zs[0]; + }).goal_millis(140); + if (dep == 99999999) { + std::cout << "data dependence"; + } +} \ No newline at end of file From a9e8f3b24d9bdca00fc28761f30b10f4a6ee524d Mon Sep 17 00:00:00 2001 From: Craig Gidney Date: Thu, 14 Mar 2024 01:48:01 -0700 Subject: [PATCH 2/5] Make `stim.Tableau.from_stabilizers` faster - Measured 20x faster (140ms -> 6ms) on a 144 qubit case - Measured 100x faster (15s -> 0.13s) on a 432 qubit case --- src/stim/stabilizers/conversions.h | 4 +- src/stim/stabilizers/conversions.inl | 79 ++++++++++++------------ src/stim/stabilizers/conversions.perf.cc | 4 +- 3 files changed, 44 insertions(+), 43 deletions(-) diff --git a/src/stim/stabilizers/conversions.h b/src/stim/stabilizers/conversions.h index 4349693fd..37166d70d 100644 --- a/src/stim/stabilizers/conversions.h +++ b/src/stim/stabilizers/conversions.h @@ -83,11 +83,13 @@ Circuit stabilizer_state_vector_to_circuit( /// ignore_noise: If the circuit contains noise channels, ignore them instead of raising an exception. /// ignore_measurement: If the circuit contains measurements, ignore them instead of raising an exception. /// ignore_reset: If the circuit contains resets, ignore them instead of raising an exception. +/// inverse: The last step of the implementation is to invert the tableau. Setting this argument +/// to true will skip this inversion, saving time but returning the inverse tableau. /// /// Returns: /// A tableau encoding the given circuit's Clifford operation. template -Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset); +Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset, bool inverse = false); /// Simulates the given circuit and outputs a state vector. /// diff --git a/src/stim/stabilizers/conversions.inl b/src/stim/stabilizers/conversions.inl index b6a8edec3..7c47097b1 100644 --- a/src/stim/stabilizers/conversions.inl +++ b/src/stim/stabilizers/conversions.inl @@ -149,7 +149,7 @@ std::vector>> tableau_to_unitary(const Tableau -Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset) { +Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ignore_measurement, bool ignore_reset, bool inverse) { Tableau result(circuit.count_qubits()); TableauSimulator sim(std::mt19937_64(0), circuit.count_qubits()); @@ -185,7 +185,10 @@ Tableau circuit_to_tableau(const Circuit &circuit, bool ignore_noise, bool ig } }); - return sim.inv_state.inverse(); + if (!inverse) { + return sim.inv_state.inverse(); + } + return sim.inv_state; } template @@ -556,7 +559,7 @@ Tableau stabilizers_to_tableau( } for (size_t k1 = 0; k1 < stabilizers.size(); k1++) { - for (size_t k2 = 0; k2 < stabilizers.size(); k2++) { + for (size_t k2 = k1 + 1; k2 < stabilizers.size(); k2++) { if (!stabilizers[k1].ref().commutes(stabilizers[k2])) { std::stringstream ss; ss << "Some of the given stabilizers anticommute.\n"; @@ -568,44 +571,31 @@ Tableau stabilizers_to_tableau( } } } - Tableau inverted(num_qubits); - - PauliString cur(num_qubits); - std::vector targets; - while (targets.size() < num_qubits) { - targets.push_back(targets.size()); - } - auto overwrite_cur_apply_recorded = [&](const PauliString &e) { - PauliStringRef cur_ref = cur.ref(); - cur.xs.clear(); - cur.zs.clear(); - cur.xs.word_range_ref(0, e.xs.num_simd_words) = e.xs; - cur.zs.word_range_ref(0, e.xs.num_simd_words) = e.zs; - cur.sign = e.sign; - inverted.apply_within(cur_ref, targets); - }; + Circuit elimination_instructions; + PauliString buf(num_qubits); size_t used = 0; for (const auto &e : stabilizers) { - overwrite_cur_apply_recorded(e); + buf = e; + buf.ref().do_circuit(elimination_instructions); // Find a non-identity term in the Pauli string past the region used by other stabilizers. size_t pivot; for (pivot = used; pivot < num_qubits; pivot++) { - if (cur.xs[pivot] || cur.zs[pivot]) { + if (buf.xs[pivot] || buf.zs[pivot]) { break; } } // Check for incompatible / redundant stabilizers. if (pivot == num_qubits) { - if (cur.xs.not_zero()) { + if (buf.xs.not_zero()) { throw std::invalid_argument("Some of the given stabilizers anticommute."); } - if (cur.sign) { + if (buf.sign) { throw std::invalid_argument("Some of the given stabilizers contradict each other."); } - if (!allow_redundant && cur.zs.not_zero()) { + if (!allow_redundant && buf.zs.not_zero()) { throw std::invalid_argument( "Didn't specify allow_redundant=True but one of the given stabilizers is a product of the others. " "To allow redundant stabilizers, pass the argument allow_redundant=True."); @@ -614,32 +604,36 @@ Tableau stabilizers_to_tableau( } // Change pivot basis to the Z axis. - if (cur.xs[pivot]) { - std::string name = cur.zs[pivot] ? "H_YZ" : "H_XZ"; - inverted.inplace_scatter_append(GATE_DATA.at(name).tableau(), {pivot}); + if (buf.xs[pivot]) { + GateType g = buf.zs[pivot] ? GateType::H_YZ : GateType::H; + GateTarget t = GateTarget::qubit(pivot); + CircuitInstruction instruction{g, {}, &t}; + elimination_instructions.safe_append(instruction); + buf.ref().do_instruction(instruction); } // Cancel other terms in Pauli string. for (size_t q = 0; q < num_qubits; q++) { - int p = cur.xs[q] + cur.zs[q] * 2; + int p = buf.xs[q] + buf.zs[q] * 2; if (p && q != pivot) { - inverted.inplace_scatter_append( - GATE_DATA.at(p == 1 ? "XCX" - : p == 2 ? "XCZ" - : "XCY") - .tableau(), - {pivot, q}); + std::array targets{GateTarget::qubit(pivot), GateTarget::qubit(q)}; + CircuitInstruction instruction{p == 1 ? GateType::XCX : p == 2 ? GateType::XCZ : GateType::XCY, {}, targets}; + elimination_instructions.safe_append(instruction); + buf.ref().do_instruction(instruction); } } // Move pivot to diagonal. if (pivot != used) { - inverted.inplace_scatter_append(GATE_DATA.at("SWAP").tableau(), {pivot, used}); + std::array targets{GateTarget::qubit(pivot), GateTarget::qubit(used)}; + CircuitInstruction instruction{GateType::SWAP, {}, targets}; + elimination_instructions.safe_append(instruction); } // Fix sign. - overwrite_cur_apply_recorded(e); - if (cur.sign) { - inverted.inplace_scatter_append(GATE_DATA.at("X").tableau(), {used}); + if (buf.sign) { + GateTarget t = GateTarget::qubit(used); + CircuitInstruction instruction{GateType::X, {}, &t}; + elimination_instructions.safe_append(instruction); } used++; @@ -653,10 +647,15 @@ Tableau stabilizers_to_tableau( } } + if (num_qubits > 0) { + GateTarget t = GateTarget::qubit(num_qubits - 1); + elimination_instructions.safe_append(CircuitInstruction{GateType::I, {}, &t}); + } + if (invert) { - return inverted; + return circuit_to_tableau(elimination_instructions.inverse(), false, false, false, true); } - return inverted.inverse(); + return circuit_to_tableau(elimination_instructions, false, false, false, true); } } // namespace stim diff --git a/src/stim/stabilizers/conversions.perf.cc b/src/stim/stabilizers/conversions.perf.cc index abc8bc6ef..2ed27a787 100644 --- a/src/stim/stabilizers/conversions.perf.cc +++ b/src/stim/stabilizers/conversions.perf.cc @@ -110,8 +110,8 @@ BENCHMARK(stabilizers_to_tableau) { benchmark_go([&]() { Tableau<128> t = stabilizers_to_tableau(stabilizers, true, true, false); dep += t.xs[0].zs[0]; - }).goal_millis(140); + }).goal_millis(6); if (dep == 99999999) { std::cout << "data dependence"; } -} \ No newline at end of file +} From 67a7f3695c6dc32081be01c0184a597260087262 Mon Sep 17 00:00:00 2001 From: Craig Gidney Date: Thu, 14 Mar 2024 01:52:04 -0700 Subject: [PATCH 3/5] 64 --- src/stim/stabilizers/conversions.perf.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/stim/stabilizers/conversions.perf.cc b/src/stim/stabilizers/conversions.perf.cc index 2ed27a787..50ad12c50 100644 --- a/src/stim/stabilizers/conversions.perf.cc +++ b/src/stim/stabilizers/conversions.perf.cc @@ -88,12 +88,12 @@ BENCHMARK(stabilizers_to_tableau) { return (int)c.real() / 2 + c.imag() * (w / 2); }; - std::vector> stabilizers; + std::vector> stabilizers; for (size_t x = 0; x < w; x++) { for (size_t y = x % 2; y < h; y += 2) { std::complex s{x % 2 ? -1.0f : +1.0f, 0.0f}; std::complex c{(float)x, (float)y}; - stim::PauliString<128> ps(w * h / 2); + stim::PauliString<64> ps(w * h / 2); for (const auto &offset : offsets) { size_t i = q2i(c + offset * s); if (x % 2 == 0) { @@ -108,9 +108,9 @@ BENCHMARK(stabilizers_to_tableau) { size_t dep = 0; benchmark_go([&]() { - Tableau<128> t = stabilizers_to_tableau(stabilizers, true, true, false); + Tableau<64> t = stabilizers_to_tableau(stabilizers, true, true, false); dep += t.xs[0].zs[0]; - }).goal_millis(6); + }).goal_millis(5); if (dep == 99999999) { std::cout << "data dependence"; } From 7790744aae014c46c9d426dc847af5fd0bfa3bfa Mon Sep 17 00:00:00 2001 From: Craig Gidney Date: Thu, 14 Mar 2024 01:55:40 -0700 Subject: [PATCH 4/5] Size fix --- src/stim/stabilizers/conversions.inl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stim/stabilizers/conversions.inl b/src/stim/stabilizers/conversions.inl index 7c47097b1..0983ca3d2 100644 --- a/src/stim/stabilizers/conversions.inl +++ b/src/stim/stabilizers/conversions.inl @@ -648,8 +648,10 @@ Tableau stabilizers_to_tableau( } if (num_qubits > 0) { + // Force size of resulting tableau to be correct. GateTarget t = GateTarget::qubit(num_qubits - 1); - elimination_instructions.safe_append(CircuitInstruction{GateType::I, {}, &t}); + elimination_instructions.safe_append(CircuitInstruction{GateType::X, {}, &t}); + elimination_instructions.safe_append(CircuitInstruction{GateType::X, {}, &t}); } if (invert) { From cc7029d6076914c651ef5698e03303349caa936b Mon Sep 17 00:00:00 2001 From: Craig Gidney Date: Thu, 14 Mar 2024 02:12:42 -0700 Subject: [PATCH 5/5] size fix --- src/stim/stabilizers/conversions.inl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/stim/stabilizers/conversions.inl b/src/stim/stabilizers/conversions.inl index 0983ca3d2..6ab5f2b17 100644 --- a/src/stim/stabilizers/conversions.inl +++ b/src/stim/stabilizers/conversions.inl @@ -576,7 +576,15 @@ Tableau stabilizers_to_tableau( size_t used = 0; for (const auto &e : stabilizers) { - buf = e; + if (e.num_qubits == num_qubits) { + buf = e; + } else { + buf.xs.clear(); + buf.zs.clear(); + memcpy(buf.xs.u8, e.xs.u8, e.xs.num_u8_padded()); + memcpy(buf.zs.u8, e.zs.u8, e.zs.num_u8_padded()); + buf.sign = e.sign; + } buf.ref().do_circuit(elimination_instructions); // Find a non-identity term in the Pauli string past the region used by other stabilizers.