Department of Materials Science and Engineering, University of Michigan, Ann
Arbor, Michigan 48109, USA

# Tucker tensor approach for accelerating exchange computations in a real-
space finite-element discretization of generalized Kohn-Sham density
functional theory

Vishal Subramanian  Department of Materials Science and Engineering,
University of Michigan, Ann Arbor, Michigan 48109, USA  Sambit Das  Department
of Mechanical Engineering, University of Michigan, Ann Arbor, Michigan 48109,
USA  Vikram Gavini  Department of Mechanical Engineering, University of
Michigan, Ann Arbor, Michigan 48109, USA<EMAIL_ADDRESS>
###### Abstract

The evaluation of Fock exchange is often the computationally most expensive
part of hybrid functional density functional theory calculations in a
systematically improvable, complete basis. In this work, we employ a Tucker
tensor based approach that substantially accelerates the evaluation of the
action of Fock exchange by transforming 3-dimensional convolutional integrals
into a tensor product of 1-dimensional convolution integrals. Our numerical
implementation uses a parallelization strategy that balances the memory and
communication bottlenecks, alongside overalapping compute and communication
operations to enhance computational efficiency and parallel scalability. The
accuracy and computational efficiency is demonstrated on various systems,
including Pt clusters of various sizes and a $\text{TiO}_{\text{2}}$ cluster
with 3,684 electrons.

Figure 1: Graphical abstract

##  1 Introduction

Kohn-Sham Density functional theory 1, 2 (DFT) is widely used to predict
various ground state properties of materials 3. The wide adoption of DFT can
be attributed to its balance between computational complexity and the accuracy
obtained 4. In principle, DFT is exact, however, it involves an unknown
exchange-correlation (XC) functional of the electron density that encapsulates
the many-body quantum-mechanical interactions. While the existence of the
universal XC functional is known, its functional form is an open question to
date 5. The development of accurate and transferable exchange-correlation
functionals continues to be an active research area 6. The progress in the
past many decades has provided model XC functionals that have successfully
predicted phase diagrams, mechanical, electronic, and magnetic properties of
materials 3, 7, 8. Based on the nature of approximation used in the models,
the XC functionals can be classified into different classes such as local
density approximation (LDA) (dependent on local electron density)9,
generalized gradient approximation (GGA) (dependent on local electron density
and its gradient) 10 , hybrid functionals (also has a dependence on single-
electron wavefunctions), etc. In general, studies have shown that the accuracy
of the hybrid XC functionals, such as B3LYP 11, PBE012, HSE06 13, 14, is
better than that of LDA and GGA functionals when compared with experimental
observations 15, 16, 17. In particular, the accuracy of the band gap
calculation is significantly improved by hybrid functionals 18, 19, 20, 21.
This improved accuracy of hybrid functionals is attributed to including a
fraction of the single-electron wavefunction-dependent Fock exchange energy.
This inclusion of the wavefunction-dependent quantities has been formally
justified under the generalized Kohn-Sham framework 22.

The Fock exchange is a non-local operator. Its explicit construction scales as
$\mathcal{O}({N_{e}}^{4})$, where $N_{e}$ is the number of the electrons, as
compared to the cubic-scaling cost of DFT calculations using semi-local
functionals. Instead of the explicit construction, the action of exchange
operator on an orbital can be computed by performing convolution integrals.
The number of convolution integrals required at each iteration scales
quadratically with system size, making it computationally expensive
nevertheless. Many algorithms have been proposed to accelerate the exact
exchange calculations. Broadly, these algorithms can be classified into two
categories: (i) algorithms that reduce the number of times the action of the
exchange operator needs to be computed; (ii) algorithms that reduce the
computational complexity of the action of the exchange operator itself.
Adaptively Compressed Exchange (ACE) 23 is an example of the former, wherein a
low-rank approximation for the exchange operator, which is exact in the space
spanned by the occupied orbitals, is constructed. ACE approach has been shown
to drastically reduce the cost of hybrid DFT calculations 23. However, even
with ACE, the cost of a hybrid calculation can be substantially higher than
using GGA, and the computational cost associated with the action of the
exchange operator still remains a bottleneck. Hence, various reduced-order
scaling strategies such as the interpolative separable density fitting (ISDF)
approach 24, 25, and linear scaling methods like integral screening approaches
and integral evaluation schemes 26, 27, 28 have been proposed. ISDF approach
exploits the rank deficiency in the product of the wavefunctions that appear
in the convolution integrals, while linear-scaling methods utilize the
locality of the wavefunctions to reduce the number of convolution integrals
required. Though these approaches can reduce the computational cost, they do
suffer from some shortcomings. Mainly, the robustness and systematic
convergence of these algorithms are not well established. Further, the
localization strategies employed in linear-scaling methods may not be
generically applicable, especially for metallic systems. Various linear-
scaling approaches also require the explicit construction of the exchange
operator and hence are not amenable to complete basis-sets like plane-waves
and finite-elements.

Over the years, many codes have been developed that solve the eigenvalue
problem arising in the Kohn-Sham formalism, with atomic orbital type basis 29,
30, 31 or the plane-wave basis 32, 33, 34 being popular discretization
choices. The atomic orbital basis is well-suited for systems such as molecules
and atomic clusters. However, systematic numerical convergence with this
basis, especially for metallic systems, is a concern. On the other hand, DFT
codes based on the systematically convergent plane-wave basis exhibit limited
parallel scalability due to the global nature of the basis. These limitations
have restricted typical DFT calculations to systems containing a few thousand
electrons. However, large-scale DFT simulations are useful in a wide range of
scientific studies, such as understanding the energetics of defects,
understanding the electronic structure of large biomolecules, and studying
diffusivity in solids at interfaces, to name a few. This need for large-scale
DFT calculations spurred the development of numerical methods that can handle
large-scale problems efficiently. In particular, the recent use of a higher-
order finite-element (FE) basis in electronic structure calculations 35, 36
has shown significant promise. The FE basis is a piecewise polynomial basis
with desirable features such as locality of the basis which is suitable for
good parallel scalability, spatial adaptivity to efficiently handle non-
periodic systems and all-electron calculations, and its ability to handle
general boundary conditions (periodic, non-periodic and semi-periodic). DFT-FE
37, 38—a massively parallel, open-source, real-space DFT code using adaptive
higher-order FE discretization—has been shown to handle large-scale systems
with $\mathcal{O}(10^{5})$ electrons, exhibiting good parallel scalability on
many-core and hybrid CPU-GPU architectures 39, 40, and has recently been
employed in studying a range of problems, including the electronic structure
of DNA molecules 41, dislocations in crystalline materials 39, 42, phase
transformations in doped nano-films 43, computing spin Hamiltonian parameters
of systems with defects 44, 45, and tackling the inverse DFT problem 46, 47,
48.

In this work, we extend the capability of DFT-FE, presently limited to semi-
local exchange-correlation functionals, to efficiently handle hybrid
functionals. To this end, we build on the ideas presented in 49, 50, 51 to
implement a tensor-based algorithm to accelerate the computation of the action
of the exchange operator. In particular, the product of single-particle
wavefunctions appearing in the exchange operator is approximated using a low-
rank Tucker tensor decomposition with systematic error control. Further, the
$1/r$ kernel appearing in the exchange operator is approximated using a sum of
Gaussians 52. These systematically controllable approximations transform the
3-D convolution integrals into a tensor product of 1-D convolution integrals.
This work presents many key numerical aspects of the implementation that
enable an efficient and scalable treatment of the action of the exchange
operator. We use an auxiliary tensor-structured finite-element mesh employing
higher-order finite-elements for a systematically convergent and efficient
evaluation of the convolution integrals in the application of the exchange
operator, thereby significantly reducing the computational cost. In order to
efficiently handle large-scale systems and obtain good parallel scalability,
we have employed communication patterns that strike a balance between the
memory required per MPI task and the amount of data communicated. In
particular, the communication patterns have been designed to minimize data
transfer and allow asynchronous communication and compute. Our numerical
implementation also developed a fast interface for the transfer of fields
between the unstructured finite-element mesh used to solve the Kohn-Sham
eigenvalue problem and the tensor-structured finite-element mesh used for
evaluating the action of the exchange operator.

We demonstrate the accuracy, efficiency, and parallel scalability of the
implementation by considering benchmark pseudopotential calculations on
platinum (Pt) nanoclusters of increasing sizes, up to 1,584 electrons and a
large $\text{TiO}_{\text{2}}$ cluster containing 3,648 electrons. The accuracy
of computed ground-state energies has been ascertained by comparison with
Quantum Espresso 53, a state-of-the-art open-source plane-wave code. In terms
of computational efficiency, the present implementation is up to $\sim
10\times$ faster than Quantum Espresso for the application of the exchange
operator, which in turn results in considerable speed-ups for the solution
time for the full ground-state problem. The scalability of the implementation
is demonstrated on a 38 atom Pt nanocluster, with the code exhibiting good
parallel scalability even at 128 nodes for a small system.

The remaining paper is structured as follows. Section 2 presents the
mathematical background relevant to the DFT formulation for hybrid functionals
and efficient solution procedures for the governing equations. Section 3
presents the algorithmic details and numerical implementation of the action of
exchange operator using the Tucker tensor based approach. Section 4 presents
the results demonstrating the accuracy and performance of the algorithm and
numerical implementation. Finally, we conclude with an outlook in Section 5.

##  2 Background

In density functional theory 1, 54, the total electronic energy can be written
as the sum of a universal functional of electron density and the energy due to
interaction with an external potential. The ground state energy is determined
by performing a minimization of the functional over the N-representable
densities:

$E_{\textrm{ground-state}}=\min_{\rho(\mathbf{r})}\Big{\\{}\int
V_{ext}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}+F[\rho]\Big{\\}}\,.$ (1)

The universal functional F[$\rho$] is obtained by minimising the expectation
of $\hat{T}+\hat{V}_{ee}$—where $\hat{T}$ is the many-body kinetic energy
operator and $\hat{V}_{ee}$ is the many-body electron-electron interaction
operator—over all anti-symmetric many-body wavefunctions that yield the
density $\rho(\mathbf{r})$:

$F[\rho]=\min_{\Psi\to\rho(\mathbf{r})}\\{\braket{\Psi}{\hat{T}+\hat{V}_{ee}}{\Psi}\\}\,.$
(2)

In the generalised Kohn-Sham (GKS) 22, 55 approach, the many-body interacting
system is replaced by an auxiliary system that can be represented by a single
Slater determinant and has the same ground state density. Since the auxiliary
system is represented by a single Slater determinant, it is described by
single particle equations. By including energy functional $S$ that depends on
the single Slater determinant($\mathbf{\\{\phi\\}}$), the auxiliary system can
partly account for the quantum electron-electron interactions. The
corresponding density functional ($F_{S}[\rho]$) is then given by

$F_{S}[\rho]=\min_{\mathbf{\\{\phi\\}}\to\rho}S[\mathbf{\\{\phi\\}}]\,.$ (3)

The minimisation in Eq. 3 is performed over all single Slater determinants
that yield the density $\rho$. Thus, the functional $F_{S}$ depends on the
choice of $S[\mathbf{\\{\phi\\}}]$. The difference between $F[\rho]$ (defined
in Eq. 2) and $F_{S}[\rho]$ (defined in Eq. 3) is denoted as the remainder
functional $R_{S}[\rho]$. The exact form of $R_{S}[\rho]$ is unknown and is
approximated in practice.

The $E_{\textrm{ground-state}}$ (Eq. 1) can thus be expressed as

$E_{\textrm{ground-state}}=\min_{\rho(\mathbf{r})\to N}\Big{\\{}\int
V_{ext}(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}+F_{S}[\rho]+R_{S}[\rho]\Big{\\}}\,.$
(4)

Using the implicit dependence on the Slater determinant, the minimization in
the above equation, reduces to a minimization over admissible Slater
determinants leading to a non-linear eigenvalue problem commonly referred to
as the generalized Kohn-Sham equation:

$\Big{(}\hat{O}[\mathbf{\\{\phi\\}}]+V_{ext}(\mathbf{r})+V_{R}(\mathbf{r})\Big{)}\psi_{i}(\mathbf{r})=\epsilon_{i}\psi_{i}(\mathbf{r})\,,$
(5)

where $\hat{O}[\mathbf{\\{\phi\\}}]$ is a non-multiplicative potential
obtained from the functional derivatives of $S[\mathbf{\\{\phi\\}}]$ with
respect to the orbitals, and $V_{R}$ is a multiplicative potential obtained
from the functional derivative of $R_{S}[\rho]$ with respect to $\rho$. In the
case of hybrid functionals, $S[\mathbf{\\{\phi\\}}]$ includes the non-
interacting kinetic energy of electrons and a fraction of the exact exchange
given by

