An Equation-By-Equation Method for Solving the Multidimensional Moment Constrained Maximum Entropy Problem

TL;DR

The paper introduces an equation-by-equation (EBE) method combining homotopy continuation and Newton's methods to solve complex multidimensional maximum entropy problems, demonstrating robustness and improved accuracy over existing methods.

Contribution

A novel EBE approach for solving nonlinear systems in maximum entropy problems, with proven local convergence and adaptive constraint selection for higher-dimensional cases.

Findings

The EBE method successfully solves high-dimensional entropy problems with explicit solutions.

It outperforms classical Newton and BFGS methods in accuracy.

The method adapts to problems where solutions may not exist, selecting relevant constraints.

Abstract

An equation-by-equation (EBE) method is proposed to solve a system of nonlinear equations arising from the moment constrained maximum entropy problem of multidimensional variables. The design of the EBE method combines ideas from homotopy continuation and Newton's iterative methods. Theoretically, we establish the local convergence under appropriate conditions and show that the proposed method, geometrically, finds the solution by searching along the surface corresponding to one component of the nonlinear problem. We will demonstrate the robustness of the method on various numerical examples, including: (1) A six-moment one-dimensional entropy problem with an explicit solution that contains components of order $10^0-10^3$ in magnitude; (2) Four-moment multidimensional entropy problems with explicit solutions where the resulting systems to be solved ranging from $70-310$ equations; (3) Four- to eight-moment of a two-dimensional entropy problem, which solutions correspond to the densities of the two leading EOFs of the wind stress-driven large-scale oceanic model. In this case, we find that the EBE method is more accurate compared to the classical Newton's method, the MATLAB generic solver, and the previously developed BFGS-based method, which was also tested on this problem. (4) Four-moment constrained of up to five-dimensional entropy problems which solutions correspond to multidimensional densities of the components of the solutions of the Kuramoto-Sivashinsky equation. For the higher dimensional cases of this example, the EBE method is superior because it automatically selects a subset of the prescribed moment constraints from which the maximum entropy solution can be estimated within the desired tolerance. This selection feature is particularly important since the moment constrained maximum entropy problems do not necessarily have solutions in general.

Tables

Table 1. Table 1 . The coordinate of the solutions of Example  1 for each iteration, starting from ( 0 , 0 , 0 ) 0 0 0 (0,0,0) . For each outer loop i 𝑖 i , the EBE takes few iterates ( m 𝑚 m ) to find the i − limit-from 𝑖 i- dimensional solution, fixing λ j = α j = 0 subscript 𝜆 𝑗 subscript 𝛼 𝑗 0 \lambda_{j}=\alpha_{j}=0 for j > i 𝑗 𝑖 j>i .

	1	2	3
0	(0,0,0)	(2.30,0,0)	(1.58,1.43,0)
1	(1.76,0,0)	(2.23,0.22,0)	(1.52,1.38,0.26)
2	(2.23,0,0)	(1.87,0.57,0)	(1.12,1.09,0.76)
3	(2.30,0,0)	(1.67,1.21,0)	(1,1,1)
4		(1.58,1.43,0)

Table 2. Table 2 . Summary of solutions for Example  4 : Moment errors for different algorithms with different grids.

Methods	order
	4	6	8
BFGS algorithm with uniform grid	$4.07\times {10}^{-2}$	$1.45\times {10}^{-4}$	$1.14\times {10}^{-2}$
EBE algorithm with uniform grid	$1.27\times {10}^{-11}$	$9.84\times {10}^{-15}$	$7.75\times {10}^{-13}$
EBE algorithm with sparse grid	$7.54\times {10}^{-12}$	$8.12\times {10}^{-15}$	$2.43\times {10}^{-13}$
Matlab \asciifamilyfsolve.m with sparse grid	$4.70\times {10}^{-7}$	$1.19\times {10}^{-4}$	$1.74\times {10}^{-4}$
Newton with sparse grid	$5.12\times {10}^{-11}$	diverge	diverge

Table 3. Table 3 . Summary of solutions for Example  5 .

d	EBE method	fsolve	Newton
2	$1.098\times {10}^{-15}$	$9.779\times {10}^{-7}$	$8.128\times {10}^{-14}$
3	$4.29\times {10}^{-13}$	$3.150\times {10}^{-2}$	diverge
4	$1.19\times {10}^{-14}$	$0.021$	diverge
5	$2.47\times {10}^{-11}$	$0.018$	diverge

Equations

$$S(\rho ):=-{\int }_{\Omega }\log (\rho (\mathbf{x}))\rho (\mathbf{x})d\mathbf{x},$$

$${\mathcal{F}}_{\mathbf{j}}:={\int }_{\Omega }{c}_{\mathbf{j}}(\mathbf{x})\rho (\mathbf{x})d\mathbf{x}=f_{\mathbf{j}},|\mathbf{j}|=0,1,2,\cdots ,p.$$

$$n=\sum _{j=1}^p{C}_{d-1}^{j+d-1},$$

$$\mathcal{L}(\rho (\mathbf{x}),{\lambda }_0,\cdots {\lambda }_n)=S(\rho )+\sum _{j=0}^n{\lambda }_j({\mathcal{F}}_j-f_j).$$

$$\displaystyle\rho(\mathbf{x})=\frac{1}{Z}\exp\Big{(}{\sum_{j=1}^{n}\lambda_{j}c_{j}(\mathbf{x})\Big{)}},$$

$$\displaystyle Z(\lambda_{1},\ldots,\lambda_{n})=\int_{\Omega}\exp\Big{(}{\sum_{j=1}^{n}\lambda_{j}c_{j}(\mathbf{x})}\Big{)}\,d\mathbf{x},$$

$$F_j({\lambda }_1,\cdots ,{\lambda }_n)$$

$${\int }_{\Omega }f(\mathbf{x})d\mathbf{x}\approx \sum _{i}f({\mathbf{x}}_i)w_i,$$

$${\mathbf{F}}_n({\mathbf{\lambda }}_n)=0,$$

$$\displaystyle\mathbf{F}_{n}(\bm{\lambda}_{n}):=\Big{(}F_{1}(\bm{\lambda}_{n}),\ldots,F_{n}(\bm{\lambda}_{n})\Big{)},$$

$${\mathbf{F}}_i({\mathbf{\lambda }}_i,{\alpha }_{i+1},\dots ,{\alpha }_n)=0,i=1,\dots ,n,$$

$${\mathbf{F}}_1({\mathbf{\lambda }}_1,{\alpha }_2,\dots ,{\alpha }_n)=0,$$

$${\lambda }_i^{m+1}$$

$${\lambda }_i^0$$

$$\displaystyle\frac{\partial F_{i}}{\partial\lambda_{i}}(\bm{\lambda}_{i-1}^{m},\lambda_{i}^{m})=\int_{\Omega}(c_{i}(\mathbf{x})-f_{i})c_{i}(\mathbf{x})\exp\Big{(}{\sum_{j=1}^{i-1}\lambda_{j}^{m}c_{j}(\mathbf{x})+\lambda_{i}^{m}c_{i}(\mathbf{x})}\Big{)}\,d\mathbf{x},$$

$${\mathbf{F}}_{i-1}({\mathbf{\lambda }}_{i-1}^{m+1},{\lambda }_i^{m+1})={\mathbf{F}}_{i-1}({\mathbf{\lambda }}_{i-1}^m+\Delta \mathbf{\lambda },{\lambda }_i^m+({\lambda }_i^{m+1}-{\lambda }_i^m))=0$$

