# Efficient first-order algorithms for large-scale, non-smooth maximum entropy
models with application to wildfire science

Gabriel P. Langlois Courant Institute of Mathematical Sciences, New York
University, New York City, NY 10012, USA<EMAIL_ADDRESS>Jatan Buch Department of Earth and Environmental Engineering, Columbia
University, New York City, NY 10027, USA<EMAIL_ADDRESS>Jérôme Darbon
Division of Applied Mathematics, Brown University, Providence, RI 02912, USA
<EMAIL_ADDRESS>
###### Abstract

Maximum entropy (Maxent) models are a class of statistical models that use the
maximum entropy principle to estimate probability distributions from data. Due
to the size of modern data sets, Maxent models need efficient optimization
algorithms to scale well for big data applications. State-of-the-art
algorithms for Maxent models, however, were not originally designed to handle
big data sets; these algorithms either rely on technical devices that may
yield unreliable numerical results, scale poorly, or require smoothness
assumptions that many practical Maxent models lack. In this paper, we present
novel optimization algorithms that overcome the shortcomings of state-of-the-
art algorithms for training large-scale, non-smooth Maxent models. Our
proposed first-order algorithms leverage the Kullback-Leibler divergence to
train large-scale and non-smooth Maxent models efficiently. For Maxent models
with discrete probability distribution of $n$ elements built from samples,
each containing $m$ features, the stepsize parameters estimation and
iterations in our algorithms scale on the order of $O(mn)$ operations and can
be trivially parallelized. Moreover, the strong $\ell_{1}$ convexity of the
Kullback–Leibler divergence allows for larger stepsize parameters, thereby
speeding up the convergence rate of our algorithms. To illustrate the
efficiency of our novel algorithms, we consider the problem of estimating
probabilities of fire occurrences as a function of ecological features in the
Western US MTBS-Interagency wildfire data set. Our numerical results show that
our algorithms outperform the state of the arts by one order of magnitude and
yield results that agree with physical models of wildfire occurrence and
previous statistical analyses of wildfire drivers.

###### keywords:

Primal-dual method, Maximum entropy estimation, Kullback–Leibler divergence,
wildfire science

90C30, 90C06, 90C90, 62P12

##  1 Introduction

Maximum entropy (Maxent) models are a class of density estimation methods that
use the maximum entropy principle [43] to reproduce key statistics of
datasets. Historically used in physics [43, 44], engineering [47, 38, 6] and
statistics [85] applications, Maxent models are now frequently used for large-
scale machine learning problems in natural language processing [7, 13, 21, 57,
74, 82], social science [41, 51, 63], neuroscience [37, 81, 75], ecological
modeling [24, 25, 30, 46, 60, 70, 69, 68, 77, 76], and wildfire science [67,
12, 87].

Maxent models in machine learning must often estimate probability
distributions from big datasets comprising hundreds of thousands to billions
of samples or features or both [56]. Due to this enormous number, large-scale
Maxent models need efficient and robust algorithms to perform well. State-of-
the-art algorithms for Maxent models, however, were not originally designed to
handle big datasets; these algorithms either rely on technical devices that
may yield unreliable numerical results [32], scale poorly in size [19, 21], or
require smoothness assumptions that many Maxent models lack in practice [57].
These limitations make it essentially impossible to scale Maxent models to big
data applications without adequate and costly computational resources [22,
80]. This constraint on computational resources, in particular, has been
recently identified as a crucial challenge to overcome for using large-scale
Maxent models in climate change and wildfire studies within a reasonable
amount of run time [27, 77, 76, 12].

### Contributions of this paper

In this paper, we present novel optimization algorithms that overcome the
shortcomings of state-of-the-art algorithms used for training large-scale,
non-smooth Maxent models. Our proposed algorithms are first-order accelerated
nonlinear primal-dual hybrid gradient (NPDHG) algorithms, whose theory was
provided by two authors of this paper in [18], based on the Kullback–Leibler
divergence. Using the Kullback–Leibler divergence over the classical proximal
operator makes it possible to train large-scale and non-smooth Maxent models
much more efficiently than the state of the arts. In particular, we show that
for a Maxent model with discrete probability distribution of $n$ elements
built from samples each containing $m$ features, the stepsize parameters
estimation and iterations in our algorithms all scale on the order of $O(mn)$
operations. This significantly improves on the known complexity bound of
$O(\min(m^{2}n,mn^{2}))$ operations for computing the optimal stepsize
parameters of classical first-order optimization methods, such as the linear
PDHG or forward-backward splitting methods. We also show, as a consequence,
that for a given tolerance level $\epsilon>0$, our algorithms provably compute
solutions using on the order of $O(mn/\sqrt{\epsilon})$ or
$O(mn/\log(1/\epsilon))$ operations, the order depending on the smoothness of
the Maxent model and which are optimal with respect to the Nesterov class of
optimal first-order methods [64]. Moreover, the computational bottleneck
consists of matrix-vector multiplications, which can be trivially
parallelized, and so our proposed algorithms exhibit scalable parallelism.
Finally, we show that the strong convexity of the Kullback–Leibler divergence
with respect to the $\ell_{1}$ norm naturally allows for significantly larger
stepsize parameters, thereby substantially speeding up the convergence rate of
our algorithms.

To illustrate the efficiency of our novel algorithms on large-scale problems,
we present an application to wildfire science. Specifically, we consider the
problem of estimating probabilities of fire occurrences across the western US
as a function of ecological features. To do so, we fit elastic net, non-
overlapping group lasso and $\ell_{\infty}$ Maxent models to a large number of
hyperparameters using our proposed algorithms and the state-of-the-art
forward-backward splitting and coordinate descent STRUCTMAXENT2 algorithms [5,
62]. Our numerical results show that our algorithms outperform the two latter
algorithms by at least one order of magnitude, and yield results that are in
good agreement with physical models of wildfire occurrence [66, 61] as well as
previous statistical analyses of wildfire drivers [67, 12, 87, 10].

### Organization of this paper

This paper is organized as follows. In Section 2, we present our setup, we
describe how Maxent models work, and we review three popular non-smooth Maxent
models that are challenging to train with big data sets: the elastic net,
group lasso and $\ell_{\infty}$ regularized Maxent models. In Section 3, we
explain why large-scale, non-smooth Maxent models are particularly challenging
to train from big datasets, and we describe the state-of-the-art algorithms
for training these Maxent models and their limitations. In Section 4, we
describe our approach for estimating Maxent models with NPDHG optimization
methods, derive explicit algorithms for regularized Maxent models, and explain
why our algorithms overcome the limitations of state-of-the-art algorithms. In
Section 5, we present an application of our algorithms to wildfire science,
where we train a large number of Maxent models on the Western US MTBS-
Interagency (WUMI) wildfire data set [45] to estimate probabilities of fire
occurrences as a function of ecological features. Finally, we review our
results and outline directions for future work in Section 6.

##  2 Preliminaries

###  2.1 Setup

Suppose we receive $l$ independent, identically distributed samples
$\\{v_{1},...,v_{l}\\}\subset\mathcal{I}$ from an unknown distribution
$\mathcal{D}$. We assume throughout, without of loss of generality, that the
input space $\mathcal{I}\equiv\\{1,\dots,n\\}$. In addition, suppose we have a
prior probability distribution $p_{\text{prior}}$ on $\mathcal{I}$ that
encapsulates some prior knowledge about the samples or unknown distribution.
Finally, suppose we have access to a set of features from the samples via a
bounded feature map $\bm{\Phi}\colon\mathcal{I}\to\mathbb{R}^{m}$. Then how do
we estimate the unknown distribution $\mathcal{D}$ from the prior distribution
$p_{\text{prior}}$, the samples $\\{v_{1},...,v_{l}\\}$ and the feature map
$\bm{\Phi}$?

The maximum entropy principle offers a systematic way to answer this question.
It states that the distribution that best estimates the unknown distribution
$\mathcal{D}$ is the one that remains as close as possible to the prior
probability $p_{\text{prior}}$ while matching the features
$\\{\bm{\Phi}(v_{1}),\dots,\bm{\Phi}(v_{l})\\}$ exactly or as closely as
possible, in some suitable sense. Specifically, let $\Delta_{n}$ denote the
$n$-dimensional probability simplex and suppose
$p_{\text{prior}}\in\text{int}~{}{\Delta_{n}}$ (i.e., $p_{\text{prior}}(j)>0$
for every $j\in\\{1,\dots,n\\}$). We measure closeness of a probability
distribution $p\in\Delta_{n}$ to the prior probability $p_{\text{prior}}$
using the Kullback–Leibler divergence
$D_{\text{KL}}\colon\Delta_{n}\times\text{int}~{}\Delta_{n}\to\mathbb{R}$:

(1)
$D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})=\sum_{j=1}^{n}p(j)\log\left(p(j)/p_{\text{prior}}(j)\right).$

Next, let $\hat{\mathcal{D}}$ denote the empirical distribution induced by
samples $\\{v_{1},...,v_{l}\\}$:

(2) $\hat{\mathcal{D}}(j)=\frac{1}{l}\left|\\{1\leqslant i\leqslant l\mid
v_{i}=j\\}\right|.$

In addition, let $\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]$ denote the
empirical average of the features induced by the samples and let
$\mathbb{E}_{p}[\bm{\Phi}]$ denote the model average induced by a probability
distribution $p\in\Delta_{n}$:

$\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]=\sum_{j=1}^{n}\hat{\mathcal{D}}(j)\bm{\Phi}(j)\quad\text{and}\quad\mathbb{E}_{p}[\bm{\Phi}]=\sum_{j=1}^{n}p(j)\bm{\Phi}(j).$

We measure how the averages $\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]$ and
$\mathbb{E}_{p}[\bm{\Phi}]$ match using a function
$H^{*}\colon\mathbb{R}^{m}\to\mathbb{R}\cup\\{+\infty\\}$ called the potential
function. Maxent models combine the two measures described above and seek to
minimize the sum

(3)
$\inf_{p\in\Delta_{n}}f(p;t)=\inf_{p\in\Delta_{n}}\left\\{D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})+tH^{*}\left(\frac{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p}[\bm{\Phi}]}{t}\right)\right\\}$

where $t>0$ is a hyperparameter selected using data-driven methods, e.g.,
cross-validation.

Standing assumptions: The potential function is proper, lower semicontinuous,
convex, and bounded from below by zero with $H^{*}(\bm{0})$ = 0. Under the
standing assumptions above on the potential function, which are fairly
general, the probability distribution $\hat{\mathcal{D}}$ is a feasible point
of Eq. 3. Moreover, the Maxent estimation problem Eq. 3 admits a unique global
solution (see [28, Page 35, Proposition 1.2]; uniqueness follows from the
strong convexity of the Kullback–Leibler divergence with respect to the
$\ell_{1}$ norm [4, 16, 49, 50, 71]).

###  2.2 Dual formulation and optimality conditions

The generalized Maxent problem Eq. 3 admits a dual problem corresponding to
regularized maximum a posteriori estimation [2]. To derive the dual problem,
we first write the second term on the right hand side of Eq. 3 in terms of its
convex conjugate,

$tH^{*}\left(\frac{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p}[\bm{\Phi}]}{t}\right)=\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\left\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p}[\bm{\Phi}]\right\rangle-
tH(\bm{w})\right\\},$

where we abuse notation and write the convex conjugate of $H^{*}$ as $H$. This
lets us express problem Eq. 3 in saddle-point form:

(4)
$\inf_{p\in\Delta_{n}}\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\left\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p}[\bm{\Phi}]\right\rangle-
tH(\bm{w})+D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})\right\\}.$

Due to our assumptions on the potential function $H^{*}$, we can swap the
infimum and supremum [28, Statement (4.1) on page 61]. Using the convex
conjugate formula

$\inf_{p\in\Delta_{n}}\left\\{D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})-\left\langle\bm{w},\mathbb{E}_{p}[\bm{\Phi}]\right\rangle\right\\}=-\log\left[\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\left\langle\bm{w},\bm{\Phi}(j)\right\rangle}\right],$