$S[\mathbf{\\{\phi\\}}]=\sum_{i}\int\psi_{i}(\mathbf{r})\Big{(}-\frac{1}{2}{\nabla_{i}}^{2}\Big{)}\psi_{i}(\mathbf{r})d\mathbf{r}-\frac{\alpha}{2}\sum_{i}\sum_{j}\int\int\frac{\psi_{i}(\mathbf{r})\psi_{j}(\mathbf{r})\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}d\mathbf{r^{\prime}}d\mathbf{r}\,,$
(6)

where $\alpha$ is generally chosen to be between $0$ and $1$. The above
expressions are for a non-periodic system, which is the focus of this work,
where the single-electron wavefunctions are real-valued. From the above,
$\hat{O}[\mathbf{\\{\phi\\}}]$ is obtained to be

$\hat{O}[\mathbf{\\{\phi\\}}]=-\frac{1}{2}{\nabla_{i}}^{2}+\alpha V_{\rm
X}\,,$ (7)

where $V_{\rm X}$ denotes the exchange operator obtained from the functional
derivative of the exact exchange. The exact form of remainder functional is
unknown and an approximation often considered is given by $R_{S}[\rho]=E_{\rm
H}[\rho]+(1-\alpha)E_{\rm X,sl}[\rho]+E_{\rm c,sl}[\rho]$, where $E_{\rm H}$
denotes the Hartree energy, $E_{\rm X,sl}[\rho]$ and $E_{\rm c,sl}[\rho]$
denote the semi-local exchange and correlation parts, respectively, of model
$\rho$-dependent exchange-correlation (XC) functionals. In hybrid XC
functionals such as PBE012, semi-local functionals (PBE 10) are combined with
a fraction of the exact exchange. The functional derivative of the above form
of $R_{S}[\rho]$ is given by

$V_{R}[\rho](\mathbf{r})=V_{\rm H}[\rho](\mathbf{r})+(1-\alpha)V_{\rm
X,sl}[\rho](\mathbf{r})+V_{\rm c,sl}[\rho](\mathbf{r})\,.$ (8)

Finally, substituting the above terms in Eq. 5, the following form of the GKS
equation employed in hybrid DFT calculations is obtained:

$\Big{(}-\frac{1}{2}{\nabla_{i}}^{2}+V_{ext}(\mathbf{r})+V_{H}[\rho](\mathbf{r})+\alpha
V_{X}+(1-\alpha)V_{X,SL}[\rho](\mathbf{r})+V_{c,SL}[\rho](\mathbf{r})\Big{)}\psi_{i}(\mathbf{r})=\epsilon_{i}\psi_{i}(\mathbf{r})\,.$
(9)

###  2.1 Exchange operator

The exchange operator is a non-local operator which requires the single-
particle wavefunctions explicitly:

$(V_{\rm
X})(\mathbf{r},\mathbf{r^{\prime}})=-\sum_{j=1}^{N_{e}}\frac{\psi_{j}(\mathbf{r})\psi_{j}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}\,.$
(10)

The explicit evaluation of the exchange operator can be computationally
expensive as the complexity scales as $\mathcal{O}(N_{e}^{4})$, where $N_{e}$
is the number of electrons in the system. Algorithms such as the resolution of
identity 56, 57 have been proposed to accelerate this computation. However,
these algorithms are not compatible with a complete basis as they are still
prohibitively expensive. The explicit construction of the exchange operator
can be avoided by employing iterative algorithms to compute the action of the
exchange operator on an occupied orbital $\psi_{i}$, and can be expressed as

$(V_{\rm
X})\psi_{i}(\mathbf{r})=-\sum_{j=1}^{N_{e}}\psi_{j}(\mathbf{r})\int\frac{\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}d\mathbf{r^{\prime}}\,.$
(11)

Conventionally, numerical implementations based on systematically improvable
and complete basis sets evaluate the convolution integrals arising in the
action of the exchange operator by either solving a Poisson equation (method
of choice for real-space basis sets) or using fast Fourier transforms (method
of choice for plane-wave basis). The exchange operator requires occupied
orbitals, which are obtained as the solution of the non-linear eigenvalue
problem (cf. (9)). Thus, these equations are often solved self-consistently,
using a nested two-level self-consistent field (SCF) procedure as discussed
below.

###  2.2 Two-level nested SCF

Two-level nested SCF refers to a nested fixed point iteration scheme involving
inner and outer fixed point iterations. In the outer fixed point iteration,
the wavefunction dependent parts of the Hamiltonian (exchange operator) are
updated. In the inner SCF, the wavefunction dependent parts of Hamiltonian are
held fixed and the electron-density dependent non-linear eigenvalue problem is
solved until convergence. With such a split, robust mixing strategies such as
Anderson mixing 58 can be used to accelerate the convergence of the electron-
density in the inner SCF iterations. Within DFT-FE, we use a Chebyshev-
filtered subspace iteration technique 59, 35 to solve for the lowest occupied
eigensubspace and the ground-state electron density. The Kohn-Sham (KS)
orbitals obtained at the end of the inner SCF iterations are used to update
the exchange operator. The simulation is considered converged when the
difference in exact-exchange energy between successive outer SCF iterations is
below a pre-determined tolerance. Such a two-level nested SCF scheme is
implemented in conjunction with Adaptively Compressed Exchange (ACE) 23, 60
operator to reduce the computational cost further, which is discussed below.

###  2.3 ACE implementation

The ACE implementation aims to construct a low-rank representation of the
exchange operator 23. The exchange operator is generally full-rank, and a low-
rank approximation can lead to aproximation errors. However, a low-rank
approximation of the exchange operator that is exact for the occupied KS
orbitals can be constructed and is the main idea behind the ACE
implementation. The ACE operator, denoted by $V^{ACE}_{\rm X}$, is constructed
such that it satisfies the Hermiticity and negative semi-definite nature of
$V_{\rm X}$

$V^{ACE}_{\rm
X}(\mathbf{r},\mathbf{r^{\prime}})=-\sum_{k=1}^{N_{e}}\xi_{k}(\mathbf{r})\xi_{k}(\mathbf{r^{\prime}})\,,$
(12)

where $\mathbf{\xi}_{k}$ is the projection vector defined as

$\xi_{k}(\mathbf{r})=\sum_{i=1}^{N_{e}}W_{i}(\mathbf{r})\big{(}\mathbf{L}^{-1}\big{)}_{ki}\,.$
(13)

In the above, $W_{i}$ is the vector obtained from the action of the exchange
operator on the $i^{th}$ KS orbital and is given by

$W_{i}(\mathbf{r})=V_{X}[\psi(\mathbf{r})]\psi_{i}(\mathbf{r})\,,$ (14)

and the matrix $\mathbf{L}$ is obtained by performing a Cholesky decomposition
on matrix $\mathbf{M}$:

$\mathbf{M}=-\mathbf{L}\mathbf{L}^{T}\,,\quad\quad
M_{ij}=\int\psi_{i}(\mathbf{r})W_{j}(\mathbf{r})d\mathbf{r}\,.$ (15)

The flowchart of the entire algorithm for solving the GKS equations in hybrid
DFT is depicted in Fig. 2. The construction of the ACE operator, as described
above, is performed in updateExchange. The advantage of using ACE operator is
that, once constructed, it significantly reduces the cost of action of the
exchange operator on trial KS orbitals, thereby reducing the walltimes of the
inner SCF iterations. The aim of this work is to reduce the cost of the
calculations performed in updateExchange using Tucker-tensor decomposition.

Figure 2: Nested SCF loop algorithm for hybrid DFT

###  2.4 Tucker-tensor decomposition

Tensor-structured techniques, in particular the Tucker-tensor decomposition,
have been used to construct low-rank approximations of data and fields in
three and higher dimensions—we refer to 61 for a comprehensive review of the
subject. The Tucker tensor decomposition of a field in 3-D is given by

$\phi(\mathbf{r})\approx\sum_{\alpha=1}^{r1}\sum_{\beta=1}^{r2}\sum_{\gamma=1}^{r3}\mathcal{G}_{\alpha\beta\gamma}{f^{\alpha}}({x}){g^{\beta}}({y}){h^{\gamma}}({z})\,.$
(16)

In the above, $r=\max\\{r_{1},r_{2},r_{2}\\}$ denotes the rank of the Tucker
decomposition, $\mathcal{G}_{\alpha\beta\gamma}$ denotes the core-tensor and
$f^{\alpha}$, $g^{\beta}$, $h^{\gamma}$ are 1-D functions. Notably, the
approximation error for regular fields in the Tucker decomposition decays
exponentially with increasing Tucker rank—also demonstrated in Sec. 3.4. For
electronic structure, including single-particle wavefunctions and electron
density, a Tucker rank of $\mathcal{O}(10)$ provides a good approximation to
use in practical calculations 62, 63, 64.

In the context of the present work, which seeks to accelerate the evaluation
of the action of the Fock exchange operator on single-particle KS
wavefunctions, the Tucker decomposition of the KS wavefunctions presents a
path to reduce the convolution integrals in 3D to convolutions in 1D as
discussed in Sec. 3.4. A proof-of-concept demonstration of this idea was
presented in 62, and this idea was employed in evaluating the convolution
integrals arising in the computation of Hartree potential in a Tucker tensor
basis 65, 64. In the present work, we build on these earlier efforts and
develop an efficient and scalable implementation for hybrid DFT calculations
in a systematically convergent and complete finite-element basis—a piecewise
polynomial basis in real-space. In the following section we discuss the
various aspects of the computational framework, the details of the numerical
implementation, including the strategies for achieving an efficient and
scalable implementation.

##  3 Algorithmic implementation

In this section, we present the various aspects of our Tucker-tensor algorithm
used to accelerate the action of the Exchange operator on the Kohn-Sham
orbitals.

###  3.1 Overview of the algorithm

We provide a brief overview of the various steps in the Tucker-tensor
algorithm before discussing the implementation details. To begin, we define a
finite-dimensional subspace ${V}^{M}_{1}$ spanned by an unstructured finite-
element basis onto which we project the KS Hamiltonian for solving the inner
fixed-point iteration on electron-density. Next, the KS wavefunctions
($\psi_{i}$) from the solution of the inner SCF iterations, belonging to
${V}^{M}_{1}$, are transferred to another finite-element basis which has a
tensor structure, as this is required to perform Tucker-tensor decomposition.
We refer to this subspace as ${V}^{M}_{2}$. Subsequently, we compute the
Tucker decomposition and the convolution integrals in ${V}^{M}_{2}$, using the
tensor structured basis. Finally, we compute action of the exchange operator
$W_{i}=V_{X}\psi_{i}$ in ${V}^{M}_{2}$, and this is transferred back to
${V}^{M}_{1}$ as $\tilde{W_{i}}\in{V}^{M}_{1}$. These steps are mathematically
represented as:

$\psi^{{V}^{M}_{2}}_{i}=T_{{V}^{M}_{1}\to{V}^{M}_{2}}\ \
\psi^{{V}^{M}_{1}}_{i}\quad\forall\,i$ (17a)
$\varphi_{ji}^{{V}^{M}_{2}}(\mathbf{r})=\int\frac{{\psi_{j}}^{{V}^{M}_{2}}(\mathbf{r^{\prime}}){\psi_{i}}^{{V}^{M}_{2}}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}d\mathbf{r^{\prime}}\quad\forall\,(i,j)$
(17b)
${W_{i}}^{{V}^{M}_{2}}=\sum_{j}{\psi_{j}}^{{V}^{M}_{2}}\varphi_{ji}^{{V}^{M}_{2}}\quad\forall\,i$
(17c)
${\tilde{W_{i}}}^{{V}^{M}_{1}}=T_{{V}^{M}_{2}\to{V}^{M}_{1}}{W_{i}}^{{V}^{M}_{2}}\quad\forall\,i$
(17d)

where $\psi_{i}$ refers to the wavefunctions, $T$ represents the transfer
operation from one subspace to another, and $\varphi_{ji}$ represents the
output of the convolution. The implementation details are discussed in the
subsequent subsections.

As the subspaces ${V}^{M}_{1}$ and ${V}^{M}_{2}$ can be different in general,
this can lead to a loss in accuracy during these transfers. To mitigate this
issue, a structured finite-element (FE) basis with a commensurate mesh
refinement as that of the unstructured FE basis is used. Using an FE basis for
the Tucker-tensor decomposition provides a simplified, efficient, and accurate
numerical framework to perform the interpolations required for the transfer
and the evaluation of convolution integrals.

###  3.2 Creation of Structured FE basis

Tensor structured operations require a structured mesh, whereas the FE basis
for the Kohn-Sham problem can work in general with an unstructured mesh
(henceforth denoting the corresponding mesh and basis as $\mathcal{U}-{\rm
FE}$ mesh and basis). Hence, we create a structured mesh—denoting the
corresponding mesh and basis as $\mathcal{S}-{\rm FE}$ mesh and basis—to aid
the tensor-structured computations. To this end, we construct three 1-D meshes
(one along each direction), and the tensor product of these 1-D meshes
provides the desired tensor-structured 3-D mesh. Figure 3 shows the schematic
of the unstructured and structured meshes, the tensor structured nature of the
$\mathcal{S}-{\rm FE}$ mesh, and the various aspects of FE discretization
including FE cells, FE nodes and FE basis. We refer the readers to 66 for a
more comprehensive discussion on the approximation properties of finite
element discretization, and to 36, 35, 37, 67, 68 for the use of finite
elements in electronic structure calculations.

