Source code for pennylane.resource.first_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 estimating the number of logical qubits and
non-Clifford gates for quantum algorithms in first quantization using a plane-wave basis.
"""
# pylint: disable=no-self-use disable=too-many-arguments disable=too-many-instance-attributes
import numpy as np
from scipy import integrate
from pennylane.operation import AnyWires, Operation
[docs]class FirstQuantization(Operation):
r"""Estimate the number of non-Clifford gates and logical qubits for a quantum phase estimation
algorithm in first quantization using a plane-wave basis.
To estimate the gate and qubit costs for implementing this method, the number of plane waves,
the number of electrons and the lattice vectors need to be defined. The costs can then be
computed using the functions :func:`~.pennylane.resource.FirstQuantization.gate_cost` and
:func:`~.pennylane.resource.FirstQuantization.qubit_cost` with a target error that has the default
value of 0.0016 Ha (chemical accuracy). Atomic units are used throughout the class.
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
charge (int): total electric charge of the system
br (int): number of bits for ancilla qubit rotation
vectors (array[float]): lattice vectors
**Example**
>>> n = 100000
>>> eta = 156
>>> vectors = np.array([[10.46219511, 0.00000000, 0.00000000],
... [ 0.00000000, 10.46219511, 0.00000000],
... [ 0.00000000, 0.00000000, 10.46219511]])
>>> algo = FirstQuantization(n, eta, vectors=vectors)
>>> print(algo.lamb, # the 1-Norm of the Hamiltonian
>>> algo.gates, # estimated number of non-Clifford gates
>>> algo.qubits # estimated number of logical qubits
>>> )
649912.4804278888 1.1e+13 4416
.. details::
:title: Theory
Following `PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_
, the target algorithm error, :math:`\epsilon`, is distributed among four different sources
of error using Eq. (131)
.. math::
\epsilon^2 \geq \epsilon_{qpe}^2 + (\epsilon_{\mathcal{M}} + \epsilon_R + \epsilon_T)^2,
where :math:`\epsilon_{qpe}` is the quantum phase estimation error and
:math:`\epsilon_{\mathcal{M}}`, :math:`\epsilon_R`, and :math:`\epsilon_T` are defined in
Eqs. (132-134).
Here, we fix :math:`\epsilon_{\mathcal{M}} = \epsilon_R = \epsilon_T = \alpha \epsilon` with
a default value of :math:`\alpha = 0.01` and obtain
.. math::
\epsilon_{qpe} = \sqrt{\epsilon^2 [1 - (3 \alpha)^2]}.
Note that the user only needs to define the target algorithm error :math:`\epsilon`. The
error distribution takes place inside the functions.
"""
num_wires = AnyWires
grad_method = None
def __init__(
self,
n,
eta,
omega=None,
error=0.0016,
charge=0,
br=7,
vectors=None,
):
self.n = n
self.eta = eta
self.omega = omega
self.error = error
self.charge = charge
self.br = br
self.vectors = vectors
self.cubic = True
if omega is None and vectors is None:
raise ValueError("The lattice vectors must be provided.")
if omega is not None and vectors is not None:
raise ValueError(
"The lattice vectors and the unit cell volume should not be provided"
" together. Please provide only one of them."
)
if self.vectors is not None:
self.omega = np.abs(np.sum((np.cross(vectors[0], vectors[1]) * vectors[2])))
recip_vectors = ( # taken from Eq. (35) of arXiv:2302.07981v1 (2023)
2
* np.pi
/ self.omega
* np.array([np.cross(vectors[i], vectors[j]) for i, j in [(1, 2), (2, 0), (0, 1)]])
)
bbt = np.matrix(recip_vectors) @ np.matrix(recip_vectors).T
self.cubic = np.linalg.norm(bbt - (recip_vectors**2).max() * np.identity(3)) < 1e-6
self.lamb = self.norm(
self.n, self.eta, self.omega, self.error, self.br, self.charge, self.cubic, self.vectors
)
self.gates = self.gate_cost(
self.n, self.eta, self.omega, self.error, self.br, self.charge, self.cubic, self.vectors
)
self.qubits = self.qubit_cost(
self.n, self.eta, self.omega, self.error, self.br, self.charge, self.cubic, self.vectors
)
super().__init__(wires=range(self.qubits))
def _flatten(self):
return (self.n, self.eta), (
("omega", self.omega),
("error", self.error),
("charge", self.charge),
("br", self.br),
("vectors", self.vectors),
)
# pylint: disable=arguments-differ
@classmethod
def _primitive_bind_call(cls, n, eta, **kwargs):
return cls._primitive.bind(n, eta, **kwargs)
@classmethod
def _unflatten(cls, data, metadata):
return cls(*data, **dict(metadata))
[docs] @staticmethod
def success_prob(n, br):
r"""Return the probability of success for state preparation.
The expression for computing the probability of success is taken from Eqs. (59, 60) of
[`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_].
Args:
n (int): number of basis states to create an equal superposition for state preparation
br (int): number of bits for ancilla qubit rotation
Returns:
float: probability of success for state preparation
**Example**
>>> n = 3
>>> br = 8
>>> success_prob(n, br)
0.9999928850303523
"""
if n <= 0:
raise ValueError("The number of plane waves must be a positive number.")
if br <= 0 or not isinstance(br, int):
raise ValueError("br must be a positive integer.")
c = n / 2 ** np.ceil(np.log2(n))
d = 2 * np.pi / 2**br
theta = d * np.round((1 / d) * np.arcsin(np.sqrt(1 / (4 * c))))
p = c * ((1 + (2 - 4 * c) * np.sin(theta) ** 2) ** 2 + np.sin(2 * theta) ** 2)
return p
[docs] @staticmethod
def norm(n, eta, omega, error, br=7, charge=0, cubic=True, vectors=None):
r"""Return the 1-norm of a first-quantized Hamiltonian in the plane-wave basis.
The expressions needed for computing the norm are taken from
[`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_].
The norm is computed assuming that amplitude amplification is performed.
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
cubic (bool): True if the unit cell is cubic
vectors (array[float]): lattice vectors
Returns:
float: 1-norm of a first-quantized Hamiltonian in the plane-wave basis
**Example**
>>> n = 10000
>>> eta = 156
>>> omega = 1145.166
>>> error = 0.001
>>> norm(n, eta, omega, error)
281053.75612801575
.. details::
:title: Theory
To compute the norm, for numerical convenience, we use the following modified
expressions to obtain parameters that contain a sum over
:math:`\frac{1}{\left \| \nu \right \|^k}` where :math:`\nu` denotes an element of the
set of reciprocal lattice vectors, :math:`G_0`, and
:math:`k \in \left \{ 1, 2 \right \}`.
For :math:`\lambda_{\nu}` defined in Eq. (25) of
`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_
as
.. math::
\lambda_{\nu} = \sum_{\nu \in G_0} \frac{1}{\left \| \nu \right \|^2},
we follow Eq. (F6) of
`PRX 8, 011044 (2018) <https://journals.aps.org/prx/abstract/10.1103/PhysRevX.8.011044>`_
and use
.. math::
\lambda_{\nu} = 4\pi \left ( \frac{\sqrt{3}}{2} N^{1/3} - 1 \right) + 3 - \frac{3}{N^{1/3}}
+ 3 \int_{x=1}^{N^{1/3}} \int_{y=1}^{N^{1/3}} \frac{1}{x^2 + y^2} dydx.
We also need to compute :math:`\lambda^{\alpha}_{\nu}` defined in Eq. (123) of
`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_
.. math::
\lambda^{\alpha}_{\nu} = \alpha \sum_{\nu \in G_0} \frac{\left \lceil
\mathcal{M}(2^{\mu - 2}) / \left \| \nu \right \|^2 \right \rceil}{\mathcal{M}
2^{2\mu - 4}},
which we compute here, for :math:`\alpha = 1`, as
.. math::
\lambda^{1}_{\nu} = \lambda_{\nu} + \epsilon_l,
where :math:`\epsilon_l` is simply defined as the difference of
:math:`\lambda^{1}_{\nu}` and :math:`\lambda_{\nu}`. We follow Eq. (113) of
`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_
to derive an upper bound for its absolute value:
.. math::
|\epsilon_l| \le \frac{4}{2^{n_m}} (7 \times 2^{n_p + 1} + 9 n_p - 11 - 3 \times 2^{-n_p}),
where :math:`\mathcal{M} = 2^{n_m}` and :math:`n_m` is defined in Eq. (132) of
`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_.
Finally, for :math:`p_{\nu}` defined in Eq. (128) of
`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_
.. math::
p_{\nu} = \sum_{\mu = 2}^{n_p + 1} \sum_{\nu \in B_{\mu}} \frac{\left \lceil M(2^{\mu-2}
/ \left \| \nu \right \|)^2 \right \rceil}{M 2^{2\mu} 2^{n_{\mu} + 1}},
we use the upper bound from Eq. (29) in
`arXiv:1807.09802v2 <https://arxiv.org/abs/1807.09802v2>`_ which gives
:math:`p_{\nu} = 0.2398`.
"""
if n <= 0:
raise ValueError("The number of plane waves must be a positive number.")
if eta <= 0 or not isinstance(eta, (int, np.integer)):
raise ValueError("The number of electrons must be a positive integer.")
if omega <= 0:
raise ValueError("The unit cell volume must be a positive number.")
if error <= 0.0:
raise ValueError("The target error must be greater than zero.")
if br <= 0 or not isinstance(br, int):
raise ValueError("br must be a positive integer.")
if not isinstance(charge, int):
raise ValueError("system charge must be an integer.")
if not cubic:
return FirstQuantization._norm_noncubic(n, eta, error, br, charge, vectors)[0]
l_z = eta + charge
# target error in the qubitization of U+V which we set to be 0.01 of the algorithm error
error_uv = 0.01 * error
# taken from Eq. (22) of PRX Quantum 2, 040332 (2021)
n_p = int(np.ceil(np.log2(n ** (1 / 3) + 1)))
n0 = n ** (1 / 3)
lambda_nu = ( # expression is taken from Eq. (F6) of PRX 8, 011044 (2018)
4 * np.pi * (np.sqrt(3) * n ** (1 / 3) / 2 - 1)
+ 3
- 3 / n ** (1 / 3)
+ 3 * integrate.nquad(lambda x, y: 1 / (x**2 + y**2), [[1, n0], [1, n0]])[0]
)
n_m = int(
np.log2( # taken from Eq. (132) of PRX Quantum 2, 040332 (2021)
(2 * eta)
/ (error_uv * np.pi * omega ** (1 / 3))
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
# computed using Eq. (113) of PRX Quantum 2, 040332 (2021)
lambda_nu_1 = lambda_nu + 4 / 2**n_m * (
7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p)
)
p_nu = 0.2398 # approximation from Eq. (29) in arxiv:1807.09802
p_nu_amp = (
np.sin(3 * np.arcsin(np.sqrt(p_nu))) ** 2
) # Eq. (129), PRX Quantum 2, 040332 (2021)
# lambda_u and lambda_v are taken from Eq. (25) of PRX Quantum 2, 040332 (2021)
lambda_u = eta * l_z * lambda_nu / (np.pi * omega ** (1 / 3))
lambda_v = eta * (eta - 1) * lambda_nu / (2 * np.pi * omega ** (1 / 3))
# taken from Eq. (71) of PRX Quantum 2, 040332 (2021)
lambda_t_p = 6 * eta * np.pi**2 / omega ** (2 / 3) * 2 ** (2 * n_p - 2)
# lambda_u_1 and lambda_v_1 are taken from Eq. (124) of PRX Quantum 2, 040332 (2021)
lambda_u_1 = lambda_u * lambda_nu_1 / lambda_nu
lambda_v_1 = lambda_v * lambda_nu_1 / lambda_nu
# taken from Eq. (63) of PRX Quantum 2, 040332 (2021)
p_eq = (
FirstQuantization.success_prob(3, 8)
* FirstQuantization.success_prob(3 * eta + 2 * charge, br)
* FirstQuantization.success_prob(eta, br) ** 2
)
# final lambda value is computed from Eq. (126) of PRX Quantum 2, 040332 (2021)
lambda_a = lambda_t_p + lambda_u_1 + lambda_v_1
lambda_b = (lambda_u_1 + lambda_v_1 / (1 - 1 / eta)) / p_nu_amp
return np.maximum(lambda_a, lambda_b) / p_eq
@staticmethod
def _cost_qrom(lz):
r"""Return the minimum number of Toffoli gates needed for erasing the output of a QROM.
Args:
lz (int): sum of the atomic numbers
Returns:
int: the minimum cost of erasing the output of a QROM
**Example**
>>> lz = 100
>>> _cost_qrom(lz)
21
"""
if lz <= 0 or not isinstance(lz, (int, np.integer)):
raise ValueError("The sum of the atomic numbers must be a positive integer.")
k_f = np.floor(np.log2(lz) / 2)
k_c = np.ceil(np.log2(lz) / 2)
cost_f = int(2**k_f + np.ceil(2 ** (-1 * k_f) * lz))
cost_c = int(2**k_c + np.ceil(2 ** (-1 * k_c) * lz))
return min(cost_f, cost_c)
[docs] @staticmethod
def unitary_cost(n, eta, omega, error, br=7, charge=0):
r"""Return the number of Toffoli gates needed to implement the qubitization unitary
operator.
The expression for computing the cost is taken from Eq. (125) of
[`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_].
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
Returns:
int: the number of Toffoli gates needed to implement the qubitization unitary operator
**Example**
>>> n = 100000
>>> eta = 156
>>> omega = 169.69608
>>> error = 0.01
>>> unitary_cost(n, eta, omega, error)
17033
"""
if n <= 0:
raise ValueError("The number of plane waves must be a positive number.")
if eta <= 0 or not isinstance(eta, (int, np.integer)):
raise ValueError("The number of electrons must be a positive integer.")
if omega <= 0:
raise ValueError("The unit cell volume must be a positive number.")
if error <= 0.0:
raise ValueError("The target error must be greater than zero.")
if br <= 0 or not isinstance(br, int):
raise ValueError("br must be a positive integer.")
if not isinstance(charge, int):
raise ValueError("system charge must be an integer.")
lamb = FirstQuantization.norm(n, eta, omega, error, br, charge)
alpha = 0.01
l_z = eta + charge
l_nu = 2 * np.pi * n ** (2 / 3)
# defined in the third and second paragraphs of page 15 of PRX Quantum 2, 040332 (2021)
n_eta = np.ceil(np.log2(eta))
n_etaz = np.ceil(np.log2(eta + 2 * l_z))
# n_p is taken from Eq. (22)
n_p = int(np.ceil(np.log2(n ** (1 / 3) + 1)))
# errors in Eqs. (132-134) are set to be 0.01 of the algorithm error
error_t = alpha * error
error_r = alpha * error
error_m = alpha * error
# parameters taken from Eqs. (132-134) of PRX Quantum 2, 040332 (2021)
n_t = int(np.log2(np.pi * lamb / error_t)) # Eq. (134)
n_r = int(np.log2((eta * l_z * l_nu) / (error_r * omega ** (1 / 3)))) # Eq. (133)
n_m = int(
np.log2( # Eq. (132)
(2 * eta)
/ (error_m * np.pi * omega ** (1 / 3))
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
e_r = FirstQuantization._cost_qrom(l_z)
# taken from Eq. (125)
cost = 2 * (n_t + 4 * n_etaz + 2 * br - 12) + 14 * n_eta + 8 * br - 36
cost += 3 * (3 * n_p**2 + 15 * n_p - 7 + 4 * n_m * (n_p + 1))
cost += l_z + e_r + 2 * (2 * n_p + 2 * br - 7) + 12 * eta * n_p
cost += 5 * (n_p - 1) + 2 + 24 * n_p + 6 * n_p * n_r + 18
cost += n_etaz + 2 * n_eta + 6 * n_p + n_m + 16
return int(np.ceil(cost))
[docs] @staticmethod
def estimation_cost(n, eta, omega, error, br=7, charge=0, cubic=True, vectors=None):
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. (125) of
[`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_].
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
cubic (bool): True if the unit cell is cubic
vectors (array[float]): lattice vectors
Returns:
int: number of calls to unitary
**Example**
>>> n = 100000
>>> eta = 156
>>> omega = 1145.166
>>> error = 0.01
>>> estimation_cost(n, eta, omega, error)
102133985
"""
if error <= 0.0:
raise ValueError("The target error must be greater than zero.")
lamb = FirstQuantization.norm(
n, eta, omega, error, br=br, charge=charge, cubic=cubic, vectors=vectors
)
if cubic:
alpha = 0.01
# qpe_error obtained to satisfy inequality (131)
error_qpe = np.sqrt(error**2 * (1 - (3 * alpha) ** 2))
else:
alpha = 0.0248759298
# qpe_error obtained to satisfy inequality (131)
error_qpe = np.sqrt(error**2 * (1 - (4 * alpha) ** 2))
return int(np.ceil(np.pi * lamb / (2 * error_qpe)))
[docs] @staticmethod
def gate_cost(n, eta, omega, error, br=7, charge=0, cubic=True, vectors=None):
r"""Return the total number of Toffoli gates needed to implement the first quantization
algorithm.
The expression for computing the cost is taken from Eq. (125) of
[`PRX Quantum 2, 040332 (2021) <https://link.aps.org/doi/10.1103/PRXQuantum.2.040332>`_].
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
cubic (bool): True if the unit cell is cubic
vectors (array[float]): lattice vectors
Returns:
int: the number of Toffoli gates needed to implement the first quantization algorithm
**Example**
>>> n = 100000
>>> eta = 156
>>> omega = 169.69608
>>> error = 0.01
>>> gate_cost(n, eta, omega, error)
3676557345574
"""
if n <= 0:
raise ValueError("The number of plane waves must be a positive number.")
if eta <= 0 or not isinstance(eta, (int, np.integer)):
raise ValueError("The number of electrons must be a positive integer.")
if omega <= 0:
raise ValueError("The unit cell volume must be a positive number.")
if error <= 0.0:
raise ValueError("The target error must be greater than zero.")
if not isinstance(charge, int):
raise ValueError("system charge must be an integer.")
if br <= 0 or not isinstance(br, int):
raise ValueError("br must be a positive integer.")
e_cost = FirstQuantization.estimation_cost(
n, eta, omega, error, br=br, charge=charge, cubic=cubic, vectors=vectors
)
if cubic:
u_cost = FirstQuantization.unitary_cost(n, eta, omega, error, br, charge)
else:
u_cost = FirstQuantization._unitary_cost_noncubic(n, eta, error, br, charge, vectors)
return e_cost * u_cost
[docs] @staticmethod
def qubit_cost(n, eta, omega, error, br=7, charge=0, cubic=True, vectors=None):
r"""Return the number of logical qubits needed to implement the first quantization
algorithm.
The expression for computing the cost is taken from Eq. (101) of
[`arXiv:2204.11890v1 <https://arxiv.org/abs/2204.11890v1>`_].
Args:
n (int): number of plane waves
eta (int): number of electrons
omega (float): unit cell volume
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
cubic (bool): True if the unit cell is cubic
vectors (array[float]): lattice vectors
Returns:
int: number of logical qubits needed to implement the first quantization algorithm
**Example**
>>> n = 100000
>>> eta = 156
>>> omega = 169.69608
>>> error = 0.01
>>> qubit_cost(n, eta, omega, error)
4377
"""
if n <= 0:
raise ValueError("The number of plane waves must be a positive number.")
if eta <= 0 or not isinstance(eta, (int, np.integer)):
raise ValueError("The number of electrons must be a positive integer.")
if omega <= 0:
raise ValueError("The unit cell volume must be a positive number.")
if error <= 0.0:
raise ValueError("The target error must be greater than zero.")
if not isinstance(charge, int):
raise ValueError("system charge must be an integer.")
if not cubic:
return FirstQuantization._qubit_cost_noncubic(n, eta, error, br, charge, vectors)
lamb = FirstQuantization.norm(n, eta, omega, error, br=br, charge=charge)
alpha = 0.01
l_z = eta + charge
l_nu = 2 * np.pi * n ** (2 / 3)
# n_p is taken from Eq. (22) of PRX Quantum 2, 040332 (2021)
n_p = np.ceil(np.log2(n ** (1 / 3) + 1))
# errors in Eqs. (132-134) of PRX Quantum 2, 040332 (2021),
# set to 0.01 of the algorithm error
error_t = alpha * error
error_r = alpha * error
error_m = alpha * error
# parameters taken from Eqs. (132-134) of PRX Quantum 2, 040332 (2021)
n_t = int(np.log2(np.pi * lamb / error_t)) # Eq. (134)
n_r = int(np.log2((eta * l_z * l_nu) / (error_r * omega ** (1 / 3)))) # Eq. (133)
n_m = int(
np.log2( # Eq. (132)
(2 * eta)
/ (error_m * np.pi * omega ** (1 / 3))
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
alpha = 0.01
# qpe_error obtained to satisfy inequality (131) of PRX Quantum 2, 040332 (2021)
error_qpe = np.sqrt(error**2 * (1 - (3 * alpha) ** 2))
# the expression for computing the cost is taken from Eq. (101) of arXiv:2204.11890v1
qubits = 3 * eta * n_p + 4 * n_m * n_p + 12 * n_p
qubits += 2 * (np.ceil(np.log2(np.ceil(np.pi * lamb / (2 * error_qpe))))) + 5 * n_m
qubits += 2 * np.ceil(np.log2(eta)) + 3 * n_p**2 + np.ceil(np.log2(eta + 2 * l_z))
qubits += np.maximum(5 * n_p + 1, 5 * n_r - 4) + np.maximum(n_t, n_r + 1) + 33
return int(np.ceil(qubits))
@staticmethod
def _norm_noncubic(
n,
eta,
error,
br,
charge,
vectors,
):
r"""Return the 1-norm of a first-quantized Hamiltonian in the plane-wave basis
for non-cubic systems.
Args:
n (int): number of plane waves
eta (int): number of electrons
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
vectors (array[float]): lattice vectors
Returns:
float: 1-norm of a first-quantized Hamiltonian in the plane-wave basis
"""
omega = np.abs(np.sum((np.cross(vectors[0], vectors[1]) * vectors[2])))
recip_vectors = ( # taken from Eq. (35) of arXiv:2302.07981v1 (2023)
2
* np.pi
/ omega
* np.array([np.cross(vectors[i], vectors[j]) for i, j in [(1, 2), (2, 0), (0, 1)]])
)
bbt = np.matrix(recip_vectors) @ np.matrix(recip_vectors).T
orthogonal = (
np.linalg.norm(bbt - np.array([np.max(b**2) for b in recip_vectors]) * np.identity(3))
< 1e-6
)
l_z = eta + charge
alpha = 0.0248759298 # optimal value for lower resource estimates
p_th = 0.95 # optimal value for lower resource estimates
error_uv = alpha * error
# taken from Eq. (22) of PRX Quantum 2, 040332 (2021)
n_p = int(np.ceil(np.log2(n ** (1 / 3) + 1)))
n0 = n ** (1 / 3)
# defined in Eq. (F3) of arXiv:2302.07981v1 (2023)
bmin = np.min(np.linalg.svd(recip_vectors)[1])
n_m = int(
np.ceil(
np.log2( # taken from Eq. (132) of PRX Quantum 2, 040332 (2021) with
# modifications taken from arXiv:2302.07981v1 (2023)
(8 * np.pi * eta)
/ (error_uv * omega * bmin**2)
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
)
lambda_nu = ( # expression is taken from Eq. (F6) of PRX 8, 011044 (2018)
# bmin**2 is needed for non-cubic systems, from Eq. (F3) of arXiv:2302.07981v1 (2023)
4 * np.pi * (np.sqrt(3) * n ** (1 / 3) / 2 - 1)
+ 3
- 3 / n ** (1 / 3)
+ 3 * integrate.nquad(lambda x, y: 1 / (x**2 + y**2), [[1, n0], [1, n0]])[0]
) / bmin**2
# computed using error term derived in Eq. (113) of PRX Quantum 2, 040332 (2021)
# bmin**2 is needed for non-cubic systems, from Eq. (F3) of arXiv:2302.07981v1 (2023)
lambda_nu_1 = lambda_nu + 4 / (2**n_m * bmin**2) * (
7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p)
)
# taken from Eq. (I7) of arXiv:2302.07981v1 (2023)
p_nu = lambda_nu_1 * bmin**2 / 2 ** (n_p + 6)
# amplitude amplification steps, taken from Eq. (117) of PRX Quantum 2, 040332 (2021)
p_nu_amp = 0
aa_steps = 0
n_steps = 30 # maximum number of amplitude amplification steps
for i in reversed(range(n_steps)):
probability = (np.sin((2 * i + 1) * np.arcsin(np.sqrt(p_nu)))) ** 2
if probability > p_th:
aa_steps = i
p_nu_amp = probability
# lambda_u and lambda_v are taken from Eq. (25) of PRX Quantum 2, 040332 (2021) with
# modifications taken from arXiv:2302.07981v1 (2023)
lambda_u = 4 * np.pi * eta * l_z * lambda_nu / omega
lambda_v = 2 * np.pi * eta * (eta - 1) * lambda_nu / omega
# lambda_u_1 and lambda_v_1 are taken from Eq. (124) of PRX Quantum 2, 040332 (2021)
lambda_u_1 = lambda_u * lambda_nu_1 / lambda_nu
lambda_v_1 = lambda_v * lambda_nu_1 / lambda_nu
b_mat = np.matrix(recip_vectors)
abs_sum = np.abs(b_mat @ b_mat.T).flatten().sum()
# taken from Appendix I.1 of arXiv:2302.07981v1 (2023)
if orthogonal:
lambda_t_p = abs_sum * eta * 2 ** (2 * n_p - 2) / 4
else:
lambda_t_p = abs_sum * eta * 2 ** (2 * n_p - 2) / 2
# taken from Eq. (63) of PRX Quantum 2, 040332 (2021) with
# modifications taken from arXiv:2302.07981v1 (2023)
p_eq = (
FirstQuantization.success_prob(3 * eta + 2 * charge, br)
* FirstQuantization.success_prob(eta, br) ** 2
)
# final lambda value is computed from case checking derived in Section E and Appendix D of
# PRX Quantum 2, 040332 (2021)
if p_nu * lambda_t_p >= (1 - p_nu) * (lambda_u_1 + lambda_v_1):
raise ValueError("The computed 1-norm is zero.") # return 0.0, 0.0
if p_nu_amp * lambda_t_p >= (1 - p_nu_amp) * (lambda_u_1 + lambda_v_1):
return (lambda_t_p + lambda_u_1 + lambda_v_1) / p_eq, aa_steps
return ((lambda_u_1 + lambda_v_1 / (1 - 1 / eta)) / p_nu_amp) / p_eq, aa_steps
@staticmethod
def _qubit_cost_noncubic(n, eta, error, br, charge, vectors):
r"""Return the number of logical qubits needed to implement the first quantization
algorithm for non-cubic systems.
Args:
n (int): number of plane waves
eta (int): number of electrons
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
vectors (array[float]): lattice vectors
Returns:
int: number of logical qubits needed to implement the first quantization algorithm
"""
lambda_total, _ = FirstQuantization._norm_noncubic(
n,
eta,
error,
br,
charge,
vectors,
)
omega = np.abs(np.sum((np.cross(vectors[0], vectors[1]) * vectors[2])))
recip_vectors = ( # taken from Eq. (35) of arXiv:2302.07981v1 (2023)
2
* np.pi
/ omega
* np.array([np.cross(vectors[i], vectors[j]) for i, j in [(1, 2), (2, 0), (0, 1)]])
)
l_z = eta + charge
l_nu = 2 * np.pi * n ** (2 / 3)
# taken from Eq. (22) of PRX Quantum 2, 040332 (2021)
n_p = np.ceil(np.log2(n ** (1 / 3) + 1))
# defined in Eq. (F3) of arXiv:2302.07981v1 (2023)
bmin = np.min(np.linalg.svd(recip_vectors)[1])
alpha = 0.0248759298 # optimal value for lower resource estimates
# errors derived according to the strategy in Appendix M.1 of arXiv:2302.07981v1 (2023)
error_t, error_r, error_m, error_b = [alpha * error] * 4
# parameters taken from PRX Quantum 2, 040332 (2021)
n_t = int(np.ceil(np.log2(np.pi * lambda_total / error_t))) # Eq. (134)
n_r = int(np.ceil(np.log2((eta * l_z * l_nu) / (error_r * omega ** (1 / 3))))) # Eq. (133)
n_m = int( # taken from Eq. (J13) arXiv:2302.07981v1 (2023)
np.ceil(
np.log2(
(8 * np.pi * eta)
/ (error_m * omega * bmin**2)
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
)
# qpe error obtained to satisfy a modification of
# inequality Eq. (131) of PRX Quantum 2, 040332 (2021) including error_b
n_errors = 4
error_qpe = np.sqrt(error**2 * (1 - (n_errors * alpha) ** 2))
# adapted from equations (L1, L2) in Appendix L of arXiv:2302.07981v1 (2023)
clean_temp_H_cost = max([5 * n_r - 4, 5 * n_p + 1]) + max([5, n_m + 3 * n_p])
reflection_cost = (
np.ceil(np.log2(eta + 2 * l_z)) + 2 * np.ceil(np.log2(eta)) + 6 * n_p + n_m + 16 + 3
)
clean_temp_cost = max([clean_temp_H_cost, reflection_cost])
# the expression for computing the cost is taken from Appendix C of
# PRX Quantum 2, 040332 (2021) and adapting to non-cubic using Appendix L of
# arXiv:2302.07981v1 (2023)
clean_cost = 3 * eta * n_p
clean_cost += np.ceil(np.log2(np.ceil(np.pi * lambda_total / (2 * error_qpe))))
clean_cost += 1 + 1 + np.ceil(np.log2(eta + 2 * l_z)) + 3 + 3
clean_cost += 2 * np.ceil(np.log2(eta)) + 5 + 3 * (n_p + 1)
clean_cost += n_p + n_m + 3 * n_p + 2 + 2 * n_p + 1 + 1 + 2 + 2 * n_p + 6 + 1
clean_cost += clean_temp_cost
# taken from Eq. (J7) of arXiv:2302.07981v1 (2023)
n_b = np.ceil(
np.log2(
4
* np.pi
* eta
* 2 ** (2 * n_p - 2)
* np.abs(np.matrix(recip_vectors) @ np.matrix(recip_vectors).T).flatten().sum()
/ error_b
)
)
clean_cost += np.max([n_r + 1, n_t, n_b]) + 6 + n_m + 1
return int(np.ceil(clean_cost))
@staticmethod
def _unitary_cost_noncubic(n, eta, error, br, charge, vectors):
r"""Return the number of Toffoli gates needed to implement the qubitization unitary
operator for non-cubic systems.
Args:
n (int): number of plane waves
eta (int): number of electrons
error (float): target error in the algorithm
br (int): number of bits for ancilla qubit rotation
charge (int): total electric charge of the system
vectors (array[float]): lattice vectors
Returns:
int: the number of Toffoli gates needed to implement the qubitization unitary operator
"""
lambda_total, aa_steps = FirstQuantization._norm_noncubic(
n,
eta,
error,
br,
charge,
vectors,
)
omega = np.abs(np.sum((np.cross(vectors[0], vectors[1]) * vectors[2])))
recip_vectors = ( # taken from Eq. (35) of arXiv:2302.07981v1 (2023)
2
* np.pi
/ omega
* np.array([np.cross(vectors[i], vectors[j]) for i, j in [(1, 2), (2, 0), (0, 1)]])
)
n_tof = 500 # optimal value for lower resource estimates
alpha = 0.0248759298 # optimal value for lower resource estimates
l_z = eta + charge
l_nu = 2 * np.pi * n ** (2 / 3)
# defined in the third and second paragraphs of page 15 of PRX Quantum 2, 040332 (2021)
n_eta = np.ceil(np.log2(eta))
n_etaz = np.ceil(np.log2(eta + 2 * l_z))
# taken from Eq. (22) of PRX Quantum 2, 040332 (2021)
n_p = int(np.ceil(np.log2(n ** (1 / 3) + 1)))
# errors in Eqs. (132-134) of PRX Quantum 2, 040332 (2021)
error_t, error_r, error_m = [alpha * error] * 3
# parameters taken from PRX Quantum 2, 040332 (2021)
n_t = int(np.ceil(np.log2(np.pi * lambda_total / error_t))) # Eq. (134)
n_r = int(np.ceil(np.log2((eta * l_z * l_nu) / (error_r * omega ** (1 / 3))))) # Eq. (133)
# defined in Eq. (F3) of arXiv:2302.07981v1 (2023)
bmin = np.min(np.linalg.svd(recip_vectors)[1])
# equivalent to Eq. (J13) of arXiv:2302.07981v1 (2023)
n_m = int(
np.ceil(
np.log2(
(8 * np.pi * eta)
/ (error_m * omega * bmin**2)
* (eta - 1 + 2 * l_z)
* (7 * 2 ** (n_p + 1) - 9 * n_p - 11 - 3 * 2 ** (-1 * n_p))
)
)
)
e_r = FirstQuantization._cost_qrom(l_z)
# taken from Eq. (125) of PRX Quantum 2, 040332 (2021)
cost = 2 * (n_t + 4 * n_etaz + 2 * br - 12) + 14 * n_eta + 8 * br - 36
cost += (2 * aa_steps + 1) * (3 * n_p**2 + 15 * n_p - 7 + 4 * n_m * (n_p + 1))
cost += l_z + e_r + 2 * (2 * n_p + 2 * br - 7) + 12 * eta * n_p
cost += 5 * (n_p - 1) + 2 + 24 * n_p + 6 * n_p * n_r + 18
cost += n_etaz + 2 * n_eta + 6 * n_p + n_m + 16
# taken from Eq. (J7) of arXiv:2302.07981v1 (2023)
error_b = alpha * error
n_b = np.ceil(
np.log2(
2
* np.pi
* eta
* 2 ** (2 * n_p - 2)
* np.abs(np.matrix(recip_vectors) @ np.matrix(recip_vectors).T).flatten().sum()
/ error_b
)
)
n_dirty = FirstQuantization._qubit_cost_noncubic(n, eta, error, br, charge, vectors)
ms_cost = FirstQuantization._momentum_state_qrom(n_p, n_m, n_dirty, n_tof, kappa=1)[0]
# adapted Eq. (125) of PRX Quantum 2, 040332 (2021) to the noncubic case using Appendix K of
# arXiv:2302.07981v1 (2023)
cost -= 2 * (3 * 2 + 2 * br - 9)
cost += 2 * (2 * (2 * (2 ** (4 + 1) - 1) + (n_b - 3) * 4 + 2**4 + (n_p - 2)))
cost += 8
cost -= (2 * aa_steps + 1) * (3 * n_p**2 + 15 * n_p - 7 + 4 * n_m * (n_p + 1))
cost += (2 * aa_steps + 1) * ms_cost
return int(np.ceil(cost))
@staticmethod
def _momentum_state_qrom(n_p, n_m, n_dirty, n_tof, kappa):
r"""Returns the Toffoli cost for preparing the momentum state superposition.
Derived from Section D.1 item (6) and Appendix K.1.f of arXiv:2302.07981v1 (2023)"""
x = 2 ** (3 * n_p)
beta_dirty = max([np.floor(n_dirty / n_m), 1])
beta_parallel = max([np.floor(n_tof / kappa), 1])
if n_tof == 1:
beta_gate = max([np.floor(np.sqrt(2 * x / (3 * n_m))), 1])
beta = np.min([beta_dirty, beta_gate])
ms_cost_qrom = 2 * np.ceil(x / beta) + 3 * n_m * beta
else:
beta_gate = max([np.floor(2 * x / (3 * n_m / kappa) * np.log(2)), 1])
beta = np.min([beta_dirty, beta_gate, beta_parallel])
ms_cost_qrom = 2 * np.ceil(x / beta) + 3 * np.ceil(n_m / kappa) * np.ceil(np.log2(beta))
ms_cost = 2 * ms_cost_qrom + n_m + 8 * (n_p - 1) + 6 * n_p + 2 + 2 * n_p + n_m + 2
return ms_cost, beta
_modules/pennylane/resource/first_quantization
Download Python script
Download Notebook
View on GitHub