# Source code for pennylane.qchem.factorization

# Copyright 2018-2022 Xanadu Quantum Technologies Inc.

# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
"""
This module contains the functions needed for two-electron tensor factorization.
"""
import pennylane as qml
from pennylane import numpy as np

[docs]def factorize(two_electron, tol_factor=1.0e-5, tol_eigval=1.0e-5):
r"""Return the double-factorized form of a two-electron integral tensor in spatial basis.

The two-electron tensor :math:V, in
chemist notation <http://vergil.chemistry.gatech.edu/notes/permsymm/permsymm.pdf>_, is first
factorized in terms of symmetric matrices :math:L^{(r)} such that
:math:V_{ijkl} = \sum_r^R L_{ij}^{(r)} L_{kl}^{(r) T}. The rank :math:R is determined by a
threshold error. Then, each matrix :math:L^{(r)} is diagonalized and its eigenvalues (and
corresponding eigenvectors) are truncated at a threshold error.

Args:
two_electron (array[array[float]]): two-electron integral tensor in the molecular orbital
basis arranged in chemist notation
tol_factor (float): threshold error value for discarding the negligible factors
tol_eigval (float): threshold error value for discarding the negligible factor eigenvalues

Returns:
tuple(array[array[float]], list[array[float]], list[array[float]]): tuple containing
symmetric matrices (factors) approximating the two-electron integral tensor, truncated
eigenvalues of the generated factors, and truncated eigenvectors of the generated factors

**Example**

>>> symbols  = ['H', 'H']
>>> geometry = np.array([[0.0, 0.0, 0.0],
>>> mol = qml.qchem.Molecule(symbols, geometry)
>>> core, one, two = qml.qchem.electron_integrals(mol)()
>>> two = np.swapaxes(two, 1, 3) # convert to chemist notation
>>> factors, eigvals, eigvecs = factorize(two, 1e-5, 1e-5)
>>> print(factors)
[[[ 1.06723440e-01  9.73575768e-15]
[ 8.36288956e-15 -1.04898533e-01]]
[[-2.20945401e-13 -4.25688222e-01]
[-4.25688222e-01 -2.98228790e-13]]
[[-8.14472856e-01  5.01669019e-13]
[ 5.01689072e-13 -8.28642140e-01]]]

.. details::
:title: Theory

The second quantized electronic Hamiltonian is constructed in terms of fermionic creation,
:math:a^{\dagger} , and annihilation, :math:a, operators as
[arXiv:1902.02134 <https://arxiv.org/abs/1902.02134>_]

.. math::

H = \sum_{\alpha \in \{\uparrow, \downarrow \} } \sum_{pq} h_{pq} a_{p,\alpha}^{\dagger}
a_{q, \alpha} + \frac{1}{2} \sum_{\alpha, \beta \in \{\uparrow, \downarrow \} } \sum_{pqrs}
h_{pqrs} a_{p, \alpha}^{\dagger} a_{q, \beta}^{\dagger} a_{r, \beta} a_{s, \alpha},

where :math:h_{pq} and :math:h_{pqrs} are the one- and two-electron integrals computed
as

.. math::

h_{pq} = \int \phi_p(r)^* \left ( -\frac{\nabla_r^2}{2} - \sum_i \frac{Z_i}{|r-R_i|} \right)
\phi_q(r) dr,

and

.. math::

h_{pqrs} = \int \frac{\phi_p(r_1)^* \phi_q(r_2)^* \phi_r(r_2) \phi_s(r_1)}{|r_1 - r_2|}
dr_1 dr_2.

The two-electron integrals can be rearranged in the so-called chemist notation which gives

.. math::

V_{pqrs} = \int \frac{\phi_p(r_1)^* \phi_q(r_1)^* \phi_r(r_2) \phi_s(r_2)}{|r_1 - r_2|}
dr_1 dr_2,

and the molecular Hamiltonian can be rewritten as

.. math::

H = \sum_{\alpha \in \{\uparrow, \downarrow \} } \sum_{pq} T_{pq} a_{p,\alpha}^{\dagger}
a_{q, \alpha} + \frac{1}{2} \sum_{\alpha, \beta \in \{\uparrow, \downarrow \} } \sum_{pqrs}
V_{pqrs} a_{p, \alpha}^{\dagger} a_{q, \alpha} a_{r, \beta}^{\dagger} a_{s, \beta},

with

.. math::

T_{pq} = h_{pq} - \frac{1}{2} \sum_s h_{pssq}.

This notation allows a low-rank factorization of the two-electron integral. The objective of
the factorization is to find a set of symmetric matrices, :math:L^{(r)}, such that

.. 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 factorization algorithm has the following steps
[arXiv:1902.02134 <https://arxiv.org/abs/1902.02134>_]:

- Reshape the :math:n \times n \times n \times n two-electron tensor to a
:math:n^2 \times n^2 matrix where :math:n is the number of orbitals.

- Diagonalize the resulting matrix and keep the :math:r eigenvectors that have
corresponding eigenvalues larger than a threshold.

- Multiply the eigenvectors by the square root of the eigenvalues to obtain
matrices :math:L^{(r)}.

- Reshape the selected eigenvectors to :math:n \times n matrices.

- Diagonalize the :math:n \times n matrices and for each matrix keep the eigenvalues (and
their corresponding eigenvectors) that are larger than a threshold.
"""
shape = two_electron.shape