Figure 3: Schematic describing various aspects of the finite element meshes
used in this work. (a)-(b) An unstructured finite element mesh
($\mathcal{U}-{\rm FE}$ mesh) and a structured finite element mesh
($\mathcal{S}-{\rm FE}$ mesh). In (a)-(b), the green triangle represents an
atom. Selective refinement around the atom results in a non-uniform,
adaptively refined mesh. The mesh parameters—Base Mesh Size, Atom Ball Radius,
and Mesh Size Near Atom—are used to construct the spatially adaptive meshes.
The purple lines mark the boundary of a finite element cell; (c) The
structured mesh represented as a tensor product of 1-D meshes; (d) Finite
element cell, with ‘x’ markers denoting the finite element nodes. The finite
element basis is the tensor product of the 1D Lagrange polynomials constructed
from the finite element nodes.

The $\mathcal{S}-{\rm FE}$ mesh is generated by using similar mesh parameters
as that of the $\mathcal{U}-{\rm FE}$ mesh, i.e. similar ‘base mesh-size’
($bms$; the initial uniform mesh size before any adaptive refinement) and
similar ‘mesh-size near atom’ ($msna$; the smaller mesh-size near the nuclei,
obtained by adaptive refinement). This is achieved by starting with a 1-D mesh
set to the base mesh-size ($bms$), followed by selective refinement until the
finite-element cells near an atom are smaller than the mesh-size near the atom
($msna$). The selective refinement uses an adaptive algorithm that refines
finite-element cells based on their proximity to an atom. The cells are then
populated with FE nodes.

In order to reduce the peak memory requirements, the coordinates of the FE
nodes in the $\mathcal{S}-{\rm FE}$ mesh are not stored. Instead, the
coordinates of the FE nodes in the 1-D meshes are stored, from which the
coordinates of the 3-D mesh are computed on the fly. We employ a lexicographic
numbering scheme that simplifies the map from nodal coordinates in the 3-D
mesh to the nodal coordinates in the three 1-D meshes, and vice versa.

###  3.3 Transfer of wavefunctions between Unstructured and Structured FE
basis

The first step in the Tucker-tensor algorithm to accelerate the action of the
exchange operator is the transfer of the wavefunctions in the
$\mathcal{U}-{\rm FE}$ basis to the $\mathcal{S}-{\rm FE}$ basis (cf. Eq.
17a), which is performed by interpolation. In the FE basis, interpolation to
an arbitrary point requires the evaluation of the shape functions—basis
functions local to a given FE cell—at the point of interest. To this end,
corresponding to each $\mathcal{U}-{\rm FE}$ cell, we find all the FE nodes of
the $\mathcal{S}-{\rm FE}$ mesh that lie in the cell. Subsequently, the
relevant $\mathcal{U}-{\rm FE}$ shape function values at the $\mathcal{S}-{\rm
FE}$ node locations are evaluated. This map and the shape function values,
used for the interpolation, are computed once and stored _a priori_.

Using the stored shaped function values, we now discuss a cache-efficient
implementation to perform the transfer operation. In the $\mathcal{U}-{\rm
FE}$ basis, we store the values of all KS wavefunctions for all nodes in a
contiguous memory layout with the wavefunction index for a given FE node being
the fastest index in this storage. Since the shape functions are independent
of the wavefunctions, the interpolation from one cell of the $\mathcal{U}-{\rm
FE}$ to the relevant nodes of the $\mathcal{S}-{\rm FE}$ is performed
simultaneously for all the wavefunctions. This can be recast as the following
matrix-matrix multiplication for each $\mathcal{U}-{\rm FE}$ cell

${\psi^{\mathcal{S}-FE}_{i,a}}=\sum_{b}N_{ab}{\psi^{\mathcal{U}-FE}_{i,b}}\,,$
(18)

where $a$ and $b$ refers to the nodal index in the $\mathcal{S}-{\rm FE}$ and
$\mathcal{U}-{\rm FE}$ basis, respectively. $N_{ab}$ refers to the shape
function value of the $b^{th}$ node of the $\mathcal{U}-{\rm FE}$ basis at the
spatial location of the $a^{th}$ nodal point of $\mathcal{S}-{\rm FE}$ basis.

The transfer of fields computed in the $\mathcal{S}-{\rm FE}$ basis (cf. Eqs.
(17b)) and (17c)) to the $\mathcal{U}-{\rm FE}$ basis proceeds along similar
lines.

###  3.4 Convolution operation

In order to compute

$\varphi_{ji}(\mathbf{r})\equiv\int\frac{\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})}{|\mathbf{r}-\mathbf{r^{\prime}}|}d\mathbf{r^{\prime}}\,,$
(19)

we compute $\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})$ by
trivially multiplying the fields at the FE nodes of $\mathcal{S}-{\rm FE}$
mesh. Subsequently, we perform a Tucker decomposition for every pair of
$i$-$j$ wavefunctions using the STHOSVD algorithm implemented in TuckerMPI 69,
providing the low-rank approximation

$\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})\approx\sum_{\alpha}\sum_{\beta}\sum_{\gamma}\mathcal{G}_{\alpha\beta\gamma}f_{ij}^{\alpha}(x^{\prime})g_{ij}^{\beta}(y^{\prime})h_{ij}^{\gamma}(z^{\prime})\,.$
(20)

The error in the Tucker decomposition decreases exponentially with the Tucker
rank as depicted in Fig. 4 for a 14-atom Pt nanocluster, a benchmark system
used in this study. Notably, very low approximation errors can be achieved
with modest values of rank around 50.

Figure 4: Approximation error in the Tucker decomposition of
$\psi_{j}(\mathbf{r^{\prime}})\psi_{i}(\mathbf{r^{\prime}})$, measured in
terms of the relative error in Frobenius norm on the FE nodes as a function of
Tucker rank. Benchmark system: Pt-14 cluster.

Further, to enable a 1-D separable form for fast evaluation of convolution
integrals, we approximate the $1/r$ kernel using a canonical tensor
decomposition with 1-D Gaussian functions 52 as

$\frac{1}{|\mathbf{r}-\mathbf{r^{\prime}}|}\approx\sum_{t=1}^{T}\mu_{t}exp(-\theta_{t}(x-x^{\prime})^{2})exp(-\theta_{t}(y-y^{\prime})^{2})exp(-\theta_{t}(z-z^{\prime})^{2})\,.$
(21)

Figure 5: Accuracy of the Gaussian approximation for the $1/r$ kernel with
T=35.

The accuracy of this approximation for the $1/r$ kernel is depicted in Fig. 5.
Though the approximation fails to capture the asymptotical behavior of the
kernel accurately at small $r$, the accuracy of the convolution integral is
negligibly affected as long as the error in the approximation is controlled
within a region $[{r_{1}},{r_{2}}]$. The smallest spacing between the 1-D FE
nodes and the domain size determines the bounds ${r_{1}}$ and ${r_{2}}$,
respectively. Finally, using Eqs. 20 and 21 in Eq. 19, the 3-D convolution
integrals can be well approximated as a tensor product of 1-D convolution
integrals as

$\begin{split}\varphi_{ji}(\mathbf{r})\approx&\sum_{\alpha\beta\gamma}\int\int\int\sum_{t}\Big{\\{}\mathcal{G}_{\alpha\beta\gamma}f_{ij}^{\alpha}(x^{\prime})g_{ij}^{\beta}(y^{\prime})h_{ij}^{\gamma}(z^{\prime})\\\
&\mu_{t}exp(-\theta_{t}(x-x^{\prime})^{2})exp(-\theta_{t}(y-y^{\prime})^{2})exp(-\theta_{t}(z-z^{\prime})^{2})\Big{\\}}dx^{\prime}dy^{\prime}dz^{\prime}\,,\\\
=&\sum_{t}\sum_{\alpha\beta\gamma}\mathcal{G}_{\alpha\beta\gamma}\mu_{t}\Big{[}\int
f_{ij}^{\alpha}(x^{\prime})exp(-\theta_{t}(x-x^{\prime})^{2})dx^{\prime}\Big{]}\times\\\
&\Big{[}\int
g_{ij}^{\beta}(y^{\prime})exp(-\theta_{t}(y-y^{\prime})^{2})dy^{\prime}\Big{]}\times\Big{[}\int
h_{ij}^{\gamma}(z^{\prime})exp(-\theta_{t}(z-z^{\prime})^{2})dz^{\prime}\Big{]}\,.\end{split}$
(22)

The evaluation of Eq. 22 requires a numerical evaluation of the 1-D integrals,
which are computed using the 1-D FE meshes that generate the $\mathcal{S}-{\rm
FE}$ mesh (cf. Figure 3(c)). Numerical integration using the 1-D FE basis
involves interpolating to a set of quadrature points and summing with
appropriate quadrature weights. The number of quadrature points required to
integrate the 1-D functions in Eq.(22) accurately can be high due to the
sharply varying exponentials in the integrand. However, the number and
location of quadrature points and the powers of the exponentials are all known
_a priori_. Since these values are independent of the wavefunction values,
they can be precomputed. Once these values are precomputed, the numerical
integration of the 1-D integrals can be recast as a summation over FE nodal
values of wavefunction specific quantities and the precomputed values. As the
number of nodal points are much lower than the number of quadrature points, a
significant reduction in computation cost can be obtained. Mathematically, the
recasting of the compute operations can be expressed as follows: the
interpolation of the field $f_{ij}^{\alpha}(x^{\prime})$ to the quadrature
points is a weighted sum of the nodal values given by

$f_{ij,b}^{\alpha}=\sum_{a}\xi_{ab}f_{ij,a}^{\alpha}\,,$ (23)

where $\xi_{ab}$ is the 1-D shape function of the $a^{th}$ node at the
$b^{th}$ quadrature point. Further, the 1-D convolution integrals numerically
evaluated at any location $x$ are given by

$\int
f_{ij}^{\alpha}(x^{\prime})exp(-\theta_{t}(x-x^{\prime})^{2})dx^{\prime}=\sum_{b}f_{ij,b}^{\alpha}exp(-\theta_{t}(x-x_{b})^{2})Jw_{b}\,,$
(24)

where the summation is over the quadrature points $b$, $x$ represents the
spatial location at which the convolution integrals are evaluated, and
$Jw_{b}$ represents the determinant of the 1-D finite-element Jacobian times
the quadrature weight. Using the expression for $f_{ij,b}^{\alpha}$ from Eq.
23 in Eq. 24, we get

$\int
f_{ij}^{\alpha}(x^{\prime})exp(-\theta_{t}(x-x^{\prime})^{2})dx^{\prime}=\sum_{b}\Big{(}\sum_{a}\xi_{ab}f_{ij,a}^{\alpha}\Big{)}exp(-\theta_{t}(x-x_{b})^{2})Jw_{b}\,.$
(25)

Rearranging the summation, we obtain

$\int
f_{ij}^{\alpha}(x^{\prime})exp(-\theta_{t}(x-x^{\prime})^{2})dx^{\prime}=\sum_{a}\Big{(}\sum_{b}exp(-\theta_{t}(x-x_{b})^{2})Jw_{b}\xi_{ab}\Big{)}f_{ij,a}^{\alpha}\,,$
(26)

where the summation over $b$ (quadrature points) can be precomputed once for
the 1-D FE basis and stored. The summation over $a$ involves wavefunction
specific values ($f_{ij,a}^{\alpha}$), which has to be performed repeatedly,
but the computational cost of this summation is cheap as the number of FE
nodes in 1-D mesh are $\mathcal{O}(10)-\mathcal{O}(100)$.

The convolution integrals are evaluated at the nodal locations of the 1-D
mesh, i.e $x$ is set to the spatial location of FE nodes in 1-D FE mesh.
Hence, the above summation has to be performed for all pairs of $i$-$j$, at
all nodes FE nodes in the 1-D mesh, for each $t$ along each direction. These
computations are performed as matrix-matrix products by choosing appropriate
memory storage for all the above variables. Once the convolution of every
$i$-$j$ pair is computed, we reconstruct the full 3-D tensor and multiply with
$\psi_{j}$ (cf. Eq. 17c). Finally, the sum along the $j$ index provides the
action of the exchange operator on $\psi_{i}$, computed in $\mathcal{S}-FE$
basis. In our implementation, this step of convolution operations—the
evaluation of the 1-D convolution integrals—together with reconstruction,
forms the most computationally expensive part of the calculation.

###  3.5 MPI implementation

