# An $\ell_{0}\ell_{2}$-norm regularized regression model for construction of
robust cluster expansions in multicomponent systems

Peichen Zhong  Department of Materials Science and Engineering, University of
California, Berkeley, California 94720, United States  Materials Sciences
Division, Lawrence Berkeley National Laboratory, California 94720, United
States  Tina Chen  Department of Materials Science and Engineering, University
of California, Berkeley, California 94720, United States  Materials Sciences
Division, Lawrence Berkeley National Laboratory, California 94720, United
States  Luis Barroso-Luque  Department of Materials Science and Engineering,
University of California, Berkeley, California 94720, United States  Materials
Sciences Division, Lawrence Berkeley National Laboratory, California 94720,
United States  Fengyu Xie  Department of Materials Science and Engineering,
University of California, Berkeley, California 94720, United States  Materials
Sciences Division, Lawrence Berkeley National Laboratory, California 94720,
United States  Gerbrand Ceder<EMAIL_ADDRESS>Department of Materials
Science and Engineering, University of California, Berkeley, California 94720,
United States  Materials Sciences Division, Lawrence Berkeley National
Laboratory, California 94720, United States

###### Abstract

We introduce $\ell_{0}\ell_{2}$-norm regularization and hierarchy constraints
into linear regression for the construction of cluster expansions to describe
configurational disorder in materials. The approach is implemented through
mixed integer quadratic programming (MIQP). The $\ell_{2}$-norm regularization
is used to suppress intrinsic data noise, while the $\ell_{0}$-norm is used to
penalize the number of non-zero elements in the solution. The hierarchy
relation between clusters imposes relevant physics and is naturally included
by the MIQP paradigm. As such, sparseness and cluster hierarchy can be well
optimized to obtain a robust, converged set of effective cluster interactions
with improved physical meaning. We demonstrate the effectiveness of
$\ell_{0}\ell_{2}$-norm regularization in two high-component disordered
rocksalt cathode material systems, where we compare the cross-validation,
convergence speed, and the reproduction of phase diagrams, voltage profiles,
and Li-occupancy energies with those of the conventional $\ell_{1}$-norm
regularized cluster expansion models.

###### pacs:

##  I Introduction

First-principles density functional theory (DFT) calculations have been
demonstrated as a reliable tool in computational materials science. Despite
the increase in computing power and accuracy of DFT methods, the scaling with
the number of atoms $\sim O(n^{3})$ intrinsically prohibits large-scale
calculations (over $10^{3}$ atoms) or sampling of a high-dimensional occupancy
space (millions of structures) Jain  _et al._ (2013). This is particularly
relevant in systems with configurational degrees of freedom that need to be
sampled at non-zero temperature to equilibrate states of partial order, and
their associated entropy and free energy. The state of configurational order
determines many materials properties, especially in systems composed of many
species (high-entropy systems). It has also recently been shown to be relevant
for mechanical properties in metallic alloy systems Zhang  _et al._ (2020);
Yin  _et al._ (2021) and the energy density of complex electrode materials for
energy storage application Ji  _et al._ (2019); Ouyang  _et al._ (2020); Lun
_et al._ (2021); Zhong  _et al._ (2020).

The cluster expansion (CE) method has been well developed to describe such
configurational energetics for metallic alloys Sanchez  _et al._ (1984);
Sanchez (2019), as well as for ionic systems Tepesch  _et al._ (1995); Seko
_et al._ (2009). The CE method expands any property (e.g., formation energy,
volume) in terms of the distribution of atoms on a set of predefined sites.
When the quantity being expanded is the energy, the expansion coefficients are
referred to as Effective Cluster Interactions (ECI). For example, in a
multicomponent system, the energy is expanded as

$E(\boldsymbol{\boldsymbol{\sigma}})=\sum_{\boldsymbol{\beta}}m_{\boldsymbol{\beta}}J_{\boldsymbol{\beta}}\left\langle\Phi_{\boldsymbol{\alpha}\in\boldsymbol{\beta}}\right\rangle_{\boldsymbol{\beta}},~{}\Phi_{\boldsymbol{\alpha}}=\prod_{i=1}^{N}\phi_{\alpha_{i}}(\sigma_{i}).$
(1)

A configuration $\boldsymbol{\sigma}$ represents a specific occupancy on all
the sites of the system, where $\sigma_{i}$ describes which species sits on
the $i$-th site of the structure. The site basis function
$\phi_{\alpha_{i}}(\sigma_{i})$ transforms the occupancy variable $\sigma_{i}$
into a scalar value. There are typically as many (non-constant) cluster basis
functions as possible occupancies on a site minus one. The cluster basis
function label $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3},\ldots)$
indicates a group of sites, each with a specific basis function on it, where
each entry $\alpha_{i}$ labels the corresponding site basis function
$\phi_{\alpha_{i}}$. Thus, the cluster basis function
$\Phi_{\boldsymbol{\alpha}}=\prod_{i=1}^{N}\phi_{\alpha_{i}}(\sigma_{i})$ can
be obtained by taking the product of site basis functions.

For example, the cation sublattice of a LiMnO2 rocksalt oxide is a binary
system, where Li and Mn share the octahedral interstitial of the FCC anion
framework. In such a system, Li can be encoded by $\sigma^{\mathrm{Li}}=0$ and
Mn by $\sigma^{\mathrm{Mn}}=1$. The parameter $\alpha_{i}$ takes a value from
$[0,1,...M-1]$, where $M$ is the number of allowed species defined on the
sublattice (e.g., $M=2$ for Li-Mn). While many forms of site basis function
can be usedSanchez  _et al._ (1984); Zhang and Sluiter (2016); Barroso-Luque
_et al._ (2021), a sinusoid (orthogonal) basis function is applied here to
transform the occupancy variable ($\sigma^{\mathrm{Li}},\sigma^{\mathrm{Mn}}$)
into a value van de Walle (2009), where

$\begin{split}\phi_{j}(\sigma_{i})=\begin{cases}1&\text{if }j=0\\\
-\cos\left(\frac{\pi(j+1)\sigma_{i}}{M}\right)&\text{if }j\text{ is odd}\\\
-\sin\left(\frac{\pi j\sigma_{i}}{M}\right)&\text{if }j\text{ is
even}\end{cases}.\end{split}$ (2)

The $j$ indicates $\alpha_{i}$ in Eq.(1) and can take a value of 0 or 1. Thus,
we have $\phi_{j=0}\equiv 1$, $\phi_{j=1}(\sigma^{\mathrm{Li}}=0)=-1$, and
$\phi_{j=1}(\sigma^{\mathrm{Mn}}=1)=1$. This situation corresponds to the spin
variables used in a generalized Ising model Ceder (1993); Huang  _et al._
(2016). For systems with species number $M>2$, the basis functions take values
beyond those of spin variables $\left\\{-1,1\right\\}$ typically used in
binary CE. Some examples of other types of site-basis functions also developed
for the CE method are the Chebyshev polynomials Sanchez  _et al._ (1984) and
the indicator function (point delta function) Zhang and Sluiter (2016).