we obtain the dual problem of Eq. 3:

(5)
$\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\left\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]\right\rangle-
tH(\bm{w})-\log\left[\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\left\langle\bm{w},\bm{\Phi}(j)\right\rangle}\right]\right\\}.$

The dual problem Eq. 5 is a regularized maximum likelihood estimation problem
over the family of Gibbs distributions [2, 62]. It has at least one global
solution [28, Proposition (1.2) on page 35]. Moreover, if $(p^{s},\bm{w}^{s})$
denotes a pair of solutions to the Maxent estimation problem Eq. 3 and its
dual Eq. 5, then this pair satisfies the optimality conditions

(6)
$\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p^{s}}[\bm{\Phi}]\in
t\partial H(\bm{w}^{s})\quad\text{and}\quad
p^{s}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{w}^{s},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{w}^{s},\bm{\Phi}(j)\rangle}},$

for every $j\in\\{1,\dots,n\\}$ where $\partial H(\bm{w}^{s})$ denotes the
subdifferential of $H$ at $\bm{w}^{s}$.

###  2.3 Examples of non-smooth Maxent models

Maxent models differ in the choice of the prior distribution
$p_{\text{prior}}$, hyperparameter $t$ and the potential function $H^{*}$. In
practice, the prior distribution is often chosen to be uniform and the
hyperparameter is chosen via data-driven methods. The choice of potential
function depends on the application. The classical Maxent model uses the
indicator function

$\bm{u}\mapsto H^{*}(\bm{u})=\begin{dcases}0,\;&\text{if}\;\bm{u}=\bm{0},\\\
+\infty\;&\text{otherwise}.\end{dcases}$

This forces the model average $\mathbb{E}_{p}[\bm{\Phi}]$ to be equal to the
empirical average $\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]$. However, this
is often too restrictive because the averages are expected to be close, and
not equal, with high probability. Forcing the model and empirical averages to
be equal can lead to severe over-fitting issues with big data sets, and more
flexibility is therefore desired.

We review below three different non-smooth Maxent models that offer this
flexibility: the elastic net, group lasso and $\ell_{\infty}$-regularized
Maxent models. We chose these models because they are used extensively in
statistics and machine learning but are also challenging to train on big data
sets due to their non-smoothness [25, 68, 23, 77]. These models will also be
used for our numerical experiments and results in Section 5.

#### Elastic net regularized Maxent models

These models use for potential function the relaxation of the convex conjugate
of the $\ell_{1}$-norm:

$\bm{u}\mapsto
H^{*}(\bm{u})\equiv\left(\frac{1-\alpha}{2}\left\|{\cdot}\right\|_{2}^{2}+\alpha\left\|{\cdot}\right\|_{1}\right)^{*}(\bm{u})=\frac{1}{2(1-\alpha)}\sum_{i=1}^{m}\max(0,|u_{i}|-\alpha)^{2}$

where $\alpha\in(0,1)$. This yields the elastic net regularized maximum
entropy estimation problem

(7)
$\min_{p\in\Delta_{n}}\left\\{D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})+\frac{1}{2t(1-\alpha)}\sum_{i=1}^{m}\max\left(0,\left|[\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]]_{i}-[\mathbb{E}_{p}[\bm{\Phi}]]_{i}\right|-t\alpha\right)^{2}\right\\}.$

This type of regularization is frequently used for feature selection in
machine learning [91, 79, 82, 25]. The corresponding dual problem is

(8)
$\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\left\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]\right\rangle-t\left(\frac{(1-\alpha)}{2}\left\|{\bm{w}}\right\|_{2}^{2}+\alpha\left\|{\bm{w}}\right\|_{1}\right)-\log\left[\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\left\langle\bm{w},\bm{\Phi}(j)\right\rangle}\right]\right\\}.$

Thanks to the elastic net penalty
$\frac{(1-\alpha)}{2}\left\|{\bm{w}}\right\|_{2}^{2}+\alpha\left\|{\bm{w}}\right\|_{1}$,
problem Eq. 8 is strongly concave and has a unique solution.

An important aspect of the elastic net penalty is that it promotes sparsity;
that is, it leads to a solution with many entries equal to zero [33, 29, 89],
the number depending on the hyperparameter $t$ and parameter $\alpha$.
Sparsity is useful when seeking accurate solutions based only on a few
features or when performing feature selection.

#### Group lasso regularized Maxent models

These models are specified starting from the dual problem. Consider a
partition $G$ of the components of a vector $\bm{w}\in\mathbb{R}^{m}$ into $g$
different and possibly overlapping groups $[\bm{w}_{1},\dots,\bm{w}_{g}]$ with
$\bm{w}_{g}\in\mathbb{R}^{m_{g}}$ and $\cup_{g=1}^{G}\bm{w}_{g}=\bm{w}$. Let
$\bm{w}\mapsto\left\|{\bm{w}_{g}}\right\|_{2,g}$ denote the $\ell_{2}$ norm
over the respective components of the $g^{\text{th}}$ group. Then the dual
version of the group lasso regularized maximum entropy estimation problem is

(9)
$\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]\rangle-t\sum_{g=1}^{G}\sqrt{m_{g}}\left\|{\bm{w}_{g}}\right\|_{2,g}-\log\left[\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\left\langle\bm{w},\bm{\Phi}(j)\right\rangle}\right]\right\\}.$

Thanks to the group lasso penalty
$\sum_{g=1}^{G}\sqrt{m_{g}}\left\|{\bm{w}_{g}}\right\|_{2,g}$, problem Eq. 10
is strictly concave and so has at least one global solution. The corresponding
primal problem follows from the convex conjugate formula

$\bm{u}\mapsto
H^{*}(\bm{u})\equiv\left(\sum_{g=1}^{G}\sqrt{m_{g}}\left\|{\cdot}\right\|_{2,g}\right)^{*}(\bm{u})=\begin{dcases}0,\;\text{if}\;\left\|{\bm{u}_{g}}\right\|_{2,g}\leqslant\sqrt{m_{g}}\;\text{for}\;g=1,\dots,G,\\\
+\infty,\;\text{otherwise,}\end{dcases}$

and it is therefore given by

(10)
$\min_{p\in\Delta_{n}}D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})\quad\text{s.t.}\quad\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}_{g}]-\mathbb{E}_{p}[\bm{\Phi}_{g}]}\right\|_{2,g}\leqslant
t\sqrt{m_{g}}\quad\text{for}\quad g=1,\dots,G.$

Similarly to the elastic net penalty, the group lasso penalty promotes
sparsity in the global solutions of Eq. 9, but it differs in that it sets
different blocks of components $[\bm{w}_{1},\dots,\bm{w}_{g}]$ to zero at the
same time.

####  $\ell_{\infty}$-regularized Maxent models

These models use for potential function the convex conjugate of the
$\ell_{\infty}$ norm,

$\bm{u}\mapsto
H^{*}(\bm{u})\equiv\left\|{\bm{u}}\right\|_{\infty}^{*}=\left\\{\hat{\bm{u}}\in\mathbb{R}^{m}\mid\left\|{\bm{u}-\hat{\bm{u}}}\right\|_{1}\leqslant
1\right\\},$

which is the characteristic set of the $\ell_{1}$ ball with unit radius. This
yields the $\ell_{\infty}$-regularized maximum entropy estimation problem

(11)
$\min_{p\in\Delta_{n}}\left\\{D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})+t\left\|{\frac{\mathbb{E}_{\hat{\mathcal{D}}}{[\bm{\Phi}]}-\mathbb{E}_{p}{[\bm{\Phi}]}}{t}}\right\|_{\infty}^{*}\right\\}.$

This type of regularization is used for certain machine learning problems,
including matrix factor decomposition, and in robust statistics applications
(see, e.g., [52, 23]). The corresponding dual problem is

(12)
$\sup_{\bm{w}\in\mathbb{R}^{m}}\left\\{\langle\bm{w},\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]\rangle-t\left\|{\bm{w}}\right\|_{\infty}-\log\left[\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\left\langle\bm{w},\bm{\Phi}(j)\right\rangle}\right]\right\\}$

Thanks to the $\ell_{\infty}$ norm, problem Eq. 12 is strictly concave and has
at least one global solution.

##  3 Related work

###  3.1 Large-scale sparse Maxent models: Computational challenges

Estimating a probability distribution from the Maxent model Eq. 3 can be
computationally prohibitive for big data sets. To illustrate this point,
suppose the Maxent model Eq. 3 with hyperparameter $t\geqslant 0$ has a global
solution $p^{s}(t)$. Let $p^{\epsilon}(t)\in\Delta_{n}$ with $\epsilon>0$
denote an $\epsilon$-approximate solution to the global solution $p^{s}(t)$,
that is, the objective function $f$ in Eq. 3 satisfies

$f(p^{\epsilon}(t);t)-f(p^{s}(t);t)<\epsilon.$

If the potential function $H^{*}$ is smooth (equivalently, its conjugate $H$
is strongly convex), the best achievable rate of convergence for computing
$p^{\epsilon}(t)$ with $t>0$ in the Nesterov class of optimal first-order
methods is linear $O(\log(1/\epsilon))$ in the number of operations [64]. If
$H$ is not smooth, then the optimal convergence rate is sublinear
$O(1/\sqrt{\epsilon})$ in the number of operations.

These rates, while optimal, require carefully fine-tuned stepsize parameters.
In classical first-order optimization algorithms, these stepsize parameters
are fine-tuned using a precise estimate of the largest singular value of the
linear operator $\bm{A}\colon\Delta_{n}\to\mathbb{R}^{m}$ defined by

(13) $\bm{A}p=\sum_{j=1}^{n}p(j)\bm{\Phi}(j)=\mathbb{E}_{p}[\bm{\Phi}].$

Unfortunately, computing the largest singular value of the linear operator
$\bm{A}$ in Eq. 13 accurately is often computationally expensive for large
matrices due to its prohibitive cubic cost of $O(\min{(m^{2}n,mn^{2})})$
operations [40]. In this situation, line search methods and other heuristics
can sometime be employed to bypass this issue, but they typically slow down
the convergence speed. Another approach is to compute a crude estimate of the
largest singular value of the matrix $\bm{A}$, but doing so significantly
reduces convergence speed as well. In fact, even if the largest singular value
can be computed quickly, the resulting stepsize parameters may be much smaller
than what is permissible to maintain convergence, that is, the largest
singular value itself may be a poor estimate for determining how large the
stepsize parameters are allowed to be to maintain convergence. This point will
discussed in more detail in Section 4.1 when we describe our proposed NPDHG
algorithms.

This issue makes solving the Maxent model Eq. 3 difficult and laborious. Even
worse, in some applications, the appropriate value of the hyperparameter $t$
in Eq. 3 is difficult to guess and must be selected by repeatedly solving Eq.
3 from a large pool of values of hyperparameters, a process that can become
particularly time-consuming and resource intensive for big data sets. This
issue has driven much research in developing robust and efficient algorithms
to minimize computational costs and maximize model performance.

###  3.2 State-of-the-art methods for large-scale, non-smooth Maxent models

State-of-the-art methods for computing solutions to large-scale, non-smooth
Maxent models are based on coordinate descent algorithms [34, 35, 39, 15] and
first-order optimization algorithms such as forward-backward splitting [5, 11,
20]. We’ll discuss these methods below, but note that other types of methods
have been developed for Maxent models; see, e.g., [26, 57, 58, 62] for surveys
and comparisons of different algorithms. In particular, we will not consider
second-order based methods suitable for smooth Maxent models such as limited-
memory BFGS algorithms [57, 3] since this work focuses on non-smooth Maxent
models.

Coordinate descent methods. The state-of-the-art is a coordinate descent
algorithm based on a technical device called infinitely weighted logistic
regression (IWLR) [32, 68]. The IWLR method approximates a Maxent model as a
logistic regression model and then fits the approximate logistic regression
model using an existing, efficient optimization algorithm for logistic
regression. This indirect approach was proposed because efficient and scalable
coordinate descent algorithms were already available for fitting logistic
regression models, which in [32] was an earlier version of the GLMNET software
package [39]. The IWLR method remains the state-of-the-art for this reason. It
is implemented, for example, in the Maxent package available in the R
programming language [39, 68].