In order to enable large-scale computations involving thousands of electrons,
a scalable numerical implementation of the Tucker-tensor algorithm is crucial.
To this end, we now discuss our task/memory parallelization and MPI
communication strategy. Eq. 19 requires the computation of $N_{e}^{2}$
convolutions, corresponding to each $i$-$j$ pair. A high level of task
parallelization can be achieved by distributing these $N_{e}^{2}$ convolutions
equitably among the tasks. For the purpose of computing the convolution
integrals, it is efficient to have the KS wavefunction fields whose
convolution are being computed to be local to the task. Thus, an efficient
layout for storing KS wavefunctions in the $\mathcal{S}-{\rm FE}$ basis is
band parallelism over the processor grid as shown in Fig. 6. To elaborate,
each task computing a subset of the convolution integrals has the relevant KS
wavefunctions in the $\mathcal{S}-{\rm FE}$ basis local to the processor. The
column index range for each MPI task determines the subset of $j$-indexed
wavefunctions assigned to the task. Similarly, the task’s row index range for
each MPI task determines the subset of $i$-indexed wavefunctions assigned to
it. Upon computing the convolution integrals, summing the output obtained from
the tasks along a row of the processor grid provides the final output of the
action of the exchange operator on a wavefunction $\psi_{i}$
($V_{X}\psi_{i}$). We note that such a memory parallelization layout for KS
wavefunctions in $\mathcal{S}-{\rm FE}$ basis is efficient from a
communication standpoint, as no transfer of wavefunctions between tasks is
necessary while computing the convolution integrals. However, this strategy
increases the memory footprint as duplicate copies of KS wavefunctions are
stored—for instance, in the layout in Fig. 6, $\psi_{0}$ is stored in Proc 0,
Proc 1, Proc 2 and Proc 4. While this does not pose a limitation on large
number of tasks, the increased memory footprint can be an issue when running
the simulation on a limited number of nodes. To address this, we implemented a
strategy that reduces the overall memory footprint. This is based on the
observation that all tasks along a column of the processor grid (cf. Fig. 6)
require the same subgroup of $j$-indexed KS wavefunctions. Leveraging this, we
implement a round-robin algorithm to transfer the j-indexed KS wavefunction
data among the tasks in a column, as shown in Fig. 7. To elaborate, $i$-$j$
pairs are assigned to a task based on its position in the grid. Once the task
completes the necessary calculations corresponding to the $i$-$j$ pairs for a
given $j$, the wavefunction $j$ temporarily assigned to this task is sent to
the next task in the column using the round-robin algorithm, reusing the pre-
allocated memory.

Figure 6: MPI task parallelization for KS wavefunctions in the structured FE
mesh, shown for six wavefunctions distributed over six tasks. Summing the
output from the tasks along a row provides the action of the exchange
operator. Figure 7: Representation of memory parallelization for different
subgroups of MPI tasks using round-robin algorithm. At the end of each step,
the wavefunction is transferred to the next task using MPI point-to-point
communication calls.

While band parallelization is adopted for storing the KS wavefunctions in
$\mathcal{S}-{\rm FE}$ basis, where the complete data for any wavefunction is
local to a task, the parallelization employed for the $\mathcal{U}-{\rm FE}$
basis is based on domain decomposition, which is more ideal for the solution
of the eigenvalue problem in the inner SCF iteration on electron density. To
elaborate, in domain decomposition, the FE mesh is distributed across various
MPI tasks, with each task holding the information for all KS wavefunctions in
the finite-element cells local to the task. We note that DFT-FE framework,
which handles all aspects of the computation in the inner SCF, is based on
domain decomposition parallelization with good parallel scalability exhibited
on both many-core CPU as well as hybrid CPU-GPU architectures 37, 38, 39, 40.
Thus, the transfer of fields between $\mathcal{S}-{\rm FE}$ basis and
$\mathcal{U}-{\rm FE}$ basis has to navigate the different parallelization
strategies for $\mathcal{S}-{\rm FE}$ and $\mathcal{U}-{\rm FE}$ meshes.
Importantly, the same set of processors are used to store the fields on both
$\mathcal{S}-{\rm FE}$ and $\mathcal{U}-{\rm FE}$ basis.

The main operation in the transfer of fields between $\mathcal{S}-{\rm FE}$
basis and $\mathcal{U}-{\rm FE}$ basis is an interpolation operation (cf. Eq.
18). In order to transfer the KS wavefunctions from $\mathcal{U}-{\rm FE}$
basis to $\mathcal{S}-{\rm FE}$ basis, we first create MPI communication
pattern that maps domain decomposition parallelization in $\mathcal{U}-{\rm
FE}$ mesh to the band parallelization in $\mathcal{S}-{\rm FE}$ mesh. This
entails mapping FE nodes of $\mathcal{S}-{\rm FE}$ to FE cells containing
these spatial locations in the $\mathcal{U}-{\rm FE}$. This creates a
communication pattern, where each task handling a set of KS wavefunctions in
the $\mathcal{S}-{\rm FE}$ basis will gather the interpolated fields from the
$\mathcal{U}-{\rm FE}$ mesh distributed over all tasks and collates the data.
We note that the communication pattern only needs to be established at the
beginning of calculations, and the interpolation and transfer of fields uses
this communication pattern in every call of updateExchange. To improve the
efficiency of the transfer operation, we also implemented overlapping compute
and communication. While the communication to task $P$ of already computed
interpolated values is happening, the computation of interpolated values to be
sent to task $P+1$ can commence. We achieve this using asynchronous point-to-
point communication calls ($\texttt{MPI}\\_\texttt{Isend()}$,
$\texttt{MPI}\\_\texttt{Irecv()}$). A similar implementation is used for
transferring fields from $\mathcal{S}-{\rm FE}$ basis to $\mathcal{U}-{\rm
FE}$ basis.

##  4 Results and Discussion

We now present the accuracy, computational efficiency and parallel scalability
of our implementation of the action of the exchange operator using Tucker
tensor decomposition in the DFT-FE framework. We demonstrate the accuracy in
ground-state energy on a range of small- to medium-sized benchmark materials
systems involving both insulating and metallic systems. Further, we
demonstrate the computational efficiency on Pt clusters of increasing sizes up
till 1,584 electrons, and a large $\text{TiO}_{\text{2}}$ cluster with 3,648
electrons. We use ONCV SG 15 pseudo potentials 70 and the PBE0 12 exchange-
correlation functional for all these calculations. We employ Anderson mixing
58 for the inner fixed-point iteration in electron density. We do not use
mixing in the outer fixed-point iteration involving the Kohn-Sham orbitals in
the construction of the exchange operator. In all the calculations, we use a
Fermi-Dirac smearing corresponding to 500K.

###  4.1 Accuracy benchmarks

The ground state energy obtained from our implementation in DFT-FE is compared
against Quantum Espresso (QE) 32 to validate the accuracy of the approach and
implementation. We consider the hybrid functional DFT calculation in DFT-FE to
be converged when the difference between the Fock exchange energy of
successive outer SCF iterations is less than a pre-determined tolerance. We
note that this is slightly different from the procedure followed by QE, in
which the simulation runs until the difference in the total energy between
successive outer SCFs is below a pre-determined tolerance. Hence, during the
performance comparison, the simulation in both the codes is considered
converged when the difference between the total energy between successive
outer SCFs is below $10^{-5}$ Ha/atom. The discretization parameters in both
QE and DFT-FE are set such that the discretization error is below $10^{-4}$
Ha/atom for the ground state energy obtained for PBE exchange-correlation
functional. The parameters used to generate the $\mathcal{S}-{\rm FE}$ mesh
are also set to be commensurate with discretization errors below $10^{-4}$
Ha/atom. The benchmark systems we use to demonstrate the accuracy include: i)
Au-Pt dimer; ii) Benzamide molecule; iii) 19-atom Pt nanocluster (Pt-19); and
iv) 38-atom Pt nanocluster (Pt-38). All QE studies are performed using a
vacuum just sufficient for the error in ground-state energy to be below
$10^{-4}$ Ha/atom. In DFT-FE large domains for non-periodic calculations are
easily accessible by leveraging the spatial adaptivity in the finite-element
basis. The computed ground-state energies for the various benchmark systems
are reported in Table 1. Notably, the ground-state energies from DFT-FE and QE
agree to within the discretization errors of $10^{-4}$ Ha/atom. Also, the
number of outer iterations for convergence, which involves updating the
exchange operator, are similar for DFT-FE and QE. The specific details of the
simulation are elaborated below.

i) Au-Pt dimer: Au-Pt dimer with bond length of 4.753 Bohr in considered,
which has 37 electrons. A plane-wave cutoff of 70 Ry is used for the KS
wavefunctions in the QE simulation with a cell size of 40 Bohr. In the DFT-FE
calculation, the $\mathcal{U}-{\rm FE}$ mesh is constructed using a mesh-size
near atom ($msna$) of 0.4, a base mesh-size ($bms$) of 6.66, and a finite-
element polynomial order of 4. The $\mathcal{S}-{\rm FE}$ mesh is constructed
using a mesh-size near atom ($msna$) of 0.625, a base mesh-size ($bms$) of
10.0, and a finite-element polynomial order of 4.

ii) Benzamide molecule: Benzamide molecule
($\text{C}_{7}\text{H}_{7}\text{NO}$) comprises of 46 electrons. A cutoff of
100 Ry is used for KS wavfunctions in QE, with a cell size of 90 Bohr. In the
DFT-FE calculation, both $\mathcal{U}-FE$ and $\mathcal{S}-FE$ meshes are
constructed using $bms$ of 10.0 and $msna$ of 1.2. The polynomial order for
$\mathcal{U}-{\rm FE}$ is chosen to be 6, while the polynomial order for
$\mathcal{S}-{\rm FE}$ is 5.

iii) Pt -19 atom cluster: This system comprises of 342 electrons. The QE
simulation is performed using a cutoff of 70 Ry for KS wavefunctions, with a
cell size of 40 Bohr. The DFT-FE calculation is performed using a
$\mathcal{U}-{\rm FE}$ with $bms$ 8.4, $msna$ 1.0, and polynomial order 6. The
$\mathcal{S}-{\rm FE}$ is constructed with $bms$ 7.8, $msna$ 1.5, and
polynomial order 6.

iv) Pt - 38 atom cluster: This system comprises of 684 electrons. The QE
simulation is performed using a cutoff of 70 Ry for KS wavefunctions and a
cell size of 50 Bohr. The DFT-FE calculation is performed using a
$\mathcal{U}-{\rm FE}$ with $bms$ 7.5, $msna$ 1.0 and polynomial order 6. The
$\mathcal{S}-{\rm FE}$ is constructed with $bms$ 11.0, $msna$ 1.5 and
polynomial order 6.

Table 1: Comparison of ground-state energies with PBE0 exchange-correlation functional for the various benchmark systems computed using DFT-FE and QE. $E_{g}$ denotes the ground-state energy in Ha. The difference in the computed ground-sate energies per atom using DFT-FE and QE is also tabulated. System | QE | DFT-FE | Diff in Energy
---|---|---|---
|  $E_{g}$ (Ha) | # outer |  $E_{g}$ (Ha) | # outer | (Ha/atom)
|  | iterations |  | iterations |
Pt- Au dimer | -258.3823 | 10 | -258.3826 | 6 | 1.89E-04
Benzamide | -69.9217 | 8 | -69.9208 | 6 | 5.625E-05
Pt - 19 | -2303.1141 | 14 | -2303.1144 | 11 | 1.38E-05
Pt - 38 | -4606.269 | 27 | -4606.2655 | 22 | 9.22E-05

###  4.2 Performance

The performance of our numerical implementation was measured on Stampede2’s
KNL nodes. Each Stampeded2 node has 96 GB of memory and 68 cores with four
hardware threads. In our simulations, 24-36 MPI tasks per node were used to
ensure sufficient memory is available per task. These simulations are
performed using discretization parameters providing $\mathcal{O}(10^{-4})$
Ha/atom accuracy in ground-state energy. In order to demonstrate the
computational efficiency, we choose three Pt nanoclusters of increasing sizes,
containing 38 atoms (684 electrons), 63 atoms (1,134 electrons), and 88 atoms
(1,584 electrons). Further, we consider a large $\text{TiO}_{\text{2}}$
cluster with 3,648 electrons.

We perform the full ground-state calculations for the benchmark Pt
nanoclusters. The walltimes are measured using a tolerance of $10^{-5}$
Ha/atom in the convergence of outer SCF iteration for both QE and DFT-FE
calculations. The performance of QE is heavily dependent on the band group
parallelization. However, the extent of the band group parallelization was
limited by the available memory. The comparison of walltimes for the
construction of one ACE operator are shown in Fig. 9(a). Figure 9(b) shows the
walltime comparison for the full ground-state calculations. We note that both
DFT-FE and QE simulations are performed on the same number of nodes, thus
walltimes are also representative of the relative computational efficiency.
The Pt-38 simulations are run on 64 nodes, whereas the Pt-63 and Pt-88
simulations are run on 128 nodes. The speedup in DFT-FE relative to QE for the
construction of ACE operator increases with increasing system size, reaching
7.6$\times$ for Pt-88. This translates to substantial improvement in the
computational efficiency of the full ground-state calculation as well,
resulting in $>$3$\times$ speedup for Pt-88.