$${\mathbf{F}}_{i-1}({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m)+{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m)\Delta \mathbf{\lambda }+{\mathbf{F}}_{i-1,{\lambda }_i}({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m)({\lambda }_i^{m+1}-{\lambda }_i^m)=0,$$

$$\Delta \mathbf{\lambda }=-{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}^{-1}({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m){\mathbf{F}}_{i-1,{\lambda }_i}({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m)({\lambda }_i^{m+1}-{\lambda }_i^m),$$

$${\tilde{\mathbf{\lambda }}}_{i-1}^{m+1}$$

$$0={\mathbf{F}}_{i-1}({\mathbf{\lambda }}_{i-1}^{m+1},{\lambda }_i^{m+1})={\mathbf{F}}_{i-1}({\tilde{\mathbf{\lambda }}}_{i-1}^{m+1},{\lambda }_i^{m+1})+{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}({\tilde{\mathbf{\lambda }}}_{i-1}^{m+1},{\lambda }_i^{m+1})\Delta \tilde{\mathbf{\lambda }},$$

$${\mathbf{\lambda }}_{i-1}^{m+1}={\tilde{\mathbf{\lambda }}}_{i-1}^{m+1}-{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}({\tilde{\mathbf{\lambda }}}_{i-1}^{m+1},{\lambda }_i^{m+1}{)}^{-1}{\mathbf{F}}_{i-1}({\tilde{\mathbf{\lambda }}}_{i-1}^{m+1},{\lambda }_i^{m+1}).$$

$$({\mathbf{\lambda }}_{i-1}^{m+1},{\lambda }_{i+1}^{m+1})={\mathbf{H}}_i({\mathbf{\lambda }}_{i-1}^m,{\lambda }_i^m),$$

$${\mathbf{H}}_i({\mathbf{\lambda }}_{i-1},{\lambda }_i):=\left(\begin{array}{c}{\mathbf{g}}_i-{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}({\mathbf{g}}_i,H_{i,2}{)}^{-1}{\mathbf{F}}_{i-1}({\mathbf{g}}_i,H_{i,2})\\ {\lambda }_i-(\frac{\partial F_i}{\partial {\lambda }_i}({\mathbf{\lambda }}_{i-1},{\lambda }_i){)}^{-1}{F}_i({\mathbf{\lambda }}_{i-1},{\lambda }_i)\end{array}\right).$$

$${\mathbf{g}}_i:={\mathbf{\lambda }}_{i-1}-{\mathbf{F}}_{i-1,{\mathbf{\lambda }}_{i-1}}({\mathbf{\lambda }}_{i-1},{\lambda }_i{)}^{-1}{\mathbf{F}}_{i-1,{\lambda }_i}({\mathbf{\lambda }}_{i-1},{\lambda }_i)(H_{i,2}-{\lambda }_i)$$

$$\displaystyle\mu_{i}^{(i)}=\mu_{i}^{(i)}-\Big{(}\frac{\partial F^{*}_{i}}{\partial\lambda_{i}}\Big{)}^{-1}F^{*}_{i},$$

$${\mathbf{g}}_i^{\ast }:={\mathbf{g}}_i({\mathbf{\mu }}_{i-1}^{(i)},{\mu }_i^{(i)})={\mathbf{\mu }}_{i-1}^{(i)}.$$

$${\mathbf{F}}_i({\mathbf{\lambda }}_{i-1},{\lambda }_i,{\alpha }_{i+1},\dots ,{\alpha }_n)=0,$$

$$\frac{\partial {\mathbf{H}}_{i,1}^{\ast }}{\partial {\lambda }_j}$$

$$\frac{\partial H_{i,2}^{\ast }}{\partial {\lambda }_j}$$

$$J$$

$$\mathbf{v}$$

$$\mathbf{c}$$

Full text

The research of W.H. was partially supported by the American Heart Association under Grant 17SDG33660722 and the Institute for CyberScience Seed Grant.

The research of J.H. was partially supported by the Office of Naval Research Grants N00014-16-1-2888, MURI N00014-12-1-0912 and the National Science Foundation Grants DMS-1317919, DMS-1619661.

1. Introduction

The maximum entropy principle provides a natural criterion for estimating the least biased density function subjected to the given moments [14]. This density estimation approach has a wide range of applications, such as, the harmonic solid and quantum spin systems [20], econometrics [26], and geophysical applications [3, 13]. In a nutshell, this moment constrained method is a parametric estimation technique where the resulting density function is in the form of an exponential of polynomials. This is a consequence of maximizing the Shannon entropy subjected to the polynomial moment constraints, which is usually transformed into an unconstrained minimization problem of a Lagrangian function [27]. Standard approaches for solving this unconstrained minimization problem are based on Newton’s iterative method [4, 27] or quasi-Newton’s based method such as the BFGS method [2, 5].

In the last two papers [2, 5], where the BFGS-based method was introduced and reviewed, they considered minimization problems that involve 44-83 equations, resulting from a 2D problem with moment constraints of up to order-eight, a 3D problem with moment constraints of up to order-six, and a 4D problem with moment constraints of up to order-four. In this paper, we introduce a novel equation solver that can be used to find density function of moderately high dimensional problems (e.g., systems of 70-310 equations resulting from moments up to order-four of 4-7 dimensional density functions) provided that the solutions exist. The proposed method, which we called the Equation-By-Equation (EBE) method, is an iterative method that solves a one-dimensional problem at the first iterate, a two-dimensional problem at the second iterate, a three-dimensional problem at the third iterate, and eventually, solves the full system of nonlinear equations corresponding to the maximum entropy problem at the last iterate. Technically, this method combines Newton’s method with ideas from homotopy continuation. We will show that the EBE method is locally convergent under appropriate conditions. Furthermore, we will provide sufficient conditions for global convergence. Through the convergence analysis, we will show that, geometrically, the proposed method finds the solution of the nonlinear system of equations by tracking along the surface corresponding to one component of the system of nonlinear equations. The EBE method automatically selects a subset of the prescribed constraints from which the maximum entropy solution can be estimated within the desired tolerance. This is an important feature since the maximum entropy problems do not necessarily have solutions for general set of moment constraints.

We shall find that the EBE method produces more accurate solutions (smaller error in the moments) compared to the classical Newton’s method, the MATLAB built-in \asciifamilyfsolve.m, and BFGS method on the test problem in [2, 5] and on test problems based on the solutions of the Kuramoto-Shivashinski equation. Numerically, we will demonstrate that the EBE method is able to solve problems where the true solutions consist of components of order $10^{0}-10^{3}$. We shall also see that the EBE method can solve a system of hundreds of equations in various examples, including those with explicit solutions as well as those with densities estimated based on solutions of complex spatially extended dynamical systems.

The remaining part of the paper is organized as follows. In Section 2, we give a brief overview of the multidimensional maximum entropy problem. In Section 3, we introduce the EBE algorithm. In Section 4, we provide the local convergence analysis. In Section 5, we discuss the practical issues with the proposed method and provide remedies. In Section 6, we demonstrate the robustness of the EBE method on various numerical examples. In Section 7, we conclude the paper with a brief summary and discussion. We include an Appendix to show some computational details that are left out in the main text. Interested readers and users can access the EBE codes (written in MATLAB) at [10].

2. An overview of the maximum entropy problem

We consider the Haussdorf moment-constrained maximum entropy problem [4, 5, 8]. That is, find the optimal probability density $\rho^{*}({\bf x})$ which maximizes the Shannon entropy,

[TABLE]

where $\mathbf{x}\in\Omega=[-1,1]^{d}$ satisfies the following linear constraints,

[TABLE]

In applications, one usually computes the statistics $f_{\bf j}$ from samples of data. For arbitrary finite domain, one can rescale the data to the domain $\Omega$.

While $c_{\bf j}(\mathbf{x})$ can be arbitrary functions in $L^{1}(\Omega,\rho)$, we will focus on the usual uncentered statistical moments with monomial basis functions, $c_{\bf j}(\mathbf{x})={\bf x}^{\bf j}$ in this article, where we have adopted the notations ${\bf x}=(x_{1},\ldots,x_{d})\in\Omega$, ${\bf j}=(j_{1},\ldots,j_{d})\in\mathbb{Z}_{+}^{d}$ with $\mathbb{Z}_{+}=\{0,1,2,\ldots\}$, and ${\bf x}^{\bf j}=\prod_{i=1}^{d}x_{i}^{j_{i}}$. In (2), the quantities $f_{\bf j}$ are the given ${\bf j}$-th moments that can be computed from the data. Since the total number of monomials ${\bf x}^{\bf j}$ where $|{\bf j}|=j$ is $C^{j+d-1}_{d-1}$, then the total number of constraints in (2) for moments up to order-$p$ is,

[TABLE]

excluding the normalization factor corresponding to $c_{\bf 0}({\bf x})=1$. For example, in two-dimensional problem, the total number of moments up to order $p=4$ is $n=14$. To simplify the notation below, we will use a single index notation and understood that the total number of constraints to be satisfied is $n$, excluding the zeroth moment. The exclusion of the zeroth moment will be clear as we discuss below.

By introducing Lagrange multipliers, the above constrained optimization problem can be transformed into the following unconstrained problem:

[TABLE]

In order to find a solution of (3), we set $\frac{\partial\mathcal{L}}{\partial\rho}=0$, which gives,

[TABLE]

where we have defined $Z=\exp(1-\lambda_{0})$. Since $\int_{\Omega}\rho(\mathbf{x})d\mathbf{x}=1$, we have

[TABLE]

which indicates that $Z$ (or implicitly $\lambda_{0}$) is a function of $\lambda_{1},\ldots,\lambda_{n}$. Therefore, the normalization factor $Z$ can be computed via (5) once $\lambda_{1},\ldots,\lambda_{n}$ are estimated. Therefore, we can just concentrate on finding the Lagrange multipliers $\lambda_{1},\ldots,\lambda_{n}$ which satisfy $n$ constraints in (2), excluding the case $c_{\bf{0}}({\bf x})=1$. In particular, the constrained maximum entropy problem is to solve the following nonlinear system of integral equations,

[TABLE]

for $\lambda_{1},\ldots,\lambda_{n}$.

In our numerical implementation, the integral in system (6) will be approximated with a nested sparse grid quadrature rule [9],

[TABLE]

where $\mathbf{x}_{i}$ are the nested sparse grid nodes, and $w_{i}$ are the corresponding weights based on the nested Clenshaw-Curtis quadrature rule [25]. The number of nodes depends on the dimension of the problem $d$ and the number of the nested set (based on the Smolyak construction [23]), is denoted with the parameter $\ell$ (referred as level). In the numerical implementation, we need to specify the parameter $\ell$.

3. An equation-by-equation algorithm

In this section, we describe the new Equation-By-Equation (EBE) technique to solve the system of equations in (6),

[TABLE]

where we have defined,

[TABLE]

and $\bm{\lambda}_{n}=(\lambda_{i},\ldots,\lambda_{n})$. In the following iterative scheme, we start the iteration with an initial condition $(\alpha_{1},\ldots,\alpha_{n})\in\mathbb{R}^{n}$. We define $\bm{\mu}^{(i)}\in\mathbb{R}^{i}$ as the exact solution to the following $i$-dimensional system,

[TABLE]

where we have fixed the last $n-i$ coefficients, $\lambda_{i+1}=\alpha_{i+1}$, …, $\lambda_{n}=\alpha_{n}$. With this notation, the exact solution for (7) is $\bm{\mu}^{(n)}\in\mathbb{R}^{n}$. We also define $\bm{\hat{\mu}}^{(i)}$ to be the numerical estimate of $\bm{\mu}^{(i)}$. With these notations, we now describe the algorithm.

Generally speaking, at each iteration-$i$, where $i=1,\ldots,n$, the EBE algorithm solves a system of $i$-dimensional system in (8). At each step-$i$, given the numerical solution at the previous step, $\bm{\hat{\mu}}^{(i-1)}\in\mathbb{R}^{i-1}$ and initial condition $\alpha_{i}$, we apply idea from homotopy continuation to find the solution $\bm{\mu}^{(i)}\in\mathbb{R}^{i}$ that solves the $i$-dimensional system of equations (8). Notice that we do not only add a new equation $F_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})=0$ but we also estimate the $i$th variable in the previous $(i-1)$ equations, $\mathbf{F}_{i-1}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})={\bf 0}$. The scheme proceeds by solving the larger systems one-by-one until $i=n$ so we eventually solve (7).

