From 7169e087367f052b0e563dbc37d9634f4b6ff2e2 Mon Sep 17 00:00:00 2001 From: Adrian Mak Date: Mon, 25 Mar 2024 09:58:05 +0800 Subject: [PATCH] refactor basis rotation --- examples/br_example.py | 48 ++++--- src/qibochem/ansatz/__init__.py | 6 +- src/qibochem/ansatz/basis_rotation.py | 184 ++++++-------------------- tests/test_basis_rotation.py | 165 +++++++++++------------ 4 files changed, 142 insertions(+), 261 deletions(-) diff --git a/examples/br_example.py b/examples/br_example.py index 36ba772..ea6226a 100644 --- a/examples/br_example.py +++ b/examples/br_example.py @@ -4,10 +4,11 @@ """ import numpy as np +from qibo import Circuit, gates, models from qibo.optimizers import optimize from scipy.optimize import minimize -from qibochem.ansatz.basis_rotation import br_circuit +from qibochem.ansatz import basis_rotation from qibochem.ansatz.hf_reference import hf_circuit from qibochem.driver.molecule import Molecule from qibochem.measurement.expectation import expectation @@ -19,7 +20,6 @@ except ModuleNotFoundError: mol.run_psi4() - # Diagonalize H_core to use as the guess wave function # First, symmetric orthogonalization of overlap matrix S using np: @@ -47,35 +47,33 @@ ) -def electronic_energy(parameters): - """ - Loss function (Electronic energy) for the basis rotation ansatz - """ - circuit = hf_circuit(mol.nso, mol.nelec) - circuit += br_circuit(mol.nso, parameters, mol.nelec) +def basis_rotation_circuit(mol, parameters=0.0): - return expectation(circuit, hamiltonian) + nqubits = mol.nso + occ = range(0, mol.nelec) + vir = range(mol.nelec, mol.nso) + U, theta = basis_rotation.unitary(occ, vir, parameters=parameters) -# Build a basis rotation_circuit -params = np.random.rand(mol.nelec * (mol.nso - mol.nelec)) # n_occ * n_virt + gate_angles, final_U = basis_rotation.givens_qr_decompose(U) + gate_layout = basis_rotation.basis_rotation_layout(nqubits) + gate_list, ordered_angles = basis_rotation.basis_rotation_gates(gate_layout, gate_angles, theta) -best, params, extra = optimize(electronic_energy, params) + circuit = Circuit(nqubits) + for _i in range(mol.nelec): + circuit.add(gates.X(_i)) + circuit.add(gate_list) -print("\nResults using Qibo optimize:") -print(f" HF energy: {mol.e_hf:.8f}") -print(f"VQE energy: {best:.8f} (Basis rotation ansatz)") -# print() -# print("Optimized parameters:", params) + return circuit, gate_angles +br_circuit, qubit_parameters = basis_rotation_circuit(mol, parameters=0.1) +vqe = models.VQE(br_circuit, hamiltonian) +vqe_result = vqe.minimize(qubit_parameters) -params = np.random.rand(mol.nelec * (mol.nso - mol.nelec)) # n_occ * n_virt +print(f" HF energy: {mol.e_hf:.8f}") +print(f"VQE energy: {vqe_result[0]:.8f} (Basis rotation ansatz)") +print() +print("Optimized qubit parameters:\n", vqe_result[1]) +print("Optimizer message:\n", vqe_result[2]) -result = minimize(electronic_energy, params) -best, params = result.fun, result.x -print("\nResults using scipy.optimize:") -print(f" HF energy: {mol.e_hf:.8f}") -print(f"VQE energy: {best:.8f} (Basis rotation ansatz)") -# print() -# print("Optimized parameters:", params) diff --git a/src/qibochem/ansatz/__init__.py b/src/qibochem/ansatz/__init__.py index 6ed1a80..a94b14c 100644 --- a/src/qibochem/ansatz/__init__.py +++ b/src/qibochem/ansatz/__init__.py @@ -1,12 +1,8 @@ from qibochem.ansatz.basis_rotation import ( basis_rotation_gates, basis_rotation_layout, - br_circuit, givens_qr_decompose, - givens_rotation_gate, - givens_rotation_parameters, - swap_matrices, - unitary_rot_matrix, + unitary, ) from qibochem.ansatz.hardware_efficient import hea from qibochem.ansatz.hf_reference import bk_matrix, bk_matrix_power2, hf_circuit diff --git a/src/qibochem/ansatz/basis_rotation.py b/src/qibochem/ansatz/basis_rotation.py index 0759b51..6c7841b 100644 --- a/src/qibochem/ansatz/basis_rotation.py +++ b/src/qibochem/ansatz/basis_rotation.py @@ -3,167 +3,59 @@ """ import numpy as np -from openfermion.linalg.givens_rotations import givens_decomposition +#from openfermion.linalg.givens_rotations import givens_decomposition from qibo import gates, models from scipy.linalg import expm - +from qibochem.driver import hamiltonian from qibochem.ansatz import ucc # Helper functions -def unitary_rot_matrix(parameters, occ_orbitals, virt_orbitals, orbital_pairs=None, conserve_spin=False): - r""" - Returns the unitary rotation matrix U = exp(\kappa) mixing the occupied and virtual orbitals, - where \kappa is an anti-Hermitian matrix - +def unitary(occ_orbitals, virt_orbitals, parameters=None): + """ + Returns the unitary rotation matrix U = exp(\kappa) mixing the occupied and virtual orbitals. Args: - parameters: List/array of rotation parameters. dim = len(occ_orbitals)*len(virt_orbitals) + occ_orbitals: Iterable of occupied orbitals virt_orbitals: Iterable of virtual orbitals - orbital_pairs: (optional) Iterable of occupied-virtual orbital pairs, sorting is optional - conserve_spin: (optional) in the case where orbital pairs are unspecified, they are generated - from ucc.generate_excitations for singles excitations. If True, only singles - excitations with spin conserved are considered + parameters: List/array of rotation parameters. dim = len(occ_orbitals)*len(virt_orbitals) Returns: exp(k): Unitary matrix of Givens rotations, obtained by matrix exponential of skew-symmetric kappa matrix """ + + # conserve_spin has to be true for SCF/basis_rotation cases, else expm(k) is not unitary + ov_pairs = ucc.generate_excitations(1, occ_orbitals, virt_orbitals, conserve_spin=True) + # print('ov_pairs presort', ov_pairs) + ov_pairs = ucc.sort_excitations(ov_pairs) + #print('ov_pairs sorted ', ov_pairs) + n_theta = len(ov_pairs) + if parameters is None: + print('basis rotation: parameters not specified') + print('basis rotation: using default occ-virt rotation parameter value = 0.0') + parameters = np.zeros(n_theta) + elif isinstance(parameters, float): + print('basis rotation: using uniform value of', parameters, 'for each parameter value') + parameters = np.full(n_theta, parameters) + else: + if len(parameters) != n_theta: + raise IndexError('parameter array specified has bad size or type') + print('basis rotation: loading parameters from input') + n_orbitals = len(occ_orbitals) + len(virt_orbitals) kappa = np.zeros((n_orbitals, n_orbitals)) - if orbital_pairs == None: - # Get all possible (occ, virt) pairs - # orbital_pairs = [(_o, _v) for _o in occ_orbitals for _v in virt_orbitals] - orbital_pairs = ucc.generate_excitations(1, occ_orbitals, virt_orbitals, conserve_spin=conserve_spin) - # occ/occ and virt/virt orbital mixing doesn't change the energy, so can ignore - for _i, (_occ, _virt) in enumerate(orbital_pairs): + + for _i, (_occ, _virt) in enumerate(ov_pairs): kappa[_occ, _virt] = parameters[_i] kappa[_virt, _occ] = -parameters[_i] - return expm(kappa) - - -def swap_matrices(permutations, n_qubits): - r""" - Build matrices that permute (swap) the columns of a given matrix - - Args: - permutations [List(tuple, ), ]: List of lists of 2-tuples permuting two consecutive numbers - e.g. [[(0, 1), (2, 3)], [(1, 2)], ...] - n_qubits: Dimension of matrix (\exp(\kappa) matrix) to be operated on i.e. # of qubits - - Returns: - List of matrices that carry out \exp(swap) for each pair in permutations - """ - exp_swaps = [] - _temp = np.array([[-1.0, 1.0], [1.0, -1.0]]) - for pairs in permutations: - gen = np.zeros((n_qubits, n_qubits)) - for pair in pairs: - # Add zeros to pad out the (2, 2) matrix to (n_qubits, n_qubits) with 0. - gen += np.pad(_temp, ((pair[0], n_qubits - pair[1] - 1),), "constant") - # Take matrix exponential, remove all entries close to zero, and convert to real matrix - exp_mat = expm(-0.5j * np.pi * gen) - exp_mat[np.abs(exp_mat) < 1e-10] = 0.0 - exp_mat = np.real_if_close(exp_mat) - # Add exp_mat to exp_swaps once cleaned up - exp_swaps.append(exp_mat) - return exp_swaps - - -def givens_rotation_parameters(n_qubits, exp_k, n_occ): - r""" - Get the parameters of the Givens rotation matrix obtained after decomposing \exp(\kappa) - - Args: - n_qubits: Number of qubits (i.e. spin-orbitals) - exp_k: Unitary rotation matrix - n_occ: Number of occupied orbitals - """ - # List of 2-tuples to link all qubits via swapping - swap_pairs = [ - [(_i, _i + 1) for _i in range(_d % 2, n_qubits - 1, 2)] - # Only 2 possible ways of swapping: [(0, 1), ... ] or [(1, 2), ...] - for _d in range(2) - ] - # Get their corresponding matrices - swap_unitaries = swap_matrices(swap_pairs, n_qubits) - - # Create a copy of exp_k for rotating - exp_k_shift = np.copy(exp_k) - # Apply the column-swapping transformations - for u_rot in swap_unitaries: - exp_k_shift = u_rot @ exp_k_shift - # Only want the n_occ by n_qubits (n_orb) sub-matrix - qmatrix = exp_k_shift.T[:n_occ, :] - - # Apply Givens decomposition of the submatrix to get a list of individual Givens rotations - # (orb1, orb2, theta, phi), where orb1 and orb2 are rotated by theta(real) and phi(imag) - qdecomp, _, _ = givens_decomposition(qmatrix) - # Lastly, reverse the list of Givens rotations (because U = G_k ... G_2 G_1) - return list(reversed(qdecomp)) - - -def givens_rotation_gate(n_qubits, orb1, orb2, theta): - r""" - Circuit corresponding to a Givens rotation between two qubits (spin-orbitals) - Args: - n_qubits: Number of qubits in the quantum circuit - orb1, orb2: Orbitals used for the Givens rotation - theta: Rotation parameter - - Returns: - circuit: Qibo Circuit object representing a Givens rotation between orb1 and orb2 - """ - # Define 2x2 sqrt(iSwap) matrix - iswap_mat = np.array([[1.0, 1.0j], [1.0j, 1.0]]) / np.sqrt(2.0) - # Build and add the gates to circuit - circuit = models.Circuit(n_qubits) - circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False)) - # circuit.add(gates.SiSWAP(orb1, orb2)) - circuit.add(gates.RZ(orb1, -theta)) - circuit.add(gates.RZ(orb2, theta + np.pi)) - circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False)) - # circuit.add(gates.SiSWAP(orb1, orb2)) - circuit.add(gates.RZ(orb1, np.pi, trainable=False)) - return circuit - - -def br_circuit(n_qubits, parameters, n_occ): - r""" - Google's basis rotation circuit, applied between the occupied/virtual orbitals. Forms the exp(kappa) matrix, decomposes - it into Givens rotations, and sets the circuit parameters based on the Givens rotation decomposition. Note: Supposed - to be used with the JW fermion-to-qubit mapping - - Args: - n_qubits: Number of qubits in the quantum circuit - parameters: Rotation parameters for exp(kappa); Must have (n_occ * n_virt) parameters - n_occ: Number of occupied orbitals - - Returns: - Qibo ``Circuit``: Circuit corresponding to the basis rotation ansatz between the occupied and virtual orbitals - """ - assert len(parameters) == (n_occ * (n_qubits - n_occ)), "Need len(parameters) == (n_occ * n_virt)" - # Unitary rotation matrix \exp(\kappa) - exp_k = unitary_rot_matrix(parameters, range(n_occ), range(n_occ, n_qubits)) - # List of Givens rotation parameters - g_rotation_parameters = givens_rotation_parameters(n_qubits, exp_k, n_occ) - - # Start building the circuit - circuit = models.Circuit(n_qubits) - # Add the Givens rotation gates - for g_rot_parameter in g_rotation_parameters: - for orb1, orb2, theta, _phi in g_rot_parameter: - assert np.allclose(_phi, 0.0), "Unitary rotation is not real!" - circuit += givens_rotation_gate(n_qubits, orb1, orb2, theta) - return circuit + return expm(kappa), parameters # clements qr routines - - def givens_qr_decompose(U): - r""" + """ Clements scheme QR decompose a unitary matrix U using Givens rotations see arxiv:1603.08788 Args: @@ -254,7 +146,7 @@ def col_op(U, row, col): def basis_rotation_layout(N): - r""" + """ generates the layout of the basis rotation circuit for Clements scheme QR decomposition Args: N: number of qubits/modes @@ -350,7 +242,7 @@ def assign_element(A, row, col, k): def basis_rotation_gates(A, z_array, parameters): - r""" + """ places the basis rotation gates on circuit in the order of Clements scheme QR decomposition Args: A: @@ -366,17 +258,23 @@ def basis_rotation_gates(A, z_array, parameters): Outputs: gate_list: list of gates which implement the basis rotation using Clements scheme QR decomposition + ordered_angles: + list of angles ordered by sequence of singles excitation gates added to circuit """ N = len(A[0]) gate_list = [] - + ordered_angles = [] # for j in range(N): for i in range(N): if A[i][j] == 0: - print("CRY", i, i + 1, j, z_array[A[i + 1][j] - 1]) + #print("CRY", i, i + 1, A[i+1][j]-1, z_array[A[i + 1][j] - 1]) gate_list.append(gates.CNOT(i + 1, i)) gate_list.append(gates.CRY(i, i + 1, z_array[A[i + 1][j] - 1])) + ordered_angles.append(z_array[A[i + 1][j] - 1]) gate_list.append(gates.CNOT(i + 1, i)) - return gate_list + return gate_list, ordered_angles + + + diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py index 911a152..94a97af 100644 --- a/tests/test_basis_rotation.py +++ b/tests/test_basis_rotation.py @@ -5,97 +5,86 @@ import numpy as np import pytest from qibo import Circuit, gates -from qibo.optimizers import optimize +#from qibo.optimizers import optimize from qibochem.ansatz import basis_rotation from qibochem.driver.molecule import Molecule from qibochem.measurement.expectation import expectation -def test_givens_rotation_gate(): - n_qubits = 2 - orb1 = 0 - orb2 = 1 - theta = -0.1 - circuit = basis_rotation.givens_rotation_gate(n_qubits, orb1, orb2, theta) - ref_u = np.array( - [ - [ - -1.0, - 0.0, - 0.0, - 0.0, - ], - [0.0, 0.99500417, 0.09983342, 0.0], - [0.0, -0.09983342, 0.99500417, 0.0], - [0.0, 0.0, 0.0, -1.0], - ] - ) - - assert np.allclose(circuit.unitary(), ref_u) - - -def test_givens_rotation_parameters(): - u = np.array( - [[0.99001666, 0.099667, 0.099667], [-0.099667, 0.99500833, -0.00499167], [-0.099667, -0.00499167, 0.99500833]] - ) - n_occ = 1 - n_qubits = 3 - params = basis_rotation.givens_rotation_parameters(n_qubits, u, n_occ) - ref_params = [((0, 1, -1.4709635780470989, 0.0),), ((1, 2, 1.4704623293305714, 0.0),)] - - assert np.allclose(params, ref_params) - - -def test_swap_matrices(): - """Test for swap_matrices""" - permutations = [[(0, 1), (2, 3)]] - n_qubits = 4 - smat = basis_rotation.swap_matrices(permutations, n_qubits) - ref_smat = np.array([[0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0]]) - - assert np.allclose(smat, ref_smat) - - -def test_unitary_rot_matrix(): - """Test for unitary rotation matrix""" - occ = [0] - vir = [1, 2] - parameters = np.zeros(len(occ) * len(vir)) - parameters += 0.1 - u = basis_rotation.unitary_rot_matrix(parameters, occ, vir, orbital_pairs=None, conserve_spin=False) - ref_u = np.array( - [[0.99001666, 0.099667, 0.099667], [-0.099667, 0.99500833, -0.00499167], [-0.099667, -0.00499167, 0.99500833]] - ) - - assert np.allclose(u, ref_u) - - -def test_br_ansatz(): - """Test of basis rotation ansatz against hardcoded HF energies""" - h2_ref_energy = -1.117349035 - - mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))]) - mol.run_pyscf(max_scf_cycles=1) - # Use an un-converged wave function to build the Hamiltonian - mol_sym_ham = mol.hamiltonian() - - # Define quantum circuit - circuit = Circuit(mol.nso) - circuit.add(gates.X(_i) for _i in range(mol.nelec)) - - # Add basis rotation ansatz - # Initialize all at zero - parameters = np.zeros(mol.nelec * (mol.nso - mol.nelec)) # n_occ * n_virt - circuit += basis_rotation.br_circuit(mol.nso, parameters, mol.nelec) - - def electronic_energy(parameters): - """ - Loss function (Electronic energy) for the basis rotation ansatz - """ - circuit.set_parameters(parameters) - return expectation(circuit, mol_sym_ham) - - hf_energy, parameters, _extra = optimize(electronic_energy, parameters) - - assert hf_energy == pytest.approx(h2_ref_energy) +def test_unitary(): + """ + Test for basis_rotation.unitary() + """ + N = 6 + occ = range(0, 2) + vir = range(2, 6) + + preset_params = [-0.1, -0.2, -0.3, -0.4] + + U1, theta1 = basis_rotation.unitary(occ, vir) + U2, theta2 = basis_rotation.unitary(occ, vir, parameters=0.1) + U3, theta3 = basis_rotation.unitary(occ, vir, parameters=preset_params) + + ref_U2 = np.array( + [[ 0.99001666, 0., 0.099667, 0., 0.099667, 0. ], + [ 0., 0.99001666, 0., 0.099667, 0., 0.099667 ], + [-0.099667, 0., 0.99500833, 0., -0.00499167, 0. ], + [ 0., -0.099667, 0., 0.99500833, 0., -0.00499167], + [-0.099667, 0., -0.00499167, 0., 0.99500833, 0. ], + [ 0., -0.099667, 0., -0.00499167, 0., 0.99500833]]) + + ref_U3 = np.array([[ 0.95041528, 0., -0.09834165, 0., -0.29502494, 0. ], + [ 0., 0.9016556, 0., -0.19339968, 0., -0.38679937], + [ 0.09834165, 0., 0.99504153, 0., -0.01487542, 0. ], + [ 0., 0.19339968, 0., 0.98033112, 0., -0.03933776], + [ 0.29502494, 0., -0.01487542, 0., 0.95537375, 0. ], + [ 0., 0.38679937, 0., -0.03933776, 0., 0.92132448]]) + + identity = np.eye(6) + + assert np.allclose(U1@U1.T, identity) + assert np.allclose(U2@U2.T, identity) + assert np.allclose(U3@U3.T, identity) + assert np.allclose(U1, identity) + assert np.allclose(U2, ref_U2) + assert np.allclose(U3, ref_U3) + + +def test_givens_qr_decompose(): + """ + Test for basis_rotation.givens_qr_decompose() + """ + N = 6 + occ = range(0,2) + vir = range(2,6) + + U2, theta2 = basis_rotation.unitary(occ, vir, parameters=0.1) + z_angles, final_U2 = basis_rotation.givens_qr_decompose(U2) + + ref_z = np.array([-3.141592653589793, -1.5707963267948966, -2.356194490192345, 0.0, -1.5707963267948966, -1.5207546393123066, -1.5707963267948966, -1.5707963267948954, -3.000171297352484, -2.356194490192345, 0.0, -1.5707963267948966, -1.5707963267948954, -0.09995829685982476, -1.5207546393123068]) + + assert np.allclose(z_angles, ref_z) + assert np.allclose(np.eye(6), final_U2) + + +def test_basis_rotation_layout(): + + N = 10 + ref_A = + np.array( + [[ 0, -1, 0, -1, 0, -1, 0, -1, 0, -1], + [ 1, 0, 6, 0, 15, 0, 28, 0, 45, 0], + [ 0, 5, 0, 14, 0, 27, 0, 44, 0, 29], + [ 4, 0, 13, 0, 26, 0, 43, 0, 30, 0], + [ 0, 12, 0, 25, 0, 42, 0, 31, 0, 16], + [11, 0, 24, 0, 41, 0, 32, 0, 17, 0], + [ 0, 23, 0, 40, 0, 33, 0, 18, 0, 7], + [22, 0, 39, 0, 34, 0, 19, 0, 8, 0], + [ 0, 38, 0, 35, 0, 20, 0, 9, 0, 2], + [37, -1, 36, -1, 21, -1, 10, -1, 3, -1]]) + + A = basis_rotation.basis_rotation_layout(N) + assert np.allclose(A, ref_A) + +