code: 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],
      ...                      [1.398397361, 0.0, 0.0]], requires_grad=False)
      >>> 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[0]
      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
  
imports: import numpy
  import pennylane
label_desc: Use image node_id 1 for calling a global function with example usage: factorize(two_electron, tol_factor, tol_eigval) and returns: factors,  eigvals_m,  eigvecs_m