Now let us describe how to numerically estimate $\bm{\mu}^{(i)}$ at every step-$i$. For the first step $i=1$, we solve the one-dimensional problem,

[TABLE]

for $\bm{\lambda}_{1}$ with Newton’s method. For the steps $i=2,\ldots,n$, we have $\bm{\hat{\mu}}^{(i-1)}$ which are the numerical estimates of ${\bf F}_{i-1}(\bm{\lambda}_{i-1},\alpha_{i},\ldots,\alpha_{n})={\bf 0}$. To simplify the expression below, let us use $F_{i}(\bm{\lambda}_{i-1},\lambda_{i})$ as a short hand notation for $F_{i}(\bm{\lambda}_{i-1},\lambda_{i},\alpha_{i},\ldots,\alpha_{n})$ to emphasize the independent variables.

We proceed to estimate $\lambda_{i}$ using Newton’s method with $Tol_{1}$ on the $i$-th equation. That is, we iterate

[TABLE]

assuming that $\frac{\partial F_{i}}{\partial\lambda_{i}}(\bm{\lambda}_{i-1}^{m},\lambda_{i}^{m})\neq 0$. Here, the partial derivative of $F_{i}$ with respect to $\lambda_{i}$ evaluated at $\lambda_{i}^{m}$ is defined as,

[TABLE]

where we have denoted $\bm{\lambda}_{i-1}^{m}=(\lambda_{i}^{m},\ldots,\lambda_{i-1}^{m})$. Notice that to proceed the iteration in (9), we need to update $\bm{\lambda}_{i-1}^{m}$ for $m>0$. We propose to follow the homotopy continuation method for this update. In particular, we are looking for $\bm{\lambda}_{i-1}^{m+1}$ that solves $\mathbf{F}_{i-1}(\bm{\lambda}_{i-1}^{m+1},\lambda_{i}^{m+1})={\bf 0}$, given the current estimate $\lambda_{i}^{m+1}$ from (9) as well as $\mathbf{F}_{i-1}(\bm{\lambda}_{i-1}^{m},\lambda_{i}^{m})={\bf 0}$. At $m=0$, this last constraint is numerically estimated by $\mathbf{F}_{i-1}(\bm{\hat{\mu}}^{(i-1)},\alpha_{i})\approx{\bf 0}$.