if len(shape) != 4 or len(set(shape)) != 1:
raise ValueError("The two-electron repulsion tensor must have a (N x N x N x N) shape.")

n = shape
two = two_electron.reshape(n * n, n * n)

eigvals_r, eigvecs_r = np.linalg.eigh(two)
eigvals_r = np.array([val for val in eigvals_r if abs(val) > tol_factor])

eigvecs_r = eigvecs_r[:, -len(eigvals_r) :]

if eigvals_r.size == 0:
raise ValueError(
"All factors are discarded. Consider decreasing the first threshold error."
)

vectors = eigvecs_r @ np.diag(np.sqrt(eigvals_r))

r = len(eigvals_r)
factors = np.array([vectors.reshape(n, n, r)[:, :, k] for k in range(r)])

eigvals, eigvecs = np.linalg.eigh(factors)
eigvals_m = []
eigvecs_m = []
for n, eigval in enumerate(eigvals):
idx = [i for i, v in enumerate(eigval) if abs(v) > tol_eigval]
eigvals_m.append(eigval[idx])
eigvecs_m.append(eigvecs[n][idx])

if np.sum([len(v) for v in eigvecs_m]) == 0:
raise ValueError(
"All eigenvectors are discarded. Consider decreasing the second threshold error."
)

return factors, eigvals_m, eigvecs_m

