refactor basis rotation

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)

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

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
+
+
+