One way to solve this problem is through the following predictor-corrector step which is usually used in homotopy continuation method [7, 24]. In particular, we apply Taylor’s expansion to

[TABLE]

at $(\bm{\lambda}_{i-1}^{m},\lambda_{i}^{m})$, which gives,

[TABLE]

which means that

[TABLE]

assuming that ${\bf F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\lambda}_{i-1}^{m},\lambda_{i}^{m})$ is invertible. Based on this linear prediction, $\bm{\lambda}_{i-1}^{m+1}$ is approximated by,

[TABLE]

Subsequently, when $\|\mathbf{F}_{i}(\bm{\tilde{\lambda}}_{i-1}^{m+1},\lambda_{i}^{m+1})\|\geq Tol_{2}$ , apply a correction using Newton’s method by expanding,

[TABLE]

assuming that $\bm{\lambda}_{i-1}^{m+1}=\bm{\tilde{\lambda}}_{i-1}^{m+1}+\Delta\bm{\tilde{\lambda}}$, to find that,

[TABLE]

This expression assumes that $\mathbf{F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\tilde{\lambda}}_{i-1}^{m+1},\lambda_{i}^{m+1})$ is invertible.

In summary, at each step-$i$, we iterate (9), (11), (12). So, the outer loop $i$ corresponds to adding one equation to the system at the time and for each $i$, we apply an inner loop, indexed with $m$, to find the solution $\bm{\mu}^{(i)}$ for ${\bf F}_{i}(\bf{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})={\bf 0}$. We denote the approximate solution as $\bm{\hat{\mu}}^{(i)}$. An adaptive tolerance technique is employed to compute the initial guess of ${\bf F}_{i}$ by using Newton’s method. In particular, when the current tolerance $Tol_{2}$ is not satisfied after executing (12), then we divide $Tol_{1}$ by ten until $Tol_{2}$ is met.

Recall that the standard Newton’s method assumes that the Jacobian ${\bf F}_{n,\bm{\lambda}_{n}}\in\mathbb{R}^{n\times n}$ is nonsingular at the root of the full system in (6) to guarantee the local convergence. In the next section, we will show that the EBE method requires the following conditions for local convergence:

Assumption 1**.**

Let $\bm{\mu}^{(i)}\in\mathbb{R}^{i}$ be a solution of ${\bf F}_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})={\bf 0}$, for each $i=1,\ldots,n$. The EBE method assumes the following conditions:

1. (1)

$\frac{\partial F_{i}}{\partial\lambda_{i}}(\bm{\mu}^{(i)},\alpha_{i+1},\ldots,\alpha_{n})\neq 0$. 2. (2)

$\mathbf{F}_{i,\bm{\lambda}_{i}}(\bm{\mu}^{(i)},\alpha_{i+1},\ldots,\alpha_{n})$* are nonsingular.* 3. (3)

Each component of ${\bf F_{i}}$ is twice differentiable in a close region whose interior contains the solution $\bm{\mu}^{(i)}$.

These conditions are similar to the standard Newton’s assumptions on each system of $i$ equations. The smoothness condition will be used in the proof of the local convergence in the next section. Of course if one can specify initial conditions that are sufficiently close to the true solution, then one can simply apply Newton’s method directly. With the EBE method, we can start with any arbitrary initial condition. Theoretically, this will require an additional condition beyond the Assumption 1 for global convergence as we shall discuss in Section 4. In Section 5, we will provide several remedies when the initial condition is not close to the solution. In fact, we will always set the initial condition to zero in our numerical implementation in Section 6, $\alpha_{i}=0,\forall i=1,\ldots,n$, and demonstrate that the EBE method is numerically accurate in the test problems with solutions that are far away from zero.

4. Convergence analysis

In this section, we study the convergence of this method. First, let’s concentrate on the convergence of the iteration (9), (11), (12) for solving the $i$-dimensional system, ${\bf F}_{i}(\bm{\lambda}_{i-1},\lambda_{i},\alpha_{i+1},\ldots,\alpha_{n}):={\bf F}_{i}(\bm{\lambda}_{i-1},\lambda_{i})={\bf 0}$ for $\bm{\lambda}_{i-1}$ and $\lambda_{i}$. In compact form, these three steps can be written as an iterative map,

[TABLE]

where the map ${\bf H}_{i}:\mathbb{R}^{i}\to\mathbb{R}^{i}$ is defined as,

[TABLE]

In (14), the notation $H_{i,2}$ denotes the second component of (14) and

[TABLE]

is defined exactly as in (11).

For notational convenience in the discussion below, we let the components of the exact solution of (8) be defined as $\bm{\mu}^{(i)}:=(\bm{\mu}_{i-1}^{(i)},\mu_{i}^{(i)})\in\mathbb{R}^{i}$. Here, we denote the first $i-1$ components as $\bm{\mu}_{i-1}^{(i)}=(\mu_{1}^{(i)},\ldots,\mu_{i-1}^{(i)})\in\mathbb{R}^{i-1}$. Similarly, we also denote ${\bf H}_{i}=({\bf H}_{i,1},{H}_{i,2})$. First, we can deduce that,

Theorem 4.1**.**

Let $\bm{\mu}^{(i)}\in\mathbb{R}^{i}$ be a fixed point of (13). Assume that ${\bf F}^{*}_{i-1,\bm{\lambda}_{i-1}}:={\bf F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\mu}^{(i)})$ is nonsingular and $\frac{\partial F^{*}_{i}}{\partial\lambda_{i}}:=\frac{\partial F_{i}}{\partial_{\lambda_{i}}}(\bm{\mu}^{(i)})\neq 0$, then ${\bf F}^{*}_{i}:={\bf F}_{i}(\bm{\mu}^{(i)})={\bf 0}$.

Proof.

Evaluating the second equation in (14) at the fixed point, we obtain

[TABLE]

which means that $F^{*}_{i}:=F_{i}(\bm{\mu}^{(i)})=0$. This also implies that $H^{*}_{i,2}=\mu_{i}^{(i)}$, where $H^{*}_{i,2}$ denotes the second component of (14) evaluated at the fixed point. Subsequently,

[TABLE]

Substituting $H^{*}_{i,2}=\mu_{i}^{(i)}$ and $\bm{g}^{*}_{i}=\bm{\mu}_{i-1}^{(i)}$ into $\bm{\mu}_{i-1}^{(i)}=\bm{H}_{i,1}^{*}$, where $\bm{H}^{*}_{i,1}$ denotes the first equation in (14) evaluated at the fixed point $\bm{\mu}^{(i)}$, we immediately obtain $\mathbf{F}^{*}_{i-1}:=\mathbf{F}_{i-1}(\bm{\mu}^{(i)})={\bf 0}$ and the proof is completed. ∎

This theorem says that the fixed points of (13) are indeed the solutions of

[TABLE]

which is what we intend to solve on each iteration $i=2,\ldots,n$. Next, we will establish the condition for the fixed point to be locally attracting. This condition will ensure that if we iterate the map in (14) with an initial condition that is close to the solution, then we will obtain the solution.

For local convergence, we want to show that eigenvalues of the Jacobian matrix $D{\bf H}_{i}^{*}:=D{\bf H}_{i}(\bm{\mu}^{(i)})$ are in the interior of the unit ball of the complex plane. One can verify that the components of the Jacobian matrix $D{\bf H}_{i}^{*}$ are given by,

[TABLE]

for $j=1,\ldots,i$, where we have used all the three conditions in the Assumption 1 (see Appendix A for the detailed derivation). Here, $\delta_{j,i}$ is one only if $j=i$ and zero otherwise. To simplify the discussion below, let’s define the following notations,