The IWLR method, however, is an approximate technical device that is not
guaranteed to work, a fact acknowledged by the authors who proposed the method
[32], and therefore may not produce reliable numerical results. Coordinate
descent algorithms themselves are generally non-parallelizable and often lack
robustness and good convergence properties. For example, the aforementioned
GLMNET software package approximates the logarithm term in logistic regression
models with a quadratic term to fit the models efficiently. Without costly
step-size optimization, which GLMNET avoids to improve performance, the GLMNET
implementation may not converge [35, 53]. Case in point, [88] provides two
numerical experiments in which GLMNET does not converge.

Other coordinate descent algorithms have been developed to compute solutions
to non-smooth Maxent models directly (see, e.g., [86, 15]) but they are also
generally non-parallelizable and often lack robustness and good convergence
properties. Finally, another issue is that many coordinate descent algorithms,
including GLMNET, were designed for sparse Maxent models (e.g., the elastic
net and group lasso regularized Maxent models described in Section 2.3) and
those algorithms depend on the sparsity of the model to converge quickly [91].
It would be desirable to have fast optimization methods for when the feature
mapping Eq. 13 is dense, as this often occurs in practice.

First-order methods. First-order optimization algorithms such as the forward-
backward splitting algorithm are popular because they are robust and can
provably compute $\epsilon$-approximate solutions of Eq. 3 with an optimal
rate of convergence. However, as discussed in detail at the end of the
previous subsection, achieving this rate of convergence requires fine-tuning
the stepsize parameters using an accurate estimate of the largest singular
value of the feature mapping Eq. 13, and this estimate is typically
computationally expensive for large matrices. This problem makes classical
first-order optimization algorithms inefficient and impractical in estimating
probability densities from large-scale Maxent models.

Summary. The current state-of-the-art algorithms for estimating probability
densities from Maxent models either produce unreliable numerical results, lack
scalable parallelism or scale poorly in size. These shortcomings in terms of
robustness and efficiency make it challenging to use non-smooth, large-scale
Maxent models in applications without access to adequate and costly
computational resources. The next section presents novel, efficient and robust
accelerated NPDHG optimization methods that address these shortcomings.

##  4 Main results

We describe here our approach for computing solutions to the generalized
Maxent problem Eq. 3 using accelerated NPDHG optimization methods [18]. In
addition to the standing assumptions on the potential function $H^{*}$, we
assume we can compute efficiently the proximal operator of the convex
conjugate of the potential function

$\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{\lambda
H(\bm{w})+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2}^{2}\right\\},$

for every $\lambda>0$ and $\hat{\bm{w}}\in\mathbb{R}^{m}$. These assumptions
are satisfied for most potential functions used in practice, including for the
potential functions described in Section 2.3.

###  4.1 Methodology

We start with the saddle-point formulation Eq. 4 of the generalized Maxent
estimation problem Eq. 3. Based on the work two of the authors provided in
[18], we propose to split the infimum and supremum in the saddle-point problem
Eq. 4 using an iterative NPDHG scheme that alternates between a nonlinear
proximal ascent step based on the Kullback–Leibler divergence and a linear
proximal descent step.

More precisely, let $\tau_{0}>0$ and $\sigma_{0}>0$ be two stepsize parameters
satisfying the inequality

(14)
$\tau_{0}\sigma_{0}\left(\max_{j\in\\{1,\dots,n\\}}\left\|{\bm{\Phi}(j)}\right\|_{2}^{2}\right)\leqslant
1,$

let $\theta_{0}=0$, let
$\bm{w}_{-1}=\bm{w}_{0}\in\mathrm{int}{~{}\mathrm{dom~{}}H}$, let
$\bm{z}_{0}\in\mathbb{R}^{m}$ and define the initial probability distribution
$p_{0}\in\Delta_{n}$ through the constant $\bm{z}_{0}$ via

(15)
$p_{0}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}}$

for each $j\in\\{1,\dots,n\\}$. Our proposed NPDHG iterative scheme computes
the iterates

(16) $\displaystyle p_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{p\in\Delta_{n}}\left\\{\tau_{k}D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})-\tau_{k}\left\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\mathbb{E}_{p}[\bm{\Phi}]\right\rangle+D_{\mathrm{KL}}(p\;\|\;p_{k})\right\\}$
$\displaystyle\hat{\bm{w}}_{k+1}$
$\displaystyle=\bm{w}_{k}+\sigma_{k}(\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k+1}}[\bm{\Phi}]),$
$\displaystyle\bm{w}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{t\sigma_{k}H(\bm{w})+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}_{k+1}}\right\|_{2}^{2}\right\\},$
$\displaystyle\theta_{k+1}$
$\displaystyle=1/\sqrt{1+\tau_{k}},\quad\tau_{k+1}=\theta_{k+1}\tau_{k}\quad\text{and}\quad\sigma_{k+1}=\sigma_{k}/\theta_{k+1}.$

According to [18, Proposition 4.1], the sequence of iterates $p_{k}$ converges
strongly to the unique solution of the generalized Maxent problem Eq. 3.
Moreover, for any $t\geqslant 0$ and a given tolerance level $\epsilon>0$, the
scheme Eq. 16 provably computes an $\epsilon$-approximate solution
$p^{\epsilon}(t)$ of the generalized Maxent model Eq. 3 in
$O(1/\sqrt{\epsilon})$ time. This rate of convergence is, without further
smoothness assumptions on the potential function $H^{*}$, the best achievable
rate of convergence with respect to the Nesterov class of optimal first-order
methods [64].

The key element in this scheme is the choice of the Kullback–Leibler
divergence as a nonlinear proximal step in the first line of Eq. 16. We use it
for two reasons: First, because the Kullback–Leibler divergence already
appears in the saddle-point problem Eq. 4. This allows us to compute $p_{k+1}$
explicitly. Indeed, thanks to the choice of initial probability distribution
Eq. 15, we have

(17)
$p_{k+1}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}\;\text{with}\;\bm{z}_{k+1}=\frac{1}{1+\tau_{k}}\left(\bm{z}_{k}+\tau_{k}(\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}))\right)$

for each $j\in\\{1,\dots,n\\}$. See Appendix A for the derivation of Eq. 17.

Second, because the Kullback–Leibler divergence is $1$-strongly convex with
respect to the $\ell_{1}$ norm, that is,

$D_{\mathrm{KL}}(p\;\|\;p_{k})\geqslant\frac{1}{2}\left\|{p-p_{k}}\right\|_{1}^{2}.$

This fact follows from a fundamental result in information theory known as
Pinsker’s inequality [4, 16, 49, 50, 71]. In particular, this means that the
scheme Eq. 16 alternate between solving a strongly convex problem over the
space $(\mathbb{R}^{n},\left\|{\cdot}\right\|_{1})$ and a concave problem over
the space $(\mathbb{R}^{m},\left\|{\cdot}\right\|_{2})$. The choice of these
spaces is significant, for the induced operator norm
$\left\|{\cdot}\right\|_{\mathrm{op}}$ of the linear operator $\bm{A}$ defined
in Eq. 13 becomes

$\left\|{\bm{A}}\right\|_{\mathrm{op}}=\sup_{\left\|{p}\right\|_{1}=1}\left\|{\bm{A}p}\right\|_{2}=\max_{j\in\\{1,\dots,n\\}}\left\|{\bm{\Phi}(j)}\right\|_{2}.$

This operator norm offers two crucial advantages: First, it can be computed in
optimal $\Theta(mn)$ time, or better if the features
$\\{\bm{\Phi}(j)\\}_{j=1}^{n}$ have structure, e.g., if they are sparse. This
means that the stepsize parameters $\tau_{0}$ and $\sigma_{0}$ of the NPDHG
scheme can be computed in Eq. 14, with equality, in optimal $\Theta(mn)$ time.
This is important because in classical first-order optimization methods, we
typically require a precise estimate estimate of the largest singular value of
the feature mapping Eq. 13, namely the number

$\left\|{\bm{A}}\right\|_{2}=\sup_{\left\|{\bm{x}}\right\|_{2}=1}\left\|{\bm{A}p}\right\|_{2},$

to fine-tune the stepsize parameters to gain an optimal convergence rate.
However, as discussed in detail in Section 3.1, this estimate is
computationally expensive for large matrices due to its prohibitive cubic cost
of $O(\min{(m^{2}n,mn^{2})})$ operations. In contrast, our NPDHG scheme Eq. 16
does not suffer from this computational bottleneck and therefore scales well.
Second, the operator norm $\left\|{\bm{A}}\right\|_{\mathrm{op}}$ can be
significantly smaller than the estimate $\left\|{\bm{A}}\right\|_{2}$, hence
allowing for bigger stepsize parameters to further speed up convergence111An
easy calculation yields
$\max_{\left\|{\bm{p}}\right\|_{2}=1}\left\|{\bm{A}\bm{p}}\right\|_{2}\geqslant\max_{j\in\\{1,\dots,n\\}}\left\|{\bm{\Phi}(j)}\right\|_{2}$..

Summary. To solve the generalized Maxent estimation problem Eq. 3 and its dual
problem Eq. 5, let $\tau_{0}>0$ and $\sigma_{0}>0$ be two stepsize parameters
satisfying inequality Eq. 14, let $\theta_{0}=0$, let
$\bm{w}_{-1}=\bm{w}_{0}\in\mathrm{int}{~{}\mathrm{dom~{}}H}$, let
$\bm{z}_{0}\in\mathbb{R}^{m}$ and define an initial probability distribution
$p_{0}$ through $\bm{z}_{0}$ via Eq. 15. Then compute the iterates

(18) $\displaystyle\bm{z}_{k+1}$
$\displaystyle=\left(\bm{z}_{k}+\tau_{k}(\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}))\right)/(1+\tau_{k}),$
$\displaystyle p_{k+1}(j)$
$\displaystyle=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}\quad\text{for}\;j\in\\{1,\dots,n\\},$
$\displaystyle\hat{\bm{w}}_{k+1}$
$\displaystyle=\bm{w}_{k}+\sigma_{k}(\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k+1}}[\bm{\Phi}]),$
$\displaystyle\bm{w}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{t\sigma_{k}H(\bm{w})+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}_{k+1}}\right\|_{2}^{2}\right\\},$
$\displaystyle\theta_{k+1}$
$\displaystyle=1/\sqrt{1+\tau_{k}},\quad\tau_{k+1}=\theta_{k+1}\tau_{k}\quad\text{and}\quad\sigma_{k+1}=\sigma_{k}/\theta_{k+1},$

until convergence is achieved. As all parameters and updates can be computed
in $O(mn)$ time, for any $t\geqslant 0$ and a given tolerance level
$\epsilon>0$ the overall complexity for computing an $\epsilon$-approximate
solution $p^{\epsilon}(t)$ is $O(mn/\sqrt{\epsilon})$.

###  4.2 Algorithm for smooth potential functions

The iterative scheme Eq. 18 does not require the potential function $H^{*}$ to
be smooth. If, however, the potential function $H^{*}$ is
$\gamma_{H^{*}}$-smooth (equivalently, $H$ is
$\frac{1}{\gamma_{H^{*}}}$-strongly convex) for some $\gamma_{H^{*}}>0$, then
we can modify the NPDHG iterative scheme Eq. 18 to achieve a linear rate of
convergence. More precisely, let $t>0$ and introduce the stepsize parameters

(19)
$\theta=1-\frac{t}{2\gamma_{H^{*}}\left\|{\bm{A}}\right\|_{\mathrm{op}}^{2}}\left(\sqrt{1+\frac{4\gamma_{H^{*}}\left\|{\bm{A}}\right\|_{\mathrm{op}}^{2}}{t}}-1\right),\quad\tau=\frac{1-\theta}{\theta}\quad\mathrm{and}\quad\sigma=\frac{\gamma_{H^{*}}\tau}{t}.$