(a) Schematic of Benzamide system

(b) Schematic of Pt-19 atoms system; the bonds are shown for schematic
purposed

(c) Schematic of $\text{TiO}_{\text{2}}$ cluster

Figure 8: Schematics of different benchmark systems used in the study

(a) Wall time of updateExchange for different Pt clusters

(b) Wall time of ground state calculation for different Pt clusters

Figure 9: Wall time measurements for different Pt clusters

The largest system for which the performance is measured is the
$\text{TiO}_{\text{2}}$ system containing 480 atoms (depicted in Fig. 8(c),
3,648 electrons and 1,900 wavefunctions). Both QE and DFT-FE calculations are
performed on 192 nodes. As this is an expensive calculation, we only
considered the comparison of the walltimes for the construction of the ACE
operator. The results of the simulation are reported in Table 2, which shows a
remarkable 11.3$\times$ speed-up over QE.

Table 2: Comparison of walltimes for ACE construction for a 480 atom $\text{TiO}_{\text{2}}$ cluster. System | No of Electrons | QE | DFT-FE | Speed-up
---|---|---|---|---
|  | $(s)$ | $(s)$ |
$\text{TiO}_{\text{2}}$ | 3,648 | 21,390 | 1,892 | 11.3

###  4.3 Parallel scalability

####  4.3.1 Strong scaling

We use the Pt-38 nanocluster for demonstrating the parallel scalability of our
implementation. The strong scaling is measured on nodes ranging from 16 nodes
to 128 nodes. Figure 10(a) shows the strong scaling of updateExchange.
Notably, the walltime for this operation reduced from 380 seconds on 16 nodes
to 56.8 seconds on 128 nodes, translating to a parallel efficiency of 84%. We
note that the number of finite-element nodes in $\mathcal{U}-{\rm FE}$ mesh in
some processors is as low as 330 at 120 nodes, which demonstrates the
excellent parallel scaling of the implementation.

(a) Strong scaling study of updateExchange

(b) Performance breakdown for updateExchange for Pt 38 atom cluster on
different number of nodes

Figure 10: Scaling analysis of updateExchange for Pt 38 atom cluster

####  4.3.2 Performance breakdown

We also measured the time taken for various steps in the updateExchange
operation. Since our implementation includes asynchronous compute by
overlapping some aspects of compute and communication, the exact time taken by
communication and computation can not be measured. Hence we measure the time
taken by processor 0 during the computation steps and subtract it from the
total time to estimate the time taken in communication that is not being
overlapped. This data is presented in Fig. 10(b). We find that communication,
Tucker decomposition, and convolution constitute the major fraction of
walltime, and this is consistent across various number of nodes. Notably, the
fraction of walltime for these steps does not change significantly as the
number of nodes are increased, also demonstrating efficient parallel scaling
across the various steps.

##  5 Conclusion

In this work, we employ a Tucker-tensor based approach to accelerate the
evaluation of Fock exchange in a real-space finite-element basis, which is a
systematically improvable and complete basis. The use of Tucker-tensor
decomposition techniques allows the transformation of 3-D convolution
integrals in the Fock exchange into a tensor product of 1-D convolution
integrals, which results in a significant reduction in computational cost.
This is achieved through a low-rank Tucker-tensor decomposition of the
$\psi_{i}\psi_{j}$ products and a canonical decomposition of the $\frac{1}{r}$
kernel appearing in the Fock exchange. In addition, we have implemented a
task/memory parallelization strategy over the $\psi_{i}\psi_{j}$ pairs that
balance the peak memory and communication costs, in conjunction with
strategies of overlapping compute-and-communication to reduce the
communication cost and maintain good parallel scalability. While the proposed
approach is compatible with any tensor structured basis, we have implemented
the approach in DFT-FE—a massively parallel real-space DFT code based on
adaptive higher order finite-element discretization—extending the capability
of DFT-FE to treat hybrid XC functionals.

The accuracy and efficiency of the implementation have been demonstrated by
comparing it against Quantum Espresso. The ground-state energies between DFT-
FE and Quantum Espresso agree to $\sim 10^{-4}$ Ha/atom for hybrid XC
functional calculations using PBE0 on various benchmark systems. In terms of
computational efficiency, the action of exchange operator on single-particle
wavefunctions demonstrates a speedup of up to 11$\times$ in comparison to
Quantum Espresso for a $\text{TiO}_{\text{2}}$ cluster comprising of 3,684
electrons. Further, the implementation exhibits good parallel scalability,
even in the extreme parallel scaling regime. In the future, we aim to extend
the present implementation to take advantage of GPU acceleration, alongside
extending the formulation to periodic systems.

## Acknowledgments

We gratefully acknowledge the support of Toyota Research Institute under the
auspices of which this work was performed. This work used Stampede2 from the
ACCESS program, which is supported by National Science Foundation grants
#2138259, #2138286, #2138307, #2137603, and #2138296.