[TABLE]

such that,

[TABLE]

We can now obtain the following result:

Theorem 4.2**.**

Let $\bm{\mu}^{(i)}\in\mathbb{R}^{i}$ be a fixed point of (13) such that the conditions in the Assumption 1 are satisfied. Let’s $\sigma_{j}({\bm{F}}_{i-1,{\bm{\lambda}}_{i-1}}^{*})$ be the eigenvalues of $\mathbf{F}_{i-1,{\bm{\lambda}}_{i-1}}^{*}$ and assume that they satisfy the following order $|\sigma_{1}|\geq|\sigma_{2}|\geq\ldots|\sigma_{i-1}|$. If

[TABLE]

then $\bm{\mu}^{(i)}$ is locally attracting.

Proof.

From (19), we only need to analyze the eigenvalues of $J^{-1}{\bf v}{\bf c}^{\top}$. From basic matrix theory, recall that the magnitude of the largest eigenvalue can be bounded above as follows,

[TABLE]

where $\|\cdot\|_{2}$ denotes the matrix $\ell_{2}$-norm. For the fixed point to be locally attracting, all of the eigenvalues of $J^{-1}{\bf v}{\bf c}^{\top}$ have to be in the interior of the unit ball in the complex plane. This means that we only need to show that $\|J^{-1}\|_{2}\|{\bf v}{\bf c}^{\top}\|_{2}<1$ or $\|{\bf v}{\bf c}^{\top}\|_{2}<|\sigma_{i-1}(J)|$, where $\sigma_{i-1}(J)$ denotes the smallest eigenvalue of the $(i-1)\times(i-1)$ matrix $J$ following the ordering in the hypothesis.

Since $\mbox{Tr}({\bf v}{\bf c}^{\top})=\sum_{j=1}^{i}\sigma_{j}({\bf v}{\bf c}^{\top})$ and ${\bf v}{\bf c}^{\top}$ is a rank-one matrix, then its nontrivial eigenvalue is given by,

[TABLE]

where we have used the definitions in (18) and the second component in (17). From the assumption in (20), we have

[TABLE]

and the proof is completed. ∎

This theorem provides the conditions for local convergence on each iteration-$i$. In particular, if the hypothesis in Theorem 4.2 is satisfied, we will find the solutions to (8) by iterating (13) provided that we start with a sufficiently close initial condition. Notice also that this condition suggests that in practice the local convergence will be difficult to satisfy if the Jacobian matrix ${F}_{i-1,\bm{\lambda}_{i-1}}$ is close to singular. With these two theorems, we can now establish

Theorem 4.3**.**

Let $\bm{\mu}^{(n)}\in\mathbb{R}^{n}$ be the solution of the n-dimensional system of equations in (7). We assume the hypothesis in Theorem 4.2, then the EBE method is locally convergent.

Proof.

Choose an initial condition, $(\alpha_{1},\ldots,\alpha_{n})$, that is sufficiently close to the solution $\bm{\mu}^{(n)}$ of ${\bf F}_{n}(\bm{\lambda}_{n})={\bf 0}$. First, let us define the surface $F_{1}(\lambda_{1},\ldots,\lambda_{n})=0$ as $\mathcal{M}_{n}$; here, the dimension of $\mathcal{M}_{n}$ is at most $n-1$. Subsequently, we define the surfaces ${\bf F}_{2}(\bm{\lambda}_{n})=\bm{0}$ as $\mathcal{M}_{n-1}$, ${\bf F}_{3}(\bm{\lambda}_{n})=\bm{0}$ as $\mathcal{M}_{n-2}$, and so on. The dimension of $\mathcal{M}_{j}$ is at most $j-1$. We assume that ${\bf F}_{n}(\bm{\lambda}_{n})=\bm{0}$ has at least a solution, then $\mathcal{M}_{1}$ contains the solution $\bm{\mu}^{(n)}$. It is clear that $\mathcal{M}_{n}\supset\mathcal{M}_{n-1}\supset\ldots\supset\mathcal{M}_{1}$.

For $i=1$, we solve $F_{1}(\lambda_{1},\alpha_{2},\ldots,\alpha_{n})=0$ for $\lambda_{1}$. Geometrically, we look for the first coordinate on the surface $\mathcal{M}_{n}$. From the Assumption 1.2, we have the local convergence of the usual Newton’s iteration. If $\alpha_{1}$ is sufficiently close to the solution $\bm{\mu}^{(1)}=\mu_{1}^{(1)}\in\mathbb{R}$, as $m\to\infty$ we obtain the solution $(\mu_{1}^{(1)},\alpha_{2},\ldots,\alpha_{n})\in\mathcal{M}_{n}$. By the smoothness assumption, $(\mu_{1}^{(1)},\alpha_{2},\ldots,\alpha_{n})$ is also close to $\bm{\mu}^{(n)}$.

Continuing with $i>1$, we want to solve ${\bf F}_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})={\bf 0}$ for $\bm{\lambda}_{i}$. Numerically, we will apply the iterative map ${\bf H}_{i}$ in (13) starting from $(\bm{\mu}^{(i-1)},\alpha_{i},\ldots,\alpha_{n})\in\mathcal{M}_{n-i+2}$. By Assumption 1.2, the Jacobian, ${\bf F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\mu}^{(i-1)},\alpha_{i},\ldots,\alpha_{n})$ is nonsingular so by implicit function theorem, for any local neighborhood $V$ of $\bm{\mu}^{(i-1)}$, there exists a neighborhood $U$ of $\alpha_{i}$ and a $C^{1}$ function $\bm{h}_{i-1}:U\to V$ such that $\bm{\mu}^{(i-1)}=\bm{h}_{i-1}(\alpha_{i})$ and ${\bf F}_{i-1}(\bm{h}_{i-1}(\lambda_{i}),\lambda_{i},\alpha_{i+1},\ldots,\alpha_{n})=0$ for all $\lambda_{i}\in U$. Since the initial condition $\alpha_{i}$ is close to $\mu_{i}^{(n)}$, by the smoothness assumption it is also close to $\mu_{i}^{(i)}$ that solves ${\bf F}_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})=0$. The continuity of $\bm{h}_{i-1}$ on $U$ means that $(\bm{\mu}^{(i)}_{i-1},\mu_{i}^{(i)})\in V\times U$. Geometrically, this means the surface $F_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})=0$ intersects with the curve $\bm{\lambda}_{i-1}=\bm{h}_{i-1}(\lambda_{i})$ at $\bm{\mu}^{(i)}=(\bm{\mu}^{(i)}_{i-1},\mu_{i}^{(i)})$. Therefore, we can find the solution for this $i$-dimensional system by tracking along the curve $\bm{\lambda}_{i-1}=\bm{h}_{i-1}(\lambda_{i})$ where we consider $\lambda_{i}$ as an independent parameter. The iterative map ${\bf H}_{i}$ in (14) is to facilitate this tracking and the conditions in Theorem 4.2 guarantee convergence to the solution. Notice that during this iteration, the solution remains on $\mathcal{M}_{n-i+2}$. The solution for this i-dimensional problem is $(\bm{\mu}^{(i)},\alpha_{i+1},\ldots,\alpha_{n})\in\mathcal{M}_{n-(i+1)+2}\subset\mathcal{M}_{n-i+2}\subset\ldots\subset\mathcal{M}_{n}$. Continuing with the same argument, we find that for $i=n$, $\bm{\mu}^{(n)}\in\mathcal{M}_{1}\subset\mathcal{M}_{n}$. ∎