[docs]def basis_rotation(one_electron, two_electron, tol_factor=1.0e-5):
r"""Return the grouped coefficients and observables of a molecular Hamiltonian and the basis
rotation unitaries obtained with the basis rotation grouping method.

Args:
one_electron (array[float]): one-electron integral matrix in the molecular orbital basis
two_electron (array[array[float]]): two-electron integral tensor in the molecular orbital
basis arranged in chemist notation
tol_factor (float): threshold error value for discarding the negligible factors

Returns:
tuple(list[array[float]], list[list[Observable]], list[array[float]]): tuple containing
grouped coefficients, grouped observables and basis rotation transformation matrices

**Example**

>>> symbols  = ['H', 'H']
>>> geometry = np.array([[0.0, 0.0, 0.0],
>>> mol = qml.qchem.Molecule(symbols, geometry)
>>> core, one, two = qml.qchem.electron_integrals(mol)()
>>> coeffs, ops, unitaries = basis_rotation(one, two, tol_factor=1.0e-5)
>>> print(coeffs)
[array([ 0.84064649, -2.59579282,  0.84064649,  0.45724992,  0.45724992]),
array([ 9.57150297e-05,  5.60006390e-03,  9.57150297e-05,  2.75092558e-03,
-9.73801723e-05, -2.79878310e-03, -9.73801723e-05, -2.79878310e-03,
-2.79878310e-03, -2.79878310e-03,  2.84747318e-03]),
array([ 0.04530262, -0.04530262, -0.04530262, -0.04530262, -0.04530262,
0.09060523,  0.04530262]),
array([-0.66913628,  1.6874169 , -0.66913628,  0.16584151, -0.68077716,
0.16872663, -0.68077716,  0.16872663,  0.16872663,  0.16872663,
0.17166195])]

.. details::
:title: Theory

A second-quantized molecular Hamiltonian can be constructed in the
chemist notation <http://vergil.chemistry.gatech.edu/notes/permsymm/permsymm.pdf>_ format
following Eq. (1) of
[PRX Quantum 2, 030305, 2021 <https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.030305>_]
as

.. math::

H = \sum_{\alpha \in \{\uparrow, \downarrow \} } \sum_{pq} T_{pq} a_{p,\alpha}^{\dagger}
a_{q, \alpha} + \frac{1}{2} \sum_{\alpha, \beta \in \{\uparrow, \downarrow \} } \sum_{pqrs}
V_{pqrs} a_{p, \alpha}^{\dagger} a_{q, \alpha} a_{r, \beta}^{\dagger} a_{s, \beta},

where :math:V_{pqrs} denotes a two-electron integral in the chemist notation and
:math:T_{pq} is obtained from the one- and two-electron integrals, :math:h_{pq} and
:math:h_{pqrs}, as

.. math::

T_{pq} = h_{pq} - \frac{1}{2} \sum_s h_{pssq}.

The tensor :math:V can be converted to a matrix which is indexed by the indices :math:pq
and :math:rs and eigendecomposed up to a rank :math:R to give

.. math::

V_{pqrs} = \sum_r^R L_{pq}^{(r)} L_{rs}^{(r) T},

where :math:L denotes the matrix of eigenvectors of the matrix :math:V. The molecular
Hamiltonian can then be rewritten following Eq. (7) of
[Phys. Rev. Research 3, 033055, 2021 <https://journals.aps.org/prresearch/abstract/10.1103/PhysRevResearch.3.033055>_]
as

.. math::

H = \sum_{\alpha \in \{\uparrow, \downarrow \} } \sum_{pq} T_{pq} a_{p,\alpha}^{\dagger}
a_{q, \alpha} + \frac{1}{2} \sum_r^R \left ( \sum_{\alpha \in \{\uparrow, \downarrow \} } \sum_{pq}
L_{pq}^{(r)} a_{p, \alpha}^{\dagger} a_{q, \alpha} \right )^2.

The orbital basis can be rotated such that each :math:T and :math:L^{(r)} matrix is
diagonal. The Hamiltonian can then be written following Eq. (2) of
[npj Quantum Information, 7, 23 (2021) <https://www.nature.com/articles/s41534-020-00341-7>_]
as

.. math::

H = U_0 \left ( \sum_p d_p n_p \right ) U_0^{\dagger} + \sum_r^R U_r \left ( \sum_{pq}
d_{pq}^{(r)} n_p n_q \right ) U_r^{\dagger},

where the coefficients :math:d are obtained by diagonalizing the :math:T and
:math:L^{(r)} matrices. The number operators :math:n_p = a_p^{\dagger} a_p can be
converted to qubit operators using

.. math::

n_p = \frac{1-Z_p}{2},

where :math:Z_p is the Pauli :math:Z operator applied to qubit :math:p. This gives
the qubit Hamiltonian

.. math::

H = U_0 \left ( \sum_p O_p^{(0)} \right ) U_0^{\dagger} + \sum_r^R U_r \left ( \sum_{q} O_q^{(r)} \right ) U_r^{\dagger},

where :math:O = \sum_i c_i P_i is a linear combination of Pauli words :math:P_i that are
a tensor product of Pauli :math:Z and Identity operators. This allows all the Pauli words
in each of the :math:O terms to be measured simultaneously. This function returns the
coefficients and the Pauli words grouped for each of the :math:O terms as well as the
basis rotation transformation matrices that are constructed from the eigenvectors of the
:math:T and :math:L^{(r)} matrices. Each column of the transformation matrix is an
eigenvector of the corresponding :math:T or :math:L^{(r)} matrix.
"""

num_orbitals = one_electron.shape * 2
one_body_tensor, chemist_two_body_tensor = _chemist_transform(one_electron, two_electron)
chemist_one_body_tensor = np.kron(one_body_tensor, np.eye(2))  # account for spin
t_eigvals, t_eigvecs = np.linalg.eigh(chemist_one_body_tensor)

factors, _, _ = factorize(chemist_two_body_tensor, tol_factor=tol_factor)
factors = [np.kron(factor, np.eye(2)) for factor in factors]  # account for spin

v_coeffs, v_unitaries = np.linalg.eigh(factors)
indices = [np.argsort(v_coeff)[::-1] for v_coeff in v_coeffs]
v_coeffs = [v_coeff[indices[idx]] for idx, v_coeff in enumerate(v_coeffs)]
v_unitaries = [v_unitary[:, indices[idx]] for idx, v_unitary in enumerate(v_unitaries)]

ops_t = 0.0
for p in range(num_orbitals):
ops_t += 0.5 * t_eigvals[p] * (qml.Identity(p) - qml.PauliZ(p))

ops_l = []
for idx in range(len(factors)):
ops_l_ = 0.0
for p in range(num_orbitals):
for q in range(num_orbitals):
ops_l_ += (
v_coeffs[idx][p]
* v_coeffs[idx][q]
* 0.25
* (
qml.Identity(p)
- qml.PauliZ(p)
- qml.PauliZ(q)
+ qml.pauli.pauli_mult_with_phase(qml.PauliZ(p), qml.PauliZ(q))
)
)
ops_l.append(ops_l_)