Figure 1: (a) The general flowchart of constructing a CE model, including
initialization of input structures, DFT calculations, fitting and convergence
check, and cluster expansion Monte Carlo (CEMC) for sampling. (b) An
illustration of $\ell_{2}$, $\ell_{1}$ and $\ell_{0}$-norm regularization in a
two-parameter space $\boldsymbol{J}$ = ($J_{1}$, $J_{2}$) . The blue circles
represent the contours of the data term
$||\boldsymbol{E}_{\mathrm{DFT},S}-\boldsymbol{\Pi}_{S}\boldsymbol{J}||^{2}_{2}$
in cost function. The red regions represent the constraints of parameters
(e.g. $J_{1}^{2}+J_{2}^{2}\leq s$ for $\ell_{2}$-norm, $|J_{1}|+|J_{2}|\leq s$
for $\ell_{1}$-norm.) The dark red point is the intersection of data term and
regularization of parameters, which jointly determines the estimation of
$\boldsymbol{J}$.

In Eq.(1), the correlation function
$\left\langle\Phi_{\boldsymbol{\alpha}}\right\rangle_{\boldsymbol{\beta}}$ is
calculated by

$\langle\Phi_{\boldsymbol{\alpha}}(\boldsymbol{\sigma})\rangle_{\boldsymbol{\beta}}=\frac{1}{N_{\boldsymbol{\sigma}}m_{\boldsymbol{\beta}}}\sum_{\boldsymbol{\alpha}\in\boldsymbol{\beta}}\Phi_{\boldsymbol{\alpha}}(\boldsymbol{\sigma}),$
(3)

where $\boldsymbol{\beta}$ is an orbit representing all symmetrically
equivalent cluster basis functions $\boldsymbol{\alpha}$, and
$m_{\boldsymbol{\beta}}$ is the corresponding multiplicity.
$N_{\boldsymbol{\sigma}}$ is the size of the supercell of configuration
$\boldsymbol{\sigma}$; thus, the correlation function is well normalized with
respect to the primitive cell. $J_{\boldsymbol{\alpha}}$ is the effective
cluster interaction (ECI). We refer readers to Ref.De Fontaine (1979); Sanchez
_et al._ (1984); van de Walle (2009); Ångqvist  _et al._ (2019) for a more
extended description of the CE. From Eq.(1), the CE energy is linearly
dependent on the ECIs $\boldsymbol{J}$ when the configuration
$\boldsymbol{\sigma}$ is fixed. We can thus write

$E_{\mathrm{CE}}(\boldsymbol{\boldsymbol{\sigma}})=\boldsymbol{\Pi}(\boldsymbol{\sigma})\cdot\boldsymbol{J},$
(4)

where $\boldsymbol{\Pi}(\boldsymbol{\sigma})$ is a row vector of correlation
functions and $\boldsymbol{J}$ is the column vector of ECIs.

Fig.1(a) presents a brief illustration of how to iteratively construct a CE
Hamiltonian. In practice, the CE model is initially fitted on a small set of
DFT calculations. Then, a simple CE is fitted that can be used in a Monte
Carlo simulation to sample new structures. DFT calculations will be applied to
a sample of the MC-obtained structures, and a new CE will be fitted. This
procedure will be performed iteratively until the model is converged (i.e.,
the cross-validation error remains low and stable, the model reproduces DFT
ground states well, etc.) Huang  _et al._ (2017). In such a process, it is
always desirable to achieve fewer training iterations, as DFT calculations are
costly in terms of CPU time. On the other hand, fewer structures may also
result in a worse fitting due to insufficient sampling of the configuration
space.

Obtaining reliable ECIs $\boldsymbol{J}$ from the DFT energy of a set of
configurations is the central problem of CE fitting. Given a set of input
occupancy configurations $S$, the set of correlation vectors forms a feature
matrix $\boldsymbol{\Pi}_{S}=[\boldsymbol{\Pi}_{1},\boldsymbol{\Pi}_{2},...]$,
and the corresponding DFT energies are used to construct the target vector
$\boldsymbol{E}_{\mathrm{DFT},S}$. Determining the ECIs is an inverse problem
of Eq.(4), also called linear regression. Generally, the problem can be solved
by minimizing the cost function

$\min_{\boldsymbol{J}}~{}||\boldsymbol{E}_{\mathrm{DFT},S}-\boldsymbol{\Pi}_{S}\boldsymbol{J}||^{2}_{2}+\mu||\boldsymbol{J}||_{p}~{},~{}||\boldsymbol{x}||_{p}=\left(\sum_{i}|x_{i}|^{p}\right)^{\frac{1}{p}},$
(5)

where the $p$-norm of $\boldsymbol{J}$ is added to regularize the fit and
suppress over-fitting, and $\mu$ controls the degree of regularization.
Fig.1(b) shows the comparison of $\ell_{2}$, $\ell_{1}$ and $\ell_{0}$-norm
regularization in a two-parameter space $\boldsymbol{J}$ = ($J_{1}$, $J_{2}$).
The blue circles are the contours of the data term error
$||\boldsymbol{E}_{\mathrm{DFT},S}-\boldsymbol{\Pi}_{S}\boldsymbol{J}||^{2}_{2}$.
The red regions represent the regularization constraints on the parameters
($||\boldsymbol{J}||_{p}\leq s$), which can be transformed to a Lagrangian
form $\mu||\boldsymbol{J}||_{p}$ as shown in Eq.(5). The dark red point is the
regularized estimation of $\boldsymbol{J}$, which is the intersection between
the data error term and the regularization term. The $\ell_{1}$-norm tends to
generate sparser solutions compared with the $\ell_{2}$-norm, because the
intersection is likely to be located on the axis. The $\ell_{0}$-norm counts
the non-zero elements of $\boldsymbol{J}$, where the intersection is exactly
located on the axis and thus the $\ell_{0}$-norm imposes an exact sparsity
constraint on $\boldsymbol{J}$.

Conventionally, $\ell_{2}$-norm (ridge regression, $p=2$) regularization can
be applied when the problem is over-determined (i.e., the number of training
structures is larger than the dimension of $\boldsymbol{J}$). The
$\ell_{2}$-norm regularized regression reduces the over-fitting caused by
intrinsic noise in the training data. This can be achieved solely by
introducing the $\ell_{2}$ regularization function and additionally using the
mixed-basis expansion Laks  _et al._ (1992); Blum and Zunger (2004); Mueller
and Ceder (2009). Bayesian approaches have also been successfully applied to
estimate the ECIs with a prior distribution in several binary systems Nelson
_et al._ (2013a); Cockayne and van de Walle (2010); Mueller and Ceder (2009).
However, the number of ECIs increases combinatorially with the number of
species, scaling approximately as $\prod_{k}(M_{k}-1)^{n_{k}}$, where $M_{k}$
is the number of species on the $k$-th sublattice, and $n_{k}$ is the number
of cluster sites in the same sublattice $k$. The explosion in the number of
basis functions when many species can occupy a site makes it difficult to
predefine which cluster basis functions contribute to the expansion for high
dimensional multi-component systems (i.e., which cluster basis function has a
non-zero element in the solution of $\boldsymbol{J}$). Therefore, a sparse
solver for ECIs selection is required.

