Performance Analysis of Trust Region Subproblem Solvers for Limited-Memory Distributed BFGS Optimization Method

Abstract

The limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimization method performs very efficiently for large-scale problems. A trust region search method generally performs more efficiently and robustly than a line search method, especially when the gradient of the objective function cannot be accurately evaluated. The computational cost of an L-BFGS trust region subproblem (TRS) solver depend mainly on the number of unknown variables (n) and the number of variable shift vectors and gradient change vectors (m) used for Hessian updating, with m << n for large-scale problems. In this paper, we analyze the performances of different methods to solve the L-BFGS TRS. The first method is the direct method using the Newton-Raphson (DNR) method and Cholesky factorization of a dense n × n matrix, the second one is the direct method based on an inverse quadratic (DIQ) interpolation, and the third one is a new method that combines the matrix inversion lemma (MIL) with an approach to update associated matrices and vectors. The MIL approach is applied to reduce the dimension of the original problem with n variables to a new problem with m variables. Instead of directly using expensive matrix-matrix and matrix-vector multiplications to solve the L-BFGS TRS, a more efficient approach is employed to update matrices and vectors iteratively. The L-BFGS TRS solver using the MIL method performs more efficiently than using the DNR method or DIQ method. Testing on a representative suite of problems indicates that the new method can converge to optimal solutions comparable to those obtained using the DNR or DIQ method. Its computational cost represents only a modest overhead over the well-known L-BFGS line-search method but delivers improved stability in the presence of inaccurate gradients. When compared to the solver using the DNR or DIQ method, the new TRS solver can reduce computational cost by a factor proportional to for large-scale problems.

Introduction

Decision-making tools based on optimization procedures have been successfully applied to solve practical problems in a wide range of areas. An optimization problem is generally defined as minimizing (or maximizing) an objective function within a user defined search domain , and subject to some linear or nonlinear constraints, where x is an n-dimensional vector that contains all controllable variables.

In the oil and gas industry, an optimal business development plan requires robust production optimization because of the considerable uncertainty of subsurface reservoir properties and volatile oil prices. Many papers have been published on the topic of robust optimization and their applications to cyclic CO₂ flooding through the Gas-Assisted Gravity Drainage process [1], well placement optimization in geologically complex reservoirs [2], optimal production schedules of smart wells for water flooding [3] and in naturally fractured reservoirs [4], just mentioning a few of them as examples.

Some researchers formulated the optimization problem under uncertainty as a single objective optimization problem, e.g., only maximizing the mean of net present value (NPV) for simplicity by neglecting the associated risk. However, it is recommended to formulate robust optimization as a bi-objective optimization problem for consistency and completeness. A bi-objective optimization problem is generally defined as minimizing (or maximizing) two different objective functions, and , within a user-defined search domain , and subject to some linear or nonlinear constraints. For example, we may maximize the mean value, denoted by , and minimize the standard deviation, denoted by , of NPV. For a bi-objective optimization problem, the optimal solutions are defined as the Pareto optimal solutions (or Pareto front). It is a very challenging task to find multiple optimal solutions on the Pareto front [5, 6].

It is also a very challenging task to properly characterize the uncertainty of reservoir properties (e.g., porosity and permeability in each grid-block) and reliably quantify the ensuing uncertainty of production forecasts (e.g., production rates of oil, gas, and water phases) by conditioning to historical production data and 4D seismic data [7, 8], which requires generating multiple conditional realizations by minimizing a properly defined objective function, e.g., using the randomized maximum likelihood (RML) method [9].

When an adjoint-based gradient of the objective function is available [10–12], we may apply a gradient-based optimization method, e.g., the Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimization method [13, 14]. However, the adjoint-based gradient is not available for many commercial simulators and not available for integer type or discrete variables (e.g., well locations defined by grid-block indices). In such a case, we must apply a derivative-free optimization (DFO) method [15–18]. For problems with smooth objective functions and continuous variables, model-based DFO methods using either radial-basis function [19] or quadratic model [16, 20, 21] performs more efficiently than other DFO methods such as direct pattern search methods [22] and stochastic search methods [23].

Traditional optimization methods only locate a single optimal solution, and they are referred to as single-thread optimization methods in this paper. It is unacceptably expensive to use a single-thread optimization method to locate multiple optimal solutions on the Pareto front or generate multiple RML samples [24]. To overcome the limitations of single-thread optimization methods, Gao et al. [17] developed a local-search distributed Gauss-Newton (DGN) DFO method to find multiple best matches concurrently. Later, Chen et al. [25] modified the local-search DGN optimization method and generalized the DGN optimizer for global search. They also integrated the global-search DGN optimization method with the RML method to generate multiple RML samples in parallel.