## References

  * Hohenberg and Kohn 1964 Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. _Phys. Rev._ 1964, _136_ , B864–B871
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. _Phys. Rev._ 1965, _140_ , A1133–A1138
  * Martin 2020 Martin, R. M. _Electronic structure: basic theory and practical methods_ ; Cambridge university press, 2020
  * Hafner et al. 2006 Hafner, J.; Wolverton, C.; Ceder, G. Toward computational materials design: the impact of density functional theory on materials research. _MRS bulletin_ 2006, _31_ , 659–668
  * Becke 2014 Becke, A. D. Perspective: Fifty years of density-functional theory in chemical physics. _The Journal of chemical physics_ 2014, _140_ , 18A301
  * Mardirossian and Head-Gordon 2017 Mardirossian, N.; Head-Gordon, M. Thirty years of density functional theory in computational chemistry: an overview and extensive assessment of 200 density functionals. _Molecular Physics_ 2017, _115_ , 2315–2372
  * Kim et al. 2006 Kim, S.; Lee, J. K.; Kang, S. O.; Ko, J.; Yum, J.-H.; Fantacci, S.; De Angelis, F.; Di Censo, D.; Nazeeruddin, M. K.; Grätzel, M. Molecular engineering of organic sensitizers for solar cell applications. _Journal of the American Chemical Society_ 2006, _128_ , 16701–16707
  * Jhi et al. 2000 Jhi, S.-H.; Louie, S. G.; Cohen, M. L. Electronic properties of oxidized carbon nanotubes. _Physical Review Letters_ 2000, _85_ , 1710
  * Kohn and Sham 1965 Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. _Physical review_ 1965, _140_ , A1133
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868
  * Becke 1993 Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. _The Journal of Chemical Physics_ 1993, _98_ , 5648–5652
  * Perdew et al. 1996 Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. _The Journal of Chemical Physics_ 1996, _105_ , 9982–9985
  * Heyd et al. 2003 Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. _The Journal of Chemical Physics_ 2003, _118_ , 8207–8215
  * Krukau et al. 2006 Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. _The Journal of chemical physics_ 2006, _125_
  * Finazzi et al. 2008 Finazzi, E.; Di Valentin, C.; Pacchioni, G.; Selloni, A. Excess electron states in reduced bulk anatase $TiO_{2}$: comparison of standard GGA, GGA+ U, and hybrid DFT calculations. _The Journal of chemical physics_ 2008, _129_ , 154113
  * Da Silva et al. 2007 Da Silva, J. L.; Ganduglia-Pirovano, M. V.; Sauer, J.; Bayer, V.; Kresse, G. Hybrid functionals applied to rare-earth oxides: The example of ceria. _Physical Review B_ 2007, _75_ , 045121
  * Garza and Scuseria 2016 Garza, A. J.; Scuseria, G. E. Predicting band gaps with hybrid density functionals. _The journal of physical chemistry letters_ 2016, _7_ , 4165–4170
  * Heyd et al. 2005 Heyd, J.; Peralta, J. E.; Scuseria, G. E.; Martin, R. L. Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional. _The Journal of chemical physics_ 2005, _123_ , 174101
  * Muscat et al. 2001 Muscat, J.; Wander, A.; Harrison, N. On the prediction of band gaps from hybrid functional theory. _Chemical Physics Letters_ 2001, _342_ , 397–401
  * Salzner et al. 1997 Salzner, U.; Lagowski, J.; Pickup, P.; Poirier, R. Design of low band gap polymers employing density functional theory—hybrid functionals ameliorate band gap problem. _Journal of computational chemistry_ 1997, _18_ , 1943–1953
  * Xiao et al. 2011 Xiao, H.; Tahir-Kheli, J.; Goddard III, W. A. Accurate band gaps for semiconductors from density functional theory. _The Journal of Physical Chemistry Letters_ 2011, _2_ , 212–217
  * Seidl et al. 1996 Seidl, A.; Görling, A.; Vogl, P.; Majewski, J. A.; Levy, M. Generalized Kohn-Sham schemes and the band-gap problem. _Physical Review B_ 1996, _53_ , 3764
  * Lin 2016 Lin, L. Adaptively compressed exchange operator. _Journal of chemical theory and computation_ 2016, _12_ , 2242–2249
  * Hu et al. 2020 Hu, W.; Liu, J.; Li, Y.; Ding, Z.; Yang, C.; Yang, J. Accelerating excitation energy computation in molecules and solids within linear-response time-dependent density functional theory via interpolative separable density fitting decomposition. _Journal of Chemical Theory and Computation_ 2020, _16_ , 964–973
  * Hu et al. 2017 Hu, W.; Lin, L.; Yang, C. Interpolative Separable Density Fitting Decomposition for Accelerating Hybrid Density Functional Calculations with Applications to Defects in Silicon. _Journal of Chemical Theory and Computation_ 2017, _13_ , 5420–5431, PMID: 28960982
  * Ko et al. 2020 Ko, H.-Y.; Jia, J.; Santra, B.; Wu, X.; Car, R.; DiStasio Jr, R. A. Enabling large-scale condensed-phase hybrid density functional theory based ab initio molecular dynamics. 1. theory, algorithm, and performance. _Journal of Chemical Theory and Computation_ 2020, _16_ , 3757–3785
  * Ko et al. 2021 Ko, H.-Y.; Santra, B.; DiStasio Jr, R. A. Enabling Large-Scale Condensed-Phase Hybrid Density Functional Theory-Based Ab Initio Molecular Dynamics II: Extensions to the Isobaric–Isoenthalpic and Isobaric–Isothermal Ensembles. _Journal of Chemical Theory and Computation_ 2021, _17_ , 7789–7813
  * 28 Rettig, A.; Lee, J.; Head-Gordon, M. Even Faster Exact Exchange for Solids via Tensor Hypercontraction. _Journal of Chemical Theory and Computation_ PMID: 37586065
  * Apra et al. 2020 Apra, E.; Bylaska, E. J.; De Jong, W. A.; Govind, N.; Kowalski, K.; Straatsma, T. P.; Valiev, M.; van Dam, H. J.; Alexeev, Y.; Anchell, J.; others NWChem: Past, present, and future. _The Journal of chemical physics_ 2020, _152_ , 184102
  * Blum et al. 2009 Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab initio molecular simulations with numeric atom-centered orbitals. _Computer Physics Communications_ 2009, _180_ , 2175–2196
  * Frisch 2009 Frisch, A. gaussian 09W Reference. _Wallingford, USA, 25p_ 2009,
  * Giannozzi et al. 2009 Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; others QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _Journal of physics: Condensed matter_ 2009, _21_ , 395502
  * Gonze et al. 2009 Gonze, X.; Amadon, B.; Anglade, P.-M.; Beuken, J.-M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Côté, M.; Deutsch, T.; Genovese, L.; Ghosez, P.; Giantomassi, M.; Goedecker, S.; Hamann, D.; Hermet, P.; Jollet, F.; Jomard, G.; Leroux, S.; Mancini, M.; Mazevet, S.; Oliveira, M.; Onida, G.; Pouillon, Y.; Rangel, T.; Rignanese, G.-M.; Sangalli, D.; Shaltaf, R.; Torrent, M.; Verstraete, M.; Zerah, G.; Zwanziger, J. ABINIT: First-principles approach to material and nanosystem properties. _Computer Physics Communications_ 2009, _180_ , 2582–2615, cited By 1869
  * Hafner 2008 Hafner, J. Ab-initio simulations of materials using VASP: Density-functional theory and beyond. _Journal of computational chemistry_ 2008, _29_ , 2044–2078
  * Motamarri et al. 2013 Motamarri, P.; Nowak, M. R.; Leiter, K.; Knap, J.; Gavini, V. Higher-order adaptive finite-element methods for Kohn–Sham density functional theory. _Journal of Computational Physics_ 2013, _253_ , 308–343
  * Pask and Sterne 2005 Pask, J.; Sterne, P. Finite element methods in ab initio electronic structure calculations. _Modelling and Simulation in Materials Science and Engineering_ 2005, _13_ , R71
  * Motamarri et al. 2020 Motamarri, P.; Das, S.; Rudraraju, S.; Ghosh, K.; Davydov, D.; Gavini, V. DFT-FE–A massively parallel adaptive finite-element code for large-scale density functional theory calculations. _Computer Physics Communications_ 2020, _246_ , 106853
  * Das et al. 2022 Das, S.; Motamarri, P.; Subramanian, V.; Rogers, D. M.; Gavini, V. DFT-FE 1.0: A massively parallel hybrid CPU-GPU density functional theory code using finite-element discretization. _Computer Physics Communications_ 2022, _280_ , 108473
  * Das et al. 2019 Das, S.; Motamarri, P.; Gavini, V.; Turcksin, B.; Li, Y. W.; Leback, B. Fast, scalable and accurate finite-element based ab initio calculations using mixed precision computing: 46 PFLOPS simulation of a metallic dislocation system. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’19). 2019; pp 1–11
  * Das et al. 2023 Das, S.; Bikash, K.; Subramanian, V.; Panigrahi, G.; Motamarri, P.; Rogers, D.; Zimmerman, P. M.; Gavini, V. Large-scale materials modeling at quantum accuracy: Ab initio simulations of quasicrystals and interacting extended defects in metallic alloys. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC’23). 2023; pp 1–12
  * Zhuravel et al. 2020 Zhuravel, R.; Huang, H.; Polycarpou, G.; Polydorides, S.; Motamarri, P.; Katrivas, L.; Rotem, D.; Sperling, J.; Zotti, L. A.; Kotlyar, A. B.; others Backbone charge transport in double-stranded DNA. _Nature Nanotechnology_ 2020, _15_ , 836–840
  * Kumar et al. 2023 Kumar, P.; Ludhwani, M. M.; Das, S.; Gavini, V.; Kanjarla, A.; Adlakha, I. Effect of hydrogen on plasticity of $\alpha$-Fe: A multi-scale assessment. _International Journal of Plasticity_ 2023, _165_ , 103613
  * Yao et al. 2022 Yao, L.; Das, S.; Liu, X.; Wu, K.; Cheng, Y.; Gavini, V.; Xiao, B. Modulating the microscopic lattice distortions through the Al-rich layers for boosting the ferroelectricity in Al: HfO2 nanofilms. _Journal of Physics D: Applied Physics_ 2022, _55_ , 455501
  * Ghosh et al. 2019 Ghosh, K.; Ma, H.; Gavini, V.; Galli, G. All-electron density functional calculations for electron and nuclear spin interactions in molecules and solids. _Physical Review Materials_ 2019, _3_ , 043801
  * Krishnendu et al. 2021 Krishnendu, G.; He, M.; Mykyta, O.; Vikram, G.; Giulia, G. Spin–spin interactions in defects in solids from mixed all-electron and pseudopotential first-principles calculations. _NPJ Computational Materials_ 2021, _7_
  * Kanungo et al. 2019 Kanungo, B.; Zimmerman, P. M.; Gavini, V. Exact exchange-correlation potentials from ground-state electron densities. _Nature Communications_ 2019, _10_ , 1–9
  * Kanungo et al. 2021 Kanungo, B.; Zimmerman, P. M.; Gavini, V. A comparison of exact and model exchange–correlation potentials for molecules. _The Journal of Physical Chemistry Letters_ 2021, _12_ , 12012–12019
  * Kanungo et al. 2023 Kanungo, B.; Hatch, J.; Zimmerman, P. M.; Gavini, V. Exact and model exchange–correlation potentials for open shell systems. _The Journal of Physical Chemistry Letters_ 2023, _14_ , 10039–10045
  * Khoromskaia and Khoromskij 2015 Khoromskaia, V.; Khoromskij, B. N. Tensor numerical methods in quantum chemistry: from Hartree–Fock to excitation energies. _Physical Chemistry Chemical Physics_ 2015, _17_ , 31491–31509
  * Khoromskij and Khoromskaia 2007 Khoromskij, B.; Khoromskaia, V. Low rank Tucker-type tensor approximation to classical potentials. _Open Mathematics_ 2007, _5_ , 523–550
  * Khoromskij et al. 2009 Khoromskij, B. N.; Khoromskaia, V.; Chinnamsetty, S. R.; Flad, H.-J. Tensor decomposition in electronic structure calculations on 3D Cartesian grids. _Journal of computational physics_ 2009, _228_ , 5749–5762
  * Braess and Hackbusch 2009 Braess, D.; Hackbusch, W. _Multiscale, Nonlinear and Adaptive Approximation: Dedicated to Wolfgang Dahmen on the Occasion of his 60th Birthday_ ; Springer, 2009; pp 39–74
  * Giannozzi et al. 2017 Giannozzi, P.; Andreussi, O.; Brumme, T.; Bunau, O.; Nardelli, M. B.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Cococcioni, M.; others Advanced capabilities for materials modelling with Quantum ESPRESSO. _Journal of physics: Condensed matter_ 2017, _29_ , 465901
  * Levy 1979 Levy, M. Universal variational functionals of electron densities, first-order density matrices, and natural spin-orbitals and solution of the v-representability problem. _Proceedings of the National Academy of Sciences_ 1979, _76_ , 6062–6065
  * Garrick et al. 2020 Garrick, R.; Natan, A.; Gould, T.; Kronik, L. Exact Generalized Kohn-Sham Theory for Hybrid Functionals. _Physical Review X_ 2020, _10_ , 021040
  * Ren et al. 2012 Ren, X.; Rinke, P.; Blum, V.; Wieferink, J.; Tkatchenko, A.; Sanfilippo, A.; Reuter, K.; Scheffler, M. Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions. _New Journal of Physics_ 2012, _14_ , 053020
  * Levchenko et al. 2015 Levchenko, S. V.; Ren, X.; Wieferink, J.; Johanni, R.; Rinke, P.; Blum, V.; Scheffler, M. Hybrid functionals for large periodic systems in an all-electron, numeric atom-centered basis framework. _Computer Physics Communications_ 2015, _192_ , 60–69
  * Anderson 1965 Anderson, D. G. Iterative procedures for nonlinear integral equations. _Journal of the ACM (JACM)_ 1965, _12_ , 547–560
  * Zhou et al. 2006 Zhou, Y.; Saad, Y.; Tiago, M. L.; Chelikowsky, J. R. Self-consistent-field calculations using Chebyshev-filtered subspace iteration. _Journal of Computational Physics_ 2006, _219_ , 172–184
  * Hu et al. 2017 Hu, W.; Lin, L.; Banerjee, A. S.; Vecharynski, E.; Yang, C. Adaptively compressed exchange operator for large-scale hybrid density functional calculations with applications to the adsorption of water on silicene. _Journal of chemical theory and computation_ 2017, _13_ , 1188–1198
  * Kolda and Bader 2009 Kolda, T. G.; Bader, B. W. Tensor decompositions and applications. _SIAM review_ 2009, _51_ , 455–500
  * Khoromskij et al. 2009 Khoromskij, B.; Khoromskaia, V.; Chinnamsetty, S.; Flad, H.-J. Tensor decomposition in electronic structure calculations on 3D Cartesian grids. _Journal of Computational Physics_ 2009, _228_ , 5749–5762
  * Blesgen et al. 2012 Blesgen, T.; Gavini, V.; Khoromskaia, V. Approximation of the electron density of Aluminium clusters in tensor-product format. _Journal of Computational Physics_ 2012, _231_ , 2551–2564
  * Lin et al. 2021 Lin, C.-C.; Motamarri, P.; Gavini, V. Tensor-structured algorithm for reduced-order scaling large-scale Kohn–Sham density functional theory calculations. _npj Computational Materials_ 2021, _7_ , 50
  * Motamarri et al. 2016 Motamarri, P.; Gavini, V.; Blesgen, T. Tucker-tensor algorithm for large-scale Kohn-Sham density functional theory calculations. _Phys. Rev. B_ 2016, _93_ , 125104
  * Hughes 2012 Hughes, T. J. _The finite element method: linear static and dynamic finite element analysis_ ; Courier Corporation, 2012
  * Tsuchida and Tsukada 1995 Tsuchida, E.; Tsukada, M. Electronic-structure calculations based on the finite-element method. _Physical Review B_ 1995, _52_ , 5573
  * Tsuchida and Tsukada 1996 Tsuchida, E.; Tsukada, M. Adaptive finite-element method for electronic-structure calculations. _Physical Review B_ 1996, _54_ , 7602
  * Ballard et al. 2020 Ballard, G.; Klinvex, A.; Kolda, T. G. TuckerMPI: A parallel C++/MPI software package for large-scale data compression via the tucker tensor decomposition. _ACM Transactions on Mathematical Software (TOMS)_ 2020, _46_ , 1–31
  * Hamann 2013 Hamann, D. R. Optimized norm-conserving Vanderbilt pseudopotentials. _Phys. Rev. B_ 2013, _88_ , 085117


|  $l\neq s$: $\downarrow$ |  $l\neq s$: $\downarrow$ if $(\mu_{s}-\bar{\mu})(\mu_{l}-\bar{\mu})>0$
$\beta_{s}\downarrow$ |  $l=s$: $\uparrow$ if $\mu_{s}<\bar{\mu}$ | no simple rule
|  $l\neq s$: $\downarrow$ if $\mu_{s}+\mu_{l}>\bar{\mu}$ |  $\downarrow$ if $[\sigma_{s}^{2}+\mu_{s}(\mu_{s}-\bar{\mu})](\mu_{l}-\bar{\mu})>0$
Table EC.1: Sensitivity of results for bank $l$ in response to a decrease in
parameter $\mu_{s}$, $\alpha_{s}$, or $\beta_{s}$ for bank $s$. Sensitivities
shown are for predicted loss $\hat{Y}_{l}(x)$ (top), mean predicted loss
$\mathsf{E}[\hat{Y}(X_{l})]$ (middle), and the bias
$\mathsf{E}[\hat{Y}(X_{l})-Y_{l}]$.

##  Appendix EC.2 Sensitivity Analysis

This section provides supporting details for Section EC.1, particularly the
conclusions summarized in the middle and bottom panels of Table EC.1. We begin
with an analysis of forecast bias that is of independent interest.

###  EC.2.1 Forecast Bias

If losses at different banks are described by different models, then forecast
bias becomes inevitable when we apply a single model to all banks. But the
distribution of bias across banks may differ under different choices of the
single model.

Let $\hat{Y}_{s}$ be any of the forecasts for bank $s$ in Table 4.1, and, as
in (1), let $Y_{s}$ denote the actual loss rate for bank $s$. Both
$\hat{Y}_{s}$ and $Y_{s}$ are evaluated at $X_{s}$. Define the forecast bias
for bank $s$ to be

$\mathsf{bias}(s)=\mathsf{E}[\hat{Y}_{s}-Y_{s}].$ (70)

The expectation integrates over the distribution of the error $\epsilon_{s}$
in (1) and the features $X_{s}$.

######  Proposition EC.2.1.

For each forecast in Table 4.1, the bias is as follows.

  * (i)

Pooled:
$\mathsf{bias}(s)=\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]+\beta_{Pool}^{\top}(\mu_{s}-\bar{\mu})$;

  * (ii)

PTF in (18): $\mathsf{bias}(s)=\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]$;

  * (iii)

Conditional expectation:
$\mathsf{bias}(s)=\mathsf{E}[\hat{Y}_{C}(X_{s})]-\mathsf{E}[Y_{s}]$;

  * (iv)

FEO:
$\mathsf{bias}(s)=\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]+\beta_{F}^{\top}(\mu_{s}-\bar{\mu})$;

  * (v)

SEO: $\mathsf{bias}(s)=\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]$.

###### Proof.

For (i), we have, using the definition of $\alpha_{Pool}$ in (12),

$\displaystyle\mathsf{E}[\hat{Y}_{s}-Y_{s}]$ $\displaystyle=$
$\displaystyle\mathsf{E}[\alpha_{Pool}+\beta^{\top}_{Pool}X_{s}-Y_{s}]$
$\displaystyle=$
$\displaystyle(\mathsf{E}[Y_{S}]-\beta^{\top}_{Pool}\bar{\mu})+\beta^{\top}_{Pool}\mu_{s}-\mathsf{E}[Y_{s}]$
$\displaystyle=$
$\displaystyle\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]+\beta^{\top}_{Pool}(\mu_{s}-\bar{\mu}).$

For the PTF forecast, (17) and (18) yield

$\mathsf{E}[\hat{Y}_{s}]=\bar{\alpha}^{o}=\sum_{s}p_{s}\alpha_{s}^{o}=\sum_{s}p_{s}\mathsf{E}[Y_{s}]=\mathsf{E}[Y_{S}],$