Rigorously, the exact sparse solution of Eq.(5) is obtained with
$\ell_{0}$-norm regularization of $\boldsymbol{J}$. However, it is hard to
compute the $||\boldsymbol{J}||_{0}$ in the cost function as it is an NP-hard
problem. In the compressive sensing paradigm, the $\ell_{0}$-norm can be
transformed to an $\ell_{1}$-norm when the feature matrix
$\boldsymbol{\Pi}_{S}$ satisfies the restricted isometry property (RIP)
condition Candes and Wakin (2008). To satisfy the RIP condition in CE, Nelson
_et al._ (2013b) proposed to generate a training set in which each row is an
identical independent distributed (i.i.d.) random vector. However, in
practical cases, the configurations in the training set $S$ are correlated,
because structures are not randomly sampled, but are mostly part of an
ensemble of configurations with low energy. Such correlations fail to satisfy
the i.i.d. condition. Moreover, generating structures from a specific
correlation vector is also an NP-hard problem. Though the strict compressive
sensing cluster expansion is not easy to construct in practice, the
$\ell_{1}$-norm (lasso, $p=1$) regularization is widely used as feature
selection, which has shown success in various alloy and ionic systems Nelson
_et al._ (2013b); Seko  _et al._ (2014); Leong and Tan (2019); Barroso-Luque
_et al._ (2021).

In this paper, we propose an $\ell_{0}\ell_{2}$-norm regularization approach
that incorporates hierarchy constraints to generate more robust and predictive
CE models. First, we introduce the $\ell_{0}\ell_{2}$ penalty term and
hierarchy constraints in the paradigm of mixed-integer quadratic programming
(MIQP). Second, we compare the sparseness and convergence rate of ECIs with
those of the conventional $\ell_{1}$ method in the Li-Mn-V-Ti-O-F disordered
rocksalt system. Finally, we demonstrate that an
$\ell_{0}\ell_{2}$-regularized CE better reproduces the correct physical
interactions by comparing with $\ell_{1}$-CE in terms of computed phase
diagrams, voltage profiles, and related physical quantities in the Li-Mn-Ti-O
system.

##  II Methods

###  II.1 The $\ell_{0}$-norm regularization

In Eq.(5), $p=0$ manifests itself as a pseudo-norm that counts the number of
non-zero elements of $\boldsymbol{J}$:

$||\boldsymbol{J}||_{0}=\sum_{i}\mathrm{Ind}(J_{i}),~{}\mathrm{Ind}(J_{i})=\left\\{\begin{aligned}
&0,J_{i}=0\\\ &1,J_{i}\neq 0\end{aligned}\right.$ (6)

Adding the $\ell_{0}$ term into the cost function directly penalizes the
number of non-zero ECIs, yielding better sparseness in its solution. However,
optimizing a cost function with an $\ell_{0}$ term is an NP-hard problem and
is difficult to present in a direct way Candes and Wakin (2008); Elad (2010).
Previously, Huang  _et al._ has approached the problem by rewriting $\ell_{0}$
optimization as a mixed-integer programming problem, such that

$\displaystyle\min||\boldsymbol{J}||_{0}\quad\Leftrightarrow\quad$
$\displaystyle\min~{}\sum_{c\in\boldsymbol{C}}z_{0,c}$ (7)
$\displaystyle\begin{array}[]{r@{~}l@{}l@{\quad}l}\boldsymbol{s.t.}~{}&Mz_{0,c}\geq
J_{c},~{}\forall c\in\boldsymbol{C}\\\ ~{}&Mz_{0,c}\geq-J_{c},~{}\forall
c\in\boldsymbol{C}\\\ ~{}&z_{0,c}\in\\{0,1\\},~{}\forall
c\in\boldsymbol{C}\end{array}$ (11)

where $M$ is a sufficiently large number (larger than the maximum possible
absolute value of any ECI), and $z_{0,c}$ is a slack variable (binary integer)
indicating whether the ECI of orbit $c$ is zero or not. $J_{c}$ is constrained
to $0$ when the slack variable $z_{0,c}=0$ (inactive) and to $[-M,M]$ when
$z_{0,c}=1$ (active). (For a rigorous mathematical background, refer to
Ref.Huang  _et al._ .) In practice, it is shown that one can at least obtain a
sparseness-improved near-optimal solution within a reasonable CPU time cutoff.
In our benchmark tests of the Li-Mn-V-Ti-O-F and Li-Mn-Ti-O systems, the
optimizations of ECIs were completed within 600s using the gurobi package
Gurobi Optimization, LLC (2021).

###  II.2 Hierarchy constraints

In a CE, clusters are usually enumerated in an iterative, low-to-high order
(i.e., from singlets to pairs, triplets, quadruplets, and so on). Practically,
the CE is truncated to a maximum of $n$ (e.g., quadruplet clusters with $n=4$
are a typical limit), ignoring the higher-order interactions to control the
model complexity. To differentiate the cluster orbits by different
significance, we take one of the basic assumptions of CE that $n$-body cluster
interactions become less important to the configurational energy (or other
scalar properties) as $n$ becomes larger. This assumption means that the
majority of the fitted property can be described by the lower-order
interactions and that the higher-order interactions serve as the fine-tuning
part in the fitting.

Such a physically inspired concept can be introduced in the form of hierarchy
constraints, as has been done successfully in some previous studies Sluiter
and Kawazoe (2005); Zarkevich and Johnson (2004); Walle and Ceder (2002). The
hierarchy constraint manifests itself as $J_{b}\neq 0$ if and only if
$J_{a}\neq 0$ ($a\subset b$), where $a$ and $b$ are a lower- and higher-order
cluster function orbit, respectively, and $b$ contains all the site bases of
$a$ as a subset. In the MIQP representation, the hierarchy relationship can be
easily expressed as a constraint between slack variables:

$z_{0,b}\leq z_{0,a},~{}a\subset b.$ (12)

This treatment was first proposed by Huang  _et al._ , where it was used in
the $\ell_{0}\ell_{1}$-norm regularization paradigm.

###  II.3 The $\ell_{2}$-norm regularization

We propose that combining $\ell_{2}$-norm and $\ell_{0}$-norm regularization
can impose true hierarchy constraints unlike the $\ell_{0}\ell_{1}$-norm. It
is to be noted that the inequality between slack variables does not
necessarily impose the hierarchy relation ($J_{b}\neq 0$, iff $J_{a}\neq 0$).
This is because the hierarchy constraints are defined on the magnitude of ECIs
$J_{a}$ and $J_{b}$, while the slack variables $z_{0,b},z_{0,a}$ are
intermediate to represent the presence or exclusion of the variables.

Figure 2: Illustration of hierarchy relations
($\alpha\subset\beta\subset\gamma$) between pair, triplet, and quadruplet
orbit. The different colors on the cluster sites represent the decorating
species for a given site-basis function. The equation in red shows a pseudo-
active hierarchy constraint that may appear in $\ell_{1}$ and its derivative
methods.

When implementing the hierarchy constraints in $\ell_{0}\ell_{1}$-norm
regularization, pseudo-active behavior can manifest itself when a $J=0$, but
its slack variable $z_{0}=1$ within the MIQP paradigm. $J$ can be regularized
to zero, which is still a valid solution between $[-M,M]$, even with
$z_{0}=1$. This is caused by the fact that the $\ell_{1}$-norm has feature-
selection properties that intrinsically produce a sparse solution Tibshirani
(1996). This pseudo-activeness can introduce excessive sparseness to the
solution and break the hierarchy constraints. Fig.2 presents an example of
pseudo-activeness in $\ell_{0}\ell_{1}$-norm regularization. The excessive
sparseness is introduced to the orbit $\beta$ with $J_{\beta}=0$, while all
orbits $\alpha,\beta,\gamma$ has active slack variables $z_{0}=1$. The higher-
order orbit $\gamma$ is erroneously activated while $J_{\beta}=0$. To avoid
such a situation and ensure proper function with $\ell_{0}$ under hierarchy
constraints, a norm with no feature-selection properties is required. The
$\ell_{2}$-norm is a natural choice.

With the introduction of the $\ell_{0}\ell_{2}$-norm and hierarchy
constraints, the final ECI optimization problem can be written as

$\displaystyle\min_{\boldsymbol{J}}~{}\boldsymbol{J}^{T}\boldsymbol{\Pi}_{S}^{T}\boldsymbol{\Pi}_{S}\boldsymbol{J}^{T}-2\boldsymbol{E}_{\mathrm{DFT}}^{T}\boldsymbol{\Pi}_{S}\boldsymbol{J}+\mu_{0}\sum_{c\in\boldsymbol{C}}z_{0,c}+\mu_{2}||\boldsymbol{J}||_{2}^{2}$
(13) $\displaystyle\begin{array}[]{r@{~}l<EMAIL_ADDRESS>Mz_{0,c}~{}&\geq J_{c},~{}\forall c\in\boldsymbol{C}\\\ Mz_{0,c}~{}&\geq-
J_{c},~{}\forall c\in\boldsymbol{C}\\\ z_{0,b}~{}&\leq z_{0,a},~{}\forall
a\subset b,~{}\\{a,b\\}\in\boldsymbol{C}\\\ z_{0,c}~{}&\in\\{0,1\\},~{}\forall
c\in\boldsymbol{C},\end{array}$ (18)

where $||\boldsymbol{J}||_{2}^{2}=\boldsymbol{J}^{T}\boldsymbol{J}$ penalizes
the magnitude of ECIs, thus avoiding over-fitting by regularizing sampling
noise while the $\ell_{0}$-term $\sum_{c}z_{0,c}$ optimizes the sparseness.
The hierarchy constraints ensure correct containment relationships by
manifesting lower-order ECIs first to reduce redundancy. The packages for our
implementation are available in Ref.111https://github.com/CederGroupHub/smol,
https://github.com/CederGroupHub/sparse-lm.

##  III Results

Figure 3: (a) An illustration of the rocksalt lattice structure. The cation
sites are labeled in red and can be occupied by Li+ and transition metals (TM,
including Mn2+, V3+, Ti4+ in our example) in DRX. The anion sites are labeled
in gray and can be occupied by O2- and F-. The lower panel gives some examples
of $n$-body ($n=2,3,4$) clusters used in the CE model, including intra and
inter-sublattice interactions. (b) The procedure to obtain an
$\ell_{0}\ell_{2}$-norm regularized solution, including finding the $\mu_{2}$
by minimizing the CV error in ridge regression, sparseness engineering with
$\ell_{0}$ using MIQP, and terminating if the solution is converged with good
sparseness, as well as good-reproduction-relevant physical properties.

As mentioned above, for systems with many species, the number of basis
functions grows rapidly. An example of such systems are the Li-excess
disordered rocksalts (DRX), which are multicomponent systems that can be
synthesized with a wide variety of elements Clément  _et al._ (2020).
Recently, high-entropy DRX materials have been synthesized with up to 12 metal
species Lun  _et al._ (2021). In addition, their configurational short-range
order is critical to their transport properties, warranting a detailed CE
approach.Huang  _et al._ (2021); Ouyang  _et al._ (2020) Here, we provide a
heuristic solution to study configurations in such high-dimensional DRX
systems by applying an $\ell_{0}\ell_{2}$-regularized CE model to fit the
formation energy of DRX compounds.

###  III.1 Robustness and convergence

The convergence of the CE when the $\ell_{0}\ell_{2}$-norm and hierarchy
constraints are enforced was tested on configurational disorder in the
LiF–MnO–LiVO2–Li2TiO3 composition space. The CE model contains pair
interactions up to 7.1 Å, triplets up to 4.0 Å, and quadruplets up to 4.0 Å
based on a lattice parameter $a=3$ Å for the primitive cell. Fig.3(a) presents
the rocksalt framework of a DRX structure. The framework contains a cation
sublattice (red) and anion sublattice (gray), where the cation sites can be
occupied by Li and transition metals (TM, including Mn2+, V3+, Ti4+ in this
example) and the anion sites can be occupied by O2- and F-. A species
indicator where the site basis function reads
$\phi_{j}(\sigma_{i})=\delta_{i,j}$ was used Zhang and Sluiter (2016). The
electrostatic energy (Ewald energy $E_{0}/\varepsilon_{r}$) is also included
to capture long-range electrostatic interactions ($E_{0}$ is the unscreened
electrostatic energy and $1/\varepsilon_{r}$ is fitted as one of the ECI
($1/\varepsilon_{r}\geq 0$) Seko and Tanaka (2014); Richards  _et al._ (2018).
In total, 162 ECIs (including the constant term $J_{0}$) are predefined in the
CE Hamiltonian. The dimension of the feature matrix
$\boldsymbol{\Pi}_{\text{DFT},S}$ is $487\times 162$. The performance of the
$\ell_{0}\ell_{2}$-CE is compared with the $\ell_{1}$-CE. We emphasize two
major improvements in the $\ell_{0}\ell_{2}$-CE.

Figure 4: (a) Cross-validation error (meV/atom) of the $\ell_{1}$-CE and the
$\ell_{0}\ell_{2}$-CE. The sparseness is the number of non-zero ECIs in the
fit ($||\boldsymbol{J}||_{0}$). The curves are generated by varying
hyperparameters $\mu_{0},\mu_{1}$, and $\mu_{2}$ in regularization. (b) ECIs
convergence test vs training set size. $\boldsymbol{J}$ is the ECIs fitted
with full training data, and $\boldsymbol{J}_{sub}$ is the ECIs fitted with a
subset of corresponding size.

####  III.1.1 (1) Sparseness vs. cross-validation error

Cross-validation (CV) error vs. model complexity is a general metric used to
evaluate the robustness of a CE model. The optimal trade-off between under-
fitting and over-fitting can be found with a CV test, where the optimal model
is fitted with the regularization hyper-parameter $\mu$ that minimizes the CV
error. In our test, a $k$-fold CV error is used,

$\text{CV}=\sqrt{\frac{1}{k}\sum_{j=1}^{k}\text{MSE}_{j}},~{}\text{MSE}=\frac{1}{N}\sum_{i=1}^{N}(E_{\text{DFT}}^{i}-E_{\text{CE}}^{i})^{2},$
(19)

where CV is the cross-validation error averaged over $k$ splits of the
validation dataset, and MSE is the mean-squared-error of each validation
dataset. Here, $N$ is the size of the validation dataset, and $k=5$ is the
number of folds. In our tests, the regularization hyperparameter $\mu$ is
selected from the logarithm space between $[10^{-6},10^{-1}]$. The sparseness
is defined as the number of non-zero elements of the solution
($||\boldsymbol{J}||_{0}$) and represents the model complexity.

The CV error versus sparseness is presented in Fig.4(a) for an $\ell_{1}$ and
$\ell_{0}\ell_{2}$-norm regularized CE. For the $\ell_{0}\ell_{2}$-CE, the CV
error remains low as the sparseness varies between 25 and 150 ECIs. In this
regime, the $\ell_{0}\ell_{2}$-CE shows no sign of over-fitting as the CV
error remains near the global minimum around 6 meV/atom. The $\ell_{1}$-CE
shows a similar optimal CV error as that of $\ell_{0}\ell_{2}$-CE near this
_minimum plateau_ regime from 50 to 100 in sparseness. However, as the model
complexity changes, the CV error increases at both low and high sparsity,
indicating that the $\ell_{1}$-CE is less robust against the choice of model
complexity. Therefore, we conclude that the $\ell_{0}\ell_{2}$-CE can reach
low CV error with a lower complexity, which is empirically believed to result
in models that better reproduce physics. A more sparse CE can increase the
computational speed of energy evaluations and is also less sensitive to model
complexity change as compared with the $\ell_{1}$-CE.

####  III.1.2 (2) Convergence of ECIs with a subset of training data

The second point that we want to emphasize is that the $\ell_{0}\ell_{2}$-CE
converges to its most accurate solution faster than the $\ell_{1}$-CE, which
lowers the risk of obtaining an over-fitted result when the configuration
space is under-sampled. This is an important improvement in the practical use
of CE constructions. To test this hypothesis and mimic the iterative sampling
process, we designed a numerical experiment based on a finished DFT dataset
(with 487 structures in total). Then, we evaluated the quality of fits
performed on subsets of training data of increasing size. We subsequently
compared the subset-fitted ECIs $\boldsymbol{J}_{sub}$ with the full-set
result. In such a comparison, the ground-truth (full set) solution is set as
follows: (1) For the $\ell_{1}$-CE, the regularization parameter $\mu_{1}$ is
chosen at the minimum CV error according to Fig.4(a). This solution has 99
non-zero ECIs when all 487 training structures are used in the fitting. (2)
For the $\ell_{0}\ell_{2}$-CE, to compare the convergence rate under a similar
degree of model complexity, hyperparameters are selected such that the
$\ell_{1}$-CE and $\ell_{0}\ell_{2}$-CE have similar sparsity. The resulting
$\ell_{0}\ell_{2}$-CE has 92 non-zero ECIs with all 487 training structures
included according to Fig.4(a).

After setting the hyperparameters for both models, we compared the normalized
absolute difference
$||\boldsymbol{J}_{sub}-\boldsymbol{J}||_{1}/||\boldsymbol{J}||_{1}$ between
the $\ell_{1}$-CE (blue line) and $\ell_{0}\ell_{2}$-CE (red line) in
Fig.4(b). For each subset size, ten randomly selected subsets with the same
size were evaluated and averaged. The solid square represents the average, and
the error bar represents the standard deviation resulting from different
subsets. Fig.4(b) indicates that the $\ell_{1}$-CE demonstrates higher
deviation from the ground-truth solution and converges more slowly to it than
the $\ell_{0}\ell_{2}$-CE as the training set is increased. This result
unambiguously demonstrates the robustness of $\ell_{0}\ell_{2}$-CE to work
with small input data sets.

Figure 5: (a) Phase diagram generated with DFT, $\ell_{0}\ell_{2}$-CE, and
$\ell_{1}$-CE. The DFT ground states are labeled in blue text. The incorrectly
predicted ground states are labeled with red circles and text. (b) The
simplified spinel voltage profile (blue line) generated by $\ell_{1}$-CE and
$\ell_{0}\ell_{2}$-CE for spinel orderings in LixMn2O4 is compared with the
DFT ground-truths (orange line). (c) Energy difference of Li occupation in
octahedral and tetrahedral sites in layered MnO2 (upper) and spinel MnO2
framework (lower).

###  III.2 ECIs with improved physics

From a general perspective of machine learning (ML), the predictions of
energies are made by fitting statistical models on a group of data points. The
statistical models can predict the absolute energy with high accuracy by
minimizing the cost function, which is constructed by the difference between
prediction and observation of the training data. However, in materials
science, relative energy quantities are of greater significance than the
absolute one (such as energy above the hull, phase diagram and the derivatives
of formation energy with respect to the compositional variables). Bartel  _et
al._ (2020) critically examined several ML models for energetics prediction,
and found that while the models predict the formation energy ($\Delta H_{f}$)
of materials well, they failed to predict the relative phase stability. Such a
dilemma indicates that the prediction error (CV or RMSE) is not the only thing
one should consider when constructing a statistical model for the energy.

To demonstrate that the $\ell_{0}\ell_{2}$-CE also leads to a more physically
informed solution, we studied a multicomponent system: Li-Mn-Ti-O oxide in an
fcc rocksalt framework, with Li+-Mn2+-Mn3+-Mn4+-Ti4+-vacancy disorder on the
octahedral cation sites and Li+-Mn2+-Mn3+-vacancy disorder on the interstitial
tetrahedral sites. The Li-Mn-Ti-O composition space contains a number of
battery-relevant systems Ji  _et al._ (2019); Ouyang  _et al._ (2020). These
battery systems are charged and discharged by adding or removing lithium
(i.e., lithiation or delithiation) and a charge-compensating electron, which
reduces or oxidizes a transition metal. As a result, an important physical
property to correctly model in the Li-Mn-Ti-O system is the energetics of Li
in octahedral vs. tetrahedral sites. One significant battery-relevant system
in which the effects of Li local environment preference are especially
presented is the LiMn2O4 spinel. When fully lithiated to Li2Mn2O4, Li-ions
occupy octahedral sites while the Li-ions occupy tetrahedral sites for
compositions LixMn2O4 when $x\leq 1$.

Thus, we design two additional tests to ensure that the CE models well
represent the physics of the Li octahedral vs. tetrahedral site preferences.
Specifically, we compare how well the CE model reproduces: (1) energy
differences between the Li in the tetrahedral and octahedral sites in layered
MnO2 and spinel MnO2 frameworks and (2) a simplified spinel voltage profile
against the DFT ground truths. The simplified spinel voltage profile includes
the fully lithiated rocksalt-like Li2Mn2O4, the spinel LiMn2O4, the commonly
seen Li0.5Mn2O4 ordering, and the fully delithiated Mn2O4 and is calculated by
taking the average voltage between each set of adjacent orderings. The average
voltage is calculated using DFT and the following equation Urban  _et al._
(2016); Seo  _et al._ (2015); Aydinol and Ceder (1997):

$\bar{V}(x_{1},x_{2})\approx-\frac{E_{\mathrm{Li}_{x_{1}}\mathrm{Mn}_{2}\mathrm{O}_{4}}-E_{\mathrm{Li}_{x_{2}}\mathrm{Mn}_{2}\mathrm{O}_{4}}-(x_{1}-x_{2})E_{\mathrm{Li}}}{F(x_{1}-x_{2})},$
(20)

where $x_{1}$ and $x_{2}$ are adjacent Li contents with $x_{1}>x_{2}$,
$E_{\mathrm{Li}}$ is the DFT energy of bcc Li metal, and $F$ is the Faraday
constant.

The CE was generated with pair interactions up to 7.1 Å, triplet interactions
up to 4.0 Å, and quadruplet interactions up to 3.0 Å based on a primitive cell
of the rocksalt structure with lattice parameter $a=3$ Å. A sinusoid site
basis was used (Eq.(2)). In total, 1475 ECIs (including the constant term
$J_{0}$) were predefined in the CE Hamiltonian. The dimension of the feature
matrix is $1137\times 1475$. Because of the high compositional dimensionality,
the possible number of ECIs within the interaction cutoffs is large. In
addition, there are some constraints on the occupancies in the Li-Mn-Ti-O
system, such as (1) the total number of Li, transition metals, and vacancies
is fixed between octahedral and tetrahedral cation sublattice; (2) the net
charge of the system must be neutral, etc. These relations and the inability
to sample all possible configurations with DFT reduce the rank of the feature
matrix below the dimension ($\text{rank}(\boldsymbol{\Pi}_{S})=557$), which
indicates that a sparse solution is required.

From the test results in Fig.4, we notice that when the sparseness varies, the
variation of the CV error is smaller when the CE is regularized with the
$\ell_{0}\ell_{2}$-norm than with the $\ell_{1}$-norm. This result indicates
that $\ell_{0}\ell_{2}$ has a hyperparameter space that is larger and more
tunable, whereas the $\ell_{1}$-CE is more deterministic with a small range of
optimal $\mu_{1}$ obtained by minimizing the CV error. Motivated by this
observation, the selection of ECIs for the Li–Mn–Ti–O system was completed as
follows.

The regularization strength $\mu_{1}$ in the $\ell_{1}$-CE was selected from
the stable plateau region when minimizing the CV error in lasso (e.g., the
$\mu_{1}$ associated with points between sparseness of 50 to 100 in Fig.4).
For the $\ell_{0}\ell_{2}$-norm, the $\mu_{2}$ was selected from the stable
plateau region by minimizing the CV error in ridge regression, similar to what
is done for $\ell_{1}$-CE. After obtaining the optimal $\mu_{2}$, the solution
for $\ell_{0}\ell_{2}$-CE was further determined by searching $\mu_{0}$ for a
solution with the proper sparseness (at least $||\boldsymbol{J}||_{0}<$
rank($\boldsymbol{\Pi}_{S}$), $\mu_{1},\mu_{2},\mu_{0}\in[10^{-6},10^{-1}]$).
For both $\ell_{1}$-CE and $\ell_{0}\ell_{2}$-CE, several models with low CV
error were tested for their ability to well reproduce physical properties,
such as minimal violation of DFT ground states in the phase diagram, voltage
profile comparison against DFT, as well as the Li-site energy difference
between tetrahedral and octahedral occupancy. The best performing models for
both $\ell_{1}$ and $\ell_{0}\ell_{2}$ are presented in Fig.5, respectively.

Fig.5(a) presents a comparison of ground-state phase diagrams with the
$\ell_{1}$-CE predictions, $\ell_{0}\ell_{2}$-CE predictions, and DFT
calculations. The phase diagrams were generated with in-sample training data
(all 1137 structures evaluated with DFT) for both DFT and CE models. We take
the DFT phase diagram as the ground truth. In formation-energy prediction, the
phase diagram is a key quantity that directly demonstrates the correct physics
near the ground states. As the ground states are formed variationally, they
are particularly discerning towards spurious ECIs, as the non-physical noisy
interactions often create new ground states leading one to miss the true
ground states. Thus, a well-reproduced phase diagram is desirable for a CE
model. In our tests, the $\ell_{1}$-CE creates 12 new ground states,
indicating that the correct physics in terms of cluster interactions is not
well captured. However, the $\ell_{0}\ell_{2}$-CE preserves most of the DFT
ground states, with only four spurious ”ground states” in the
$\ell_{0}\ell_{2}$-CE phase diagram.

The improvement in the physics of the predictions associated with applying the
$\ell_{0}\ell_{2}$-norm with hierarchy constraints is further demonstrated by
the voltage profile and Li-occupancy energy. In Fig.5(b), the voltage profiles
generated by prediction using the $\ell_{1}$-CE and $\ell_{0}\ell_{2}$-CE
(blue lines) are compared with those from DFT (orange lines), taken as the
ground truth. We can see that the $\ell_{1}$-CE incorrectly predicts the
voltage plateau between $x=0.5$ to $1$ in the LixMn2O4 spinel-like structure
such that the $x=0.5$ configuration is no longer stable (the voltage between
$x=0.5$ and $x=1.0$ is higher than that between $x=0.0$ and $x=0.5$). In
contrast, the $\ell_{0}\ell_{2}$-CE matches very well with the DFT-generated
voltage profiles. The erroneous predictions of the $\ell_{1}$-CE are further
confirmed by the Li-occupancy energy. In Fig.5(c), the energy difference
between Li in octahedral and tetrahedral occupancy was evaluated in the
layered-MnO2 and spinel-MnO2 frameworks. The absolute error compared with DFT
is 0.52 eV (layered) and 0.18 eV (spinel) for the $\ell_{1}$-CE, whereas that
for the $\ell_{0}\ell_{2}$-CE is 0.09 eV (layered) and 0.09 eV (spinel),
respectively. A significant reduction of prediction error is observed with the
$\ell_{0}\ell_{2}$-norm regularized CE.

##  IV Discussion & Summary

In two complex oxide systems, we showed that the $\ell_{0}\ell_{2}$-CE with
hierarchy constraints outperforms the conventional $\ell_{1}$-CE in terms of
sparseness against CV error, convergence rate with a subset of training-data,
and some critical physical quantities in Li intercalation materials. More
generally, the optimization of the ECIs is not deterministic within a single
method, and the successful construction of a CE model typically relies on two
aspects: (1) choosing a valid interaction space by truncating the clusters or
orbits and (2) applying a proper optimization algorithm to obtain the ECIs.
The results in this paper show that for the second step, the
$\ell_{0}\ell_{2}$-norm method is the superior choice for a robust and
physical solution compared to the conventional $\ell_{1}$ method.

We note that one limitation of the $\ell_{0}\ell_{2}$-norm method in the MIQP
paradigm is the computational efficiency. As solving the $\ell_{0}$-norm is an
NP-hard problem, more computational time is required to solve the MIQP when
more predefined ECIs are included. The $\ell_{0}\ell_{2}$-CE works well for
relatively small or well-predefined systems
($\mathrm{dim}(\boldsymbol{\Pi}_{S})\leq 2000)$. Therefore, the most
applicable way to use $\ell_{0}\ell_{2}$-norm regularized CE with hierarchy
constraints is likely to be as follows: (1) define a CE within a relatively
small cutoff and truncate to quadruplet or quintuplet clusters at most
(ideally staying within $\mathrm{dim}(\boldsymbol{\Pi}_{S})\leq 2000$) and (2)
follow the procedure described in Fig.3(b) to determine the optimal
hyperparameter to obtain the ECIs. However, we note that
$\mathrm{dim}(\boldsymbol{\Pi}_{S})\leq 2000$ applies to virtually all known
published CE.

To obtain a model that represents the physics of a system well, the relative
difference of energies between configurations is of greater significance than
the absolute ones. In ordinary least-squares fitting, the cost function only
focuses on the global averaged error of the training set, which leads to over-
fitting. Adding regularization of the ECIs can alleviate this issue by
constraining the optimization space of parameters, but our results show that
not all regularization creates physically meaningful solutions. We propose
that it is beneficial to include the physically inspired constraints into the
design of the cost function, such as adding hierarchy constraints with
$\ell_{0}\ell_{2}$-norm implementation. The $\ell_{0}\ell_{2}$-CE can improve
the physical meaning of the solution and break the correlation between coupled
clusters, which is achieved by directly penalizing the number of non-zero ECIs
for feature selection and enforcing hierarchy relations between ECIs via the
slack variables in the MIQP paradigm. The $\ell_{0}\ell_{2}$-CE gives an
estimation of the ECIs with reasonable physics near the ground-states, but
does not strictly enforce the preservation of ground-states. Additionally, the
ground-state preservation can be further achieved by adding inequality
constraints on the energies into Eq.(13) as shown in previous work Huang  _et
al._ (2017).

In summary, our method sheds light on how to obtain good ECIs for simulations
in complex and coupled multicomponent systems, with several proposed criteria
in ECIs optimization: (1) minimize the CV error under general regression level
(e.g., ridge regression); (2) as the sparseness describes the complexity of
and number of independent variables, the sparseness of the solution shall be
improved (reduced) with reasonable in-sample training error; and (3) check the
near-ground-state behavior and related physics for the optimal ECIs selection.

##  V Acknowledgement

This work was funded by the U.S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, Materials Sciences and Engineering Division
under Contract No. DE-AC0205CH11231 (Materials Project program KC23MP). The
work was also supported with computational resources provided by the Extreme
Science and Engineering Discovery Environment (XSEDE), supported by National
Science Foundation grant number ACI1053575; the National Energy Research
Scientific Computing Center (NERSC), a U.S. Department of Energy Office of
Science User Facility located at Lawrence Berkeley National Laboratory; the
Center for Functional Nanomaterials (CFN), which is a U.S. Department of
Energy Office of Science User Facility, at Brookhaven National Laboratory
under Contract No. DE-SC0012704; and the Lawrencium computational cluster
resource provided by the IT Division at the Lawrence Berkeley National
Laboratory.

*

##  Appendix A Appendix: DFT calculations

DFT calculations were performed with the Vienna ab initio simulation package
(VASP) using the projector-augmented wave method Kresse and Furthmüller
(1996); Kresse and Joubert (1999), a plane-wave basis set with an energy
cutoff of 520 eV, and a reciprocal space discretization of 25 k-points per Å.
All the calculations were converged to $10^{-6}$ eV in total energy for
electronic loops and 0.02 eV/Å in interatomic forces for ionic loops. In the
LiF-MnO-LiVO2-Li2TiO3 system, we used the Perdew-Burke-Ernzerhof (PBE)
generalized gradient approximation exchange-correlation functional Perdew  _et
al._ (1996) with rotationally averaged Hubbard U correction (GGA+U) to
compensate for the self-interaction error on all transition-metal atoms except
titanium Wang  _et al._ (2006). The U parameters were obtained from the
literature, where they were calibrated to transition-metal oxide formation
energies (3.9 eV for Mn and 3.1 eV for V). The GGA+U computational framework
is believed to be reliable in determining the formation enthapies of similar
compounds Jain  _et al._ (2011). In the Li-Mn-Ti-O oxide system, the strongly
constrained and appropriately normed (SCAN) meta-GGA exchange-correlation
functional was used Sun  _et al._ (2015). The SCAN functional is believed to
have better performance at capturing charge transfer due to better redox and
atomic coordination prediction Yang  _et al._ (2019); Zhang  _et al._ (2018),
which would improve the accuracy of energetics involving introducing vacancies
on octahedral and interstitial tetrahedral sublattices in the rocksalt
framework.

## References

  * Jain  _et al._ (2013) A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, _et al._ , APL materials 1, 011002 (2013).
  * Zhang  _et al._ (2020) R. Zhang, S. Zhao, J. Ding, Y. Chong, T. Jia, C. Ophus, M. Asta, R. O. Ritchie, and A. M. Minor, Nature 581, 283 (2020).
  * Yin  _et al._ (2021) S. Yin, Y. Zuo, A. Abu-Odeh, H. Zheng, X.-G. Li, J. Ding, S. P. Ong, M. Asta, and R. O. Ritchie, Nature communications 12, 1 (2021).
  * Ji  _et al._ (2019) H. Ji, A. Urban, D. A. Kitchaev, D.-H. Kwon, N. Artrith, C. Ophus, W. Huang, Z. Cai, T. Shi, J. C. Kim, H. Kim, and G. Ceder, Nature Communications 10, 592 (2019).
  * Ouyang  _et al._ (2020) B. Ouyang, N. Artrith, Z. Lun, Z. Jadidi, D. A. Kitchaev, H. Ji, A. Urban, and G. Ceder, Advanced Energy Materials 10, 1903240 (2020).
  * Lun  _et al._ (2021) Z. Lun, B. Ouyang, D.-h. Kwon, Y. Ha, E. E. Foley, T.-Y. Huang, Z. Cai, H. Kim, M. Balasubramanian, Y. Sun, J. Huang, Y. Tian, H. Kim, B. D. McCloskey, W. Yang, R. J. Clément, H. Ji, and G. Ceder, Nature Materials 20, 214 (2021).
  * Zhong  _et al._ (2020) P. Zhong, Z. Cai, Y. Zhang, R. Giovine, B. Ouyang, G. Zeng, Y. Chen, R. Clément, Z. Lun, and G. Ceder, Chemistry of Materials 32, 10728 (2020).
  * Sanchez  _et al._ (1984) J. Sanchez, F. Ducastelle, and D. Gratias, Physica A: Statistical Mechanics and its Applications 128, 334 (1984).
  * Sanchez (2019) J. M. Sanchez, Physical Review B 99, 134206 (2019).
  * Tepesch  _et al._ (1995) P. D. Tepesch, G. D. Garbulsky, and G. Ceder, Physical Review Letters 74, 2272 (1995).
  * Seko  _et al._ (2009) A. Seko, Y. Koyama, and I. Tanaka, Physical Review B 80, 165122 (2009).
  * Zhang and Sluiter (2016) X. Zhang and M. H. F. Sluiter, Journal of Phase Equilibria and Diffusion 37, 44 (2016).
  * Barroso-Luque  _et al._ (2021) L. Barroso-Luque, J. H. Yang, and G. Ceder, Physical Review B 104, 224203 (2021).
  * van de Walle (2009) A. van de Walle, Calphad 33, 266 (2009).
  * Ceder (1993) G. Ceder, Computational Materials Science 1, 144 (1993).
  * Huang  _et al._ (2016) W. Huang, D. A. Kitchaev, S. T. Dacek, Z. Rong, A. Urban, S. Cao, C. Luo, and G. Ceder, Physical Review B 94, 134424 (2016).
  * De Fontaine (1979) D. De Fontaine, Solid state physics 34, 73 (1979).
  * Ångqvist  _et al._ (2019) M. Ångqvist, W. A. Muñoz, J. M. Rahm, E. Fransson, C. Durniak, P. Rozyczko, T. H. Rod, and P. Erhart, Advanced Theory and Simulations 2, 1900015 (2019).
  * Huang  _et al._ (2017) W. Huang, A. Urban, Z. Rong, Z. Ding, C. Luo, and G. Ceder, npj Computational Materials 3, 30 (2017).
  * Laks  _et al._ (1992) D. B. Laks, L. G. Ferreira, S. Froyen, and A. Zunger, Phys. Rev. B 46, 12587 (1992).
  * Blum and Zunger (2004) V. Blum and A. Zunger, Phys. Rev. B 70, 155108 (2004).
  * Mueller and Ceder (2009) T. Mueller and G. Ceder, Physical Review B 80, 024103 (2009).
  * Nelson  _et al._ (2013a) L. J. Nelson, V. Ozoliņš, C. S. Reese, F. Zhou, and G. L. W. Hart, Physical Review B 88, 155105 (2013a).
  * Cockayne and van de Walle (2010) E. Cockayne and A. van de Walle, Phys. Rev. B 81, 012104 (2010).
  * Candes and Wakin (2008) E. J. Candes and M. B. Wakin, IEEE Signal Processing Magazine 25, 21 (2008).
  * Nelson  _et al._ (2013b) L. J. Nelson, G. L. W. Hart, F. Zhou, and V. Ozoliņš, Physical Review B 87, 035125 (2013b).
  * Seko  _et al._ (2014) A. Seko, A. Takahashi, and I. Tanaka, Phys. Rev. B 90, 024101 (2014).
  * Leong and Tan (2019) Z. Leong and T. L. Tan, Physical Review B 100, 1 (2019).
  * Elad (2010) M. Elad, in  _Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing_, edited by M. Elad (Springer, New York, NY, 2010) pp. 35–54.
  * (30) W. Huang, A. Urban, P. Xiao, Z. Rong, H. Das, T. Chen, N. Artrith, A. Toumar, and G. Ceder,  arXiv:1807.10753 .
  * Gurobi Optimization, LLC (2021) Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,”  (2021).
  * Sluiter and Kawazoe (2005) M. H. F. Sluiter and Y. Kawazoe, Phys. Rev. B 71, 212201 (2005).
  * Zarkevich and Johnson (2004) N. A. Zarkevich and D. D. Johnson, Phys. Rev. Lett. 92, 255702 (2004).
  * Walle and Ceder (2002) A. Walle and G. Ceder, Journal of Phase Equilibria 23, 348 (2002).
  * Tibshirani (1996) R. Tibshirani, Journal of the Royal Statistical Society: Series B (Methodological) 58, 267 (1996).
  * Note (1) https://github.com/CederGroupHub/smol, https://github.com/CederGroupHub/sparse-lm.
  * Clément  _et al._ (2020) R. J. Clément, Z. Lun, and G. Ceder, Energy and Environmental Science 13, 345 (2020).
  * Huang  _et al._ (2021) J. Huang, P. Zhong, Y. Ha, D.-h. Kwon, M. J. Crafton, Y. Tian, M. Balasubramanian, B. D. McCloskey, W. Yang, and G. Ceder, Nature Energy 6, 706 (2021).
  * Seko and Tanaka (2014) A. Seko and I. Tanaka, Journal of physics: Condensed matter 26, 115403 (2014).
  * Richards  _et al._ (2018) W. D. Richards, S. T. Dacek, D. A. Kitchaev, and G. Ceder, Advanced Energy Materials 8, 1701533 (2018).
  * Bartel  _et al._ (2020) C. J. Bartel, A. Trewartha, Q. Wang, A. Dunn, A. Jain, and G. Ceder, npj Computational Materials 6, 1 (2020).
  * Urban  _et al._ (2016) A. Urban, D.-H. Seo, and G. Ceder, npj Computational Materials 2, 16002 (2016).
  * Seo  _et al._ (2015) D.-H. Seo, A. Urban, and G. Ceder, Physical Review B 92, 115118 (2015).
  * Aydinol and Ceder (1997) M. Aydinol and G. Ceder, Journal of the Electrochemical Society 144, 3832 (1997).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Physical review B 59, 1758 (1999).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Wang  _et al._ (2006) L. Wang, T. Maxisch, and G. Ceder, Physical Review B 73, 195107 (2006).
  * Jain  _et al._ (2011) A. Jain, G. Hautier, S. P. Ong, C. J. Moore, C. C. Fischer, K. A. Persson, and G. Ceder, Physical Review B 84, 1 (2011).
  * Sun  _et al._ (2015) J. Sun, A. Ruzsinszky, and J. P. Perdew, Physical Review Letters 115, 036402 (2015).
  * Yang  _et al._ (2019) J. H. Yang, D. A. Kitchaev, and G. Ceder, Phys. Rev. B 100, 035132 (2019).
  * Zhang  _et al._ (2018) Y. Zhang, D. A. Kitchaev, J. Yang, T. Chen, S. T. Dacek, R. A. Sarmiento-Pérez, M. A. Marques, H. Peng, G. Ceder, J. P. Perdew, _et al._ , npj Computational Materials 4, 1 (2018).