Because multiple search points are generated in each iteration, finally, multiple optimal solutions can be found in parallel. These distributed optimization methods are referred to as multiple-thread optimization methods. Both the local- and global-search DGN optimization methods are only applicable to a specific optimization problem, i.e., history matching problem or least-squares optimization problem, of which the objective function can be expressed as a form of least-squares, , but they cannot be applied to other type or generic optimization problems where the objective function cannot be expressed as a form of least-squares. Furthermore, the DGN optimization methods may become less efficient for history matching problems when both the number of variables (n) and the number of observed data (m) are large (e.g., in the order of thousands or more).

Recently, Gao et al. [18] developed a well-paralleled distributed quasi-Newton (DQN) DFO method for generic optimization problems by introducing a generalized form of the objective function, . Here, x represents the vector of controllable variables or model parameters to be optimized (explicit variables) and y(x) denotes the vector of simulated responses (implicit variables) of a reservoir using explicit variables x. Using the generalized form of the objective function, the gradient of the objective function can be evaluated analytically by,In Eq. 1, is the transpose of the sensitivity matrix. Using the generalized expression of the objective function F(x,y), different objective functions can be evaluated using the same y(x) that is simulated from the same reservoir model, e.g., for multi-objective optimization problems. A specific case of the generalized expression is , where the transpose of the sensitivity matrix becomes the gradient of the objective function.

The DQN optimization algorithm runs N_e optimization threads in parallel. In the k-th iteration, there are non-converged threads. The DQN optimizer generates a new search points for each non-converged thread, where the search step is solved from a trust region subproblem (TRS). The simulation cases are submitted to a high-performance computing (HPC) cluster in parallel, to generate corresponding simulated responses . Then, we evaluate the objective function and its associated partial derivatives and . If the new search point improves the objective function, i.e., , then we update and .

All simulation results generated by all DQN optimization threads in previous and current iterations are recorded in a training data set, and they are shared among all DQN optimization threads. The sensitivity matrix is approximated using a modified QR-method proposed by Gao, et al. [18] by linear interpolation of training data points that are closest to x^(i,k+1). Then, we approximate the gradient using Eq. 1. Finally, the Hessian H^(i,k+1) for each thread is updated using either the BFGS method or the symmetric rank-1 (SR1) method [14]. For simplicity, we will drop the superscript “i” (the optimization thread index) in the following discussions.

Both DGN and DQN optimization methods approximate the objective function by a quadratic model, and they are designed for problems with smooth objective function and continuous variables. Although their convergence is not guaranteed for problems with integer type variables (e.g., well location optimization), if those integers can be treated as truncated continuous variables, our numerical tests indicate that these distributed optimization methods can improve the objective function significantly and locate multiple suboptimal solutions for problems with integer type variables in only a few iterations.

The number of variables for real-world problems may vary from a few to thousands or even more, depending on the problem and the parameterization techniques employed. For example, the number of variables could be in the order of millions if we do not apply any parameter reduction techniques to reduce the dimension of some history matching problems (e.g., to tune permeability and porosity in each grid block). Because reservoir properties are generally correlated with each other with long correlation lengths, we may reduce the number of variables to be tuned to only a few hundred, e.g., using the principal component analysis (PCA) or other parameter reduction techniques [26].

In this paper, our focus is on performance analysis of different methods to solve the TRS formulated with the limited-memory BFGS (L-BFGS) Hessian updating method for unconstrained optimization problems. In the future, we will further integrate the new TRS solver with our newly proposed limited-memory distributed BFGS (L-DBFGS) optimization algorithm and then apply the new optimizer to some realistic oil/gas field optimization cases and benchmark its overall performance against other distributed optimization methods such as those only using the gradient [27]. We will also continue our investigation in the future for constrained optimization problems (e.g., variables with lower- and/or upper-bounds and problems with nonlinear constraints). Gao et al. [18] applied the popular Newton-Raphson method [28] to directly solve the TRS using matrix factorization, the DNR method in short. The DNR method is quite expensive when it is applied to solve TRSs of a distributed optimization method, especially for large-scale problems with thousands or more variables [29]. This paper follows the ideas and concepts presented in the book “Matrix Computation” [30]. Flops are used to quantify the volume of work associated with a computation, a count for floating-point operations of add, subtract, multiply, or divide. Computational cost (flops) of some commonly used algebraic operations and numerical methods are summarized in Table 1 for reference.

For completeness, we discuss the compact representation of the L-BFGS Hessian updating formulation [31] and the algorithm to directly update the Hessian in the next section directly. In the third section, we present three different methods to solve the L-BFGS TRS: the DNR method, the direct method using inverse quadratic (DIQ) interpolation approach proposed by Gao et al. [32], and the technique using matrix inversion lemma (MIL) together with an efficient matrix updating algorithm. Some numerical tests and performance comparisons are discussed in the fourth section. We finally draw some conclusions in the last section.