ops = [ops_t] + ops_l
c_group = [op.coeffs for op in ops]
o_group = [op.ops for op in ops]
u_transform = list([t_eigvecs] + list(v_unitaries))  # Inverse of diagonalizing unitaries

return c_group, o_group, u_transform

def _chemist_transform(one_body_tensor=None, two_body_tensor=None, spatial_basis=True):
r"""Transforms one- and two-body terms in physicists' notation to chemists' notation <http://vergil.chemistry.gatech.edu/notes/permsymm/permsymm.pdf>_\ .

This converts the input two-body tensor :math:h_{pqrs} that constructs :math:\sum_{pqrs} h_{pqrs} a^\dagger_p a^\dagger_q a_r a_s
to a transformed two-body tensor :math:V_{pqrs} that follows the chemists' convention to construct :math:\sum_{pqrs} V_{pqrs} a^\dagger_p a_q a^\dagger_r a_s
in the spatial basis. During the tranformation, some extra one-body terms come out. These are returned as a one-body tensor :math:T_{pq} in the
chemists' notation either as is or after summation with the input one-body tensor :math:h_{pq}, if provided.

Args:
one_body_tensor (array[float]): a one-electron integral tensor giving the :math:h_{pq}.
two_body_tensor (array[float]): a two-electron integral tensor giving the :math:h_{pqrs}.
spatial_basis (bool): True if the integral tensor are passed in spatial-orbital basis. False if they are in spin basis.

Returns:
tuple(array[float], array[float]) or tuple(array[float],): transformed one-body tensor :math:T_{pq} and two-body tensor :math:V_{pqrs} for the provided terms.

**Example**

>>> symbols  = ['H', 'H']
>>> geometry = np.array([[0.0, 0.0, 0.0],
>>> mol = qml.qchem.Molecule(symbols, geometry)
>>> core, one, two = qml.qchem.electron_integrals(mol)()
>>> qml.qchem.factorization._chemist_transform(two_body_tensor=two, spatial_basis=True)
(tensor([[-0.427983, -0.      ],
tensor([[[[0.337378, 0.      ],
[0.       , 0.331856]],
[[0.      , 0.090605],
[0.090605 , 0.      ]]],
[[[0.      , 0.090605],
[0.090605 , 0.      ]],
[[0.331856, 0.      ],

.. details::
:title: Theory

The two-electron integral in physicists' notation is defined as:

.. math::

\langle pq \vert rs \rangle = h_{pqrs} = \int \frac{\chi^*_{p}(x_1) \chi^*_{q}(x_2) \chi_{r}(x_1) \chi_{s}(x_2)}{|r_1 - r_2|} dx_1 dx_2,

while in chemists' notation it is written as:

.. math::

[pq \vert rs] = V_{pqrs} = \int \frac{\chi^*_{p}(x_1) \chi_{q}(x_1) \chi^*_{r}(x_2) \chi_{s}(x_2)}{|r_1 - r_2|} dx_1 dx_2.

In the spin basis, this index reordering :math:pqrs \rightarrow psrq leads to formation of one-body terms :math:h_{prrs} that come out during
the coversion:

.. math::

h_{prrs} = \int \frac{\chi^*_{p}(x_1) \chi^*_{r}(x_2) \chi_{r}(x_1) \chi_{s}(x_2)}{|x_1 - x_2|} dx_1 dx_2,

where both :math:\chi_{r}(x_1) and :math:\chi_{r}(x_2) will have same spin functions, i.e.,
:math:\chi_{r}(x_i) = \phi(r_i)\alpha(\omega) or :math:\chi_{r}(x_i) = \phi(r_i)\beta(\omega)\ . These are added to the one-electron
integral tensor :math:h_{pq} to compute :math:T_{pq}\ .

"""

chemist_two_body_coeffs, chemist_one_body_coeffs = None, None

if one_body_tensor is not None:
chemist_one_body_coeffs = one_body_tensor.copy()

if two_body_tensor is not None:
chemist_two_body_coeffs = np.swapaxes(two_body_tensor, 1, 3)
one_body_coeffs = -np.einsum("prrs", chemist_two_body_coeffs)

if chemist_one_body_coeffs is None:
chemist_one_body_coeffs = np.zeros_like(one_body_coeffs)

if spatial_basis:
chemist_two_body_coeffs = 0.5 * chemist_two_body_coeffs
one_body_coeffs = 0.5 * one_body_coeffs

chemist_one_body_coeffs += one_body_coeffs

return (x for x in [chemist_one_body_coeffs, chemist_two_body_coeffs] if x is not None)


Using PennyLane

Development

API

Internals