Let $\bm{w}_{-1}=\bm{w}_{0}\in\mathrm{int}{~{}\mathrm{dom~{}}H}$, let
$\bm{z}_{0}\in\mathbb{R}^{m}$ and define $p_{0}$ through $\bm{z}_{0}$ via Eq.
15. Then the explicit NPDHG iterative scheme is

(20) $\displaystyle\bm{z}_{k+1}$
$\displaystyle=\left(\bm{z}_{k}+\tau(\bm{w}_{k}+\theta(\bm{w}_{k}-\bm{w}_{k-1}))\right)/(1+\tau),$
$\displaystyle p_{k+1}(j)$
$\displaystyle=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}\quad\text{for}\;j\in\\{1,\dots,n\\},$
$\displaystyle\hat{\bm{w}}_{k+1}$
$\displaystyle=\bm{w}_{k}+\sigma(\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k+1}}[\bm{\Phi}]),$
$\displaystyle\bm{w}_{k+1}$
$\displaystyle=\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{t\sigma
H(\bm{w})+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}_{k+1}}\right\|_{2}^{2}\right\\}.$

According to [18, Proposition 4.3], the sequences of iterates $p_{k}$ and
$\bm{w}_{k}$ converge strongly to the unique solution of the generalized
Maxent estimation problem Eq. 3 and its dual problem Eq. 5. Moreover, for any
$t>0$ and a given tolerance level $\epsilon>0$, this scheme provably computes
an $\epsilon$-approximate solution $p^{\epsilon}(t)$ of the generalized Maxent
estimation problem Eq. 3 in $O(\log(1/\epsilon))$. This rate of convergence is
the best achievable rate of convergence with respect to the Nesterov class of
optimal first-order methods [64]. As all parameters and updates can be
computed in $O(mn)$ time, the overall complexity for computing an
$\epsilon$-approximate solution $p^{\epsilon}(t)$ is $O(mn\log(1/\epsilon))$.

##  5 Application to wildfire science

To illustrate the efficiency of our novel algorithms on large-scale problems,
we present here an application to wildfire science. The problem at hand is to
combine fire occurrence data with ecological data in a fixed geographical
region to estimate the probability of fire occurrences as a function of
ecological features. Maxent models achieve this goal by translating fire
occurrence and ecological data into probabilities of fire occurrences and
ecological features. This approach closely mirrors how Maxent models are used
for modeling species geographic distributions [24, 25, 46, 70, 69, 68, 77,
76]. Another related goal is to identify what ecological features correlate
most with fire occurrences. This can be achieved using a sparse Maxent model,
e.g., an elastic net or group lasso regularized Maxent model, to identify
ecological features correlating significantly with fire occurrences.

For this application, we use the Western US MTBS-Interagency (WUMI) wildfire
data set [45], which we describe in Section 5.1 below. We formulate the
problem of combining the fire occurrence and ecological data from the WUMI
wildfire data set into a Maxent estimation problem in Section 5.2. Using this
data, we then fit the elastic net, (non-overlapping) group lasso and
$\ell_{\infty}$ Maxent models for a large number of hyperparameters weighting
the regularization. We detail this fitting procedure and the explicit NPDHG
algorithms for these Maxent models in Section 5.3. Following this, we compare
in Section 5.4 the running times required to fit the aforementioned Maxent
models using our NPDHG algorithms with the forward-backward splitting
algorithm [5, 11] and the STRUCTMAXENT2 coordinate descent algorithm of [62].
Finally, in Section 5.5, we interpret the results obtained from fitting the
aforementioned Maxent models to the WUMI wildfire data set.

###  5.1 WUMI wildfire data set

The Western US MTBS-Interagency wildfire data set [45] consists of all fires
($\geqslant 1\,{\rm km}^{2}$) from multiple state and federal agencies,
supplemented by satellite observations of large fires ($\geqslant 4\,{\rm
km}^{2}$) from the Monitoring Trends in Burn Severity (MTBS) program, in the
continental United States west of $103^{\circ}$ W longitude. For this
application, we extracted all wildfires from the WUMI data set that occurred
between 1984-2020 inclusive (accessed May 18, 2023). The locations of all
fires used are shown in Fig. 1.

Figure 1: Wildfire activity in the western United States from 1984 to 2020.
(Left) Fire locations of all fires (black dots) in the Western US MTBS-
Interagency (WUMI) data set; also shown are three ecological divisions
characterized by their primary vegetation type – forests (green), deserts
(yellow), and plains (gray). (Right) Prior distribution indicating mean fire
probability across all calendar months.

Next, following the procedure outlined in [10], we overlay the fire locations
on a $12\,{\rm km}\times 12\,{\rm km}$ grid to construct a data frame of
prevailing climate, vegetation, topographic, and human-related features for
each fire. Altogether we include a total of 35 potential fire-related features
from various sources; we provide a summary of all features used in our
analysis along with their sources in Appendix B. For each grid cell, we also
provide an index identifying its Environmental Protection Agency (EPA) Level
III ecoregion. Defined on the basis of roughly homogeneous climate and
vegetation patterns, ecoregions are commonly used in the wildfire science
literature to identify climate-fire relationships at a coarse spatial scale
[55, 67]. The full WUMI wildfire data set for all months from 1984-2020 is
publicly available [9] as part of the Stochastic Machine Learning for wildfire
activity in the western US (SMLFire1.0) package.

###  5.2 Data preprocessing

We construct a weakly informative prior [59, 36] for incorporating existing
knowledge of wildfire conducive environmental conditions in our Maxent model.
First, we prepare a training data set of 10000 fire occurrences and absences
chosen randomly across all months between 1984-2020 correlated with the values
of 35 fire-related features described in the previous section. The features
are averaged over each calendar month (i.e., January, February, …) between
1984 and 2020. We then apply min-max scaling to each feature, ensuring that
all features lie in the same range. Second, we construct two Random Forest
(RF) models, one each for fires in dry (May – September) and wet (October –
April) seasons respectively. We withhold 20$\%$ of the training data to tune
the model hyperparameters, such that the optimal hyperparameters ensure a
trade-off between model precision and recall. Next, we predict the fire
probability for all grid cells with fire-related features using either the
trained wet or dry season RF model depending on the calendar month. Since the
Maxent algorithm assumes a presence-only framework, that is, the absence of
fires in a grid cell does not imply a non-zero fire probability, we impute
grid cells without a fire with probability, $p_{\rm nfire}={\rm min}(p_{\rm
fire})/10$. Last, we normalized the prior probability distribution to ensure
that the fire probability across all grid cells sums up to one. For
convenience, we represent the prior distribution as the normalized mean of 12
monthly fire probabilities predicted by the RF models in Fig. 1.

To construct the empirical distribution $\hat{\mathcal{D}}$, we divide the
study region into its 18 different EPA level III ecoregions and weigh the
relative frequencies of fire among the ecoregions using the following
strategy. For each ecoregion $r\in\\{1,\dots,18\\}$, let $n_{\text{r,total}}$
denote the total number of grid cells in ecoregion $r$ and let
$n_{\text{r,fire}}$ denote the total number grid cell in ecoregion $r$ where
at least one fire was recorded. In addition, let
$Z=\sum_{r=1}^{18}\frac{n_{\text{r,fire}}}{n_{\text{r,total}}}$ denote the sum
of the proportions of grid cells where at least one fire was recorded among
all ecoregions. Then we compute the empirical probability at cell
$j\in\\{1,\dots,n\\}$, $\hat{\mathcal{D}}(j)$, as follows:

(21) $\hat{{\mathcal{D}}}(j)=\frac{1}{Z}\begin{cases}0&\;\text{if no fire was
recorded in grid cell $j$},\\\
1/n_{\text{r,total}}&\;\parbox[t]{260.17464pt}{if at least one fire was
recorded in grid $j$ and the grid cell $j$ belongs to ecoregion
$r$}\end{cases}$

This construction gives the empirical distribution more weight to ecoregions
where fires are more frequent and widespread.

###  5.3 Fitting procedure and algorithmic setup

For the analysis, we fit six different Maxent models to the wildfire data:
four elastic net Maxent models with parameters
$\alpha=\\{0.95,0.4,0.15,0.05\\}$, a non-overlapping group lasso Maxent model,
and an $\ell_{\infty}$-regularized Maxent model. For each model, we fit a
regularization path of 141 hyperparameters as follows:

$t^{(l)}=\begin{cases}(1-l/100)t^{(0)}&\;\text{from}\;l=0\,\text{to}\,50\\\
(0.5-(l-50)/200)t^{(0)}&\;\text{from}\;l=51\,\text{to}\,140,\end{cases}$

where $t^{(0)}$ depends on the Maxent model and corresponds to the smallest
hyperparameter for which the primal and dual solutions are equal to the prior
distribution and zero. More precisely, for each model we set
$(p^{(0)},\bm{w}^{(0)})=(p_{\text{prior}},\bm{0})$ and compute the sequence of
solutions $\\{(p^{(l)},\bm{w}^{(l)})\\}_{l=1}^{140}$ of the corresponding
Maxent primal and dual problems Eq. 7 and Eq. 8 with hyperparameter
$t=t^{(l)}$.

We chose these Maxent models to study the impact of the regularization on the
primal and dual solutions as a function of the hyperparameters. In particular,
for the elastic net and non-overlapping group lasso Maxent models, we are
interested in identifying the set of features that are selected or discarded
as a function of the sequence of hyperparameters $\\{t^{(l)}\\}_{l=0}^{140}$.
In the parts below, we describe the value of $t^{(0)}$ and the NPDHG algorithm
used for each model.

#### Elastic net Maxent models

For these Maxent models, the smallest hyperparameter for which the solutions
to the primal and dual problems Eq. 7 and Eq. 8 are equal to the prior
distribution and zero is

$t^{(0)}=\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{\text{prior}}}[\bm{\Phi}]}\right\|_{\infty}/\alpha.$

This follows from the optimality condition

(22)
$\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p^{s}}[\bm{\Phi}]-t(1-\alpha)\bm{w}^{s}\in
t\alpha\partial\left\|{\cdot}\right\|_{1}\bm{w}^{s}$

and computing the smallest parameter $t$ for which Eq. 22 is satisfied at
$(p^{s},\bm{w}^{s})=(p_{\text{prior}},\bm{0})$.

For the elastic net Maxent model with $\alpha=0.95$, we use the NPDHG
algorithm Eq. 18 (with sublinear convergence rate), while for the models with
$\alpha=\\{0.4,0.15,0.05\\}$, we use the NPDHG algorithm Eq. 20 (with linear
convergence rate). Starting from $l=1$, we compute the pair of solutions
$(p^{(l)},\bm{w}^{(l)})$ to the primal and dual problems Eq. 7 and Eq. 8 at
hyperparameter $t=t^{(l)}$ using the previously computed pair of solutions
$(p^{(l-1)},\bm{w}^{(l-1)})$ by setting the initial vectors
$\bm{z}_{0}=\bm{w}_{0}=\bm{w}^{(l-1)}$ in both NPDHG algorithms. For the
stepsize parameters, we set $\theta_{0}=0$, $\tau_{0}=2$ and
$\sigma_{0}=1/2\left\|{\bm{\left\|{A}\right\|}}\right\|_{\mathrm{op}}^{2}$ in
Eq. 18, and we set $\theta$, $\tau$ and $\sigma$ according to the formulas in
Eq. 19 with $\gamma_{H^{*}}=1-\alpha$ in Eq. 20. We compute the update
$\bm{w}_{k+1}$ in both NPDHG algorithms using the classical soft thresholding
operator [20, 31, 54]. Specifically, for any $\lambda>0$,
$\hat{\bm{w}}\in\mathbb{R}^{m}$ and $i\in\\{1,\dots,m\\}$ we have

