diff --git a/src/qforte/gate.cc b/src/qforte/gate.cc index 8512d53b..7b3fa134 100644 --- a/src/qforte/gate.cc +++ b/src/qforte/gate.cc @@ -27,25 +27,67 @@ size_t Gate::control() const { return control_; } const complex_4_4_mat& Gate::matrix() const { return gate_; } const SparseMatrix Gate::sparse_matrix(size_t nqubit) const { - size_t nbasis = std::pow(2, nqubit); - if (target_ != control_) { - throw std::runtime_error("Gate must be a Pauli to convert to matrix!"); - } else if (target_ >= nqubit) { + if (target_ >= nqubit) { throw std::runtime_error("Target index is too large for specified nqbits!"); } + if (control_ >= nqubit) { + throw std::runtime_error("Control index is too large for specified nqbits!"); + } + size_t nbasis = std::pow(2, nqubit); SparseMatrix Spmat = SparseMatrix(); - for (size_t i = 0; i < 2; i++) { - for (size_t j = 0; j < 2; j++) { - auto op_i_j = gate_[i][j]; - if (std::abs(op_i_j) > 1.0e-16) { - for (size_t I = 0; I < nbasis; I++) { - QubitBasis basis_I = QubitBasis(I); - if (basis_I.get_bit(target_) == j) { - QubitBasis basis_J = basis_I; - basis_J.set_bit(target_, i); - Spmat.set_element(basis_J.index(), basis_I.index(), op_i_j); + if (target_ == control_) { + // single-qubit case + // Iterate over the gate matrix elements + for (size_t i = 0; i < 2; i++) { + for (size_t j = 0; j < 2; j++) { + auto op_i_j = gate_[i][j]; + // Consider only non-zero elements (within a tolerance) + if (std::abs(op_i_j) > 1.0e-14) { + // Iterate over all basis states + for (size_t I = 0; I < nbasis; I++) { + QubitBasis basis_I = QubitBasis(I); + // Check if the target qubit is in the correct state + if (basis_I.get_bit(target_) == j) { + QubitBasis basis_J = basis_I; + // Set the target qubit to the new state + basis_J.set_bit(target_, i); + // Set the corresponding element in the sparse matrix + Spmat.set_element(basis_J.index(), basis_I.index(), op_i_j); + } + } + } + } + } + } else { + // two-qubit gate + for (size_t i = 0; i < 4; i++) { + for (size_t j = 0; j < 4; j++) { + auto op_i_j = gate_[i][j]; + if (std::abs(op_i_j) > 1.0e-14) { + for (size_t I = 0; I < nbasis; I++) { + QubitBasis basis_I = QubitBasis(I); + size_t target_bit = basis_I.get_bit(target_); + size_t control_bit = basis_I.get_bit(control_); + + // Combine the target and control bits into a 2-bit state + size_t combined_state = (control_bit << 1) | target_bit; + + // Check if the combined state matches the gate matrix column index + if (combined_state == j) { + QubitBasis basis_J = basis_I; + + // Extract the new states for control and target bits + size_t new_control_bit = (i >> 1) & 1; + size_t new_target_bit = i & 1; + + // Set the target and control qubits to the new states + basis_J.set_bit(target_, new_target_bit); + basis_J.set_bit(control_, new_control_bit); + + Spmat.set_element(basis_J.index(), basis_I.index(), op_i_j); + } } } } diff --git a/tests/test_sparse_op.py b/tests/test_sparse_op.py index ce00f69e..fb43ea7b 100644 --- a/tests/test_sparse_op.py +++ b/tests/test_sparse_op.py @@ -1,5 +1,8 @@ from qforte import build_circuit, Computer, QubitOperator import numpy as np +import qforte as qf +import random +from copy import deepcopy class TestSparseOp: @@ -102,3 +105,169 @@ def test_sparse_operator(self): print("Operator used for sparse matrix operator test: \n", qubit_op) print("||∆||: ", diff_val) assert diff_val < 1.0e-15 + + def test_sparse_gates(self): + """ + This test ensures that the sparse matrix representations of various + gates agrees with those obtained using the kronecker product with + the identity matrix + """ + + def sparse_matrix_to_numpy_array(sp_mat: dict, dim: int) -> np.array: + mat = np.zeros((dim, dim), dtype=complex) + for row in sp_mat: + for column in sp_mat[row]: + mat[row, column] = sp_mat[row][column] + return mat + + one_qubit_gate_pool = [ + "X", + "Y", + "Z", + "Rx", + "Ry", + "Rz", + "H", + "S", + "T", + "R", + "V", + "adj(V)", + ] + + two_qubit_gate_pool = [ + "CNOT", + "aCNOT", + "cY", + "cZ", + "cV", + "SWAP", + "cRz", + "cR", + "A", + "adj(cV)", + ] + + parametrized_gate_periods = { + "Rx": 4 * np.pi, + "Ry": 4 * np.pi, + "Rz": 4 * np.pi, + "R": 2 * np.pi, + "cR": 2 * np.pi, + "cRz": 4 * np.pi, + "A": 2 * np.pi, + } + + dim = 2 + identity = np.eye(2) + for gatetype in one_qubit_gate_pool: + if gatetype in parametrized_gate_periods: + period = parametrized_gate_periods[gatetype] + parameter = random.uniform(-period / 2, period / 2) + gate_0 = qf.gate(gatetype, 0, parameter) + gate_1 = qf.gate(gatetype, 1, parameter) + else: + if gatetype == "adj(V)": + gate_0 = qf.gate("V", 0) + gate_0 = gate_0.adjoint() + gate_1 = qf.gate("V", 1) + gate_1 = gate_1.adjoint() + else: + gate_0 = qf.gate(gatetype, 0) + gate_1 = qf.gate(gatetype, 1) + sparse_mat_0 = gate_0.sparse_matrix(2) + sparse_mat_1 = gate_1.sparse_matrix(2) + mat = sparse_matrix_to_numpy_array(gate_0.sparse_matrix(1).to_map(), dim) + mat_0 = sparse_matrix_to_numpy_array(sparse_mat_0.to_map(), dim * 2) + mat_1 = sparse_matrix_to_numpy_array(sparse_mat_1.to_map(), dim * 2) + mat_0_kron = np.kron(identity, mat) + mat_1_kron = np.kron(mat, identity) + assert np.all(mat_0 == mat_0_kron) + assert np.all(mat_1 == mat_1_kron) + + dim = 4 + for gatetype in two_qubit_gate_pool: + if gatetype in parametrized_gate_periods: + period = parametrized_gate_periods[gatetype] + parameter = random.uniform(-period / 2, period / 2) + gate_01 = qf.gate(gatetype, 0, 1, parameter) + gate_10 = qf.gate(gatetype, 1, 0, parameter) + gate_12 = qf.gate(gatetype, 1, 2, parameter) + gate_21 = qf.gate(gatetype, 2, 1, parameter) + else: + if gatetype == "adj(cV)": + gate_01 = qf.gate("cV", 0, 1) + gate_01 = gate_01.adjoint() + gate_10 = qf.gate("cV", 1, 0) + gate_10 = gate_10.adjoint() + gate_12 = qf.gate("cV", 1, 2) + gate_12 = gate_12.adjoint() + gate_21 = qf.gate("cV", 2, 1) + gate_21 = gate_21.adjoint() + else: + gate_01 = qf.gate(gatetype, 0, 1) + gate_10 = qf.gate(gatetype, 1, 0) + gate_12 = qf.gate(gatetype, 1, 2) + gate_21 = qf.gate(gatetype, 2, 1) + sparse_mat_01 = gate_01.sparse_matrix(3) + sparse_mat_10 = gate_10.sparse_matrix(3) + sparse_mat_12 = gate_12.sparse_matrix(3) + sparse_mat_21 = gate_21.sparse_matrix(3) + mat1 = sparse_matrix_to_numpy_array(gate_01.sparse_matrix(2).to_map(), dim) + mat2 = sparse_matrix_to_numpy_array(gate_10.sparse_matrix(2).to_map(), dim) + mat_01 = sparse_matrix_to_numpy_array(sparse_mat_01.to_map(), dim * 2) + mat_10 = sparse_matrix_to_numpy_array(sparse_mat_10.to_map(), dim * 2) + mat_12 = sparse_matrix_to_numpy_array(sparse_mat_12.to_map(), dim * 2) + mat_21 = sparse_matrix_to_numpy_array(sparse_mat_21.to_map(), dim * 2) + mat_01_kron = np.kron(identity, mat1) + mat_10_kron = np.kron(identity, mat2) + mat_12_kron = np.kron(mat1, identity) + mat_21_kron = np.kron(mat2, identity) + assert np.all(mat_01 == mat_01_kron) + assert np.all(mat_10 == mat_10_kron) + assert np.all(mat_12 == mat_12_kron) + assert np.all(mat_21 == mat_21_kron) + + def test_circuit_sparse_matrix(self): + Rhh = 1.5 + + mol = qf.system_factory( + system_type="molecule", + build_type="psi4", + basis="sto-6g", + mol_geometry=[ + ("H", (0, 0, -3 * Rhh / 2)), + ("H", (0, 0, -Rhh / 2)), + ("H", (0, 0, Rhh / 2)), + ("H", (0, 0, 3 * Rhh / 2)), + ], + symmetry="d2h", + multiplicity=1, + charge=0, + num_frozen_docc=0, + num_frozen_uocc=0, + run_mp2=0, + run_ccsd=0, + run_cisd=0, + run_fci=0, + ) + + for compact in [False, True]: + uccsd = qf.UCCNVQE( + mol, compact_excitations=compact, qubit_excitations=False + ) + uccsd.run( + pool_type="SD", opt_maxiter=200, optimizer="bfgs", opt_thresh=1.0e-5 + ) + + qc1 = qf.Computer(mol.hamiltonian.num_qubits()) + qc1.apply_circuit(uccsd.build_Uvqc()) + coeff_vec1 = deepcopy(qc1.get_coeff_vec()) + + Uvqc = uccsd.build_Uvqc() + sparse_Uvqc = Uvqc.sparse_matrix(mol.hamiltonian.num_qubits()) + qc2 = qf.Computer(mol.hamiltonian.num_qubits()) + qc2.apply_sparse_matrix(sparse_Uvqc) + coeff_vec2 = qc2.get_coeff_vec() + + assert np.linalg.norm(np.array(coeff_vec2) - np.array(coeff_vec1)) < 1.0e-14