This iterative procedure finds the solution by searching along the manifold $\mathcal{M}_{n}$ in the direction of the hypersurfaces of a single parameter at a time, which local existence is guaranteed by the Assumption 1. It is clear that after each step-$i$, the estimated solution may not necessarily be closer to the true solution since the estimates do not minimize the closest path to the true solution along the manifold $\mathcal{M}_{n}$ (or the geodesic distance). This means that, locally,

[TABLE]

for $i<n-1$ is not true.

In practice, when initial conditions are not closed to the solution, the (global) convergence of EBE requires the following additional condition: For every $i$, there exists a nonempty connected set that contains $(\bm{\mu}^{(i)},\alpha_{i+1})$ and $\bm{\mu}^{(i+1)}$ such that $\mathbf{F}_{i,\bm{\lambda}_{i}}$ evaluated at any point in this set is nonsingular. The existence of this set will allow us to build a path to connect these two points that are far apart. If this condition is not met, we need an additional treatment to overcome this issue which will be discussed in the next section.

5. Practical challenges

In this section, we will discuss several practical challenges related to our algorithm with remedies. They include non-locality of the initial condition, mistracking due to multiple solutions, non-existence of solutions within the desired numerical tolerance, and the computational complexity.

5.1. Adaptive tracking

As we mentioned in the previous section, the EBE method only converges locally, which means that it requires an adequate initial condition which is practically challenging. In our numerical simulations below, in fact, we always start from zero initial condition, $\alpha_{i}=0,\forall i=1,\ldots,n$. In this case, notice that even when we obtain an accurate solution at step-$i$, that is, ${\bf F}_{i}(\bm{\hat{\mu}}^{(i)})\approx{\bm{0}}$, as we proceed to the next iteration, $|{F}_{i+1}(\bm{\hat{\mu}}^{(i)},\alpha_{i+1})|\gg 0$, meaning that $(\bm{\hat{\mu}}^{(i)},\alpha_{i+1})$ is not close to the solution, $\bm{\mu}^{(i+1)}$. Even when $\frac{\partial F_{i+1}}{\partial\lambda_{i+1}}(\bm{\hat{\mu}}^{(i)},\alpha_{i+1})$ is not singular, according to equation (9), $\lambda_{i}^{m+1}$ could be very far away from $\lambda_{i}^{m}$. In this case, Newton’s method could fail in Eq. (12) because the initial guess could be very far from the solution.

As a remedy, we employ an adaptive tracking on $\lambda_{i}$ to guarantee that the application of Newton’s method is within its zone of convergence for each predictor-corrector step. The idea of the adaptive tracking is that we cut the tracking step, $\Delta\lambda_{i}:=\lambda_{i+1}-\lambda_{i}$, by half until the prediction-correction step in (11)-(12) converges. The detail algorithm is outlined below.

5.2. Bifurcation

In order to solve $F_{i}(\lambda_{1},\lambda_{2},\cdots,\lambda_{i})=0$, we track $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})=\mathbf{0}$ along $\lambda_{i}$ as a parameter. During this parameter tracking, we may have some bifurcation points of $\lambda_{i}$ for the nonlinear system $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})=\bm{0}$. This means that the Jacobian, $\bm{F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\lambda}_{i-1},\lambda_{i})$ is rank deficient such that $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})=\bm{0}$ has multiple solutions $\bm{\lambda}_{i-1}$ for a given $\lambda_{i}$. In this situation, $F_{i}$ has multiple realizations functions of $\lambda_{i}$ (see the illustration in Figure 1 where the bifurcation point is the intersection of the two possible realizations of $F_{i}$). In this illustration, the goal is to track along the red branch to find the root, $F_{i}(\lambda_{i})=0$. As we get closer to the bifurcation point, the Jacobian, $\bm{F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\lambda}_{i-1},\lambda_{i})$, is singular such that we can’t evaluate (11). Intuitively, the existence of multiple solutions near the bifurcation point induces a possibility of mistracking from the red curve to the green curve (as shown by the arrows) which prohibits one to find the solution.

To avoid such mistracking, we apply the deflation technique to compute the bifurcation point directly [12, 16]. Once the bifurcation point is estimated, we approximate the correct branches using the Richardson extrapolation to avoid mistracking. Denote the bifurcation point as $\lambda_{i}^{*}$, the nonlinear system $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})=\bm{0}$ is difficult to solve when $\lambda_{i}$ is close to $\lambda_{i}^{*}$ since the Jacobian of $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})$ becomes near singular. If the last attempt is $(\tilde{\bm{\lambda}}_{i-1},\tilde{\lambda}_{i})$, we compute $(\bm{\lambda}_{i-1}^{*},\lambda_{i}^{*})$ by solving the following deflated system:

[TABLE]

where $\bm{v}$ is the kernel of $\bm{F}_{i-1,\bm{\lambda}_{i-1}}(\bm{\lambda}_{i-1},\lambda_{i})$ and $\bm{\xi}$ is a random vector to guarantee that $\bm{v}$ is not a zero eigenvector. In this case, $G(\bm{\lambda}_{i-1}^{*},\lambda_{i}^{*},\bm{v})$ is well-conditioned [12, 16]. Once the bifurcation point $(\bm{\lambda}_{i-1}^{*},\lambda_{i}^{*})$ is estimated, we can avoid mistracking by setting $\lambda_{i}=2\lambda_{i}^{*}-\tilde{\lambda}_{i}$ and solve $\bm{F}_{i-1}(\bm{\lambda}_{i-1},\lambda_{i})=\bm{0}$ by using Newton’s method with an initial guess $2\bm{\lambda}_{i-1}^{*}-\tilde{\bm{\lambda}}_{i-1}$ (which is a Richardson extrapolation).

5.3. Nonexistence of solutions

In general, the moment constrained maximum entropy problems may not necessarily have solutions. Even when the solutions exist theoretically, they could be difficult to find numerically due to the noisy dataset, error in the numerical integration, etc. In this case, we simply discard the equation $F_{i}$ when the minimum is larger than the desired tolerance. This feature (discarding the constraints that give no solutions) is only feasible in the EBE algorithm. However, some theories are needed to preserve the convexity of the polynomials in the exponential term of Eq. (4) while discarding some of these constraints. In our numerical simulations below, we handle this issue by re-ordering the constraints. In particular, for a problem with moment constraints up to order-$4$, we include the constraints corresponding to $\mathbb{E}[x_{i}^{4}]$ $(i=1,\cdots,d)$ in the earlier step of the EBE iterations to avoid these constraints being discarded. Note that this method is sensitive to ordering, that is, different ordering of constraints yields different path to compute the solution. Therefore, a systematic ordering technique that simultaneously preserves the convexity of the polynomial in the exponential term of Eq. (4) is an important problem to be addressed in the future.

5.4. Computational complexity

The most expensive computational part in EBE is the numerical evaluation of (6). For a fast numerical integration, we store the monomial basis $c_{\bm{j}}(\bm{x})$ as a matrix of size $N_{\ell}\times n$, where $N_{\ell}$ is the number of sparse grid points and $n$ is number of monomial basis. In this case, the computational cost in evaluating $F_{j}$ is $(2j+1)N_{\ell}$ ($j-1$ additions, $j+1$ multiplications and 1 subtraction for each grid point), excluding the computational cost for exponential function evaluation, which is on the order of $\log^{2}m$ to obtain an error of resolution $2^{-m}$ [6]. For the $i$-th iteration of the EBE algorithm, the computational cost to evaluate the $i$-dimensional system $\bm{F}_{i}$ is $\sum_{j=1}^{i}(2j+1)N_{\ell}=\frac{i^{2}+i}{2}N_{\ell}$, excluding the exponentiation.