The Limited-Memory Hessian Updating Formulation

Let x^(k) be the best solution obtained in the current (k-th) iteration, and f^(k), g^(k), and H^(k), respectively, the objective function, its gradient and Hessian evaluated at x^(k). The Hessian H^(k+1) is updated using the BFGS formulation as follows,In Eq. 2, and . The Hessian H^(k+1) updated using Eq. 2 is guaranteed positive definite if the condition is satisfied where c₃>0 is a small positive number.

Compact Representation

To save memory usage and computational cost, the L-BFGS Hessian updating method limits the maximum history size or the maximum number of pairs (s^(k), z^(k)) used for the BFGS Hessian updating to L_M > 1. The recommended value of L_M ranges from 5 to 20, using smaller number for problems with more variables.

Let 1 ≤ l_k ≤ L_M denote the number of pairs of variable shift vectors and gradient change vectors, (s^(j), z^(j)) for j = 1,2,...l_k, used to update the Hessian using the L-BFGS method in the k-th iteration. Let be the n × l_k variable shift matrix and the n × l_k gradient change matrix. Both matrices S^(k) and Z^(k) are updated iteration by iteration. Let m=2l_k and is an n × m matrix. Let and . We decompose B^(k) into three parts: the strictly lower triangular part L^(k), the strictly upper triangular part U^(k), and the diagonal part E^(k) that contains the main diagonals of B^(k), i.e.,The Hessian can be updated using the compact form of the L-BFGS Hessian updating formulation [31, 33–35] as follows,In Eq. 4, I_n is an n × n identity matrix and W^(k) is an m×m symmetric matrix defined as,In Eq. 4 and Eq. 5, is a scaling factor. It is required that a^(k) > a_cr > 0 using Eq. 4 and Eq. 5 to update the Hessian, where a_cr is a small positive number.

The Algorithm to Directly Update the Hessian Using the L-BFGS Compact Representation

The computational cost to directly update the Hessian using the L-BFGS method as described in Algorithm-1 mainly depends on the number of variables (n) and the number of vectors used to update the Hessian (m = 2l_k). We should reemphasize that m is updated iteratively but limited to m ≤ 2L_k. The computational cost is (flops), which includes the following seven operations: 1) computing ‖s‖, ‖z‖, and in step 1 (6n flops); 2) computing A^(k) and B^(k) in step 5 (m²n flops); 3) forming the matrix in step 6 (m² flops); 4) computing W^(k) in step 7 (O(m³) flops); 5) compute U^(k) in step 8 (2m²n flops); 6) compute T^(k) in step 9 ( flops); and 7) updating in step 10 (n flops).

Solving the L-BFGS Trust Region Subproblem

The Trust Region Subproblem

In the neighborhood of the best solution x^(k), the objective function f(x) can be approximated by a quadratic function of ,If the Hessian is positive definite, then has a unique minimum s^*(k), which is the solution of ,When a line search strategy is applied, we accept the full search step by setting if it improves the objective function sufficiently (e.g., satisfying the two Wolfe conditions). Otherwise, we can find a better search step size such that improves the objective function sufficiently. If the gradient can be accurately evaluated, it has theoretically proved that a line search strategy can converge [14]. However, for most real-world optimization problems, the gradient cannot be accurately evaluated, and as discussed by Gao et al. [16], a trust region search strategy performs more robustly and more efficiently than a line search strategy.

The trust region search step for the next iteration is the global minimum of the quadratic model within a ball-shaped trust region with radius which is solved from the following trust region subproblem (TRS) [28],If the solution given by Eq. 7 satisfies , then is the solution of the TRS defined in Eq. 8. Otherwise, the solution of the TRS must lie on the boundary of , and we should solve the Lagrange multiplier together with the search step from the following nonlinear TRS equations,For the trivial case with , we have and . The solution of the TRS defined in Eq. 8 is and if or and otherwise.

We may apply the popular DNR method to solve the TRS when is a dense matrix. However, the DNR method is quite expensive for large-scale problems and we should seek and apply a more efficient method to solve the TRS.

Solving the L-BFGS TRS Directly Using the Newton-Raphson Method

The DNR method solves the TRS defined in Eq. 8 using the Cholesky decomposition. It applies the Newton-Raphson method to directly solve the nonlinear TRS equation [28], , iteratively, which requires computing the first order derivative where v is solved from together with the Cholesky factorization of or LU-decomposition.

Given the trust region size , threshold of convergence , maximum number of iterations allowed , Hessian and gradient evaluated at the current best solution , both the Lagrange multiplier and the trust region search step can be solved from the TRS defined in Eq. 8, using the DNR method as summarized in Algorithm-2.