and the bias in (ii) follows. The expression in (iii) holds by definition. The
argument for (iv) is the same as the argument for (i). The bias in (v) follows
from (iv) because we see from (38) that the SEO forecast for bank $s$
subtracts $\beta^{\top}_{F}(\mu_{s}-\bar{\mu})$ from the FEO forecast. ∎

In every case of Proposition EC.2.1, the average bias
$\sum_{s}p_{s}\mathsf{bias}(s)$ is zero, but the methods differ in how they
distribute bias across banks. We saw previously that the PTF and SEO methods
go the farthest in equalizing differences; we now see that the bias for each
of these methods is the difference $\mathsf{E}[Y_{S}]-\mathsf{E}[Y_{s}]$
between the average loss rate for all banks and the average for an individual
bank.

Using the relationship $\beta_{Pool}=\beta_{F}+\Lambda\delta$ from (30), we
see that the difference between the expressions in (i) and (iv) is

$\mathsf{bias}_{Pool}(s)-\mathsf{bias}_{FEO}(s)=\delta^{\top}\Lambda^{\top}(\mu_{s}-\bar{\mu}).$

In light of the discussion in Section 4.2, this difference is the expected
disparate impact on bank $s$ of using the pooled model.

###  EC.2.2 Improvement in Intercept $\alpha_{s}$

By taking the expectation of (66), we get

$\mathsf{E}[\hat{Y}_{F}(X_{l})]=\sum_{s}ps(\alpha_{s}+\beta_{s}\mu_{s})+\beta_{F}(\mu_{l}-\bar{\mu}),$
(71)

and the same holds for the expected pooled forecast with $\beta_{F}$ replaced
by $\beta_{Pool}$. It follows that, for any banks $s$ and $l$,

$\frac{\partial{\mathsf{E}[\hat{Y}_{F}(X_{l})]}}{\partial{\alpha_{s}}}=p_{s}>0.$

In other words, all expected forecasts decrease following a reduction in
$\alpha_{s}$.

In contrast, for the pooled model we get

$\frac{\partial{\mathsf{E}[\hat{Y}_{P}(X_{l})]}}{\partial{\alpha_{s}}}=p_{s}-\frac{{\partial\mathsf{cov}(\alpha_{S},\mu_{S})}/{\partial\alpha_{s}}}{\sum_{t}p_{t}\sigma_{t}^{2}+\mathsf{var}(\mu_{S})}\beta_{s}(\bar{\mu}-\mu_{l})=p_{s}+p_{s}\frac{(\mu_{s}-\bar{\mu})(\mu_{l}-\bar{\mu})}{\sum_{s}p_{s}\sigma_{s}^{2}+\mathsf{var}(\mu_{S})}\beta_{s}.$
(72)

Bank $s$ benefits from its reduction of $\alpha_{s}$, in the sense that the
derivative with $l=s$ is positive. For $l\neq s$, the sign of (72) does not
admit a simple description. In particular, it may be negative when $\mu_{s}$
and $\mu_{l}$ are on opposite sides of $\bar{\mu}$, meaning that one bank’s
loans are riskier than average and the other bank’s loans are less risky than
average.

For the bias under FEO we have

$\frac{\partial{\mathsf{bias}_{F}(l)}}{\partial{\alpha_{s}}}=p_{s}-\mathbf{1}\\{l=s\\}$

It is then immediate that

$\frac{\partial{\mathsf{bias}_{F}(l)}}{\partial{\alpha_{s}}}>0\text{ if }l\neq
s\quad\mbox{and}\quad\frac{\partial{\mathsf{bias}_{F}(s)}}{\partial{\alpha_{s}}}<0.$

The direction of change makes sense. If the bias for a bank is positive,
meaning that the industry model overestimates its losses, then improvements at
other banks will reduce loss forecasts and thus reduce the bias. The bank’s
own improvements will increase the bias by reducing the bank’s own losses by
more than they reduce the model’s forecasts. The situation is reversed for a
bank with a negative bias.

However, for the pooled regression method,

$\frac{\partial{\mathsf{bias}_{P}(l)}}{\partial{\alpha_{s}}}=p_{s}+p_{s}\frac{(\mu_{s}-\bar{\mu})(\mu_{l}-\bar{\mu})}{\sum_{s}p_{s}\sigma_{s}^{2}+\mathsf{var}(\mu_{S})}\beta_{s}-\mathbf{1}\\{l=s\\},$

and the direction of change is unclear.

###  EC.2.3 Improvement in Loan Quality

Now suppose bank $s$ improves the quality of its loan portfolio, resulting in
a smaller $\mu_{s}$. This has no effect on $\beta_{F}$, which makes sense —
changing one bank’s loan quality should not change the sensitivity of losses
to loan quality. However, it is evident from (11) that $\beta_{Pool}$ does
change with $\mu_{s}$.

Under FEO, the mean the mean predicted loss rate satisfies

$\frac{\partial{\mathsf{E}\hat{Y}_{F}(X_{l})}}{\partial{\mu_{s}}}=p_{s}\beta_{s}+\beta_{F}(\mathbf{1}\\{l=s\\}-p_{s}),$

which is always positive if $l=s$. This means that an improvement in bank
$l$’s loan quality (a reduction in $\mu_{l}$) reduces bank $l$’s mean
predicted losses. In the pooled model,

$\frac{\partial{\mathsf{E}\hat{Y}_{P}(X_{l})}}{\partial{\mu_{s}}}=p_{s}\beta_{s}+\beta_{Pool}(\mathbf{1}\\{l=s\\}-p_{s})+(\mu_{s}-\bar{\mu})\frac{\partial{\beta_{Pool}}}{\partial{\mu_{s}}};$

this expression could be negative, even with $l=s$, meaning that a bank could
be penalized (through a higher mean predicted loss rate) as a result of
improving its loan quality.

The sensitivity of the bias under FEO is given by

$\frac{\partial{\mathsf{bias}_{F}(l)}}{\partial{\mu_{s}}}=(\mathbf{1}\\{l=s\\}-p_{s})(\beta_{F}-\beta_{s});$

in particular, the bias for bank $l$ moves in opposite directions with respect
to changes in $\mu_{l}$ and $\mu_{s}$, $s\not=l$. Suppose industry model
overestimates bank $l$’s losses, in the sense that the bias is positive, and
suppose the industry model overestimates bank $l$’s sensitivity to loan
quality, in the sense that $\beta_{F}>\beta_{l}$. Then bank $l$ will benefit
(in the sense of reducing the bias) from improving its loan quality by
reducing $\mu_{l}$.

For the pooled regression,

$\frac{\partial{\mathsf{bias}_{P}(l)}}{\partial{\mu_{s}}}=(\beta_{Pool}-\beta_{s})(\mathbf{1}\\{l=s\\}-p_{s})+(\mu_{s}-\bar{\mu})\frac{\partial{\beta_{Pool}}}{\partial{\mu_{s}}}.$

The sign of this expression does not admit a simple condition.

###  EC.2.4 Improvement in Loan Management

Now suppose bank $s$ improves its abilities in loan management, resulting in a
reduction in $\beta_{s}$. The mean predicted loss rate under FEO satisfies

$\frac{\partial{\mathsf{E}[\hat{Y}_{F}(X_{l})]}}{\partial{\beta_{s}}}=p_{s}(\mu_{s}+\mu_{l}-\bar{\mu}),$

and is positive if $\mu_{s}+\mu_{l}>\bar{\mu}$. In the pooled model

$\frac{\partial{\mathsf{E}[\hat{Y}_{P}(X_{l})]}}{\partial{\beta_{s}}}=p_{s}\mu_{s}+\frac{p_{s}(\sigma^{2}_{s}+\mu_{s}(\mu_{s}-\bar{\mu}))}{\sum_{i}p_{i}(\sigma^{2}_{i}+\mu_{i}(\mu_{i}-\bar{\mu}))}(\mu_{l}-\bar{\mu}),$

so $[\sigma^{2}_{s}+\mu_{s}(\mu_{s}-\bar{\mu})](\mu_{l}-\bar{\mu})>0$ is a
sufficient condition for the sensitivity to be positive. Regardless of the
value of $\mu_{l}$, the sensitivity is positive for all sufficiently large
$\mu_{s}$.

For $l\not=s$, the sensitivity of the bias for bank $l$ with respect to
$\beta_{s}$ equals the sensitivity of the mean predicted loss because the
actual expected loss $\mathsf{E}[Y_{l}]$ is unaffected by $\beta_{s}$. We
therefore focus on the case $l=s$. Under FEO,

$\frac{\partial{\mathsf{bias}_{F}(s)}}{\partial{\beta_{s}}}=(p_{s}-1)\mu_{s}+p_{s}(\mu_{s}-\bar{\mu}),$

which is guaranteed to be negative if $\mu_{s}<\bar{\mu}$. Under the pooled
model, the sign of

$\frac{\partial{\mathsf{bias}_{P}(s)}}{\partial{\beta_{s}}}=(p_{s}-1)\mu_{s}+\frac{p_{s}(\sigma^{2}_{s}+\mu_{s}(\mu_{s}-\bar{\mu}))}{\sum_{i}p_{i}(\sigma^{2}_{i}+\mu_{i}(\mu_{i}-\bar{\mu}))}(\mu_{s}-\bar{\mu})$

does not admit a simple characterization.

##  Appendix EC.3 Convex Combinations of Coefficients

Equation (13) aggregates the individual scalar slopes $\beta_{s}$ into a
single value. We can generalize this perspective and ask what properties we
would like in an aggregation function, meaning a function
$f:\mathbb{R}^{\bar{S}}\to\mathbb{R}$,

$\beta_{*}=f(\beta_{1},\dots,\beta_{\bar{S}}),$

that combines bank-specific coefficients $\beta_{s}$ into an “industry”
parameter $\beta_{*}$.

We consider the following properties:

  * (i)

$f(kb_{1},\dots,kb_{\bar{S}})=kf(b_{1},\dots,b_{\bar{S}})$, for all $k$,
$b_{1},\dots,b_{\bar{S}}\in\mathbb{R}$;

  * (ii)

$f(b,\dots,b)=b$, for at least one nonzero $b\in\mathbb{R}$;

  * (iii)

$b_{s}>0$, for all $s$, implies $f(b_{1},\dots,b_{\bar{S}})\geq 0$;

  * (iv)

$f$ is differentiable at zero.

Property (i) is needed for the aggregation to perform sensibly under a change
of units in the measurement of $X_{s}$: if we divide each $X_{s}$ by $k$, each
$\beta_{s}$ increases by a factor of $k$, and it is natural to require that
$\beta_{*}$ scale accordingly. Properties (ii) and (iii) are also very modest
requirements. Property (iv) is harder to motivate but not unreasonable. These
properties constrain the aggregation function as follows:

######  Proposition EC.3.1.

If (i)–(iv) hold, then $f(\beta_{1},\dots,\beta_{\bar{S}})$ is a convex
combination of its arguments.

###### Proof.

Fix $\beta\in\mathbb{R}^{\bar{S}}$. Let $g(t)=f(t\beta)$. By condition (iv),
$g^{\prime}(0)=\beta^{\top}f^{\prime}(0)$. Condition (i) and (ii) imply
$g(t)=tf(\beta)$, so $g^{\prime}(t)=f(\beta)$ for any $t$. Thus,
$f(\beta)=g^{\prime}(0)=\beta^{\top}f^{\prime}(0)=\sum_{i=1}^{\bar{S}}f_{i}^{\prime}(0)\beta_{i}$.
Condition (ii) now implies $\sum_{i=1}^{\bar{S}}f_{i}^{\prime}(0)=1$, and
condition (iii) implies $f_{i}^{\prime}(0)\geq 0$, for all $i$. Thus,
$f(\beta)=\sum_{i=1}^{\bar{S}}f_{i}^{\prime}(0)\beta_{i}$ is a convex
combination of the components of $\beta$. ∎

The scalar FEO coefficient in (27) is a convex combination of the bank-
specific coefficients $\beta_{s}$, but the pooled coefficient (13) is
generally not. This property of the FEO model extends to the multivariate case
under additional conditions. If all the bank-specific covariance matrices
$\Sigma_{s}$, $s=1,\dots,\bar{S}$, coincide, then in (25) we get
$\beta_{F}=\mathsf{E}[\beta_{S}]=\sum_{s}p_{s}\beta_{s}$. If all $\Sigma_{s}$
are diagonal (but not necessarily identical), then the representation of the
scalar FEO coefficient in (27) applies to each coordinate of $\beta_{F}$. If
all $\Sigma_{s}$ have the same eigenvectors, then we can transform the
original features $X_{s}$ into uncorrelated features using principal
components. Using these transformed features, each coordinate of $\beta_{F}$
is a convex combination of bank-specific coefficients.

##  Appendix EC.4 Additional Information on Empirical Analysis

Table EC.1 lists the bank holding companies included in our empirical analysis
and the symbols we use to refer to them. The companies are listed in order of
size by total assets.