$\left[\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{\lambda\left\|{\bm{w}}\right\|_{1}+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2}^{2}\right\\}\right]_{i}\equiv\left[\text{shrink}_{1}(\hat{\bm{w}},\lambda)\right]_{i}=\begin{dcases}[\hat{\bm{w}}]_{i}-\lambda&\quad\mathrm{if}\,[\hat{\bm{w}}]_{i}>\lambda,\\\
0&\quad\mathrm{if}\,|[\hat{\bm{w}}]_{i}|\leqslant\lambda,\\\
[\hat{\bm{w}}]_{i}+\lambda&\quad\mathrm{if}\,[\hat{\bm{w}}]_{i}<-\lambda,\end{dcases}$

and so, for every $\lambda>0$, $\alpha\in[0,1]$,
$\hat{\bm{w}}\in\mathbb{R}^{m}$ and $i\in\\{1,\dots,m\\}$ we have

$\left[\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{\lambda\left(\alpha\left\|{\bm{w}}\right\|_{1}+\frac{(1-\alpha)}{2}\left\|{\bm{w}}\right\|_{2}^{2}\right)+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2}^{2}\right\\}\right]_{i}=\frac{[\text{shrink}_{1}(\hat{\bm{w}},\lambda\alpha)]_{i}}{1+\lambda(1-\alpha)}.$

Finally, we let the NPDHG algorithms run for at least 40 iterations before
checking for convergence. We stop the NPDHG algorithms when the optimality
condition Eq. 22 is satisfied within some tolerance $10^{-5}$:

$\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k}}[\bm{\Phi}]-t(1-\alpha)\bm{w}_{k}}\right\|_{\infty}\leqslant
t\alpha(1+10^{-5}).$

#### Non-overlapping group lasso Maxent model

For this Maxent model, we divide the features into five disjoint groups of
features, as described in Supplement Appendix B. Then, the smallest
hyperparameter for which the solutions to the primal and dual problems Eq. 7
and Eq. 8 are equal to the prior distribution and zero is

$t^{(0)}=\max_{g\in\\{1,\dots,5\\}}\left\\{\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}_{g}]-\mathbb{E}_{p}[\bm{\Phi}_{g}]}\right\|_{2,g}/\sqrt{m_{g}}\right\\},$

where $m_{g}$ is the number of features in the $g^{\text{th}}$ group. This
follows from the optimality condition

(23)
$\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p^{s}}[\bm{\Phi}]\in\bigcup_{g=1}^{5}\left\\{t\sqrt{m_{g}}\partial\left\|{\cdot}\right\|_{2,g}\right\\}\bm{w}_{g}^{s}$

and computing the smallest parameter $t$ for which Eq. 23 is satisfied at
$(p^{s},\bm{w}^{s})=(p_{\text{prior}},\bm{0})$.

For this model, we use the NPDHG algorithm Eq. 18. Starting from $l=1$, we
compute the pair of solutions $(p^{(l)},\bm{w}^{(l)})$ to the primal and dual
problems Eq. 10 and Eq. 9 at hyperparameter $t=t^{(l)}$ using the previously
computed pair of solutions $(p^{(l-1)},\bm{w}^{(l-1)})$ by setting the initial
vectors $\bm{z}_{0}=\bm{w}_{0}=\bm{w}^{(l-1)}$. For the stepsize parameters,
we set $\theta_{0}=0$, $\tau_{0}=2$ and
$\sigma_{0}=1/2\left\|{\bm{\left\|{A}\right\|}}\right\|_{\mathrm{op}}^{2}$.
Finally, we compute the update $\bm{w}_{k+1}$ in Eq. 18 using the following
proximal operator formula: for every group $g\in\\{1.\dots,5\\}$, $\lambda>0$
and $\hat{\bm{w}}\in\mathbb{R}^{m}$,

$\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{\lambda\sqrt{m_{g}}\left\|{\bm{w}}\right\|_{2,g}+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2,g}^{2}\right\\}=\max{\left(0,1-\lambda\sqrt{m_{g}}/\left\|{\hat{\bm{w}}}\right\|_{2}\right)}\hat{\bm{w}}.$

Finally, we let the NPDHG algorithms run for at least 40 iterations before
checking for convergence. We stop the NPDHG algorithms when the optimality
condition Eq. 23 is satisfied within some tolerance $10^{-5}$:

$\max_{g\in\\{1,\dots,5\\}}\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k}}[\bm{\Phi}]}\right\|_{2,g}\leqslant
t(1+10^{-5}).$

####  $\ell_{\infty}$-regularized Maxent model

For this Maxent model, the smallest hyperparameter for which the solutions to
the primal and dual problems Eq. 7 and Eq. 8 are equal to the prior
distribution and zero is

$t^{(0)}=\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{\text{prior}}}[\bm{\Phi}]}\right\|_{1}.$

This follows from the optimality condition

(24)
$\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p^{s}}[\bm{\Phi}]\in
t\partial\left\|{\cdot}\right\|_{\infty}(\bm{w}^{s})$

and computing the smallest parameter $t$ for which Eq. 24 is satisfied at
$(p^{s},\bm{w}^{s})=(p_{\text{prior}},\bm{0})$.

For this model, we use the NPDHG algorithm Eq. 18. Starting from $l=1$, we
compute the pair of solutions $(p^{(l)},\bm{w}^{(l)})$ to the primal and dual
problems Eq. 11 and Eq. 12 at hyperparameter $t=t^{(l)}$ using the previously
computed pair of solutions $(p^{(l-1)},\bm{w}^{(l-1)})$ by setting the initial
vectors $\bm{z}_{0}=\bm{w}_{0}=\bm{w}^{(l-1)}$. For the stepsize parameters,
we set $\theta_{0}=0$, $\tau_{0}=2$ and
$\sigma_{0}=1/2\left\|{\bm{\left\|{A}\right\|}}\right\|_{\mathrm{op}}^{2}$.
Finally, we compute the update $\bm{w}_{k+1}$ in Eq. 18 using Moreau’s
decomposition [42, Theorem 3.2.5]: for every $\lambda>0$ and
$\hat{\bm{w}}\in\mathbb{R}^{m}$,

$\operatorname*{arg\,min}_{\bm{w}\in\mathbb{R}^{m}}\left\\{\lambda\left\|{\bm{w}}\right\|_{\infty}+\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2}^{2}\right\\}=\hat{\bm{w}}-\operatorname*{arg\,min}_{\left\|{\bm{w}}\right\|_{1}\leqslant\lambda}\frac{1}{2}\left\|{\bm{w}-\hat{\bm{w}}}\right\|_{2}^{2}.$

The second term on the right amounts to projecting $\hat{\bm{w}}$ on the
$\ell_{1}$ ball of radius $\lambda$. There are fast algorithms for doing do;
we use Algorithm 1 described in [14]. Finally, we let the NPDHG algorithms run
for at least 40 iterations before checking for convergence. We stop the NPDHG
algorithms when the optimality condition Eq. 24 is satisfied within some
tolerance $10^{-5}$:

$\left\|{\mathbb{E}_{\hat{\mathcal{D}}}[\bm{\Phi}]-\mathbb{E}_{p_{k}}[\bm{\Phi}]}\right\|_{1}\leqslant
t(1+10^{-5}).$

###  5.4 Comparison of timings

In this section, we compare the run times of our NPDHG algorithms with two
state-of-the-art optimization algorithms for solving nonsmooth Maxent models:
the forward-backward splitting algorithm (specifically, Algorithm 5 in [11];
see also [5, 20]) and the STRUCTMAXENT2 coordinate descent algorithm from [15,
62]. All numerical experiments were performed on a single core Intel(R)
Core(TM) i7-10750H CPU @ 2.60 GHz.

We initially chose the GLMNET implementation [32, 68] in MATLAB over
STRUCTMAXENT2 for the numerical comparisons, but we found that GLMNET produced
unreliable numerical solutions when compared to both the NPDHG and forward-
backward splitting algorithms. We also tried using GLMNET’s implementation in
the R language, but to no avail. We think this problem arises because, as
discussed in Section 3.2, the GLMNET algorithm approximates a Maxent model as
a logistic regression model and then invokes a coordinate descent method
tailored to logistic regression to approximate the solution to the Maxent
model. Our observations suggest this approach does not work well for our data
set. In contrast, we found that the STRUCTMAXENT2 algorithm produced correct
numerical results.

For the forward-backward splitting algorithm, the stepsize parameters were set
to $1/\left\|{\bm{A}}\right\|_{2}$ and $1$ (corresponding to $\tau$ and
$t_{0}$ in Algorithm 5 of [11]) and for computing the pair of solutions
$(p^{(l)},\bm{w}^{(l)})$ at hyperparameter $t^{(l)}$, the initial iterate was
set to $\bm{w}^{(l-1)}$. In addition, for the elastic net Maxent models, the
acceleration quantity
$q=(1-\alpha)t^{(l)}/(\left\|{\bm{A}}\right\|_{2}+(1-\alpha)t^{(l)})$ was
employed. We used the stopping criteria of the NPDHG algorithms for the
forward-backward splitting algorithms. For the coordinate descent algorithm,
we modified the STRUCTMAXENT2 algorithm from [62, Page 305-306] to make it
applicable to the elastic net penalty. For computing the pair of solutions
$(p^{(l)},\bm{w}^{(l)})$ at hyperparameter $t^{(l)}$, the initial iterate was
set to $\bm{w}^{(l-1)}$. We did not use the STRUCTMAXENT2 algorithm for the
non-overlapping group lasso or $\ell_{\infty}$-regularized Maxent models, as
it was not designed for these Maxent models. We used the stopping criteria of
the NPDHG algorithms for the STRUCTMAXENT2 algorithm.

Table 1 shows the average timings for computing the entire regularization path
of the WUMI wildfire data set using the Coordinate descent, the forward-
backward splitting and NPDHG algorithms. All timings were averaged over five
runs, and for the Forward-backward splitting and NPDHG algorithms, they
include the time required to compute all the stepsize parameters. All
algorithms were implemented in MATLAB.

| STRUCTMAXENT2 | Forward-backward splitting | NPDHG
---|---|---|---
Elastic net ($\alpha=0.95$) | 5562.19 | 4208.01 | 365.55
Elastic net ($\alpha=0.40$) | 1018.73 | 1407.22 | 113.53
Non-overlapping group lasso | N/A | 3036.38 | 278.14
$\ell_{\infty}$-regularization | N/A | 2534.65 | 289.98
Table 1: Timings results (in seconds) for fitting the Maxent models described
in Section 5.3. All timings are averaged over five runs.

The NPDHG algorithm outperformed both the forward-backward splitting and
coordinate descent algorithms by at least one order of magnitude. In
particular, we observed that the NPDHG algorithm required far fewer iterations
to achieve convergence compared to both the forward-backward splitting and
STRUCTMAXENT2 algorithms. This difference is because the stepsize parameters
for the NPDHG algorithm were much larger compared to either the forward-
backward splitting and STRUCTMAXENT2 algorithm. Indeed, the stepsize
parameters for the NPDHG and forward-backward splitting algorithms are
inversely proportional to the norms $\left\|{A}\right\|_{\mathrm{op}}$ and
$\left\|{A}\right\|_{2}$, and for the wildfire data set these were
$\left\|{\bm{A}}\right\|_{\mathrm{op}}\approx 3.30$ and
$\left\|{A}\right\|_{2}\approx 854.08$. Thus larger stepsize parameters were
permitted thanks to the Kullback–Leibler divergence term in the NPDHG
algorithm, enabling a major speedup gain.

Figure 2: Spatial probability plot for different hyperparameter values with
elastic net penalty parameter $\alpha=\\{0.95,0.40,0.15,0.05\\}$.

###  5.5 Analysis of the Maxent regularization paths and estimated fire
probabilities

Figure 3: Same as Fig. 2 but for (left) the non-overlapping group lasso with
$\alpha=1$, and (right) the $l_{\infty}$ Maxent models respectively.

As a final validation step for the NPDHG algorithm, we use the fitted Maxent
models to compute the normalized mean fire probability in each grid cell for
all calendar months between 1984 and 2020. The probabilities are visualized in
Fig. 2 and Fig. 3. In each case, we have chosen the prior distribution as a
benchmark for our spatial plots of fire probability.