Let denote the number of iterations required for the DNR TRS solver to converge. The total computational cost (flops) to solve the TRS using Algorithm-2 is, , including the following three operations: (1) computing and Cholesky decomposition in step 2 and step 7(b) (flops); (2) solving from , from , and from in step 3 and step 7(c) ( flops); (3) computing and in step 4 and step 7(d) ( flops).

The total computational cost (flops) used to solve the L-BFGS TRS using the DNR method (Algorithm-2) together with the L-BFGS Hessian updating method (Algorithm-1) is,Our numerical results indicate that the DNR TRS solver may fail when tested on the well-known Rosenbrock function, especially for problems with large . The root cause for failure of convergence using the DNR method is the same as in the GNTRS solver using the traditional Newton-Raphson method as discussed by Gao et al. [36]. Very small value of may result in a very large search step and thus result in failure of converging. Gao et al. [36] proposed integrating the DNR method with a bisection line search to overcome this issue.

The Inverse Quadratic Model Interpolation Method to Directly Solve the L-BFGS TRS

Instead of applying the Newton-Raphson method which requires evaluating , Gao et al. [32] proposed a method to directly solve the TRS using inverse quadratic model interpolation (called the DIQ method), i.e., approximating by the following inverse quadratic function,Both coefficients of and in Eq. 12 can be determined by interpolating the values of evaluated at two different points and with .

In the first iteration, we set and accept as the solution if . Otherwise, holds and we set and holds if it is not accepted as the solution. In the following iteration, either or will be updated accordingly such that the two conditions and always hold.

Let and the following solution of will be used as the trial search point for the next iteration,We either update if or update if iteratively until convergence. We accept as the desired solution and terminate the iterative process either when , the tolerance for convergence, or when the number of iterations used to solve the TRS () reaches the user specified maximum iteration number (), i.e.,.

Given the scaling factor , trust region size , threshold of convergence , maximum number of iterations allowed , Hessian and gradient , both the Lagrange multiplier and the trust region search step can be solved from the TRS defined in Eq. 8, using the DIQ method as summarized in Algorithm-3.

Algorithm-3: Solving the L-BFGS TRS Using the DIQ Method

Let denote the number of iterations required for the DIQ TRS solver to converge. The total computational cost to solve the L-BFGS TRS using Algorithm-3 together with Algorithm-1 is,Generally, the DIQ method converges faster than the DNR method and it performs more efficiently and robustly than the DNR method, especially for large scale problems. Both the DNR method and the DIQ method become quite expensive for large-scale problems.

Using Matrix Inversion Lemma (MIL) to Solve the L-BFGS TRS

To save both memory usage and computational cost, Gao, et al. [32, 36, 37] proposed an efficient algorithm to solve the Gauss-Newton TRS (GNTRS) for large-scale history matching problems using the matrix inversion lemma (or the Woodbury matrix identity). With appropriate normalization of both parameters and residuals, the Hessian of the objective function for a history matching problem can be approximated by the well-known Gauss-Newton equation as follows,In Eq. 15, is an matrix, the transpose of the sensitivity matrix that is evaluated at the current best solution . Here, denotes the number of observed data to be matched and it is not the same as used in the L-BFGS Hessian updating Eq. 4.

The GNTRS solver proposed by Gao, et al. [32, 36, 37] using the matrix inversion lemma (MIL) has been implemented and integrated with the distributed Gauss-Newton (DGN) optimizer. Because the compact representation of the L-BFGS Hessian updating formula Eq. 4 is similar to Eq. 15, we follow the similar idea as proposed by Gao, et al. [32, 36, 37] to compute by applying the matrix inversion lemma and then solve the L-BFGS trust region search step .Let and . We first solve from,Then, we compute the trust region search step and as,From Eq. 15, we have and Eq. 19 can be rewritten as,We first try . If holds, we accept and as the solution of the TRS defined in Eq. 8. Otherwise, the solution is on the boundary defined in Eq. 10, which can be solved by the Newton-Raphson (NR) method iteratively.

Let and we have where is computed by,In Eq. 21, is solved from,Because is an symmetric and positive definite matrix, it only requires. flops to solve and ) from Eq. 17 and Eq. 22 using the Cholesky factorization.

Given the variable shift matrix and gradient change matrix , we may directly compute the following three matrices, , and , and then form the matrix as,Directly computing the three matrices, , and , requires flops. We also need to compute the vector ( flops) and ( flops). To further reduce the computational cost, the three matrices , and and the vector should be updated iteratively instead.

Updating Matrices and Vectors Used for Solving the L-BFGS TRS

We first initialize and . If the trial search point does not improve the objective function, i.e., , we set and , recompute the sensitivity matrix and the gradient when needed (e.g., using updated training data points), and shrink the trust region size with . We repeatedly generate trial search point until improves the objective function, i.e., .