Ticker | Bank Name
---|---
JPM | JPMORGAN CHASE & CO.
BAC | BANK OF AMERICA CORPORATION
C | CITIGROUP INC.
WFC | WELLS FARGO & COMPANY
GS | GOLDMAN SACHS GROUP, INC., THE
MS | MORGAN STANLEY
SCHW | CHARLES SCHWAB CORPORATION, THE
USB | U.S. BANCORP
PNC | PNC FINANCIAL SERVICES GROUP, INC., THE
TFC | TRUIST FINANCIAL CORPORATION
TD | TD GROUP US HOLDINGS LLC
BK | BANK OF NEW YORK MELLON CORPORATION, THE
COF | CAPITAL ONE FINANCIAL CORPORATION
STT | STATE STREET CORPORATION
HSBC | HSBC NORTH AMERICA HOLDINGS INC.
SIVB | SVB FINANCIAL GROUP
FITB | FIFTH THIRD BANCORP
USAA | UNITED SERVICES AUTOMOBILE ASSOCIATION
BMO | BMO FINANCIAL CORP.
CFG | CITIZENS FINANCIAL GROUP, INC.
AXP | AMERICAN EXPRESS COMPANY
KEY | KEYCORP
NTRS | NORTHERN TRUST CORPORATION
ALLY | ALLY FINANCIAL INC.
AMP | AMERIPRISE FINANCIAL, INC.
RY | RBC US GROUP HOLDINGS LLC
HBAN | HUNTINGTON BANCSHARES INCORPORATED
RF | REGIONS FINANCIAL CORPORATION
MUFG | MUFG AMERICAS HOLDINGS CORPORATION
BCS | BARCLAYS US LLC
SAN | SANTANDER HOLDINGS USA, INC.
MTB | M&T BANK CORPORATION
BNPQY | BNP PARIBAS USA, INC.
DB | DB USA CORPORATION
DFS | DISCOVER FINANCIAL SERVICES
Table EC.1: Symbols and names of included bank holding companies.

We construct the loss rates, past due rates, and allowance rates using the
entries in FR Y-9C forms outlined in Table EC.2.

Variables | Loan Types | 2007Q1 – Present | 2003Q1 – 2006Q4
---|---|---|---
Loan Amount | CC | BHCKB538 | BHCKB538
FL | BHDM5367 | BHDM5367
CRE |  | Owned: BHCKF160
---
Other: BHCKF161
BHDM1480
CI | BHCK1763 | BHCK1763
Total | BHCK2122 | BHCK2122
Charge-Offs | CC | BHCKB514 | BHCKB514
FL | BHCKC234 | BHCKC234
CRE |  | Owned: BHCKC895
---
Other: BHCKC897
BHCK3590
CI | BHCK4645 | BHCK4645
Recoveries | CC | BHCKB515 | BHCKB515
FL | BHCKC217 | BHCKC217
CRE |  | Owned: BHCKC896
---
Other: BHCKC898
BHCK3591
CI | BHCK4617 | BHCK4617
Past Due: 30-89 days and accruing | CC | BHCKB575 | BHCKB575
FL | BHCKC236 | BHCKC236
CRE |  | Owned: BHCKF178
---
Other: BHCKF179
BHCK3502
CI | BHCK1606 | BHCK1606
Past Due: 90 days and accruing | CC | BHCKB576 | BHCKB576
FL | BHCKC237 | BHCKC237
CRE |  | Owned: BHCKF180
---
Other: BHCKF181
BHCK3503
CI | BHCK1607 | BHCK1607
Past Due: non-accrual | CC | BHCKB577 | BHCKB577
FL | BHCKC229 | BHCKC229
CRE |  | Owned: BHCKF182
---
Other: BHCKF183
BHCK3504
CI | BHCK1608 | BHCK1608
Table EC.2: Loan variables and FR Y-9C form correspondence.

##  Appendix EC.5 Robustness Checks

###  EC.5.1 Pre-COVID Data and Allowance Rate

We repeat the analysis of Section 6, limiting the data to 2001–2019. This
serves two purposes. It addresses the possibility that our results are driven
by a few extreme values during the COVID period 2020–2021. It also accounts
for a change in how banks measure allowances (the Current Expected Credit
Losses methodology) that began to take effect at the end of 2019. We also
consider adding allowance rate as another proxy for portfolio risks. Tables
EC.1 and EC.2 report the results for heterogeneity tests and differences of
the parameter estimates under this setting. The evidence for heterogeneity and
its impact is generally at least as strong using the pre-COVID data as using
data through 2021.

Covariance | $\alpha$ | $\beta_{PDR}$ | $\beta_{AR}$ | $\gamma$
---|---|---|---|---
Estimation | CC | FL | CRE | CI | CC | FL | CRE | CI | CC | FL | CRE | CI | CC | FL | CRE | CI
bank clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |  |  |  |  |
time clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |  |  |  |  |
bank clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |
time clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |
bank clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |  | 0.00 | 0.00 | 0.00 | 0.00
time clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |  |  |  |  | 0.00 | 0.00 | 0.00 | 0.00
bank clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
time clustered | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
Table EC.1: Heterogeneity tests using pre-COVID data with allowance rate as an additional proxy for banks’ portfolio risks. Loan | Past Due Rate | Allowance Rate | Macro PC
---|---|---|---
Type | $\beta_{Pool}$ | $\beta_{F}$ | diff |  $p$-value | $\beta_{Pool}$ | $\beta_{F}$ | diff |  $p$-value | $\gamma_{Pool}$ | $\gamma_{F}$ | diff |  $p$-value
CC | 0.781 | 0.838 | -0.057 | 0.006∗∗∗ |  |  |  |  |  |  |  |
FL | 0.045 | 0.061 | -0.016 | 0.002∗∗∗ |  |  |  |  |  |  |  |
CRE | 0.135 | 0.133 | 0.002 | 0.449 |  |  |  |  |  |  |  |
CI | 0.207 | 0.217 | -0.011 | 0.069∗ |  |  |  |  |  |  |  |
CC | 0.794 | 0.852 | -0.058 | 0.109 | -0.029 | -0.032 | 0.003 | 0.938 |  |  |  |
FL | 0.010 | 0.021 | -0.011 | 0.306 | 0.335 | 0.307 | 0.028 | 0.424 |  |  |  |
CRE | 0.095 | 0.068 | 0.026 | 0.052∗ | 0.127 | 0.191 | -0.064 | 0.090∗ |  |  |  |
CI | 0.178 | 0.196 | -0.018 | 0.000∗∗∗ | 0.064 | 0.054 | 0.010 | 0.514 |  |  |  |
CC | 0.595 | 0.633 | -0.039 | 0.001∗∗∗ |  |  |  |  | 0.417 | 0.397 | 0.020 | 0.000∗∗∗
FL | 0.037 | 0.052 | -0.015 | 0.001∗∗∗ |  |  |  |  | 0.110 | 0.106 | 0.003 | 0.329
CRE | 0.129 | 0.126 | 0.002 | 0.477 |  |  |  |  | 0.057 | 0.057 | 0.000 | 0.601
CI | 0.155 | 0.155 | -0.001 | 0.886 |  |  |  |  | 0.190 | 0.190 | 0.000 | 0.981
CC | 0.577 | 0.618 | -0.041 | 0.046∗∗ | 0.037 | 0.032 | 0.005 | 0.846 | 0.420 | 0.399 | 0.021 | 0.000∗∗∗
FL | 0.006 | 0.018 | -0.012 | 0.146 | 0.315 | 0.268 | 0.047 | 0.041∗∗ | 0.095 | 0.094 | 0.001 | 0.768
CRE | 0.097 | 0.078 | 0.019 | 0.182 | 0.102 | 0.147 | -0.045 | 0.265 | 0.052 | 0.050 | 0.002 | 0.267
CI | 0.140 | 0.153 | -0.012 | 0.035∗∗ | 0.032 | 0.007 | 0.025 | 0.105 | 0.189 | 0.189 | -0.001 | 0.583
∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01
Table EC.2: Comparison of pooled and FEO coefficients using pre-COVID data
with allowance rate as an additional proxy for banks’ portfolio risks.

##  Appendix EC.6 Revenue Models

The Federal Reserve’s stress testing framework includes models of revenues as
well as models of losses. We have focused on loan portfolio loss models
because they fit most clearly within the Fed’s policy of equal treatment and
its preference for industry models. In this section, we show that the
heterogeneity documented for loss models in Section 6 extends to revenue
models, referred to in the Fed’s framework as models of pre-provision net
revenue or PPNR.

The Fed uses a suite of PPNR models to forecast difference sources of revenue.
These models differ from the loss models in at least two important respects:
they are typically autoregressive (AR) models, and, unlike the portfolio loss
models, they do not rule out bank fixed effects; see [9]. This feature points
to the presence of unmodeled bank heterogeneity in the revenue forecasts. Our
goal in this section is to check for heterogeneity in a simple PPNR model and
to compare coefficient estimates in the pooled and fixed-effect models.

We consider the modeling of trading revenue, which is one of the PPNR
components in the Fed’s framework. We compare AR models with fixed-effects,

$Y_{s,t}^{FE}=\alpha_{s}+\rho_{FE}Y_{s,t-1}+\beta_{FE}X_{s,t}+\gamma_{FE}\textit{VIX}_{t}+\epsilon^{FE}_{s,t}$
(73)

or pooled without fixed effects,

$Y_{s,t}^{P}=\alpha_{P}+\rho_{P}Y_{s,t-1}+\beta_{P}X_{s,t}+\gamma_{P}\textit{VIX}_{t}+\epsilon_{s,t}^{P}.$
(74)

In both models, $Y$ is trading revenue normalized by total trading assets;
this choice of normalization is consistent with [9]. For the AR term, we use
either a one-quarter lag $Y_{s,t-1}$ or a four-quarter average lag, in which
case we use $1/4\sum_{j=1}^{4}Y_{s,t-j}$ in place of $Y_{s,t-1}$ in (73) and
(74) to capture average performance over the past year.

For $X_{s,t}$ we use the size of bank $s$ in quarter $t$, as measured by the
log of total assets. For the macro variable, we use the VIXt, the market
volatility index taken from the Federal Reserve’s stress testing historical
dataset. Market volatility is expected to have a direct impact on trading
revenue, and indeed we observe a more significant effect of VIXt than MacroPCt
(from Section 6) in this setting.

Models (73) and (74) differ in their intercept terms: (73) captures banks’
fixed effects, but (74) requires the same intercept across all banks. The
fixed-effect coefficient estimates $\rho_{FE},\beta_{FE}$ and $\gamma_{FE}$
are identical to those of FEO; the methods differ in their forecasts: FEO uses
the average fixed effect, rather than bank-specific fixed effects in its
forecasts.

As in Section 6, we use Y-9C financial reporting data for the top 35 banks by
total asset size (as of year-end 2021), and we include only banks with at
least 18 years of data to ensure our panel is mostly balanced. We fit the
models using weighted least squares, weighting each observation by quarter
stress and bank asset balance. As in (53), we choose the weights so that, for
the same asset level, the quarters with the highest market volatility get
twice the weight as the quarters with the lowest market volatility. As in the
AR models in [32], we cluster standard errors by time.

Table EC.1 reports the results. The first two columns correspond to the one-
quarter AR setting, and the last two columns are the one-year-average AR
setting. The numbers in parentheses are standard errors.

As expected, both the lagged response and the VIX term are statistically
significant. The volatility term is more significant in the one-year-average
AR setting, presumably because the market environment changes less over one
quarter, and its effect is partly captured by the lagged response.

The pooled and fixed-effect methods result in different estimates of $\rho$,
and the differences in estimates in all settings are more than 20%. We tested
the hypotheses that all banks share the same (i) intercept, (ii) AR
coefficient term, (iii) VIX coefficient, or (iv) coefficients of total asset
size, following the approach used in Section 6.2; all tests strongly reject
that banks have identical model coefficients, with $p$-values less than 0.01,
extending what we found for loss models. We have also examined the forecasts
of trading interest income and trading interest expense (two other components
of PPNR), and observed similar patterns across all three components.

| Normalized Trading Revenue
---|---
AR | 0.466∗∗∗ | 0.588∗∗∗ | 0.406∗∗∗ | 0.596∗∗∗
| (0.134) | (0.111) | (0.078) | (0.094)
VIX | -0.024∗∗ | -0.022∗ | -0.027∗∗∗ | -0.022∗∗∗
| (0.011) | (0.012) | (0.008) | (0.008)
Log Assets | -0.003 | -0.007∗∗∗ | 0.000 | -0.006∗∗∗
| (0.004) | (0.001) | (0.003) | (0.002)
Bank FE | included | – | included | –
|  ∗p$<$0.1; ∗∗p$<$0.05; ∗∗∗p$<$0.01
Table EC.1: Coefficient estimates for the AR models. The first two columns use
a one-quarter lag, and the last two use a one-year average lag. Columns 1 and
3 correspond to AR models with bank fixed effects, and columns 2 and 4 are
pooled AR models.