The spatial fire probability for $\alpha=\\{0.95,0.40,0.15,0.05\\}$ are shown
in Fig. 2. The range of $\alpha$ values roughly corresponds to varying the
regularization from a purely $l_{1}$ norm ($\alpha=1$) to a purely $l_{2}$
norm ($\alpha=0$). For each value of $\alpha$, we also consider the evolution
of the spatial fire probability as we vary the hyperparameter $t$ along the
regularization path, or equivalently include additional features while fitting
the Maxent model to wildfire data. Broadly, we observe that for a fixed value
of $t/t_{\rm max}$, the ratio of predicted to prior fire probability decreases
in sharpness as $\alpha$ decreases. On the other hand, for a fixed $\alpha$,
decreasing $t/t_{\rm max}$ values enables the model to make sharper
distinctions between grid cells with high and low fire probability as
evidenced by the sharper contrast between the prior and predicted fire
probabilities. We also note a similar pattern in Fig. 3 for the Maxent models
with non-overlapping group lasso (corresponding roughly to the elastic net
case with $\alpha=1$) and the $l_{\infty}$-Maxent models, all which converge
to the empirical distribution quicker than any of the elastic net cases.

Figure 4: Number of non-zero coefficients along the regularization path plots for elastic net penalty parameter $\alpha=\\{0.95,0.40,0.15,0.05\\}$. The dashed vertical lines highlight the $t/t_{\rm max}$ value at which the first feature of the group indicated by inset text is selected. Feature group | $\alpha=0.95$ | $\alpha=0.05$
---|---|---
Fire weather | Tmax | Tmax | VPD
Prec | Prec | Tmin
Wind | Wind | ${\rm VPD}_{\rm max3}$
FM1000 | FM1000 | ${\rm VPD}_{\rm max7}$
${\rm Tmin}_{\rm max3}$ | ${\rm Tmin}_{\rm max3}$ | ${\rm Tmax}_{\rm max3}$
Lightning | Lightning | ${\rm Tmax}_{\rm max7}$
| ${\rm Tmin}_{\rm max7}$ | ${\rm SWE}_{\rm max}$
Topography | Slope | Slope | Southness
Southness |  |
Vegetation | Grassland | Biomass | Forest
| Shrub |
Human | Camp_dist | Camp_dist | Road_dist
Camp_num | Camp_num |
Urban | Urban |
Popdensity | Popdensity |
Antecedent | ${\rm AvgPrec}_{\rm 2mo}$ | ${\rm AvgPrec}_{\rm 2mo}$ | ${\rm AntPrec}_{\rm lag2}$
| ${\rm AvgVPD}_{\rm 3mo}$ | ${\rm AntPrec}_{\rm lag1}$ |
Table 2: List of non-zero features at the end of the regularization path for
two elastic net penalty parameters, organized by different feature groups. See
Appendix B for additional description of the selected features.

In Fig. 4, we show the cumulative number of non-zero coefficients at fixed
intervals along the regularization path for different $\alpha$ values. The
plot helps in visualizing the $t/t_{\rm max}$ values at which new features are
introduced in the elastic net Maxent model, with the dashed vertical lines
indicating the first time a feature from a new group is selected. Across all
$\alpha$ values we find that features appear in the same order, with fire
weather features being selected first, followed by topography, vegetation,
human, and antecedent features. We tabulate the final set of non-zero features
at the end of the regularization path for $\alpha=\\{0.95,0.05\\}$ across
various groups in Table 2. These selected features are in good agreement with
physical models of wildfire occurrence [66, 61] as well as previous
statistical analyses of wildfire drivers [67, 12, 87, 10].

##  6 Discussion

In this paper, we have introduced novel first-order NPDHG algorithms that
overcome the shortcomings of state-of-the-art algorithms for training large-
scale, non-smooth Maxent models. The crucial ingredient common to our
algorithms is the Kullback–Leibler divergence. Using it over the classical
proximal operator makes it possible to train large-scale and non-smooth Maxent
models much more efficiently than the state of the arts. In particular, all
stepsize parameters and iterations in our algorithms can be calculated on the
order $O(mn)$ operations, improving on the complexity bound of
$O(\min(m^{2}n,mn^{2}))$ operations for computing the optimal stepsize
parameters of classical first-order optimization methods, such as the linear
PDHG or forward-backward splitting methods. As a consequence, for a given
tolerance level $\epsilon>0$, our algorithms provably compute solutions using
on the order of $O(mn/\sqrt{\epsilon})$ or $O(mn/\log(1/\epsilon))$
operations, the order depending on the smoothness of the Maxent model and
which are optimal with respect to the Nesterov class of optimal first-order
methods [64]. Moreover, the computational bottleneck consists of matrix-vector
multiplications, which can be trivially parallelized, and so our algorithms
exhibit scalable parallelism.

Finally, we have shown that the strong convexity of the Kullback–Leibler
divergence with respect to the $\ell_{1}$ norm allows for significantly larger
stepsize parameters, thereby speeding up the convergence rate of our
algorithms. This was, in particular, observed in Section 5, when we applied
our algorithms to fit the WUMI wildfire data set [45] on several non-smooth
Maxent models to estimate probabilities of fire occurrences as a function of
ecological features. Our algorithms outperformed the state-of-the-art forward-
backward splitting and coordinate descent STRUCTMAXENT2 algorithms by at least
one order of magnitude. They also yielded results that are in good agreement
with physical models of wildfire occurrence [66, 61] as well as previous
statistical analyses of wildfire drivers [67, 12, 87, 10]. Future work will
also explore the scalability of our algorithms for modeling daily scale
wildfire probability [48].

We expect our algorithms to provide efficient methods for solving non-smooth
Maxent models that arise in large-scale machine learning applications beyond
the wildfire application explored in this paper. It would be interesting and
impactful to extend our algorithms to continuous regularized Maxent models.
Interesting, because the continuous version of the Maxent problem is posed
over a non-reflexive space, which makes this problem more technically
challenging. Impactful, because such an algorithm would enable a much broader
class of continuous probability distributions to be used in Maxent modelling.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Funding Information

GPL and JD are supported by ONR N00014-22-1-2667. JD is also supported by DOE-
MMICS SEA-CROGS DE-SC0023191. JB is supported by the Zegar Family Foundation
and the National Science Foundation LEAP Science and Technology Center Award
#2019625.

##  Appendix A Derivation of update (17)

We divide the derivation into two parts, first deriving an expression for the
update $p_{k+1}$ in terms of the previous iterate $p_{k}$ and then deriving
the explicit update for $p_{k+1}$ as in the third line of (16).

Part 1. First, note the minimization problem

$\min_{p\in\Delta_{n}}\left\\{\tau_{k}D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})-\tau_{k}\left\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\mathbb{E}_{p}[\bm{\Phi}]\right\rangle+D_{\mathrm{KL}}(p\;\|\;p_{k})\right\\}$

has a unique global minimum, denoted by $p_{k+1}$, in the relative interior of
the unit simplex $\Delta_{n}$ for every $k\in\mathbb{N}$. See [18, Proposition
4.1] for details and note that Proposition 4.1 applies because the
Kullback–Leibler divergence is strongly convex with respect to the $\ell_{1}$
norm in its first argument.

Having established that the minimum $p_{k+1}$ exists and is unique, we
introduce the Karush-Kuhn-Tucker multipliers $\xi\in\mathbb{R}$ and
$\bm{\mu}\in\mathbb{R}^{n}$ to write the minimization problem for $p_{k+1}$ in
terms of an unconstrained optimization problem:

$\min_{\begin{subarray}{c}p\in\mathbb{R}^{n}\\\ \xi\in\mathbb{R}\\\
\mu\in\mathbb{R}^{n}\end{subarray}}\left\\{\tau_{k}D_{\mathrm{KL}}(p\;\|\;p_{\text{prior}})-\tau_{k}\left\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\mathbb{E}_{p}[\bm{\Phi}]\right\rangle+D_{\mathrm{KL}}(p\;\|\;p_{k})+\xi\left(1-\sum_{j=1}^{n}p(j)\right)-\langle\bm{\mu},p\rangle\right\\}.$

We can invoke the Karush–Kuhn–Tucker Theorem and use the linearity constraint
qualifications to find that $p_{k+1}$ satisfies the first-order optimality
condition

(25) $\displaystyle\tau_{k}\log(p_{k+1}(j)/p_{\text{prior}}(j))+\tau_{k}$
$\displaystyle-\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle$
$\displaystyle+\log(p_{k+1}(j)/p_{k}(j))+1-\xi-[\bm{\mu}]_{j}=0$

for each $j\in\\{1,\dots,n\\}$ [8]. Moreover, the complementary slackness
condition $\langle\bm{\mu},p_{k+1}\rangle=0$ with the fact that $p_{k+1}>0$
for every $j\in\\{1,\dots,m\\}$ implies $\bm{\mu}=\bm{0}$.

Next, we use the first-order optimality condition to compute $p_{k+1}(j)$
explicitly. Rearranging (25) in terms of $\log(p_{k+1}(j))$, using
$\bm{\mu}=0$, and taking the exponential yield

$p_{k+1}(j)=\left(e^{(\xi-1)/(1+\tau_{k})-1}\right)\left((p_{\text{prior}}(j))^{\tau_{k}}p_{k}(j)e^{\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}$

Since $\sum_{j=1}^{n}p_{k+1}(j)=1$, we find that $\xi$ satisfies the relation

$e^{\tau_{k}\xi/(1+\tau_{k})-1}=\frac{1}{\sum_{j=1}^{n}\left((p_{\text{prior}}(j))^{\tau_{k}}p_{k}(j)e^{\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}}.$

Hence

$p_{k+1}(j)=\frac{\left((p_{\text{prior}}(j))^{\tau_{k}}p_{k}(j)e^{\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}}{\sum_{j=1}^{n}\left((p_{\text{prior}})^{\tau_{k}}(j)p_{k}(j)e^{\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}}$

for each $j\in\\{1,\dots,n\\}$. This yields an explicit update for $p_{k+1}$
in terms of the previous iterate $p_{k}$.

Part 2. We now induct on $k$ to prove that the explicit expression for the
update $p_{k+1}$ in the third line of (16) is correct. For $k=0$, we have

$p_{0}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}}$

and so

$\displaystyle\left(\left(p_{\text{prior}}(j)\right)^{\tau_{0}}p_{0}(j)e^{\tau_{0}\langle\bm{w}_{0}+\theta_{0}(\bm{w}_{0}-\bm{w}_{-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{0})}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad=\frac{p_{\text{prior}}(j)e^{\langle(\bm{z}_{0}+\tau_{0}(\bm{w}_{0}+\theta_{0}(\bm{w}_{0}-\bm{w}_{-1})))/(1+\tau_{0}),\bm{\Phi}(j)\rangle}}{\left(\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{0})}}.$

Writing

$\bm{z}_{1}=(\bm{z}_{0}+\tau_{0}(\bm{w}_{0}+\theta_{0}(\bm{w}_{0}-\bm{w}_{-1})))/(1+\tau_{0}),$

the term inside the exponential on the numerator simplifies to
$\langle\bm{z}_{1},\bm{\Phi}(j)\rangle$. Hence

$\left(\left(p_{\text{prior}}(j)\right)^{\tau_{0}}p_{0}(j)e^{\tau_{0}\langle\bm{w}_{0}+\theta_{0}(\bm{w}_{0}-\bm{w}_{-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{0})}=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{1},\bm{\Phi}(j)\rangle}}{\left(\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{0},\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{0})}}.$

It follows on substitution that

$p_{1}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{1},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{1},\bm{\Phi}(j)\rangle}}$

for every $j\in\\{1,\dots,n\\}$.

For arbitrary $k\in\mathbb{N}$, we use the induction hypothesis

$p_{k}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k},\bm{\Phi}(j)\rangle}}$

to find

