Source code for pennylane.resource.second_quantization

# Copyright 2018-2022 Xanadu Quantum Technologies Inc.

# 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.
"""
This module contains the functions needed for resource estimation with the double factorization
method.
"""
# pylint: disable=no-self-use disable=too-many-arguments disable=too-many-instance-attributes
import numpy as np
from pennylane.operation import AnyWires, Operation
from pennylane.qchem import factorize


[docs]class DoubleFactorization(Operation): r"""Estimate the number of non-Clifford gates and logical qubits for a quantum phase estimation algorithm in second quantization with a double-factorized Hamiltonian. Atomic units are used throughout the class. Args: one_electron (array[array[float]]): one-electron integrals two_electron (tensor_like): two-electron integrals error (float): target error in the algorithm rank_r (int): rank of the first factorization of the two-electron integral tensor rank_m (int): average rank of the second factorization of the two-electron integral tensor tol_factor (float): threshold error value for discarding the negligible factors tol_eigval (float): threshold error value for discarding the negligible factor eigenvalues br (int): number of bits for ancilla qubit rotation alpha (int): number of bits for the keep register beta (int): number of bits for the rotation angles chemist_notation (bool): if True, the two-electron integrals need to be in chemist notation **Example** >>> symbols = ['O', 'H', 'H'] >>> geometry = np.array([[0.00000000, 0.00000000, 0.28377432], >>> [0.00000000, 1.45278171, -1.00662237], >>> [0.00000000, -1.45278171, -1.00662237]], requires_grad = False) >>> mol = qml.qchem.Molecule(symbols, geometry, basis_name='sto-3g') >>> core, one, two = qml.qchem.electron_integrals(mol)() >>> algo = DoubleFactorization(one, two) >>> print(algo.lamb, # the 1-Norm of the Hamiltonian >>> algo.gates, # estimated number of non-Clifford gates >>> algo.qubits # estimated number of logical qubits >>> ) 53.62085493277858 103969925 290 .. details:: :title: Theory To estimate the gate and qubit costs for implementing this method, the Hamiltonian needs to be factorized using the :func:`~.pennylane.qchem.factorize` function following [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. The objective of the factorization is to find a set of symmetric matrices, :math:`L^{(r)}`, such that the two-electron integral tensor in `chemist notation <http://vergil.chemistry.gatech.edu/notes/permsymm/permsymm.pdf>`_, :math:`V`, can be computed as .. math:: V_{ijkl} = \sum_r^R L_{ij}^{(r)} L_{kl}^{(r) T}, with the rank :math:`R \leq n^2`, where :math:`n` is the number of molecular orbitals. The matrices :math:`L^{(r)}` are diagonalized and for each matrix the eigenvalues that are smaller than a given threshold (and their corresponding eigenvectors) are discarded. The average number of the retained eigenvalues, :math:`M`, determines the rank of the second factorization step. The 1-norm of the Hamiltonian can then be computed using the :func:`~.pennylane.resource.DoubleFactorization.norm` function from the electron integrals and the eigenvalues of the matrices :math:`L^{(r)}`. The total number of gates and qubits for implementing the quantum phase estimation algorithm for the given Hamiltonian can then be computed using the functions :func:`~.pennylane.resource.DoubleFactorization.gate_cost` and :func:`~.pennylane.resource.DoubleFactorization.qubit_cost` with a target error that has the default value of 0.0016 Ha (chemical accuracy). The costs are computed using Eqs. (C39-C40) of [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. """ num_wires = AnyWires grad_method = None def __init__( self, one_electron, two_electron, error=0.0016, rank_r=None, rank_m=None, rank_max=None, tol_factor=1.0e-5, tol_eigval=1.0e-5, br=7, alpha=10, beta=20, chemist_notation=False, ): self.one_electron = one_electron if chemist_notation: self.two_electron = two_electron else: self.two_electron = np.swapaxes(two_electron, 1, 3) self.error = error self.rank_r = rank_r self.rank_m = rank_m self.rank_max = rank_max self.tol_factor = tol_factor self.tol_eigval = tol_eigval self.br = br self.alpha = alpha self.beta = beta self.n = two_electron.shape[0] * 2 self.factors, self.eigvals, self.eigvecs = factorize( self.two_electron, self.tol_factor, self.tol_eigval ) self.lamb = self.norm(self.one_electron, self.two_electron, self.eigvals) if not rank_r: self.rank_r = len(self.factors) if not rank_m: self.rank_m = np.mean([len(v) for v in self.eigvals]) if not rank_max: self.rank_max = int(np.max([len(v) for v in self.eigvals])) self.gates = self.gate_cost( self.n, self.lamb, self.error, self.rank_r, self.rank_m, self.rank_max, self.br, self.alpha, self.beta, ) self.qubits = self.qubit_cost( self.n, self.lamb, self.error, self.rank_r, self.rank_m, self.rank_max, self.br, self.alpha, self.beta, ) super().__init__(wires=range(self.qubits)) def _flatten(self): return (self.one_electron, self.two_electron), ( ("error", self.error), ("rank_r", self.rank_r), ("rank_m", self.rank_m), ("rank_max", self.rank_max), ("tol_factor", self.tol_factor), ("tol_eigval", self.tol_eigval), ("br", self.br), ("alpha", self.alpha), ("beta", self.beta), ("chemist_notation", True), ) @classmethod def _unflatten(cls, data, metadata): return cls(*data, **dict(metadata))
[docs] @staticmethod def estimation_cost(lamb, error): r"""Return the number of calls to the unitary needed to achieve the desired error in quantum phase estimation. The expression for computing the cost is taken from Eq. (45) of [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. Args: lamb (float): 1-norm of a second-quantized Hamiltonian error (float): target error in the algorithm Returns: int: number of calls to unitary **Example** >>> lamb = 72.49779513025341 >>> error = 0.001 >>> estimation_cost(lamb, error) 113880 """ if error <= 0.0: raise ValueError("The target error must be greater than zero.") if lamb <= 0.0: raise ValueError("The 1-norm must be greater than zero.") return int(np.ceil(np.pi * lamb / (2 * error)))
@staticmethod def _qrom_cost(constants): r"""Return the number of Toffoli gates and the expansion factor needed to implement a QROM for the double factorization method. The complexity of a QROM computation in the most general form is given by (see Eq. (C39) in [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]) .. math:: \text{cost} = \left \lceil \frac{a + b}{k} \right \rceil + \left \lceil \frac{c}{k} \right \rceil + d \left ( k + e \right ), where :math:`a, b, c, d, e` are constants that depend on the nature of the QROM implementation and the expansion factor :math:`k = 2^n` minimizes the cost. This function computes the optimum :math:`k` and the minimum cost for a QROM specification. To obtain the optimum values of :math:`k`, we first assume that the cost function is continuous and use differentiation to obtain the value of :math:`k` that minimizes the cost. This value of :math:`k` is not necessarily an integer power of 2. We then obtain the value of :math:`n` as :math:`n = \log_2(k)` and compute the cost for :math:`n_{int}= \left \{\left \lceil n \right \rceil, \left \lfloor n \right \rfloor \right \}`. The value of :math:`n_{int}` that gives the smaller cost is used to compute the optimim :math:`k`. Args: constants (tuple[float]): constants specifying a QROM Returns: tuple(int, int): the cost and the expansion factor for the QROM **Example** >>> constants = (151.0, 7.0, 151.0, 30.0, -1.0) >>> _qrom_cost(constants) 168, 4 """ a, b, c, d, e = constants n = np.log2(((a + b + c) / d) ** 0.5) k = np.array([2 ** np.floor(n), 2 ** np.ceil(n)]) cost = np.ceil((a + b) / k) + np.ceil(c / k) + d * (k + e) return int(cost[np.argmin(cost)]), int(k[np.argmin(cost)])
[docs] @staticmethod def unitary_cost(n, rank_r, rank_m, rank_max, br=7, alpha=10, beta=20): r"""Return the number of Toffoli gates needed to implement the qubitization unitary operator for the double factorization algorithm. The expression for computing the cost is taken from Eq. (C39) of [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. Args: n (int): number of molecular spin-orbitals rank_r (int): rank of the first factorization of the two-electron integral tensor rank_m (float): average rank of the second factorization of the two-electron tensor rank_max (int): maximum rank of the second factorization of the two-electron tensor br (int): number of bits for ancilla qubit rotation alpha (int): number of bits for the keep register beta (int): number of bits for the rotation angles Returns: int: number of Toffoli gates to implement the qubitization unitary **Example** >>> n = 14 >>> rank_r = 26 >>> rank_m = 5.5 >>> rank_max = 7 >>> br = 7 >>> alpha = 10 >>> beta = 20 >>> unitary_cost(n, rank_r, rank_m, rank_max, br, alpha, beta) 2007 """ if n <= 0 or not isinstance(n, (int, np.integer)) or n % 2 != 0: raise ValueError("The number of spin-orbitals must be a positive even integer.") if rank_r <= 0 or not isinstance(rank_r, int): raise ValueError("The rank of the first factorization step must be a positive integer.") if rank_m <= 0: raise ValueError("The rank of the second factorization step must be a positive number.") if rank_max <= 0 or not isinstance(rank_max, int): raise ValueError( "The maximum rank of the second factorization step must be a positive integer." ) if br <= 0 or not isinstance(br, int): raise ValueError("br must be a positive integer.") if alpha <= 0 or not isinstance(alpha, int): raise ValueError("alpha must be a positive integer.") if beta <= 0 or not isinstance(beta, int): raise ValueError("beta must be a positive integer.") rank_rm = rank_r * rank_m # eta is computed based on step 1.(a) in page 030305-41 of PRX Quantum 2, 030305 (2021) eta = np.array([np.log2(n) for n in range(1, rank_r + 1) if rank_r % n == 0]) eta = int(np.max([n for n in eta if n % 1 == 0])) nxi = np.ceil(np.log2(rank_max)) # Eq. (C14) of PRX Quantum 2, 030305 (2021) nl = np.ceil(np.log2(rank_r + 1)) # Eq. (C14) of PRX Quantum 2, 030305 (2021) nlxi = np.ceil(np.log2(rank_rm + n / 2)) # Eq. (C15) of PRX Quantum 2, 030305 (2021) bp1 = nl + alpha # Eq. (C27) of PRX Quantum 2, 030305 (2021) bo = nxi + nlxi + br + 1 # Eq. (C29) of PRX Quantum 2, 030305 (2021) bp2 = nxi + alpha + 2 # Eq. (C31) of PRX Quantum 2, 030305 (2021) # cost is computed using Eq. (C39) of PRX Quantum 2, 030305 (2021) cost = ( 9 * nl - 6 * eta + 12 * br + 34 * nxi + 8 * nlxi + 9 * alpha + 3 * n * beta - 6 * n - 43 ) cost += DoubleFactorization._qrom_cost((rank_r, 1, 0, bp1, -1))[0] cost += DoubleFactorization._qrom_cost((rank_r, 1, 0, bo, -1))[0] cost += DoubleFactorization._qrom_cost((rank_r, 1, 0, 1, 0))[0] * 2 cost += DoubleFactorization._qrom_cost((rank_rm, n / 2, rank_rm, n * beta, 0))[0] cost += DoubleFactorization._qrom_cost((rank_rm, n / 2, rank_rm, 2, 0))[0] * 2 cost += DoubleFactorization._qrom_cost((rank_rm, n / 2, rank_rm, 2 * bp2, -1))[0] return int(cost)
[docs] @staticmethod def gate_cost(n, lamb, error, rank_r, rank_m, rank_max, br=7, alpha=10, beta=20): r"""Return the total number of Toffoli gates needed to implement the double factorization algorithm. The expression for computing the cost is taken from Eqs. (45) and (C39) of [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. Args: n (int): number of molecular spin-orbitals lamb (float): 1-norm of a second-quantized Hamiltonian error (float): target error in the algorithm rank_r (int): rank of the first factorization of the two-electron integral tensor rank_m (float): average rank of the second factorization of the two-electron tensor rank_max (int): maximum rank of the second factorization of the two-electron tensor br (int): number of bits for ancilla qubit rotation alpha (int): number of bits for the keep register beta (int): number of bits for the rotation angles Returns: int: the number of Toffoli gates for the double factorization method **Example** >>> n = 14 >>> lamb = 52.98761457453095 >>> error = 0.001 >>> rank_r = 26 >>> rank_m = 5.5 >>> rank_max = 7 >>> br = 7 >>> alpha = 10 >>> beta = 20 >>> gate_cost(n, lamb, error, rank_r, rank_m, rank_max, br, alpha, beta) 167048631 """ if n <= 0 or not isinstance(n, (int, np.integer)) or n % 2 != 0: raise ValueError("The number of spin-orbitals must be a positive even integer.") if error <= 0.0: raise ValueError("The target error must be greater than zero.") if lamb <= 0.0: raise ValueError("The 1-norm must be greater than zero.") if rank_r <= 0 or not isinstance(rank_r, int): raise ValueError("The rank of the first factorization step must be a positive integer.") if rank_m <= 0: raise ValueError("The rank of the second factorization step must be a positive number.") if rank_max <= 0 or not isinstance(rank_max, int): raise ValueError( "The maximum rank of the second factorization step must be a positive integer." ) if br <= 0 or not isinstance(br, int): raise ValueError("br must be a positive integer.") if alpha <= 0 or not isinstance(alpha, int): raise ValueError("alpha must be a positive integer.") if beta <= 0 or not isinstance(beta, int): raise ValueError("beta must be a positive integer.") e_cost = DoubleFactorization.estimation_cost(lamb, error) u_cost = DoubleFactorization.unitary_cost(n, rank_r, rank_m, rank_max, br, alpha, beta) return int(e_cost * u_cost)
[docs] @staticmethod def qubit_cost(n, lamb, error, rank_r, rank_m, rank_max, br=7, alpha=10, beta=20): r"""Return the number of logical qubits needed to implement the double factorization method. The expression for computing the cost is taken from Eq. (C40) of [`PRX Quantum 2, 030305 (2021) <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>`_]. Args: n (int): number of molecular spin-orbitals lamb (float): 1-norm of a second-quantized Hamiltonian error (float): target error in the algorithm rank_r (int): rank of the first factorization of the two-electron integral tensor rank_m (float): average rank of the second factorization of the two-electron tensor rank_max (int): maximum rank of the second factorization of the two-electron tensor br (int): number of bits for ancilla qubit rotation alpha (int): number of bits for the keep register beta (int): number of bits for the rotation angles Returns: int: number of logical qubits for the double factorization method **Example** >>> n = 14 >>> lamb = 52.98761457453095 >>> error = 0.001 >>> rank_r = 26 >>> rank_m = 5.5 >>> rank_max = 7 >>> br = 7 >>> alpha = 10 >>> beta = 20 >>> qubit_cost(n, lamb, error, rank_r, rank_m, rank_max, br, alpha, beta) 292 """ if n <= 0 or not isinstance(n, (int, np.integer)) or n % 2 != 0: raise ValueError("The number of spin-orbitals must be a positive even integer.") if error <= 0.0: raise ValueError("The target error must be greater than zero.") if lamb <= 0.0: raise ValueError("The 1-norm must be greater than zero.") if rank_r <= 0 or not isinstance(rank_r, int): raise ValueError("The rank of the first factorization step must be a positive integer.") if rank_m <= 0: raise ValueError("The rank of the second factorization step must be a positive number.") if rank_max <= 0 or not isinstance(rank_max, int): raise ValueError( "The maximum rank of the second factorization step must be a positive integer." ) if br <= 0 or not isinstance(br, int): raise ValueError("br must be a positive integer.") if alpha <= 0 or not isinstance(alpha, int): raise ValueError("alpha must be a positive integer.") if beta <= 0 or not isinstance(beta, int): raise ValueError("beta must be a positive integer.") rank_rm = rank_r * rank_m nxi = np.ceil(np.log2(rank_max)) # Eq. (C14) of PRX Quantum 2, 030305 (2021) nl = np.ceil(np.log2(rank_r + 1)) # Eq. (C14) of PRX Quantum 2, 030305 (2021) nlxi = np.ceil(np.log2(rank_rm + n / 2)) # Eq. (C15) of PRX Quantum 2, 030305 (2021) bo = nxi + nlxi + br + 1 # Eq. (C29) of PRX Quantum 2, 030305 (2021) bp2 = nxi + alpha + 2 # Eq. (C31) of PRX Quantum 2, 030305 (2021) # kr is taken from Eq. (C39) of PRX Quantum 2, 030305 (2021) kr = DoubleFactorization._qrom_cost((rank_rm, n / 2, rank_rm, n * beta, 0))[1] # the cost is computed using Eq. (C40) of PRX Quantum 2, 030305 (2021) e_cost = DoubleFactorization.estimation_cost(lamb, error) cost = n + 2 * nl + nxi + 3 * alpha + beta + bo + bp2 cost += kr * n * beta / 2 + 2 * np.ceil(np.log2(e_cost + 1)) + 7 return int(cost)
[docs] @staticmethod def norm(one, two, eigvals): r"""Return the 1-norm of a molecular Hamiltonian from the one- and two-electron integrals and eigenvalues of the factorized two-electron integral tensor. The 1-norm of a double-factorized molecular Hamiltonian is computed using Eqs. (15-17) of [`Phys. Rev. Research 3, 033055 (2021) <https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.3.033055>`_] .. math:: \lambda = ||T|| + \frac{1}{4} \sum_r ||L^{(r)}||^2, where the Schatten 1-norm for a given matrix :math:`T` is defined as .. math:: ||T|| = \sum_k |\text{eigvals}[T]_k|. The matrices :math:`L^{(r)}` are obtained from factorization of the two-electron integral tensor :math:`V` such that .. math:: V_{ijkl} = \sum_r L_{ij}^{(r)} L_{kl}^{(r) T}. The matrix :math:`T` is constructed from the one- and two-electron integrals as .. math:: T = h_{ij} - \frac{1}{2} \sum_l V_{illj} + \sum_l V_{llij}. Note that the two-electron integral tensor must be arranged in chemist notation. Args: one (array[array[float]]): one-electron integrals two (array[array[float]]): two-electron integrals eigvals (array[float]): eigenvalues of the matrices obtained from factorizing the two-electron integral tensor Returns: array[float]: 1-norm of the Hamiltonian **Example** >>> symbols = ['H', 'H', 'O'] >>> geometry = np.array([[0.00000000, 0.00000000, 0.28377432], >>> [0.00000000, 1.45278171, -1.00662237], >>> [0.00000000, -1.45278171, -1.00662237]], requires_grad=False) >>> mol = qml.qchem.Molecule(symbols, geometry, basis_name='sto-3g') >>> core, one, two = qml.qchem.electron_integrals(mol)() >>> two = np.swapaxes(two, 1, 3) # convert to the chemists notation >>> _, eigvals, _ = qml.qchem.factorize(two, 1e-5) >>> print(norm(one, two, eigvals)) 52.98762043980203 """ t_matrix = one - 0.5 * np.einsum("illj", two) + np.einsum("llij", two) t_eigvals, _ = np.linalg.eigh(t_matrix) lambda_one = np.sum(abs(t_eigvals)) lambda_two = 0.25 * np.sum([np.sum(abs(v)) ** 2 for v in eigvals]) return lambda_one + lambda_two