6. Numerical results

In this section, we show numerical results of the EBE method on five examples. In all of the simulations below, unless stated, we set the Newton’s tolerance $Tol_{1}={10^{-1}}$ and the predictor tolerance $Tol_{2}={10^{-10}}$. In the first test example, we will describe how the EBE method works on each iteration. The goal of the second example is to demonstrate the global convergence with solutions that are far away from initial condition, $\alpha_{j}=0$. In particular, we will test the EBE method on a problem with solutions, $\lambda_{j}$, that have magnitudes ranging from order $10^{0}-10^{3}$. In this example, we will show the robustness of the estimate as a function of the number of integration points (or the sparse grid level $\ell$). The third example is to demonstrate the performance on high dimensional problems (with $70\leq n\leq 310$ of order hundreds), induced from order-four moments of four to seven dimensional density functions. While these first three examples involve estimating densities of the form (4), in the next two examples, we also test the EBE method to estimate densities from a given data set where the maximum entropy solutions may or may not exist. The first data-driven problem is to estimate densities of the first two leading EOFS of the wind stress-driven large-scale oceanic model [2, 5]. The second data-driven problem is to estimate two- to five-dimensional densities arising from solutions of the Kuramoto-Sivashinsky equation. In these two problems, we compare our method with the classical Newton’s method, the MATLAB built-in solver \asciifamilyfsolve.m, and the previously developed BFGS-based method [2, 5].

Example 1**.**

We consider a simple example $\rho(x)\propto\exp(x+x^{2}+x^{3})$ for $x\in[-1,1]$ so that the exact solution is ${\bm{\lambda}}=(1,1,1)$. Here, the moments $f_{j}$ can be computed numerically as follows,

[TABLE]

In order to numerically integrate both the denominator and numerator, we used a regular one-dimensional sparse grid of level $\ell=7$ (the number of nodes is 65). Our goal here is to illustrate the method and to show the trajectory of the solutions after each iteration of the inner loop $m$ and outer loop $i$. In Figure 2, we show the surface of $F_{1}(\lambda_{1},\lambda_{2},\lambda_{3})=0$ (grey). For $i=1$, we solve the $F_{1}(\lambda_{1},0,0)=0$, after three iterations ($m=3$) the solution converges to $\lambda_{1}=2.3$ (see Table 1). For $i=2$, we start with this solution and introduce the second variable $\lambda_{2}$ for solving the second equation $F_{2}(\lambda_{1},\lambda_{2},0)=0$ with constraint $F_{1}(\lambda_{1},\lambda_{2},0)=0$. Here, the solution follows the path $\lambda_{1}=h_{1}(\lambda_{2})$ thanks to the implicit function theorem (black curve). Numerically, a sequence of (green) points following this path converges to a point that satisfies $F_{1}(\lambda_{1},\lambda_{2},0)=F_{2}(\lambda_{1},\lambda_{2},0)=0$ (the green point in the intersection between black and red curves in Figure 2). In the next iteration $i=3$, we introduce the third variable $\lambda_{3}$ for solving the third equation $F_{3}(\lambda_{1},\lambda_{2},\lambda_{3})=0$ with constraints $F_{1}(\lambda_{1},\lambda_{2},\lambda_{3})=F_{2}(\lambda_{1},\lambda_{2},\lambda_{3})=0$. By the implicit function theorem, we have $(\lambda_{1},\lambda_{2})=h_{2}(\lambda_{3})$ that satisfies $F_{1}(h_{2}(\lambda_{3}),\lambda_{3})=F_{2}(h_{2}(\lambda_{3}),\lambda_{3})=0$, which is shown in red curve in Figure 2. On this red curve, we have a sequence of (blue) points which converges to the solution of the full system (cyan point shown in Figure 2). The coordinate of the solution on each iteration is shown in Table 1. Notice that the solutions always lie on the surface $F_{1}(\lambda_{1},\lambda_{2},\lambda_{3})=0$.

Example 2**.**

We consider a one-dimensional example with up to order-six moment constraints with explicit solution given by,

[TABLE]

as shown in Figure 3. This example is a tough test problem since the solution, $\bm{\lambda}=(2,16,24,96,-256,1024)$, has components of order $10^{0}-10^{3}$. Similar to Example 1, we compute the moments $f_{i}$ by using a one-dimensional sparse grid of level $\ell=7$ (65 nodes). The EBE algorithm converges to the exact solution with error, $\|\bm{\lambda}-\bm{\lambda}^{*}\|=5.44\times 10^{-13}$. Since the numerical experiment is performed with an initial condition $\alpha_{j}=0$ that is far from the solution, this result demonstrates a global convergence of the EBE method.

Next, we investigate the sensitivity of the estimates to the number of sparse grid points used in approximating the integral. In our numerical experiments, we estimate the true moments $f_{i}$ using one-dimensional sparse grid of level $\ell=20$ (524,289 nodes) and feed these moment estimates into the EBE algorithm. In Figure 4, we show the error in $\lambda$ (with $\ell_{2}$ metric) for different levels of the sparse grid from 6 to 15 that are used in the EBE method. Notice that the error decreases as a function of $\mathcal{\ell}$ and the improvement becomes negligible for $\ell>8$.

Example 3**.**

In this example, we consider a $d$-dimensional example with an explicit solution,

[TABLE]

on domain $\Omega=[-1,1]^{d}$ where we will vary $d=4,\ldots,7$. For these simulations, we consider up to order-four moment constraints and fix the sparse grid level $\ell=8$ to compute the integration.

Here, the EBE method is able to estimate $\bm{\lambda}$ with $\ell_{2}$-errors of order $10^{-13}$ (the error in $\bm{\lambda}$ is $1.11\times 10^{-13}$ and moments error is $3.15\times 10^{-15}$). In this computation, the dimensions of the nonlinear system are $70$ for $d=4$, $126$ for $d=5$, $210$ for $d=6$, and $310$ for $d=7$. Here, the EBE method is able to recover the true density even if we prescribe more constraints, corresponding to $d$ larger than four.

Example 4**.**

Next, we consider estimating a two-dimensional probability density of the two leading empirical orthogonal functions of a geophysical model for wind stress-driven large-scale oceanic model [18, 19]. This is exactly the same test example in the previously developed BFGS-based method [2, 5]. In fact, the two-dimensional density that we used here was supplied by Rafail Abramov. First, we compare the EBE method with the BFGS algorithm of [2] which code can be downloaded from [1]. In this comparison, we use the same uniformly distributed grid points where the total number of nodes are $85\times 85=7,225$. We set the Newton’s tolerance of the EBE algorithm to be $10^{-10}$. In Table 2 notice that the moment errors of the EBE are much smaller compared to those of the BFGS method.