If , we set , , evaluate , and compute and . In the following iterations, we update and associated matrices and vectors used for solving the L-BFGS TRS for different cases accordingly. We use to indicate whether the new search point improves the objective function () or not ().

The L-BFGS Hessian updating formulation of Eq. 4 requires and where and are small positive numbers. If or , we simply set , , , , , , , , , and , with no updating. However, we need to compute using the new gradient . In the following cases, we assume that both conditions and are satisfied.

Case-1: Sflag^(k) = “False”

Although the best solution does not change when Sflag^(k) = “False”, the sensitivity matrix and thus the gradient of the objective function evaluated at using simulation results and training data points updated in the (current) iteration may be different from those evaluated at the same point but using simulation results and training data points obtained in the (previous) iteration. To use the right gradient, we should replace the last column in with .

Let be the submatrix of by removing its last column , the submatrix of by removing its last column, and the submatrix of by removing its last row and last column. We first compute the following two vectors and and two scalars and , with computational cost flops.

We set , , , , and , and Update as,It is straightforward to derive the following equations to update and accordingly,By definition,Thus, we have,

Case-2: Sflag^(k) = “True” and l_k−1 < L_M

We set and and , compute four -dimensional vectors and and and and five scalars, and and and and , with computational cost flops. Both and are updated by,It is straightforward to derive the following equations to update A^(k), B^(k), and C^(k),We first split the -dimensional vector into two -dimensional sub-vectors and , i.e.,where is composed of the first rows in and is composed of the last rows in .

Because and , thus we have,

Case-3: Sflag^(k) = “True” and l_k−1 = L_M

We set . Let and denote the submatrices of and by removing the first column and from and , respectively. Let and and denote the submatrices of A^(k), B^(k), and C^(k) by removing their first row and first column.

apply Eq. 28 through Eq. 32 to update S^(k), Z^(k), A^(k), B^(k), and C^(k) by simply replacing S^(k−1), Z^(k−1), A^(k−1), B^(k−1), and C^(k−1) with and and and and , respectively. We can also apply Eq. 33 to update by replacing with that is composed of the second row through the -th rows in and replacing with that is composed of the last rows in .

The Algorithm to Update Matrices and Vectors Used for Solving the L-BFGS TRS

Given L_M ≥ 1, k > 1, l_k−1 ≥ 1, the thresholds c₃ > 0 and a_cr > 0, -dimensional vectors g^(k+1), s and z, l_k−1-dimensional vector u^(k−1), n × l_k−1 matrices S^(k−1) and Z^(k−1), l_k−1 × l_k−1 matrices A^(k−1), B^(k−1) and C^(k−1), 2l_k−1 × 2l_k−1 matrices P^(k−1) and , we can update and the -dimensional vector u^(k), n × l_k matrices S^(k) and Z^(k), l_k × l_k matrices A^(k), B^(k), and C^(k), and 2l_k × 2l_k matrices and , using the algorithm as summarized in Algorithm-4.

Algorithm-4: Updating Matrices and Vectors Used for Solving the L-BFGS TRS

The computational cost to update matrices and vectors used for solving the L-BFGS TRS as described in Algorithm-4 is at most (taking Case-2 in step 5 as an example), including the following 4 operations: (1) computing , , and in step 1 ( flops); (2) computing ,,, and in step 5(b) ( flops),; (3) computing and in step 5(c) ( flops).

In this paper, we apply the same idea of reducing the dimension from to using the matrix inversion lemma as presented in the three papers [32, 36, 37]. However, the implementation is quite different. For the DGN optimization method, the matrix (or the transpose of the sensitivity matrix) in Eq. 15 is directly computed. Its dimensions (both the number of variables and the number of observed data ) are fixed, i.e., they remain the same for all iterations. Therefore, we have to directly compute the matrix and the -dimensional vector , with the computational cost flops. Using the algorithm as summarized in Algorithm-4, the computational cost can be further reduced to , by a factor of , 5 to 20 roughly.

The Algorithm to Solve the L-BFGS TRS Using the Matrix Inversion Lemma

Given the trust region size , threshold of convergence , maximum number of iterations allowed , the scaling factor , , matrix , matrices and , -dimensional vector and -dimensional gradient vector evaluated at the current best solution , both the Lagrange multiplier and the trust region search step can be solved from the TRS defined in Eq. 8 using the matrix inversion lemma (MIL) as summarized in Algorithm-5.

Algorithm-5: Solving the L-BFGS TRS Using the Matrix Inversion Lemma

Let denote the number of iterations required for the L-BFGS TRS solver to converge. The total computational cost (flops) to solve the L-BFGS TRS using Algorithm-5 is, , including the following eight operations: (1) computing in step 1 ( flops); (2) computing and Cholesky decomposition in step 4 or step 9(g) ( flops); (3) solving from , and from in step 5 or step 9(h) ( flops); (4) computing , and in step 6 or step 9(i) ( flops); (5) computing in step 9(a) ( flops); (6) Solving from and from in step 9(b) ( flops); and (7) computing , , and in step 9(c) ( flops); and (8) computing and in step 10 ().

