From 516c360d6c9ba36cd799fb893bfe90e7da45c897 Mon Sep 17 00:00:00 2001 From: omegaxyv Date: Wed, 2 Sep 2026 15:58:46 +0800 Subject: [PATCH] perf(QSEncode): optimize Walsh transform and top-k selection --- .../pyqpanda_alg/QSEncode/QSEncode.py | 94 +++++++++- .../Test_QSEncode_numeric_pipeline.py | 165 ++++++++++++++++++ 2 files changed, 255 insertions(+), 4 deletions(-) create mode 100644 test/QAlgBase/Test_QSEncode_numeric_pipeline.py diff --git a/pyqpanda-algorithm/pyqpanda_alg/QSEncode/QSEncode.py b/pyqpanda-algorithm/pyqpanda_alg/QSEncode/QSEncode.py index d1d4fde7..628cd542 100644 --- a/pyqpanda-algorithm/pyqpanda_alg/QSEncode/QSEncode.py +++ b/pyqpanda-algorithm/pyqpanda_alg/QSEncode/QSEncode.py @@ -14,10 +14,55 @@ import numpy as np from pyqpanda3.core import CPUQVM, QCircuit, QProg, Encode, H from scipy.fft import fft -from sympy import fwht from .. plugin import * + +def _fast_walsh_hadamard_transform(values): + """Return the unnormalized Walsh-Hadamard transform of a 1-D array. + + The input is zero-padded to the next power-of-two length, matching the + behavior of :func:`sympy.fwht`. The butterfly is evaluated with NumPy + ufuncs, so numeric values stay in a native numeric dtype instead of being + converted to symbolic Python objects. + + The input is never modified. Runtime is ``O(N log N)`` and auxiliary + memory is ``O(N)``. + """ + array = np.asarray(values) + if array.ndim != 1: + raise ValueError('values must be a 1D array') + if not np.issubdtype(array.dtype, np.number): + raise TypeError('values must contain numeric data') + + size = array.size + dtype = np.result_type(array.dtype, np.float64) + if size == 0: + return np.empty(0, dtype=dtype) + + padded_size = 1 << (size - 1).bit_length() + output = np.zeros(padded_size, dtype=dtype) + output[:size] = array + if padded_size == 1: + return output + + # One reusable half-size buffer avoids allocating a Python/SymPy object + # for every coefficient and avoids a fresh temporary at every stage. + scratch = np.empty(padded_size // 2, dtype=dtype) + step = 1 + while step < padded_size: + blocks = output.reshape(-1, 2 * step) + left = blocks[:, :step] + right = blocks[:, step:] + saved_left = scratch.reshape(-1, step) + np.copyto(saved_left, left) + np.add(saved_left, right, out=left) + np.subtract(saved_left, right, out=right) + step *= 2 + + return output + + class QSpare_Code: """ @@ -146,15 +191,50 @@ def select_top_n_complex_numbers(self, arr, n): ------- np.ndarray A new array with only top-n magnitudes retained, rest set to zero. + + Notes + ----- + Selection uses ``numpy.argpartition`` and therefore takes average + ``O(N)`` time. If equal magnitudes cross the selection boundary, the + larger original indices are retained to make the result deterministic. """ if type(n) != int or n <= 0: raise ValueError('n must > 0 and with class int') + arr = np.asarray(arr) + if arr.ndim != 1: + raise ValueError('arr must be a 1D array') + if arr.size == 0: + return arr.copy() + magnitudes = np.abs(arr) - top_n_indices = np.argsort(magnitudes)[-n:] + if not np.all(np.isfinite(magnitudes)): + raise ValueError('arr must contain only finite values') + if n >= arr.size: + return arr.copy() + + # argpartition finds the top-n boundary in average O(N), whereas a + # complete argsort performs O(N log N) work. Values tied at the + # boundary are resolved by index so results are deterministic. + partition_at = arr.size - n + partitioned = np.argpartition(magnitudes, partition_at) + threshold = magnitudes[partitioned[partition_at]] + top_n_indices = partitioned[partition_at:] + + # The partition result is already complete for the overwhelmingly + # common unique-boundary case. Only scan and rebuild the selection + # when equal magnitudes straddle the top-n boundary. + selected_ties = np.count_nonzero(magnitudes[top_n_indices] == threshold) + all_ties = np.count_nonzero(magnitudes == threshold) + if selected_ties != all_ties: + top_n_indices = np.flatnonzero(magnitudes > threshold) + slots_left = n - top_n_indices.size + tied_indices = np.flatnonzero(magnitudes == threshold) + top_n_indices = np.concatenate((top_n_indices, tied_indices[-slots_left:])) + result = np.zeros_like(arr, dtype=arr.dtype) result[top_n_indices] = arr[top_n_indices] - return np.array(result) + return result def Transform(self, amp): """ @@ -169,9 +249,15 @@ def Transform(self, amp): ------- np.ndarray Transformed amplitude vector in the selected basis. + + Notes + ----- + Walsh mode uses a native NumPy fast Walsh-Hadamard butterfly with + ``O(N log N)`` time and ``O(N)`` memory. """ if self.mode == 'walsh': - transform = np.array(fwht(amp)) / np.sqrt(2 ** self.qubits_num) + transform = _fast_walsh_hadamard_transform(amp) + transform /= np.sqrt(2 ** self.qubits_num) elif self.mode == 'fourier': transform = fft(amp) else: diff --git a/test/QAlgBase/Test_QSEncode_numeric_pipeline.py b/test/QAlgBase/Test_QSEncode_numeric_pipeline.py new file mode 100644 index 00000000..73ece5e8 --- /dev/null +++ b/test/QAlgBase/Test_QSEncode_numeric_pipeline.py @@ -0,0 +1,165 @@ +import numpy as np +import pytest + +from pyqpanda_alg.QSEncode.QSEncode import ( + QSpare_Code, + _fast_walsh_hadamard_transform, +) + + +def _hadamard_matrix(size): + """Build a small independent Sylvester Hadamard matrix.""" + matrix = np.array([[1.0]]) + while matrix.shape[0] < size: + matrix = np.block([[matrix, matrix], [matrix, -matrix]]) + return matrix + + +def _dense_fwht(values): + """Reference transform using an explicit matrix, not a butterfly.""" + values = np.asarray(values) + if values.size == 0: + return np.empty(0, dtype=np.result_type(values.dtype, np.float64)) + padded_size = 1 << (values.size - 1).bit_length() + padded = np.zeros(padded_size, dtype=np.result_type(values.dtype, np.float64)) + padded[:values.size] = values + return _hadamard_matrix(padded_size) @ padded + + +def _dense_sparse_walsh_probabilities(probabilities, cut): + """Independent end-to-end oracle for the Walsh sparse encoder.""" + probabilities = np.asarray(probabilities, dtype=float) + amplitudes = np.sqrt(probabilities / probabilities.sum()) + size = amplitudes.size + hadamard = _hadamard_matrix(size) / np.sqrt(size) + coefficients = hadamard @ amplitudes + + order = np.argsort(np.abs(coefficients), kind='stable') + sparse = np.zeros_like(coefficients) + sparse[order[-cut:]] = coefficients[order[-cut:]] + sparse /= np.linalg.norm(sparse) + return np.abs(hadamard @ sparse) ** 2 + + +@pytest.mark.parametrize( + 'values', + [ + np.array([], dtype=float), + np.array([3.25]), + np.array([1.0, -2.0, 4.0]), + np.array([1.0, 2.0, 3.0, 4.0]), + np.array([1.0 + 2.0j, -3.0j, 0.5 - 0.25j]), + ], +) +def test_fast_walsh_transform_matches_dense_oracle(values): + actual = _fast_walsh_hadamard_transform(values) + expected = _dense_fwht(values) + + assert actual.dtype != object + np.testing.assert_allclose(actual, expected, rtol=0.0, atol=1e-12) + + +def test_fast_walsh_transform_matches_dense_oracle_for_lengths_1_to_64(): + rng = np.random.default_rng(20260901) + + for size in range(1, 65): + values = rng.normal(size=size) + if size % 2 == 0: + values = values + 1j * rng.normal(size=size) + actual = _fast_walsh_hadamard_transform(values) + expected = _dense_fwht(values) + np.testing.assert_allclose(actual, expected, rtol=1e-13, atol=1e-12) + + +@pytest.mark.parametrize('qubits', range(13)) +def test_normalized_walsh_transform_preserves_norm_and_is_self_inverse(qubits): + size = 1 << qubits + values = np.random.default_rng(qubits).normal(size=size) + transformed = _fast_walsh_hadamard_transform(values) / np.sqrt(size) + restored = _fast_walsh_hadamard_transform(transformed) / np.sqrt(size) + + np.testing.assert_allclose(np.linalg.norm(transformed), np.linalg.norm(values), rtol=1e-13) + np.testing.assert_allclose(restored, values, rtol=1e-13, atol=1e-13) + + +def test_fast_walsh_transform_does_not_modify_input_and_promotes_numeric_dtype(): + values = np.array([1, 2, 3], dtype=np.int16) + original = values.copy() + + result = _fast_walsh_hadamard_transform(values) + + np.testing.assert_array_equal(values, original) + assert result.dtype == np.float64 + + +def test_fast_walsh_transform_rejects_non_numeric_or_non_vector_input(): + with pytest.raises(ValueError, match='1D'): + _fast_walsh_hadamard_transform(np.ones((2, 2))) + with pytest.raises(TypeError, match='numeric'): + _fast_walsh_hadamard_transform(['a', 'b']) + + +@pytest.mark.parametrize( + ('size', 'cut'), + [(8, 1), (8, 7), (257, 17), (4096, 63)], +) +def test_linear_top_k_matches_full_stable_sort_for_unique_magnitudes(size, cut): + rng = np.random.default_rng(size + cut) + phases = np.exp(2j * np.pi * rng.random(size)) + values = np.arange(1, size + 1, dtype=float) * phases + encoder = QSpare_Code([0.5, 0.5]) + + actual = encoder.select_top_n_complex_numbers(values, cut) + expected = np.zeros_like(values) + expected_indices = np.argsort(np.abs(values), kind='stable')[-cut:] + expected[expected_indices] = values[expected_indices] + + np.testing.assert_array_equal(actual, expected) + + +def test_linear_top_k_has_deterministic_boundary_ties(): + values = np.array([-3.0, 3.0, 2.0, -2.0, 1.0]) + encoder = QSpare_Code([0.5, 0.5]) + + actual = encoder.select_top_n_complex_numbers(values, 3) + + np.testing.assert_array_equal(actual, [-3.0, 3.0, 0.0, -2.0, 0.0]) + + +def test_linear_top_k_keeps_input_unchanged_and_fast_paths_full_selection(): + values = np.array([1.0 + 2.0j, -4.0j, 3.0]) + original = values.copy() + encoder = QSpare_Code([0.5, 0.5]) + + result = encoder.select_top_n_complex_numbers(values, values.size) + + np.testing.assert_array_equal(values, original) + np.testing.assert_array_equal(result, original) + assert result is not values + + +def test_linear_top_k_rejects_invalid_shape_and_non_finite_values(): + encoder = QSpare_Code([0.5, 0.5]) + with pytest.raises(ValueError, match='1D'): + encoder.select_top_n_complex_numbers(np.ones((2, 2)), 1) + with pytest.raises(ValueError, match='finite'): + encoder.select_top_n_complex_numbers(np.array([1.0, np.nan]), 1) + with pytest.raises(ValueError, match='finite'): + encoder.select_top_n_complex_numbers(np.array([1.0, np.inf]), 2) + + +@pytest.mark.parametrize( + ('probabilities', 'cut'), + [ + ([0.50, 0.30, 0.15, 0.05], 2), + ([0.31, 0.19, 0.16, 0.12, 0.09, 0.06, 0.04, 0.03], 3), + ], +) +def test_walsh_quantum_result_matches_independent_dense_pipeline(probabilities, cut): + encoder = QSpare_Code(list(map(float, probabilities)), mode='walsh', cut_length=cut) + + actual = encoder.Quantum_Res() + expected = _dense_sparse_walsh_probabilities(probabilities, cut) + + np.testing.assert_allclose(actual, expected, rtol=1e-9, atol=1e-9) + np.testing.assert_allclose(np.sum(actual), 1.0, rtol=0.0, atol=1e-12)