$\displaystyle\left(\left(p_{\text{prior}}(j)\right)^{\tau_{k}}p_{k}(j)e^{\tau_{k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad=\frac{p_{\text{prior}}(j)e^{\langle(\bm{z}_{k}+\tau_{k}(\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1})))/(1+\tau_{k}),\bm{\Phi}(j)\rangle}}{\left(\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k},\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}}.$

Writing

$\bm{z}_{k+1}=(\bm{z}_{k}+\tau_{k}(\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1})))/(1+\tau_{k}),$

the term inside the exponential on the numerator simplifies to
$\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle$. Hence

$\left(\left(p_{\text{prior}}(j)\right)^{\tau_{k}}p_{k}(j)e^{\tau_{0k}\langle\bm{w}_{k}+\theta_{k}(\bm{w}_{k}-\bm{w}_{k-1}),\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}{\left(\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k},\bm{\Phi}(j)\rangle}\right)^{1/(1+\tau_{k})}}.$

It follows on substitution that

$p_{k+1}(j)=\frac{p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}{\sum_{j=1}^{n}p_{\text{prior}}(j)e^{\langle\bm{z}_{k+1},\bm{\Phi}(j)\rangle}}$

for every $j\in\\{1,\dots,n\\}$, which proves the desired result.

##  Appendix B Summary of features extracted from the WUMI wildfire data set

We provide here a summary of all the features used in our statistical
analysis. The features are aggregated to a 12 km spatial resolution during
data preprocessing. Considering each feature’s $M$ antecedent month average
and maximum $X$-day running average components as distinct features, the total
number of features considered in our analysis adds up to 35. We organize the
features into five groups with each group reflecting a major environmental
driver of fire occurrences. These groups are also an important ingredient of
the non-overlapping group lasso Maxent model.

Table 3: Summary table for all input features organized by group, identifier, description, spatial resolution of raw data, and source of data. Feature group | Identifier |  Description | Resolution | Source
---|---|---|---|---
Fire weather | VPD |  Mean vapor pressure deficit | 5 km | Climgrid [83], PRISM [17]
${\rm VPD}^{{\rm max}X}$ |  Maximum $X$-day running average of VPD; $X\in\\{3,7\\}$ | 9 km | UCLA-ERA5 [73]
Tmax |  Daily maximum temperature | 5 km | Climgrid
${\rm Tmax}^{{\rm max}X}$ |  Maximum $X$-day running average of Tmax | 9 km | UCLA-ERA5
Tmin |  Daily minimum temperature | 5 km | Climgrid
${\rm Tmin}^{{\rm max}X}$ |  Maximum $X$-day running average of Tmin | 9 km | UCLA-ERA5
Prec |  Precipitation total | 5 km | Climgrid
${\rm SWE}_{\rm mean}$ |  Mean snow water equivalent | 500 m | NSIDC [90]
${\rm SWE}_{\rm max}$ |  Daily maximum snow water equivalent | 500 m | NSIDC
FM1000 |  1000-hour dead fuel moisture | 4 km | gridMET [1]
${\rm FFWI}^{{\rm max}7}$ |  Maximum 7-day running average of Fosberg Fire Weather Index | 9 km | UCLA-ERA5
Wind |  Mean wind speed | 9 km | UCLA-ERA5
Lightning |  Lightning strike density | 500 m | NLDN [84, 65]
Antecedent | ${\rm AntVPD}_{M{\rm mon}}$ |  Average VPD in $M$ antecedent months; $M\in\\{2,3\\}$ | 5 km | Climgrid
${\rm AntPrec}_{M{\rm mon}}$ |  Average precipitation total in $M$ antecedent months; $M\in\\{2,4\\}$ | 5 km | Climgrid
${\rm AntPrec}_{\rm lag1}$ |  Mean annual precipitation in lag year 1 | 5 km | Climgrid
${\rm AntPrec}_{\rm lag2}$ |  Mean annual precipitation in lag year 2 | 5 km | Climgrid
${\rm AvgSWE}_{3{\rm mon}}$ |  Average snow water equivalent in 3 antecedent months | 500 m | NSIDC
Vegetation | Forest |  Fraction of forest landcover | 30 m | NLCD
Grassland |  Fraction of grassland cover | 30 m | NLCD
Shrubland |  Fraction of shrubland cover | 30 m | NLCD
Biomass |  Aboveground biomass map | 300 m | Ref. [78]
Human | Urban |  Fraction of land covered by urban areas | 30 m | NLCD
Camp_num |  Mean number of camp grounds | 1km | Open source
| Camp_dist |  Mean distance from nearest camp ground | 1km | Open source
| Road_dist |  Mean distance from nearest highway | 1km | Open source
| Popdensity |  Mean population density | 1km | SILVIS [72]
| Housedensity |  Mean housing density | 1km | SILVIS
Topography | Slope |  Mean slope | 1m | USGS
Southness |  Mean south-facing degree of slope | 1m | USGS

