Efficient First-Order Algorithms for Large-Scale, Non-Smooth Maximum Entropy Models with Application to Wildfire Science

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 [38], 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 art. In particular, we show that, for a MaxEnt model with a discrete probability distribution of n elements built from samples each containing m features, the stepsize parameter estimation and iterations in our algorithms all scale on the order of $O\left(mn\right)$ operations. This significantly improves on the known complexity bound of $O\left(\min (m^{2}n,m{n}^2)\right)$ 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 $ϵ>0$, our algorithms provably compute solutions using on the order of $O(mn/\sqrt{ϵ})$ or $O(mn/\log (1/ϵ\left)\right)$ 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 [39]. 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 the 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 [40,41]. 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 [42,43] as well as previous statistical analyses of wildfire drivers [29,30,31,44].

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 data sets, 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 the probabilities of fire occurrences as a function of ecological features. Finally, we review our results and outline directions for future work in Section 6.

Suppose we receive l independent, identically distributed samples $\{v_1,\dots ,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_{\mathrm{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 $\mathbf{\Phi }:\mathcal{I}\to {\mathbb{R}}^m$. Then, how do we estimate the unknown distribution $\mathcal{D}$ from the prior distribution $p_{\mathrm{prior}}$, the samples $\{v_1,\dots ,v_l\}$, and the feature map $\mathbf{\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_{\mathrm{prior}}$ while matching the features $\{\mathbf{\Phi }\left(v_1\right),\dots ,\mathbf{\Phi }\left(v_l\right)\}$ exactly or as closely as possible, in some suitable sense. Specifically, let ${\Delta }_n$ denote the n-dimensional probability simplex and suppose $p_{\mathrm{prior}}\in \mathrm{int}{\Delta }_n$ (i.e., $p_{\mathrm{prior}}\left(j\right)>0$ for every $j\in \{1,\dots ,n\}$). We measure the closeness of a probability distribution $p\in {\Delta }_n$ to the prior probability $p_{\mathrm{prior}}$ using the Kullback–Leibler divergence $D_{\mathrm{KL}}:{\Delta }_n\times \mathrm{int}{\Delta }_n\to \mathbb{R}$: Next, let $\widehat{\mathcal{D}}$ denote the empirical distribution induced by samples $\{v_1,\dots ,v_l\}$: In addition, let ${\mathbb{E}}_{\widehat{\mathcal{D}}}[\mathbf{\Phi }]$ denote the empirical average of the features induced by the samples and let ${\mathbb{E}}_p[\mathbf{\Phi }]$ denote the model average induced by a probability distribution $p\in {\Delta }_n$: We measure how the averages ${\mathbb{E}}_{\widehat{\mathcal{D}}}[\mathbf{\Phi }]$ and ${\mathbb{E}}_p[\mathbf{\Phi }]$ match using a function $H^{*}:{\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 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^{*}\left(\mathbf{0}\right)$ = 0. Under the standing assumptions above on the potential function, which are fairly general, the probability distribution $\widehat{\mathcal{D}}$ is a feasible point of (3). Moreover, the MaxEnt estimation problem (3) admits a unique global solution (see ([46], Page 35, Proposition 1.2); uniqueness follows from the strong convexity of the Kullback–Leibler divergence with respect to the ${\ell }_1$ norm [47,48,49,50,51]).

The generalized MaxEnt problem (3) admits a dual problem corresponding to regularized maximum a posteriori estimation [52]. To derive the dual problem, we first write the second term on the right hand side of (3) in terms of its convex conjugate, where we abuse notation and write the convex conjugate of $H^{*}$ as H. This lets us express problem (3) in saddle-point form: Due to our assumptions on the potential function $H^{*}$, we can swap the infimum and supremum ([46], Statement (4.1) on page 61). Using the convex conjugate formula we obtain the dual problem of (3)

The dual problem (5) is a regularized maximum likelihood estimation problem over the family of Gibbs distributions [41,52]. It has at least one global solution ([46], Proposition (1.2) on page 35). Moreover, if $(p^s,{\mathit{w}}^s)$ denotes a pair of solutions to the MaxEnt estimation problem (3) and its dual (5), then this pair satisfies the optimality conditions for every $j\in \{1,\dots ,n\}$ where $\partial H\left({\mathit{w}}^s\right)$ denotes the subdifferential of H at ${\mathit{w}}^s$.

MaxEnt models differ in the choice of the prior distribution $p_{\mathrm{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 This forces the model average ${\mathbb{E}}_p[\mathbf{\Phi }]$ to be equal to the empirical average ${\mathbb{E}}_{\widehat{\mathcal{D}}}[\mathbf{\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 [20,26,27,53]. These models will also be used for our numerical experiments and results in Section 5.

Estimating a probability distribution from the MaxEnt model (3) can be computationally prohibitive for big data sets. To illustrate this point, suppose the MaxEnt model (3) with hyperparameter $t⩾0$ has a global solution $p^s\left(t\right)$. Let $p^{ϵ}\left(t\right)\in {\Delta }_n$ with $ϵ>0$ denote an $ϵ$-approximate solution to the global solution $p^s\left(t\right)$, that is, the objective function f in (3) satisfies If the potential function $H^{*}$ is smooth (equivalently, its conjugate H is strongly convex), the best achievable rate of convergence for computing $p^{ϵ}\left(t\right)$ with $t>0$ in the Nesterov class of optimal first-order methods is linear $O\left(\log \right(1/ϵ\left)\right)$ in the number of operations [39]. If H is not smooth, then the optimal convergence rate is sublinear $O(1/\sqrt{ϵ})$ 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 $\mathit{A}:{\Delta }_n\to {\mathbb{R}}^m$ defined by Unfortunately, computing the largest singular value of the linear operator $\mathit{A}$ in (13) accurately is often computationally expensive for large matrices due to its prohibitive cubic cost of $O\left(\min (m^{2}n,m{n}^2)\right)$ operations [60]. 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 $\mathit{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 (3) difficult and laborious. Even worse, in some applications, the appropriate value of the hyperparameter t in (3) is difficult to guess and must be selected by repeatedly solving (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.

State-of-the-art methods for computing solutions to large-scale, non-smooth MaxEnt models are based on coordinate descent algorithms [61,62,63,64] and first-order optimization algorithms such as forward–backward splitting [40,65,66]. We will discuss these methods below, but note that other types of methods have been developed for MaxEnt models; see, e.g., [10,41,67,68] 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 [10,69] 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) [26,33]. 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 [33] was an earlier version of the GLMNET software package [63]. 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 [26,63].

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 [33], 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 stepsize optimization, which GLMNET avoids to improve performance, the GLMNET implementation may not converge [62,70]. A case in point is that [71] 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., [64,72]) 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 [54]. It would be desirable to have fast optimization methods for when the feature mapping (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 $ϵ$-approximate solutions of (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 (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.

We start with the saddle-point formulation (4) of the generalized MaxEnt estimation problem (3). Based on the work two of the authors provided in [38], we propose to split the infimum and supremum in the saddle-point problem (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 let ${\theta }_0=0$, let ${\mathit{w}}_{-1}={\mathit{w}}_0\in \mathrm{int}\mathrm{dom}H$, let ${\mathit{z}}_0\in {\mathbb{R}}^m$, and define the initial probability distribution $p_0\in {\Delta }_n$ through the constant ${\mathit{z}}_0$ via for each $j\in \{1,\dots ,n\}$. Our proposed NPDHG iterative scheme computes the iterates According to ([38], Proposition 4.1), the sequence of iterates $p_k$ converges strongly to the unique solution of the generalized MaxEnt problem (3). Moreover, for any $t⩾0$ and a given tolerance level $ϵ>0$, the scheme (16) provably computes an $ϵ$-approximate solution $p^{ϵ}\left(t\right)$ of the generalized MaxEnt model (3) in $O(1/\sqrt{ϵ})$ 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 [39].

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

Second, because the Kullback–Leibler divergence is 1-strongly convex with respect to the ${\ell }_1$ norm, that is, This fact follows from a fundamental result in information theory known as Pinsker’s inequality [47,48,49,50,51]. In particular, this means that scheme (16) alternates between solving a strongly convex problem over the space $({\mathbb{R}}^n,{∥·∥}_1)$ and a concave problem over the space $({\mathbb{R}}^m,{∥·∥}_2)$. The choice of these spaces is significant, for the induced operator norm ${∥·∥}_{\mathrm{op}}$ of the linear operator $\mathit{A}$ defined in (13) becomes This operator norm offers two crucial advantages: First, it can be computed in optimal $\Theta \left(mn\right)$ time, or better if the features ${\{\mathbf{\Phi }\left(j\right)\}}_{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 (14), with equality, in optimal $\Theta \left(mn\right)$ time. This is important because, in classical first-order optimization methods, we typically require a precise estimate of the largest singular value of the feature mapping (13), namely the number 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\left(\min (m^{2}n,m{n}^2)\right)$ operations. In contrast, our NPDHG scheme (16) does not suffer from this computational bottleneck and therefore scales well. Second, the operator norm ${∥\mathit{A}∥}_{\mathrm{op}}$ can be significantly smaller than the estimate ${∥\mathit{A}∥}_2$, hence allowing for bigger stepsize parameters to further speed up convergence. (An easy calculation yields ${max}_{{∥\mathit{p}∥}_2=1}{∥\mathit{A}\mathit{p}∥}_2⩾{max}_{j\in \{1,\dots ,n\}}{∥\mathbf{\Phi }\left(j\right)∥}_2$.)

Summary. To solve the generalized MaxEnt estimation problem (3) and its dual problem (5), let ${\tau }_0>0$ and ${\sigma }_0>0$ be two stepsize parameters satisfying inequality (14), let ${\theta }_0=0$, let ${\mathit{w}}_{-1}={\mathit{w}}_0\in \mathrm{int}\mathrm{dom}H$, let ${\mathit{z}}_0\in {\mathbb{R}}^m$, and define an initial probability distribution $p_0$ through ${\mathit{z}}_0$ via (15). Then, compute the iterates until convergence is achieved. As all parameters and updates can be computed in $O\left(mn\right)$ time, for any $t⩾0$ and a given tolerance level $ϵ>0$, the overall complexity for computing an $ϵ$-approximate solution $p^{ϵ}\left(t\right)$ is $O(mn/\sqrt{ϵ})$.

The iterative scheme (18) does not require the potential function $H^{*}$ to be smooth. If, however, the potential function $H^{*}$ is ${\gamma }_{H^{*}}$-smooth (equivalently, H is ${\displaystyle \frac{1}{{\gamma }_{H^{*}}}}$-strongly convex) for some ${\gamma }_{H^{*}}>0$, then we can modify the NPDHG iterative scheme (18) to achieve a linear rate of convergence. More precisely, let $t>0$ and introduce the stepsize parameters Let ${\mathit{w}}_{-1}={\mathit{w}}_0\in \mathrm{int}\mathrm{dom}H$, let ${\mathit{z}}_0\in {\mathbb{R}}^m$, and define $p_0$ through ${\mathit{z}}_0$ via (15). Then, the explicit NPDHG iterative scheme is According to ([38], Proposition 4.3), the sequences of iterates $p_k$ and ${\mathit{w}}_k$ converge strongly to the unique solution of the generalized MaxEnt estimation problem (3) and its dual problem (5). Moreover, for any $t>0$ and a given tolerance level $ϵ>0$, this scheme provably computes an $ϵ$-approximate solution $p^{ϵ}\left(t\right)$ of the generalized MaxEnt estimation problem (3) in $O\left(\log \right(1/ϵ\left)\right)$. This rate of convergence is the best achievable rate of convergence with respect to the Nesterov class of optimal first-order methods [39]. As all parameters and updates can be computed in $O\left(mn\right)$ time, the overall complexity for computing an $ϵ$-approximate solution $p^{ϵ}\left(t\right)$ is $O\left(mn\log \right(1/ϵ\left)\right)$.

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

Next, following the procedure outlined in [44], we overlay the fire locations on a $12\mathrm{km}\times 12\mathrm{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 [29,73]. The full WUMI wildfire data set with associated fire-related features for all months from 1984 to 2020 is publicly available as part of the Stochastic Machine Learning for wildfire activity in the western US (SMLFire1.0) package [44].

We construct a weakly informative prior [74,75] for incorporating existing knowledge of wildfire conducive environmental conditions in our MaxEnt model. First, we prepare a training data set of 10,000 fire occurrences and absences chosen randomly across all months between 1984 and 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_{\mathrm{nfire}}=\min \left(p_{\mathrm{fire}}\right)/10$. Last, we normalize 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 Figure 1.

To construct the empirical distribution $\widehat{\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_{\mathrm{r},\mathrm{total}}$ denote the total number of grid cells in ecoregion r and let $n_{\mathrm{r},\mathrm{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}{\displaystyle \frac{n_{\mathrm{r},\mathrm{fire}}}{n_{\mathrm{r},\mathrm{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\}$, $\widehat{\mathcal{D}}\left(j\right)$, as follows: This construction gives the empirical distribution more weight to ecoregions where fires are more frequent and widespread.

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: where $t^{\left(0\right)}$ 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^{\left(0\right)},{\mathit{w}}^{\left(0\right)})=(p_{\mathrm{prior}},\mathbf{0})$ and compute the sequence of solutions ${\left\{(p^{\left(l\right)},{\mathit{w}}^{\left(l\right)})\right\}}_{l=1}^{140}$ of the corresponding MaxEnt primal and dual problems (7) and (8) with hyperparameter $t=t^{\left(l\right)}$.

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 ${\left\{t^{\left(l\right)}\right\}}_{l=0}^{140}$. In the parts below, we describe the value of $t^{\left(0\right)}$ and the NPDHG algorithm used for each model.

In this section, we compare the run times of our NPDHG algorithms with two state-of-the-art optimization algorithms for solving non-smooth MaxEnt models: the forward–backward splitting algorithm (specifically, Algorithm 5 in [65]; see also [40,66]) and the STRUCTMAXENT2 coordinate descent algorithm from [41,64]. 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 [26,33] in MATLAB 2023b 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/{∥\mathit{A}∥}_2$ and 1 (corresponding to $\tau$ and $t_0$ in Algorithm 5 of [65]) and, for computing the pair of solutions $(p^{\left(l\right)},{\mathit{w}}^{\left(l\right)})$ at hyperparameter $t^{\left(l\right)}$, the initial iterate was set to ${\mathit{w}}^{(l-1)}$. In addition, for the elastic net MaxEnt models, the acceleration quantity $q=(1-\alpha )t^{\left(l\right)}/({∥\mathit{A}∥}_2+(1-\alpha )t^{\left(l\right)})$ 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 ([41], pp. 305–306) to make it applicable to the elastic net penalty. For computing the pair of solutions $(p^{\left(l\right)},{\mathit{w}}^{\left(l\right)})$ at hyperparameter $t^{\left(l\right)}$, the initial iterate was set to ${\mathit{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, 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.

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 or the STRUCTMAXENT2 algorithm. Indeed, the stepsize parameters for the NPDHG and forward–backward splitting algorithms are inversely proportional to the norms ${∥A∥}_{\mathrm{op}}$ and ${∥A∥}_2$, and for the wildfire data set these were ${∥\mathit{A}∥}_{\mathrm{op}}\approx 3.30$ and ${∥A∥}_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.

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 Figure 2 and Figure 3. In each case, we have chosen the prior distribution as a benchmark for our spatial plots of fire probability.

The spatial fire probabilities for $\alpha =\{0.95,0.40,0.15,0.05\}$ are shown in Figure 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_{\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_{\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 Figure 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 of which converge to the empirical distribution quicker than any of the elastic net cases.

In Figure 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_{\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 [42,43] as well as previous statistical analyses of wildfire drivers [29,30,31,44].