diff --git a/pyqpanda-algorithm/example/QAlgBase/testeg_VQLS.py b/pyqpanda-algorithm/example/QAlgBase/testeg_VQLS.py new file mode 100644 index 0000000..fd6472c --- /dev/null +++ b/pyqpanda-algorithm/example/QAlgBase/testeg_VQLS.py @@ -0,0 +1,31 @@ +import numpy as np +from pyqpanda_alg.VQLS import VQLS + +""" +VQLS solves a linear system Ax = b for a Hermitian matrix A variationally. +A hardware-efficient RY/CZ ansatz is optimized with scipy to minimize a global +cost built from Pauli expectation values and overlap amplitudes between the +ansatz state and b. VQLS(A, b).run() returns the normalized solution vector. +final_cost near 0 indicates convergence. ansatz_layers, seed, and optimizer are +tunable. Unlike HHL, VQLS is NISQ-friendly and needs no ancilla or clock register. +""" + + +def main(): + A = np.array([[1.0, -1.0 / 3], + [-1.0 / 3, 1.0]]) + b = np.array([0.0, 1.0]) + + solver = VQLS(A, b, seed=0) + x_q = solver.run() + + x_c = np.linalg.solve(A, b) + x_c = x_c / np.linalg.norm(x_c) + + print("quantum solution :", np.round(x_q, 4)) + print("classical solution:", np.round(x_c, 4)) + print("final cost:", round(solver.final_cost, 4)) + + +if __name__ == "__main__": + main() diff --git a/pyqpanda-algorithm/pyqpanda_alg/VQLS/VQLS.py b/pyqpanda-algorithm/pyqpanda_alg/VQLS/VQLS.py new file mode 100644 index 0000000..10f3f48 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/VQLS/VQLS.py @@ -0,0 +1,231 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import warnings +import itertools + +import numpy as np +from scipy.optimize import minimize + +from pyqpanda3.core import CPUQVM, QProg, QCircuit, RY, CZ, X, Y, Z, Encode, expval_pauli_operator +from pyqpanda3.hamiltonian import PauliOperator + + +class VQLS: + """ + This class provides a framework for the Variational Quantum Linear Solver + (VQLS) algorithm for solving linear systems of equations Ax = b [1]. + + For a Hermitian matrix A it decomposes A into a weighted sum of Pauli + strings and trains a hardware-efficient RY ansatz to minimize a global cost + built from the Pauli expectation and the overlap amplitudes + . The optimized ansatz prepares a state proportional to the + solution x. Unlike phase-estimation solvers it is clock-register-free and + ancilla-free, making it suited to NISQ hardware; the real RY ansatz limits + it to solutions expressible with real amplitudes. + + Parameters + matrix : ``numpy.ndarray``\n + The Hermitian coefficient matrix A of shape (N, N), N = 2^n. + vector : ``numpy.ndarray``, ``list``\n + The right-hand-side vector b of length N. + ansatz_layers : ``int``\n + The number of entangling layers in the hardware-efficient ansatz. + optimizer : ``str``\n + The classical optimizer, any scipy.optimize.minimize method name. + maxiter : ``int``\n + The maximum number of classical optimizer iterations. + seed : ``int``\n + The random seed for the initial ansatz parameters. + + References + [1] Bravo-Prieto C, LaRose R, Cerezo M, et al. Variational quantum linear + solver[J]. Quantum, 2023, 7: 1188. + + >>> from pyqpanda_alg.VQLS import VQLS + >>> import numpy as np + >>> A = np.array([[1.0, -1.0 / 3], [-1.0 / 3, 1.0]]) + >>> b = np.array([0.0, 1.0]) + >>> x = VQLS(A, b, seed=0).run() + >>> print(np.round(x, 4)) + [0.3162 0.9487] + + """ + + def __init__(self, matrix=None, + vector=None, + ansatz_layers: int = 3, + optimizer: str = 'SLSQP', + maxiter: int = 200, + seed: int = 0, + machine_type: str = 'CPU'): + + A = np.array(matrix, dtype=complex) + b = np.array(vector, dtype=complex).flatten() + + if A.ndim != 2 or A.shape[0] != A.shape[1]: + raise ValueError('matrix must be a square 2D array') + if not np.allclose(A, A.conj().T): + raise ValueError('matrix must be Hermitian') + N = A.shape[0] + self.n_qubits = int(np.log2(N)) + if 2 ** self.n_qubits != N: + raise ValueError('matrix dimension must be a power of 2') + if b.shape[0] != N: + raise ValueError('vector length must match matrix dimension') + if np.allclose(A, 0): + raise ValueError('matrix must be nonzero') + if np.linalg.norm(b) < 1e-12: + raise ValueError('vector must be nonzero') + + self.A = A + self.b = b + self.N = N + self.ansatz_layers = ansatz_layers + self.optimizer = optimizer + self.maxiter = maxiter + self.seed = seed + self.machine_type = machine_type + if machine_type == 'CPU': + self.machine = CPUQVM() + else: + raise NameError("Support 'CPU' only") + + n = self.n_qubits + paulis = {'I': np.eye(2, dtype=complex), + 'X': np.array([[0, 1], [1, 0]], dtype=complex), + 'Y': np.array([[0, -1j], [1j, 0]], dtype=complex), + 'Z': np.array([[1, 0], [0, -1]], dtype=complex)} + + pauli_dict = {} + for combo in itertools.product('IXYZ', repeat=n): + # little-endian kron: qubit n-1 is outermost factor, qubit 0 innermost + mat = np.array([[1.0 + 0j]]) + for q in range(n - 1, -1, -1): + mat = np.kron(mat, paulis[combo[q]]) + c_l = np.trace(mat.conj().T @ A) / N + if abs(c_l) > 1e-9: + tokens = [combo[q] + str(q) for q in range(n) if combo[q] != 'I'] + label = ' '.join(tokens) if tokens else 'I' + pauli_dict[label] = float(c_l.real) + + self.pauli_dict = pauli_dict + self.A_op = PauliOperator(pauli_dict) + # pyqpanda3 pauli Y is conjugated; harmless for the real ansatz, + # = when x is real + self.A2_op = self.A_op * self.A_op + + b_unit = b / np.linalg.norm(b) + if np.allclose(b_unit.imag, 0.0): + data = [float(v.real) for v in b_unit] + else: + data = [complex(v) for v in b_unit] + enc = Encode() + enc.amplitude_encode(list(range(n)), data) + self.u_b = enc.get_circuit() + + self.n_params = n * (ansatz_layers + 1) + + self.final_cost = None + self.opt_params = None + self.n_iters = None + self.cost_history = [] + self.converged = None + self.state = None + + def __del__(self): + pass + + def _ansatz(self, params): + cir = QCircuit() + n = self.n_qubits + k = 0 + for q in range(n): + cir << RY(q, params[k]) + k += 1 + for _ in range(self.ansatz_layers): + for q in range(n - 1): + cir << CZ(q, q + 1) + # ring entangler only for n > 2, CZ(0,1)+CZ(1,0) cancels at n=2 + if n > 2: + cir << CZ(n - 1, 0) + for q in range(n): + cir << RY(q, params[k]) + k += 1 + return cir + + def _pauli_gates(self, label): + cir = QCircuit() + if label == 'I': + return cir + for tok in label.split(): + p = tok[0] + q = int(tok[1:]) + if p == 'X': + cir << X(q) + elif p == 'Y': + cir << Y(q) + elif p == 'Z': + cir << Z(q) + return cir + + def _overlap(self, params, label): + n = self.n_qubits + prog = QProg(n) << self._ansatz(params) << self._pauli_gates(label) << self.u_b.dagger() + self.machine.run(prog, 1) + # + return complex(np.array(self.machine.result().get_state_vector())[0]) + + def _cost(self, params): + n = self.n_qubits + D = expval_pauli_operator(QProg(n) << self._ansatz(params), self.A2_op, 1) + if D < 1e-12: + cost = 1.0 + else: + amp = 0.0 + 0j + for label, c_l in self.pauli_dict.items(): + amp += c_l * self._overlap(params, label) + cost = 1.0 - abs(amp) ** 2 / D + self.cost_history.append(cost) + return cost + + def run(self): + """ + Run the variational quantum linear solver. + + Returns + x : ``numpy.ndarray``\n + The normalized solution vector x of Ax = b (up to a global phase). + + """ + n = self.n_qubits + x0 = np.random.default_rng(self.seed).uniform(-np.pi, np.pi, self.n_params) + options = {'maxiter': self.maxiter} + if self.optimizer == 'SLSQP': + options['ftol'] = 1e-10 + res = minimize(self._cost, x0, method=self.optimizer, options=options) + + self.final_cost = float(res.fun) + self.opt_params = res.x + self.n_iters = int(res.get('nit', res.nfev)) + self.converged = self.final_cost < 1e-2 + if not self.converged: + warnings.warn('VQLS did not converge: final cost %.3e' % self.final_cost) + + self.machine.run(QProg(n) << self._ansatz(res.x), 1) + x = np.array(self.machine.result().get_state_vector()) + + pivot = np.argmax(np.abs(x)) + x = x * np.exp(-1j * np.angle(x[pivot])) + x = np.real_if_close(x, tol=1e6) + self.state = x + return x diff --git a/pyqpanda-algorithm/pyqpanda_alg/VQLS/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/VQLS/__init__.py new file mode 100644 index 0000000..d6c4df0 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/VQLS/__init__.py @@ -0,0 +1,4 @@ + +from .VQLS import VQLS + +__all__ = ["VQLS"] diff --git a/pyqpanda-algorithm/pyqpanda_alg/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/__init__.py index 12d6808..d866647 100644 --- a/pyqpanda-algorithm/pyqpanda_alg/__init__.py +++ b/pyqpanda-algorithm/pyqpanda_alg/__init__.py @@ -51,4 +51,5 @@ from . import Grover from . import QmRMR from . import QSEncode +from . import VQLS diff --git a/test/QAlgBase/Test_VQLS_VQLS.py b/test/QAlgBase/Test_VQLS_VQLS.py new file mode 100644 index 0000000..90ccb1c --- /dev/null +++ b/test/QAlgBase/Test_VQLS_VQLS.py @@ -0,0 +1,66 @@ +import pytest +import numpy as np +from pyqpanda_alg.VQLS import VQLS +import warnings + + +def _classical_unit_solution(A, b): + x = np.linalg.solve(A, b) + x = x / np.linalg.norm(x) + pivot = np.argmax(np.abs(x)) + return x * np.sign(x[pivot]) + + +def _fidelity(x, y): + return float(np.abs(np.vdot(x, y)) ** 2) + + +class Test_VQLS_VQLS: + + def setup_method(self): + warnings.filterwarnings("ignore") + + def test_2x2_textbook(self): + A = np.array([[1.0, -1.0 / 3], [-1.0 / 3, 1.0]]) + b = np.array([0.0, 1.0]) + x_q = np.array(VQLS(A, b, seed=0).run(), dtype=complex) + x_c = _classical_unit_solution(A, b) + assert _fidelity(x_q, x_c) > 0.99, "VQLS solution should match classical solver" + + def test_2x2_diagonal(self): + A = np.diag([2.0, 1.0]) + b = np.array([1.0, 1.0]) + x_q = np.array(VQLS(A, b, seed=0).run(), dtype=complex) + x_c = _classical_unit_solution(A, b) + assert _fidelity(x_q, x_c) > 0.99 + + def test_4x4_diagonal(self): + A = np.diag([1.0, 2.0, 3.0, 4.0]).astype(float) + b = np.array([1.0, 1.0, 1.0, 1.0]) + x_q = np.array(VQLS(A, b, seed=0).run(), dtype=complex) + x_c = _classical_unit_solution(A, b) + assert _fidelity(x_q, x_c) > 0.99 + + def test_final_cost_converged(self): + A = np.diag([2.0, 1.0]) + b = np.array([1.0, 1.0]) + solver = VQLS(A, b, seed=0) + solver.run() + assert 0.0 <= solver.final_cost < 1e-2 + assert solver.converged + + def test_non_hermitian_raises(self): + A = np.array([[1.0, 2.0], [0.0, 1.0]]) # not Hermitian + b = np.array([1.0, 0.0]) + with pytest.raises(ValueError): + VQLS(A, b, seed=0) + + def test_non_power_of_two_raises(self): + A = np.eye(3) + b = np.array([1.0, 0.0, 0.0]) + with pytest.raises(ValueError): + VQLS(A, b, seed=0) + + +if __name__ == "__main__": + pytest.main([__file__, "-v", "-s"])