## References

  * [1] J. T. Abatzoglou, Development of gridded surface meteorological data for ecological applications and modelling, International Journal of Climatology, 33 (2013), pp. 121–131, https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.3413.
  * [2] Y. Altun and A. Smola, Unifying divergence minimization and statistical inference via convex duality, in International Conference on Computational Learning Theory, Springer, 2006, pp. 139–153.
  * [3] G. Andrew and J. Gao, Scalable training of l 1-regularized log-linear models, in Proceedings of the 24th international conference on Machine learning, 2007, pp. 33–40.
  * [4] A. Beck and M. Teboulle, Mirror descent and nonlinear projected subgradient methods for convex optimization, Operations Research Letters, 31 (2003), pp. 167–175.
  * [5] A. Beck and M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM journal on imaging sciences, 2 (2009), pp. 183–202.
  * [6] P. Bereziński, B. Jasiul, and M. Szpyrka, An entropy-based network anomaly detection method, Entropy, 17 (2015), pp. 2367–2408.
  * [7] A. Berger, S. A. Della Pietra, and V. J. Della Pietra, A maximum entropy approach to natural language processing, Computational linguistics, 22 (1996), pp. 39–71.
  * [8] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization, Cambridge university press, 2004.
  * [9] J. Buch, A. P. Williams, C. Juang, W. D. Hansen, and P. Gentine, SMLFire1.0: a stochastic machine learning (SML) model for wildfire activity in the western United States, 2022, https://doi.org/10.5281/zenodo.7277980, https://doi.org/10.5281/zenodo.7277980.
  * [10] J. Buch, A. P. Williams, C. S. Juang, W. D. Hansen, and P. Gentine, Smlfire1.0: a stochastic machine learning (sml) model for wildfire activity in the western united states, EGUsphere, 2022 (2022), pp. 1–39, https://doi.org/10.5194/egusphere-2022-1148, https://egusphere.copernicus.org/preprints/2022/egusphere-2022-1148/.
  * [11] A. Chambolle and T. Bock, An introduction to continuous optimization for imaging, Acta Numer., 25 (2016), pp. 161–319.
  * [12] B. Chen, Y. Jin, E. Scaduto, M. A. Moritz, M. L. Goulden, and J. T. Randerson, Climate, fuel, and land use shaped the spatial pattern of wildfire in California’s Sierra Nevada, Journal of Geophysical Research: Biogeosciences, 126 (2021).
  * [13] S. F. Chen and R. Rosenfeld, A survey of smoothing techniques for me models, IEEE transactions on Speech and Audio Processing, 8 (2000), pp. 37–50.
  * [14] L. Condat, Fast projection onto the simplex and the l 1 ball, Mathematical Programming, 158 (2016), pp. 575–585.
  * [15] C. Cortes, V. Kuznetsov, M. Mohri, and U. Syed, Structural maxent models, in Proceedings of the 32nd International Conference on International Conference on Machine Learning - Volume 37, ICML’15, JMLR.org, 2015, p. 391–399.
  * [16] I. Csiszár, Information-type measures of difference of probability distributions and indirect observation, studia scientiarum Mathematicarum Hungarica, 2 (1967), pp. 229–318.
  * [17] C. Daly, W. Gibson, M. Doggett, J. Smith, and G. Taylor, Up-to-date monthly climate maps for the conterminous United States, Proc., 14th AMS Conf. on Applied Climatology, 84th AMS Annual Meeting Combined Preprints, (2004).
  * [18] J. Darbon and G. P. Langlois, Accelerated nonlinear primal-dual hybrid gradient algorithms with applications to machine learning, 2021, https://arxiv.org/abs/2109.12222.
  * [19] J. N. Darroch and D. Ratcliff, Generalized iterative scaling for log-linear models, The annals of mathematical statistics, (1972), pp. 1470–1480.
  * [20] I. Daubechies, M. Defrise, and C. De Mol, An iterative thresholding algorithm for linear inverse problems with a sparsity constraint, Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 57 (2004), pp. 1413–1457.
  * [21] S. Della Pietra, V. Della Pietra, and J. Lafferty, Inducing features of random fields, IEEE transactions on pattern analysis and machine intelligence, 19 (1997), pp. 380–393.
  * [22] Y. Demchenko, P. Grosso, C. De Laat, and P. Membrey, Addressing big data issues in scientific data infrastructure, in 2013 International Conference on Collaboration Technologies and Systems (CTS), IEEE, 2013, pp. 48–55.
  * [23] M. Dewaskar, C. Tosh, J. Knoblauch, and D. B. Dunson, Robustifying likelihoods by optimistically re-weighting data, arXiv preprint arXiv:2303.10525, (2023).
  * [24] M. Dudík, S. J. Phillips, and R. E. Schapire, Performance guarantees for regularized maximum entropy density estimation, in Learning Theory, J. Shawe-Taylor and Y. Singer, eds., Berlin, Heidelberg, 2004, Springer Berlin Heidelberg, pp. 472–486.
  * [25] M. Dudík, S. J. Phillips, and R. E. Schapire, Maximum entropy density estimation with generalized regularization and an application to species distribution modeling, Journal of Machine Learning Research, (2007).
  * [26] M. Dudík and R. E. Schapire, Maximum entropy distribution estimation with generalized regularization, in International Conference on Computational Learning Theory, Springer, 2006, pp. 123–138.
  * [27] M. K. Dyderski, S. Paź, L. E. Frelich, and A. M. Jagodziński, How much does climate change threaten european forest tree species distributions?, Global Change Biology, 24 (2018), pp. 1150–1163.
  * [28] I. Ekeland and R. Temam, Convex analysis and variational problems, SIAM, 1999.
  * [29] M. El Guide, K. Jbilou, C. Koukouvinos, and A. Lappa, Comparative study of l 1 regularized logistic regression methods for variable selection, Communications in Statistics-Simulation and Computation, (2020), pp. 1–16.
  * [30] J. Elith, S. J. Phillips, T. Hastie, M. Dudík, Y. E. Chee, and C. J. Yates, A statistical explanation of maxent for ecologists, Diversity and distributions, 17 (2011), pp. 43–57.
  * [31] M. A. Figueiredo and R. D. Nowak, Wavelet-based image estimation: an empirical Bayes approach using Jeffrey’s noninformative prior, IEEE Transactions on Image Processing, 10 (2001), pp. 1322–1331.
  * [32] W. Fithian and T. Hastie, Finite-sample equivalence in statistical models for presence-only data, The annals of applied statistics, 7 (2013), p. 1917.
  * [33] S. Foucart and H. Rauhut, Sparse Solutions of Underdetermined Systems, Springer New York, New York, NY, 2013, pp. 41–59, https://doi.org/10.1007/978-0-8176-4948-7_2.
  * [34] J. Friedman, T. Hastie, H. Höfling, R. Tibshirani, et al., Pathwise coordinate optimization, The annals of applied statistics, 1 (2007), pp. 302–332.
  * [35] J. Friedman, T. Hastie, and R. Tibshirani, Regularization paths for generalized linear models via coordinate descent, Journal of statistical software, 33 (2010), p. 1.
  * [36] A. Gelman, J. Carlin, H. Stern, D. Dunson, A. Vehtari, and D. Rubin, Bayesian Data Analysis, Third Edition, Chapman & Hall/CRC Texts in Statistical Science, Taylor & Francis, 2013.
  * [37] E. Granot-Atedgi, G. Tkačik, R. Segev, and E. Schneidman, Stimulus-dependent maximum entropy models of neural population codes, PLoS computational biology, 9 (2013), p. e1002922.
  * [38] Y. Gu, A. McCallum, and D. Towsley, Detecting anomalies in network traffic using maximum entropy estimation, in Proceedings of the 5th ACM SIGCOMM conference on Internet Measurement, 2005, pp. 32–32.
  * [39] T. Hastie, J. Qian, and K. Tay, An introduction to glmnet (2021).  Available at ”https://glmnet.stanford.edu/articles/glmnet.html”. Accessed on 17 February 2021., 2021.
  * [40] T. Hastie, R. Tibshirani, and J. Friedman, The elements of statistical learning, Springer Series in Statistics, Springer, New York, second ed., 2009, https://doi.org/10.1007/978-0-387-84858-7.  Data mining, inference, and prediction.
  * [41] B. Hayes and C. Wilson, A maximum entropy model of phonotactics and phonotactic learning, Linguistic inquiry, 39 (2008), pp. 379–440.
  * [42] J.-B. Hiriart-Urruty and C. Lemaréchal, Convex analysis and minimization algorithms I: Fundamentals, vol. 305 of Grundlehren Text Editions, Springer Science & Business Media, 1993.
  * [43] E. T. Jaynes, Information theory and statistical mechanics, Physical review, 106 (1957), p. 620.
  * [44] E. T. Jaynes, Information theory and statistical mechanics. II, Physical review, 108 (1957), p. 171.
  * [45] C. S. Juang, A. P. Williams, J. T. Abatzoglou, J. K. Balch, M. D. Hurteau, and M. A. Moritz, Rapid growth of large forest fires drives the exponential response of annual forest-fire area to aridity in the western united states, Geophysical Research Letters, 49 (2022), https://doi.org/https://doi.org/10.1029/2021GL097131, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021GL097131.
  * [46] C. E. Kalinski, Building Better Species Distribution Models with Machine Learning: Assessing the Role of Covariate Scale and Tuning in Maxent Models, PhD thesis, University of Southern California, 2019.
  * [47] J. N. Kapur, Maximum-entropy models in science and engineering, John Wiley & Sons, 1989.
  * [48] T. Keeping, S. P. Harrison, and I. C. Prentice, Modelling the daily probability of wildfire occurrence in the contiguous United States, Environmental Research Letters, 19 (2024), p. 024036, https://dx.doi.org/10.1088/1748-9326/ad21b0.
  * [49] J. H. Kemperman, On the optimum rate of transmitting information, in Probability and information theory, Springer, 1969, pp. 126–169.
  * [50] S. Kullback, A lower bound for discrimination information in terms of variation (corresp.), IEEE transactions on Information Theory, 13 (1967), pp. 126–127.
  * [51] T. M. Kwon, Tmc traffic data automation for mndot’s traffic monitoring program, tech. report, University of Minnesota, Twin Cities, 2004\.
  * [52] J. D. Lee, B. Recht, N. Srebro, J. Tropp, and R. R. Salakhutdinov, Practical large-scale optimization for max-norm regularization, Advances in neural information processing systems, 23 (2010).
  * [53] S.-I. Lee, H. Lee, P. Abbeel, and A. Y. Ng, Efficient l~ 1 regularized logistic regression, in Aaai, vol. 6, 2006, pp. 401–408.
  * [54] P.-L. Lions and B. Mercier, Splitting algorithms for the sum of two nonlinear operators, SIAM Journal on Numerical Analysis, 16 (1979), pp. 964–979.
  * [55] J. S. Littell, D. McKenzie, D. L. Peterson, and A. L. Westerling, Climate and wildfire area burned in western u.s. ecoprovinces, 1916–2003, Ecological Applications, 19 (2009), pp. 1003–1021, https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/07-1183.1.
  * [56] A. L’heureux, K. Grolinger, H. F. Elyamany, and M. A. Capretz, Machine learning with big data: Challenges and approaches, Ieee Access, 5 (2017), pp. 7776–7797.
  * [57] R. Malouf, A comparison of algorithms for maximum entropy parameter estimation, in COLING-02: The 6th Conference on Natural Language Learning 2002 (CoNLL-2002), 2002.
  * [58] G. Mann, R. McDonald, M. Mohri, N. Silberman, and D. Walker IV, Efficient large-scale distributed training of conditional maximum entropy models, tech. report, Google Research, 2009.
  * [59] R. McElreath, Statistical rethinking: A Bayesian course with examples in R and Stan, Chapman and Hall/CRC, 2018.
  * [60] C. Merow, M. J. Smith, and J. A. Silander Jr, A practical guide to maxent for modeling species’ distributions: what it does, and why inputs and settings matter, Ecography, 36 (2013), pp. 1058–1069.
  * [61] K. Mezuman, K. Tsigaridis, G. Faluvegi, and S. E. Bauer, The interactive global fire module pyre (v1.0), Geoscientific Model Development, 13 (2020), pp. 3091–3118, https://gmd.copernicus.org/articles/13/3091/2020/.
  * [62] M. Mohri, A. Rostamizadeh, and A. Talwalkar, Foundations of Machine Learning, MIT press, second ed., 2018.
  * [63] L. A. Muttaqin, S. H. Murti, and B. Susilo, Maxent (maximum entropy) model for predicting prehistoric cave sites in karst area of gunung sewu, gunung kidul, yogyakarta, in Sixth Geoinformation Science Symposium, vol. 11311, International Society for Optics and Photonics, 2019, p. 113110B.
  * [64] Y. Nesterov, Introductory lectures on convex optimization, vol. 87 of Applied Optimization, Kluwer Academic Publishers, Boston, MA, 2004.  A basic course.
  * [65] R. E. Orville and G. R. Huffines, Cloud-to-Ground lightning in the United States: NLDN results in the first decade, 1989–98, Monthly Weather Review, 129 (2001), pp. 1179 – 1193.
  * [66] M.-A. Parisien, D. A. Dawe, C. Miller, C. A. Stockdale, and O. B. Armitage, Applications of simulation-based burn probability modelling: a review, International Journal of Wildland Fire, 28 (2019), p. 913, http://dx.doi.org/10.1071/WF19069.
  * [67] M.-A. Parisien and M. A. Moritz, Environmental controls on the distribution of wildfire at multiple spatial scales, Ecological Monographs, 79 (2009), pp. 127–154.
  * [68] S. J. Phillips, R. P. Anderson, M. Dudík, R. E. Schapire, and M. E. Blair, Opening the black box: An open-source release of maxent, Ecography, 40 (2017), pp. 887–893.
  * [69] S. J. Phillips, R. P. Anderson, and R. E. Schapire, Maximum entropy modeling of species geographic distributions, Ecological modelling, 190 (2006), pp. 231–259.
  * [70] S. J. Phillips, M. Dudík, and R. E. Schapire, A maximum entropy approach to species distribution modeling, in Proceedings of the twenty-first international conference on Machine learning, 2004, p. 83.
  * [71] M. S. Pinsker, Information and information stability of random variables and processes, Holden-Day, 1964.
  * [72] V. C. Radeloff, R. B. Hammer, S. I. Stewart, J. S. Fried, S. S. Holcomb, and J. F. McKeefry, The Wildland-Urban Interface in the United States, Ecological Applications, 15 (2005), pp. 799–805, https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/04-1413.
  * [73] S. Rahimi, W. Krantz, Y. Lin, B. Bass, N. Goldenson, A. Hall, Z. Jebo, and J. Norris, Evaluation of a reanalysis-driven configuration of wrf4 over the western united states from 1980-2020., Journal of Geophysical Research: Atmospheres, 127 (2022), https://doi.org/10.1029/2021JD035699.
  * [74] A. Ratnaparkhi, Maximum Entropy Models for Natural Language Processing, Springer US, Boston, MA, 2017, pp. 800–805, https://doi.org/10.1007/978-1-4899-7687-1_525.
  * [75] C. Savin and G. Tkačik, Maximum entropy models as a tool for building precise neural controls, Current opinion in neurobiology, 46 (2017), pp. 120–126.
  * [76] J. L. Schnase and M. L. Carroll, Automatic variable selection in ecological niche modeling: A case study using cassin’s sparrow (peucaea cassinii), PloS one, 17 (2022), p. e0257502.
  * [77] J. L. Schnase, M. L. Carroll, R. L. Gill, G. S. Tamkin, J. Li, S. L. Strong, T. P. Maxwell, M. E. Aronne, and C. S. Spradlin, Toward a monte carlo approach to selecting climate variables in maxent, PloS one, 16 (2021), p. e0237208.
  * [78] S. A. Spawn, C. C. Sullivan, T. J. Lark, and H. K. Gibbs, Harmonized global maps of above and belowground biomass carbon density in the year 2010, Scientific Data, 7 (2020), p. 112, https://doi.org/10.1038/s41597-020-0444-4, https://www.nature.com/articles/s41597-020-0444-4.
  * [79] J. K. Tay, B. Narasimhan, and T. Hastie, Elastic net regularization paths for all generalized linear models, Journal of statistical software, 106 (2023).
  * [80] N. C. Thompson, K. Greenewald, K. Lee, and G. F. Manso, Deep learning’s diminishing returns: The cost of improvement is becoming unsustainable, IEEE Spectrum, 58 (2021), pp. 50–55.
  * [81] G. Tkačik, O. Marre, T. Mora, D. Amodei, M. J. Berry II, and W. Bialek, The simplest maximum entropy model for collective behavior in a neural network, Journal of Statistical Mechanics: Theory and Experiment, 2013 (2013), p. P03011.
  * [82] J. Tsujii and J. Kazama, Evaluation and extension of maximum entropy models with inequality constraints, in Proceedings of the 2003 conference on Empirical methods in natural language processing, 2003, pp. 137–144.
  * [83] R. Vose, S. Applequist, M. Squires, I. Durre, M. Menne, C. Williams, C. Fenimore, K. Gleason, and D. Arndt, Improved historical Temperature and Precipitation time series for U.S. climate divisions, Journal of Applied Meteorology and Climatology, 53 (2014), pp. 1232–1251.
  * [84] R. S. Wacker and R. E. Orville, Changes in measured lightning flash count and return stroke peak current after the 1994 U.S. National Lightning Detection Network upgrade: 1. Observations, Journal of Geophysical Research: Atmospheres, 104 (1999), pp. 2151–2157, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/1998JD200060.
  * [85] M. J. Wainwright and M. I. Jordan, Graphical models, exponential families, and variational inference, Now Publishers Inc, 2008.
  * [86] H.-F. Yu, F.-L. Huang, and C.-J. Lin, Dual coordinate descent methods for logistic regression and maximum entropy models, Machine Learning, 85 (2011), pp. 41–75.
  * [87] Y. Yu, S. Saatchi, B. Walters, S. Ganguly, S. Li, S. Hagen, L. Melendy, R. Nemani, G. Domke, and C. Woodall, Carbon pools across conus using the maxent model, 2005, 2010, 2015, 2016, and 2017, 2021, https://doi.org/10.3334/ORNLDAAC/1752, https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1752.
  * [88] G.-X. Yuan, K.-W. Chang, C.-J. Hsieh, and C.-J. Lin, A comparison of optimization methods and software for large-scale l1-regularized linear classification, The Journal of Machine Learning Research, 11 (2010), pp. 3183–3234.
  * [89] M. Zanon, G. Zambonin, G. A. Susto, and S. McLoone, Sparse logistic regression: Comparison of regularization and bayesian implementations, Algorithms, 13 (2020), p. 137.
  * [90] X. Zeng, P. Broxton, and N. Dawson, Snowpack change from 1982 to 2016 over conterminous United States, Geophysical Research Letters, 45 (2018), pp. 12,940–12,947, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018GL079621.
  * [91] H. Zou and T. Hastie, Regularization and variable selection via the elastic net, Journal of the royal statistical society: series B (statistical methodology), 67 (2005), pp. 301–320.