The total computational cost (flops) used to solve the L-BFGS TRS using Algorithm-5 together with Algorithm-4 is,It is straightforward to compute the normalized computational cost, or the ratio of computational cost of the DNR method (Algorithm-2) or the DIQ method (Algorithm-3) together with the algorithm that directly updates Hessian (Algorithm-1) over the MIL TRS solver using the matrix inversion lemma (Algorithm-5) and the efficient matrix updating algorithm (Algorithm-4),Because for large-scale problems, Eq. 35 can be further simplified as,For some problems the solution is accepted with . For other problems, it takes roughly iterations for a TRS solver to converge. Figure 1A,B illustrate the plots of vs. for different by fixing and , respectively. As shown in both Figure 1A and Figure 1B, the MIL TRS solver may reduce the computational cost by a factor for .

FIGURE 1

Numerical Validation

We benchmarked the MIL TRS solver against the DNR (or DIQ) TRS solver on two well-known analytic optimization problems using analytical gradients: the Rosenbrock function and the Sphere function defined as follows,The Rosenbrock function defined in Eq. 37 has one local minimum located at and for with the objective function being 4 (approximately), and one global minimum located at for with the objective function being 0. Occasionally, a local-search optimizer may converge to the local minimum.

We implemented different algorithms described in this paper in our in-house C++ optimization library (OptLib). We employ “Armadillo” as our foundational template-based C++ library for linear algebra (http://arma.sourceforge.net/). We ran all test cases reported on Tables 2–5 on a virtual machine computer equipped with an Intel Xeon Platinum Processor at 2.60 MHz and 16.0 GB of RAM. We varied the number of parameters, as powers of two and fixed twice the maximum number of snapshots, . Tables 2–5 summarize the preliminary numerical results, including the value of the objective function, the gradient norm, the number of iterations, and the elapsed CPU time (s). The last column displays the resulting speedup, namely, the ratio of CPU times using the DNR (or DIQ) method over the MIL method. We prescribe the tolerance to stop the L-BFGS TRS solver for all the numerical experiments in this subsection as and the maximum number of iterations as . We set the initial trust region size at ∆ = 0.5.

We noticed that the matrix may become ill-conditioned for a few cases. Indeed, we could observe that the first entries on the diagonal are nearly zero. The same situation occurs for different optimization steps where grows from 0, 2, 4, 6, up to ten. The Cholesky decomposition of such a matrix will fail, i.e., Armadillo throws an exception. We advise replacing the former with an LU-decomposition to handle the pivots, namely, entries on the diagonal of the lower triangular matrix, to be nearly zero. The rest of the algorithm remains the same.

Table 2 compares performance of the DNR method vs. the MIL methods for the Rosenbrock problem encompassing 8, 16, 32, and 48 parameters, respectively. Beyond those problem sizes, the DNR method could not converge to the solution. Indeed, the DNR method becomes relatively slow, i.e., renders very large lambda values that imply tiny shifts. Often, after 2,048 iterations, the objective function is still greater than one. It seems that the DNR method is not robust enough to tackle the Rosenbrock function with more than 48 parameters. It appears that the MIL method slightly outperforms the DNR method for problems with fewer parameters.

Table 3 compares performance of the DIQ method vs. the MIL method. These numerical results confirm that the MIL method outperforms the DIQ method. As expected, the speedup improves with the increase of the problem size n.

Table 4 shows similar results for the Sphere function problem. We could benchmark the DNR method against the MIL method with ranks up to 2,048 parameters. Again, the speedup increases with problem size. Finally, Table 5 benchmarks the DIQ method against the MIL method testing on the Sphere function. The speedup reduces when compared to Table 4 for the same ranks, but we still clearly see that the performance of the MIL solver exceeds its DIQ counterpart. Results shown in Table 4 and Table 5 also confirm that the DIQ method performs better than the DNR method. We believe that the DIQ TRS solver converges faster with smaller number of iterations than the DNR TRS solver (i.e., ).

The computational costs (in terms of CPU time) summarized in Tables 2–5 are the overall costs for all iterations, including the cost of computing both value and gradient of the objective function in addition to the cost of solving the L-BFGS TRS. In contrast, the theoretical analysis results of computational cost for β in Eq. 35 only count for the cost of solving the L-BFGS TRS in just one iteration. Therefore, numerical results listed in Tables 2–5 are not quantitatively in agreement with the theoretically derived β.

Conclusion

References

• 1.

  Ai-MudhafarWJRaoDNSrinivasanS. Robust Optimization of Cyclic CO2 Flooding through the Gas-Assisted Gravity Drainage Process under Geological Uncertainties. J Pet Sci Eng (2018). 166:490–509. 10.1016/j.petrol.2018.03.044

  • CrossRef
  • Google Scholar

• 2.

  AlpakFOJinLRamirezBA. Robust Optimization of Well Placement in Geologically Complex Reservoirs. in: Paper SPE-175106-MS Presented at the SPE Annual Technical Conference and Exhibition; September28–30 2015; Houston, TX (2015).

  • Google Scholar

• 3.

  van EssenGMZandvlietMJDen HofPMJBosgraOHJansenJDet alRobust Optimization of Oil Reservoir Flooding. in: Paper Presented at the IEEE International Conference on Control Applications; Oct 4-6 2006; Munich, Germany (2006). p. 4–6.

  • Google Scholar

• 4.

  LiuZForouzanfarF. Ensemble Clustering for Efficient Robust Optimization of Naturally Fractured Reservoirs. Comput Geosci (2018). 22:283–96. 10.1007/s10596-017-9689-1

  • CrossRef
  • Google Scholar

• 5.

  LiuXReynoldsAC. Augmented Lagrangian Method for Maximizing Expectation and Minimizing Risk for Optimal Well-Control Problems with Nonlinear Constraints. SPE J (2016). 21(5). 10.1007/s10596-017-9689-1

  • CrossRef
  • Google Scholar

• 6.

  LiuXReynoldsAC. Gradient-based Multi-Objective Optimization with Applications to Waterflooding Optimization. Comput Geosci (2016). 20:677–93. 10.1007/s10596-015-9523-6

  • CrossRef
  • Google Scholar

• 7.

  OliverDSReynoldsACLiuN. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge, United Kingdom: Cambridge University Press (2008). 10.1017/cbo9780511535642

  • CrossRef
  • Google Scholar

• 8.

  TarantolaA. Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia, PA: SIAM (2005). 10.1137/1.9780898717921

  • CrossRef
  • Google Scholar

• 9.

  OliverDS. On Conditional Simulation to Inaccurate Data. Math Geol (1996). 28:811–7. 10.1007/bf02066348

  • CrossRef
  • Google Scholar

• 10.

  JansenJD. Adjoint-based Optimization of Multi-phase Flow through Porous Media - A Review. Comput Fluids (2011). 46(1):40–51. 10.1016/j.compfluid.2010.09.039

  • CrossRef
  • Google Scholar

• 11.

  LiRReynoldsACOliverDS. History Matching of Three-phase Flow Production Data. SPE J (2003). 8(4):328–40. 10.2118/87336-pa

  • CrossRef
  • Google Scholar

• 12.

  SarmaPDurlofskyLAzizK. Implementation of Adjoint Solution for Optimal Control of Smart Wells. in: Paper SPE 92864 Presented at the SPE Reservoir Simulation Symposium Held in the Woodlands; Jan.-2 Feb; Woodlands, TX (2005). p. 31.

  • Google Scholar

• 13.

  GaoGReynoldsAC. An Improved Implementation of the LBFGS Algorithm for Automatic History Matching. SPE J (2006). 11(1):5–17. 10.2118/90058-pa

  • CrossRef
  • Google Scholar

• 14.

  NocedalJWrightSJ. Numerical Optimization. New York City: Springer (1999).

  • Google Scholar

• 15.

  ChenCet alAssisted History Matching Using Three Derivative-free Optimization Algorithms. in: Paper SPE-154112-MS Presented at the SPE Europec/EAGE Annual Conference Held in Copenhagen; June 4–7 2012; Copenhagen, Denmark (2012). p. 4–7.

  • Google Scholar

• 16.

  GaoGVinkJCChenCAlpakFODuK. A Parallelized and Hybrid Data-Integration Algorithm for History Matching of Geologically Complex Reservoirs. SPE J (2016). 21(6):2155–74. 10.2118/175039-pa

  • CrossRef
  • Google Scholar

• 17.

  GaoGVinkJCChenCEl KhamraYTarrahiMet alDistributed Gauss-Newton Optimization Method for History Matching Problems with Multiple Best Matches. Comput Geosciences (2017). 21(5-6):1325–42. 10.1007/s10596-017-9657-9

  • CrossRef
  • Google Scholar

• 18.

  GaoGWangYVinkJWellsTSaafFet alDistributed Quasi-Newton Derivative-free Optimization Method for Optimization Problems with Multiple Local Optima. in: 17th European Conference on the Mathematics of Oil Recovery ECMOR XVII; September 2020; Edinburgh, United Kingdom (2020). p. 14–7.

  • Google Scholar

• 19.

  WildSM. Derivative Free Optimization Algorithms for Computationally Expensive Functions. [PhD thesis]. Ithaca (new york): Cornell University (2009).

  • Google Scholar

• 20.

  PowellMJD. Least Frobenius Norm Updating of Quadratic Models that Satisfy Interpolation Conditions. Math Programming (2004). 100(1):183–215. 10.1007/s10107-003-0490-7

  • CrossRef
  • Google Scholar

• 21.

  ZhaoHLiGReynoldsACYaoJ. Large-scale History Matching with Quadratic Interpolation Models. Comput Geosci (2012). 17:117–38. 10.1007/s10596-012-9320-4

  • CrossRef
  • Google Scholar

• 22.

  AudetCDennisJEJr.Mesh Adaptive Direct Search Algorithms for Constrained Optimization. SIAM J Optim (2006). 17:188–217. 10.1137/040603371

  • CrossRef
  • Google Scholar

• 23.

  SpallJC. Introduction to Stochastic Search and Optimization. Hoboken, NJ: John Wiley & Sons (2003). 10.1002/0471722138

  • CrossRef
  • Google Scholar

• 24.

  ChenCet alEUR Assessment of Unconventional Assets Using Parallelized History Matching Workflow Together with RML Method. In: Paper Presented the 2016 Unconventional Resources Technology Conference; August 1–3 2016; San Antonio, TX(2016).

  • Google Scholar

• 25.

  ChenCGaoGLiRCaoRChenTVinkJCet alGlobal-Search Distributed-Gauss-Newton Optimization Method and its Integration with the Randomized-Maximum-Likelihood Method for Uncertainty Quantification of Reservoir Performance. SPE J (2018). 23(05):1496–517. 10.2118/182639-pa

  • CrossRef
  • Google Scholar

• 26.

  ChenCGaoGRamirezBAVinkJCGirardiAM. Assisted History Matching of Channelized Models by Use of Pluri-Principal-Component Analysis. SPE J (2016). 21(05):1793–812. 10.2118/173192-pa

  • CrossRef
  • Google Scholar

• 27.

  ArmackiAJakoveticDKrejicNKrklec JerinkicNet alDistributed Trust-Region Method with First Order Models. in Paper Presented at the IEEE EUROCON-2019, 18th International Conference on Smart TechnologiesJuly 1–14 2019Novi Sad, Serbia (2019). p. 1–4.

  • Google Scholar

• 28.

  MoréJJSorensenDC. Computing a Trust Region Step. SIAM J Sci Stat Comput (1983). 4(3):553–72. 10.1137/0904038

  • CrossRef
  • Google Scholar

• 29.

  MohamedAW. Solving Large-Scale Global Optimization Problems Using Enhanced Adaptive Differential Evolution Algorithm. Complex Intell Syst (2017). 3:205–31. 10.1007/s40747-017-0041-0

  • CrossRef
  • Google Scholar

• 30.

  GolubGHVan LoanCF. Matrix Computations. 4th ed. The Johns Hopkins University Press (2013).

  • Google Scholar

• 31.

  BrustJJLeyfferSPetraCG. Compact Representations of BFGS Matrices.Preprint ANL/MCS-P9279-0120. Lemont, IL: ARGONNE NATIONAL LABORATORY (2020).

  • Google Scholar

• 32.

  GaoGJiangHVinkJCvan HagenPPHWellsTJet alPerformance Enhancement of Gauss-Newton Trust-Region Solver for Distributed Gauss-Newton Optimization Method. Comput Geosci (2020). 24(2):837–52. 10.1007/s10596-019-09830-x

  • CrossRef
  • Google Scholar

• 33.

  BurdakovOGongLZikrinSYuanY. On Efficiently Combining Limited-Memory and Trust-Region Techniques. Math Prog Comp (2017). 9:101–34. 10.1007/s12532-016-0109-7

  • CrossRef
  • Google Scholar

• 34.

  BurkeJVWiegmannAXuL. Limited Memory BFGS Updating in a Trust-Region framework Tech. Rep. Seattle, WA: University of Washington (2008).

  • Google Scholar

• 35.

  ByrdRHNocedalJSchnabelRB. Representations of Quasi-Newton Matrices and Their Use in Limited Memory Methods. Math Programming (1994). 63:129–56. 10.1007/bf01582063

  • CrossRef
  • Google Scholar

• 36.

  GaoGJiangHvan HagenPVinkJCWellsT. A Gauss-Newton Trust Region Solver for Large Scale History Matching Problems. SPE J (2017b). 22(6):1999–2011. 10.2118/182602-pa

  • CrossRef
  • Google Scholar

• 37.

  GaoGSaafFVinkJKrymskayaMWellsTet alGauss-Newton Trust Region Search Optimization Method for Ill-Conditioned Least Squares Problem. in: ECMOR XVII 17th European Conference on the Mathematics of Oil Recovery; 14 September 2020; Edinburgh, UK (2020). p. 14–17.

  • Google Scholar