While the EBE is superior compared to BFGS, we should note that the BFGS method does not use the Hessian of ${\bf F}_{i}$ whereas the EBE does. For a fair comparison, we include results using the MATLAB built-in function \asciifamilyfsolve.m, which default algorithm is the trust-region-dogleg (see the documentation for detail [17]). In our numerical implementation, we apply \asciifamilyfsolve.m with a specified Hessian function ${\bf F}_{n}$. We also include the classical Newton’s method with a specified Hessian function ${\bf F}_{n}$. In this comparison, we use the same sparse grid of level $\ell=11$ (or 7,169 nodes) to compute the two-dimensional integral. Notice that the EBE method is still superior compared to these two schemes as reported in Table 2. In fact, Newton’s method does not converge for higher-order moment constraints. The joint two-dimensional PDFs are shown in Figure 5. The first row is the two-dimensional density function provided by R. Abramov. The second row shows the EBE estimates using up to order four-, six-, and eight-moment constraints. The third and fourth rows show the BFGS and MATLAB \asciifamilyfsolve.m estimates, respectively.

Example 5**.**

In this example, we consider estimating multidimensional densities of the solutions of the Kuramoto-Sivashinsky equation. Here, the solutions are integrated with a fourth-order time-differencing method on $128$ equally spaced grid points over a domain of $[0,32\pi]$ as in [15]. We use initial condition $u(x,0)=\cos(x/(16\xi))(1+\sin(x/16))$, where $\xi\sim U[0,1]$ and integration time step of $0.25$. The data is generated by integrating 10,000 time steps. Based on this data set, we randomly select $d$ components and estimate the $d$-dimensional joint density associated to these components. For visual comparison, we also show the results from a two-dimensional kernel density estimation method [22, 21] as a reference. Numerically, we use the MATLAB built-in function, \asciifamilyksdensity.m. Note that the BFGS algorithm [2] does not work on this data set while the classical Newton’s method only converges for the two-dimensional case. We also show the corresponding results with the MATLAB \asciifamilyfsolve.m with specified Hessian function as in the previous example. The moment errors of these three schemes are reported in Table 3.

In Figure 6, we show the two-dimensional density estimated by EBE algorithm compared to those from the \asciifamilyfsolve.m, the classical Newton’s method, and the 2D kernel density estimate. For the two-dimensional case, the resulting densities are visually identical although the corresponding moment error of the EBE method is still the smallest compared to the Newton’s and the MATLAB \asciifamilyfsolve.m (see Table 3). In Figure 7, we show the contour plot of the two-dimensional marginal densities obtained from solving the three-dimensional problem given four-moment constraints with the EBE method and the MATLAB \asciifamilyfsolve.m. For diagnostic purpose, we also provide the corresponding contour plots of the two-dimensional kernel density estimates. Notice that the MATLAB \asciifamilyfsolve.m produces completely inaccurate estimate. The EBE method produces an estimate that qualitatively agrees to the corresponding two-dimensional KDE estimates. The slight disagreement between these estimates are expected since we only provide up to order-four moments information.

In Figure 8, we show the results for the four-dimensional problem. We do not show the estimate from the MATLAB \asciifamilyfsolve.m since it is not accurate at all. Here, we include more than four-order moments. Specifically, the total number of constraints for up to order-four moments is 70 while this result is based on 87 constraints, including 17 additional higher-order moment constraints that include order-six moments, $\mathbb{E}[x_{i}^{6}],i=1,\ldots,4$. See the movie of the density estimates for each iteration in the supplementary material [11]. Notice that the marginal densities estimated by the EBE look very similar to those estimated by the two-dimensional kernel density estimation. If more constraints are included, we found that we lose the convexity of the polynomial terms in (4). As we mentioned before, we need a better criteria to preserve the convexity of the solutions.

In Figure 9, we include the result from a five-dimensional simulation. We also do not show the estimate from the MATLAB \asciifamilyfsolve.m since it is not accurate at all. In this five-dimensional case, the EBE method automatically discards 34 equations (moment constraints). In this case, we suspect that either the maximum entropy solution that accounts for all of the constraints does not exist or the EBE method cannot find the solution. Here, the EBE method just estimates the best fitted solution within the tolerance of $10^{-10}$ by solving 91 out of 125 moment constraints.

7. Summary

In this paper, we introduced a novel equation-by-equation algorithm for solving a system of nonlinear equations arising from the moment constrained maximum entropy problem. Theoretically, we have established the local convergence and provided a sufficient condition for global convergence. Through the convergence analysis, we understood that the method, geometrically, finds the solution by searching along the surface corresponding to one component of the nonlinear equations. Numerically, we have demonstrated its accuracy and efficiency on various examples. In one of the examples, we found that the EBE algorithm produces more accurate solutions compared to the previously developed BFGS-based algorithm which does not use the Hessian information [2, 5]. In this same example, we also found that the EBE is superior compared to two schemes that use the Hessian information, including the current MATLAB built-in solver which uses the trust-region-dogleg algorithm and the classical Newton’s method.

We also found that the proposed EBE algorithm is able to solve a system of 70-310 equations when the maximum entropy solution exists compared to the previously developed BFGS method which was shown to work for a system of size 44-83 equations. On the Kuramoto-Shivashinski example, the EBE method is able to reconstruct the density of a four-dimensional problem accounting up to order-four moments (or 70 constraints). In this case, we showed that the estimate is improved by accounting for 17 additional constraints of order-six moments. For the five-dimensional problem with moments up to order-four, the EBE method reconstructs the solution within the desired precision, $10^{-10}$, by automatically selecting a subset of 91 constraints from the total prescribed 125 constraints induced by moments of up to order-four.

While the automatic constraint selection is a desirable feature since the maximum entropy solutions within the tolerance may not be easily estimated (nor theoretically available), further study is required to fully take advantage of this feature. In particular, an important open problem is to develop a mathematical theory for ordering the constraints since the path of the solution is sensitive to the order of the constraints. Simultaneously, the ordering of the constraints need to preserve the convexity of the polynomials in the exponential term of (4). We should stress that the EBE method is computationally not the most efficient method since it is designed to avoid singularities by tracking along the surface corresponding to one component of the nonlinear equations. Therefore, a more efficient EBE method will be one of future directions.

Appendix A The detailed calculation of the Jacobian of the map $H_{i}$

In this Appendix, we will give the detailed computation for the Jacobian of the map $H_{i}$ in (14) evaluated at $\bm{\mu}^{(i)}$, the solution of ${\bf F}_{i}(\bm{\lambda}_{i},\alpha_{i+1},\ldots,\alpha_{n})={\bf 0}$. Recall that for $\bm{H}_{i}=(\bm{H}_{i,1},H_{i,2})$ in (14),

[TABLE]

where $\bm{g}_{i}:\mathbb{R}^{i-1}\to\mathbb{R}^{i-1}$ is defined as in (15).

To take another derivative of $\bm{H}_{i,1}$ with respect to $\lambda_{j}$. We use the fact that if ${\bf F}_{i-1,\bm{\lambda}_{i-1}}$ is a nonsingular matrix, then

[TABLE]

and the Hessian $\frac{\partial{\bf F}^{*}_{i-1,\bm{\lambda}_{i-1}}}{\partial\lambda_{j}}$ is well-defined, which are the Assumptions 1.2 and 1.3. We can deduce that for $j=1,\ldots,i$,

[TABLE]

Evaluating these two equations at $\bm{\mu}^{(i)}$ and using the fact that ${\bf F}_{i}^{*}:={\bf F}_{i}(\bm{\mu}^{(i)})={\bf 0}$, the second terms in the right-hand-side of (21)-(22) vanish and we have,

[TABLE]

where $\delta_{j,i}$ is one only if $j=i$ and zero otherwise.

Acknowledgments

We thank Rafail Abramov for supplying the two-dimensional density data set for 4. The BFGS code that we used for comparison in 4 was downloaded from [1].