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

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 ${\mathit{n}}^{\mathbf{2}}/\mathit{m}$ 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 $\mathit{f}\left(\mathit{x}\right)$ within a user defined search domain $\mathit{x}\in \mathbf{\Omega }$, 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, ${\mathit{f}}_{\mathrm{1}}\left(\mathit{x}\right)$ and ${\mathit{f}}_2\left(\mathit{x}\right)$, within a user-defined search domain $\mathit{x}\in \mathbf{\Omega }$, and subject to some linear or nonlinear constraints. For example, we may maximize the mean value, denoted by ${\mathit{f}}_1\left(\mathit{x}\right)$, and minimize the standard deviation, denoted by ${\mathit{f}}_2\left(\mathit{x}\right)$, 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, $\mathit{f}\left(\mathit{x}\right)=\frac{1}{2}\mathit{\alpha }{\mathit{x}}^{\mathit{T}}\mathit{x}+\frac{1}{2}{\mathit{y}}^{\mathit{T}}\left(\mathit{x}\right)\mathit{y}\left(\mathit{x}\right)$, 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, $\mathit{F}\left(\mathit{x},\mathit{y}\right)=\mathit{f}\left(\mathit{x}\right)$. 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,

$\mathit{g}\left(\mathit{x}\right)={\nabla }_{\mathit{x}}\mathit{F}\left(\mathit{x},\mathit{y}\right)+{\mathit{J}}^{\mathit{T}}\left(\mathit{x}\right){\nabla }_{\mathit{y}}\mathit{F}\left(\mathit{x},\mathit{y}\right)(1)$

In Eq. 1, ${\mathit{J}}^{\mathit{T}}\left(\mathit{x}\right)={\nabla }_{\mathit{x}}{\mathit{y}}^{\mathit{T}}\left(\mathit{x}\right)$ 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 $\mathit{F}\left(\mathit{x},\mathit{y}\right)=\mathit{y}=\mathit{f}\left(\mathit{x}\right)$, 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 ${\mathit{N}}_{\mathit{LC}}^{\left(\mathit{k}\right)}\le {\mathit{N}}_{\mathit{e}}$ non-converged threads. The DQN optimizer generates a new search points ${\mathit{x}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}={\mathit{x}}^{\left(\mathit{i},\mathit{k}\right)}+{\mathit{s}}^{\left(\mathit{i},\mathit{k}\right)}$ for each non-converged thread, where the search step ${\mathit{s}}^{\left(\mathit{i},\mathit{k}\right)}$ is solved from a trust region subproblem (TRS). The ${\mathit{N}}_{\mathit{LC}}^{\left(\mathit{k}\right)}$ simulation cases are submitted to a high-performance computing (HPC) cluster in parallel, to generate corresponding simulated responses ${\mathit{y}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}=\mathit{y}\left({\mathit{x}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}\right)$. Then, we evaluate the objective function ${\mathit{f}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}=\mathit{F}\left({\mathit{x}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)},{\mathit{y}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}\right)$ and its associated partial derivatives ${\nabla }_{\mathit{x}}{\mathit{F}}^{\left(\mathit{i},\mathit{k}\right)}$ and ${\nabla }_{\mathit{y}}{\mathit{F}}^{\left(\mathit{i},\mathit{k}\right)}$. If the new search point ${\mathit{x}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}$ improves the objective function, i.e., ${\mathit{f}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}<{\mathit{f}}^{\left(\mathit{i},\mathit{k}\right)}=\mathit{f}\left({\mathit{x}}^{\left(\mathit{i},\mathit{k}\right)}\right)$, then we update ${\mathit{x}}^{\left(\mathit{i},\mathit{k}+1\right)}={\mathit{x}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}$ and ${\mathit{f}}^{\left(\mathit{i},\mathit{k}+1\right)}={\mathit{f}}_{\mathit{T}}^{\left(\mathit{i},\mathit{k}\right)}$.

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 ${\mathit{J}}^{\left(\mathit{i},\mathit{k}+1\right)}=\mathit{J}\left({\mathit{x}}^{\left(\mathit{i},\mathit{k}+1\right)}\right)$ 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 ${\mathit{g}}^{\left(\mathit{i},\mathit{k}+1\right)}=\mathit{g}\left({\mathit{x}}^{\left(\mathit{i},\mathit{k}+1\right)}\right)$ 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 ${\mathit{N}}_{\mathit{e}}$ 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.

TABLE 1

TABLE 1. A summary of computational costs of commonly used algebraic operations.

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,

${\mathit{H}}^{\left(\mathit{k}\mathrm{+1}\right)}={\mathit{H}}^{\left(\mathit{k}\right)}+\frac{{\mathit{z}}^{\left(\mathit{k}\right)}{\left[{\mathit{z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}}{{\left[{\mathit{z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{s}}^{\left(\mathit{k}\right)}}-\frac{{\mathit{H}}^{\left(\mathit{k}\right)}{\mathit{s}}^{\left(\mathit{k}\right)}{\left[{\mathit{s}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{H}}^{\left(\mathit{k}\right)}}{{\left[{\mathit{s}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{H}}^{\left(\mathit{k}\right)}{\mathit{s}}^{\left(\mathit{k}\right)}}(2)$

In Eq. 2, ${\mathit{s}}^{\left(\mathit{k}\right)}={\mathit{x}}^{\left(\mathit{k}+1\right)}-{\mathit{x}}^{\left(\mathit{k}\right)}$ and ${\mathit{z}}^{\left(\mathit{k}\right)}={\mathit{g}}^{\left(\mathit{k}+1\right)}-{\mathit{g}}^{\left(\mathit{k}\right)}$. The Hessian H^(k+1) updated using Eq. 2 is guaranteed positive definite if the condition ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\left[{\mathit{z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{s}}^{\left(\mathit{k}\right)}>{\mathit{c}}_3{\mathit{z}}^{\left(\mathit{k}\right)}{\mathit{s}}^{\left(\mathit{k}\right)}$ 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 ${\mathit{S}}^{\left(\mathit{k}\right)}=\left[{\mathit{s}}^{\left(\mathrm{1}\right)},{\mathit{s}}^{\left(\mathrm{2}\right)},\dots {\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}}\right)}\right]$ be the n × l_k variable shift matrix and ${\mathit{Z}}^{\left(\mathit{k}\right)}=\left[{\mathit{z}}^{\left(1\right)},{\mathit{z}}^{\left(2\right)},\dots {\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}\right]$ the n × l_k gradient change matrix. Both matrices S^(k) and Z^(k) are updated iteration by iteration. Let m=2l_k and ${\mathit{V}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{\left(\mathit{k}\right)},{\mathit{Z}}^{\left(\mathit{k}\right)}\right]$ is an n × m matrix. Let ${\mathit{A}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{S}}^{\left(\mathit{k}\right)}$ and ${\mathit{B}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{\left(\mathit{k}\right)}$. 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.,

${\mathit{B}}^{\left(\mathit{k}\right)}={\mathit{L}}^{\left(\mathit{k}\right)}+{\mathit{E}}^{\left(\mathit{k}\right)}+{\mathit{U}}^{\left(\mathit{k}\right)}(3)$

The Hessian can be updated using the compact form of the L-BFGS Hessian updating formulation [31, 33–35] as follows,

${\mathit{H}}^{\left(\mathit{k}+1\right)}={\mathit{\alpha }}^{\left(\mathit{k}\right)}{\mathit{I}}_{\mathit{n}}-{\mathit{V}}^{\left(\mathit{k}\right)}{\mathit{W}}^{\left(\mathit{k}\right)}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}(4)$

In Eq. 4, I_n is an n × n identity matrix and W^(k) is an m×m symmetric matrix defined as,

${\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}=\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}}\left\{\begin{array}{c}{\mathit{A}}^{\left(\mathit{k}\right)} {\mathit{L}}^{\left(\mathit{k}\right)}\\ {\left[{\mathit{L}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} -{\mathit{\alpha }}^{\left(\mathit{k}\right)}{\mathit{E}}^{\left(\mathit{k}\right)}\end{array}\right\}(5)$

In Eq. 4 and Eq. 5, ${\mathit{\alpha }}^{\left(\mathit{k}\right)}=\frac{{\Vert {\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}\Vert }^2}{{\left[{\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}\right]}^{\mathit{T}}{\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}}\right)}}$ 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

Given L_M ≥ 1, k ≥ 1, l_k−1 ≥ 1, the thresholds c₃ > 0 and a_cr > 0, search step $\mathit{s}={\mathit{x}}^{\left(\mathit{k}+1\right)}-{\mathit{x}}^{\left(\mathit{k}\right)}$ and gradient change $\mathit{z}={\mathit{g}}^{\left(\mathit{k}+1\right)}-{\mathit{g}}^{\left(\mathit{k}\right)}$, the two n × l_k−1 matrices, S^(k−1) and Z^(k−1), we can update S^(k), Z^(k) and directly compute the Hessian H^(k+1) using the compact representation of the L-BFGS Hessian updating Eq. 4, as summarized in Algorithm-1.

Algorithm-1: Updating the Hessian directly using the L-BFGS compact representation

1. Compute ||s||, ||z||, ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{s}$and${\mathit{\alpha }}^{\left(\mathit{k}\right)}=\frac{\Vert {\mathit{z}}^2\Vert }{{\mathit{\epsilon }}^{\left(\mathit{k}\right)}}$ if ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}>{\mathit{c}}_3\Vert \mathit{s}\Vert \mathit{z}\Vert$;

2. If${\mathit{\epsilon }}^{\left(\mathit{k}\right)}\le {\mathit{c}}_3\Vert \mathit{s}\Vert \mathit{z}\Vert$or${\mathit{\alpha }}^{\left(\mathit{k}\right)}\le {\mathit{\alpha }}_{\mathit{cr}}$, set${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}$, ${\mathit{S}}^{\left(\mathit{k}\right)}={\mathit{S}}^{\left(\mathit{k}-1\right)}$, ${\mathit{Z}}^{\left(\mathit{k}\right)}={\mathit{Z}}^{\left(\mathit{k}-1\right)}$, and goto step 4;

3. Otherwise,

a. Set${\mathit{l}}_{\mathit{k}}={\mathit{L}}_{\mathit{M}}$and remove the first column fromS^(k−1)andZ^(k−1)if${\mathit{l}}_{\mathit{k}-1}={\mathit{L}}_{\mathit{M}}$;else set${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}+1$;

b. Update${\mathit{S}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{(\mathit{k}-1)}\mathit{ s}\right]$and${\mathit{Z}}^{\left(\mathit{k}\right)}=\left[{\mathit{Z}}^{(\mathit{k}-1)}\mathit{ z}\right]$;

4. Form${\mathit{V}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{\left(\mathit{k}\right)}\mathit{ }{\mathit{Z}}^{\left(\mathit{k}\right)}\right]$;

5. Compute${\mathit{A}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{S}}^{\left(\mathit{k}\right)}$and${\mathit{B}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{\left(\mathit{k}\right)}$;

6. Form the matrix${\mathit{W}}_{\mathit{I}}^{\left(\mathit{k}\right)}={\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}$usingEq. 5;

7. Compute${\mathit{W}}^{\left(\mathit{k}\right)}={\left[{\mathit{W}}_{\mathit{I}}^{\left(\mathit{k}\right)}\right]}^{-1}$using Cholesky factorization;

8. Compute${\mathit{U}}^{\left(\mathit{k}\right)}={\mathit{W}}^{\left(\mathit{k}\right)}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}$;

9. Compute${\mathit{T}}^{\left(\mathit{k}\right)}={\mathit{V}}^{\left(\mathit{k}\right)}{\mathit{U}}^{\left(\mathit{k}\right)}$;

10. Update${\mathit{H}}^{(\mathit{k}+1)}={\mathit{\alpha }}^{\left(\mathit{k}\right)}{\mathit{I}}_{\mathit{n}}-{\mathit{T}}^{\left(\mathit{k}\right)}$.

11. Stop.

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 ${\mathit{c}}_1=2\mathit{m}{\mathit{n}}^2+\left(7+3{\mathit{m}}^2\right)\mathit{n}+\mathit{O}\left({\mathit{m}}^3\right)$ (flops), which includes the following seven operations: 1) computing ‖s‖, ‖z‖, and ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{s}$ in step 1 (6n flops); 2) computing A^(k) and B^(k) in step 5 (m²n flops); 3) forming the matrix ${\mathit{W}}_{\mathit{I}}^{\left(\mathit{k}\right)}$ 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 ($2\mathit{m}{\mathit{n}}^2$ flops); and 7) updating ${\mathit{H}}^{(\mathit{k}+1)}$ 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 $\mathit{s}=\mathit{x}-{\mathit{x}}^{\left(\mathit{k}\right)}$,

${\mathit{q}}^{\left(\mathit{k}\right)}\left(\mathit{s}\right)={\mathit{f}}^{\left(\mathit{k}\right)}+{\left[{\mathit{g}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{s}+\frac{1}{2}{\mathit{s}}^{\mathit{T}}{\mathit{H}}^{\left(\mathit{k}\right)}\mathit{s}(6)$

If the Hessian ${\mathit{H}}^{\left(\mathit{k}\right)}$ is positive definite, then ${\mathit{q}}^{\left(\mathit{k}\right)}\left(\mathit{s}\right)$ has a unique minimum s^*(k), which is the solution of ${\mathit{H}}^{\left(\mathit{k}\right)}\mathit{s}=-{\mathit{g}}^{\left(\mathit{k}\right)}$,

${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}=-{\left[{\mathit{H}}^{\left(\mathit{k}\right)}\right]}^{-1}{\mathit{g}}^{\left(\mathit{k}\right)}(7)$

When a line search strategy is applied, we accept the full search step ${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ by setting ${\mathit{x}}^{(\mathit{k}+1)}={\mathit{x}}^{\left(\mathit{k}\right)}+{\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ if it improves the objective function sufficiently (e.g., satisfying the two Wolfe conditions). Otherwise, we can find a better search step size $0<{\mathit{\gamma }}^{\left(\mathit{k}\right)}\le 1$ such that ${\mathit{x}}^{(\mathit{k}+1)}={\mathit{x}}^{\left(\mathit{k}\right)}+{\mathit{\gamma }}^{\left(\mathit{k}\right)}{\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ improves the objective function sufficiently. If the gradient ${\mathit{g}}^{\left(\mathit{k}\right)}$ 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 ${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ for the next iteration is the global minimum of the quadratic model ${\mathit{q}}^{\left(\mathit{k}\right)}\left(\mathit{s}\right)$ within a ball-shaped trust region with radius ${\mathrm{\Delta }}^{\left(\mathit{k}\right)}\ge {\mathrm{\Delta }}_{\mathit{min}}>0$ which is solved from the following trust region subproblem (TRS) [28],

$\underset{\mathit{s}}{\min }{\mathit{q}}^{\left(\mathit{k}\right)}\left(\mathit{s}\right)={\mathit{f}}^{\left(\mathit{k}\right)}+{\left[{\mathit{g}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{s}+\frac{1}{2}{\mathit{s}}^{\mathit{T}}{\mathit{H}}^{\left(\mathit{k}\right)}\mathit{s}\text{subject}\text{to}\mathit{s}{}_2\le {\mathrm{\Delta }}^{\left(\mathit{k}\right)}(8)$

If the solution ${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ given by Eq. 7 satisfies ${\Vert {\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}\Vert }_2\le {\mathrm{\Delta }}^{\left(\mathit{k}\right)}$, then ${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ is the solution of the TRS defined in Eq. 8. Otherwise, the solution of the TRS must lie on the boundary of $\Vert \mathit{s}\Vert {}_2={\mathrm{\Delta }}^{\left(\mathit{k}\right)}$, and we should solve the Lagrange multiplier ${\mathit{\lambda }}^{\mathit{*}\left(\mathit{k}\right)}$ together with the search step ${\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ from the following nonlinear TRS equations,

$\left[{\mathit{H}}^{\left(\mathit{k}\right)}+\mathit{\lambda }{\mathit{I}}_{\mathit{n}}\right]\mathit{s}\left(\mathit{\lambda }\right)=-{\mathit{g}}^{\left(\mathit{k}\right)}(9)$

$\mathit{\pi }\left(\mathit{\lambda }\right)={\left[\mathit{s}\left(\mathit{\lambda }\right)\right]}^{\mathit{T}}\mathit{s}\left(\mathit{\lambda }\right)={\left[{\mathrm{\Delta }}^{\left(\mathit{k}\right)}\right]}^2(10)$

For the trivial case with ${\mathit{H}}^{\left(\mathit{k}\right)}=\mathit{\alpha }{\mathit{I}}_{\mathit{n}}$, we have $\mathit{s}\left(\mathit{\lambda }\right)=-\frac{{\mathit{g}}^{\left(\mathit{k}\right)}}{\mathit{\alpha }+\mathit{\lambda }}$ and $\mathit{\pi }\left(\mathit{\lambda }\right)=\frac{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}{{\left(\mathit{\alpha }+\mathit{\lambda }\right)}^2}$. The solution of the TRS defined in Eq. 8 is ${\mathit{\lambda }}^{\text{*}}=0$ and ${\mathit{s}}^{\text{*}\left(\mathit{k}\right)}=-{\mathit{g}}^{\left(\mathit{k}\right)}/\mathit{\alpha }$ if ${\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2\le \mathit{\alpha }{\mathrm{\Delta }}^{\left(\mathit{k}\right)}$ or ${\mathit{\lambda }}^{\text{*}}=\frac{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}{{\mathrm{\Delta }}^{\left(\mathit{k}\right)}}-\mathit{\alpha }$ and ${\mathit{s}}^{\text{*}\left(\mathit{k}\right)}=-{\mathit{g}}^{\left(\mathit{k}\right)}\frac{{\mathrm{\Delta }}^{\left(\mathit{k}\right)}}{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}$ otherwise.

We may apply the popular DNR method to solve the TRS when ${\mathit{H}}^{\left(\mathit{k}\right)}$ 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], $\varphi \left(\mathit{\lambda }\right)=\frac{1}{\mathit{s}\left(\mathit{\lambda }\right)}-\frac{1}{\mathit{\Delta }}=0$, iteratively, which requires computing the first order derivative ${\varphi }^{\mathit{\prime }}\left(\mathit{\lambda }\right)=-\frac{\mathit{v}{\left(\mathit{\lambda }\right)}^2}{\mathit{s}{\left(\mathit{\lambda }\right)}^3}$ where v is solved from ${\mathit{L}}^{\mathit{T}}\mathit{v}=\mathit{s}$ together with the Cholesky factorization of $\mathit{D}={\mathit{H}}^{\left(\mathit{k}\right)}+\mathit{\lambda }{\mathit{I}}_{\mathit{n}}=\mathit{L}{\mathit{L}}^{\mathit{T}}$ or LU-decomposition.

Given the trust region size $\mathbf{\Delta }={\mathrm{\Delta }}^{\left(\mathit{k}\right)}>0$, threshold of convergence ${\mathit{\delta }}_{\mathit{cr}}>0$, maximum number of iterations allowed ${\mathit{N}}_{\mathit{TRS},\mathit{max}}>0$, Hessian $\mathit{H}={\mathit{H}}^{\left(\mathit{k}\right)}$ and gradient $\mathit{g}={\mathit{g}}^{\left(\mathit{k}\right)}$ evaluated at the current best solution ${\mathit{x}}^{\left(\mathit{k}\right)}$, both the Lagrange multiplier ${\mathit{\lambda }}^{\mathit{*}}$and the trust region search step ${\mathit{s}}^{\mathit{*}}=\mathit{s}\left({\mathit{\lambda }}^{\mathit{*}}\right)$ can be solved from the TRS defined in Eq. 8, using the DNR method as summarized in Algorithm-2.

Algorithm-2: Solving the L-BFGS TRS Using the DNR Method

1. Initialize$\mathit{l}=0$, ${\mathit{\lambda }}_0=0$;

2. Compute$\mathit{D}=\mathit{H}+{\mathit{\lambda }}_0{\mathit{I}}_{\mathit{n}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

3. Solve${\mathit{u}}^{\mathit{*}}$from$\mathit{Lu}=-\mathit{g}$,${\mathit{s}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}^{\mathit{*}}$, and${\mathit{v}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{s}}^{\mathit{*}}$;

4. Compute$\Vert {\mathit{s}}^{\mathit{*}}\Vert$and$\Vert {\mathit{v}}^{\mathit{*}}\Vert$;

5. If$\Vert {\mathit{s}}^{\mathit{*}}\Vert \le \mathit{\Delta }$, then accept${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_0$and go to step 8;

6. Set $\mathit{\delta }=\left|1-\frac{\Vert {\mathit{s}}^{\mathit{*}}\Vert }{\mathit{\Delta }}\right|$;

7. Repeat steps (a) through (f) below, until convergence ($\mathit{\delta }\le {\mathit{\delta }}_{\mathit{cr}}$or$\mathit{l}>{\mathit{N}}_{\mathit{TRS},\mathit{max}}$):

a. Update${\mathit{\lambda }}_{\mathit{l}+1}={\mathit{\lambda }}_{\mathit{l}}+{\left(\frac{\Vert {\mathit{s}}^{\mathit{*}}\Vert }{\Vert {\mathit{v}}^{\mathit{*}}\Vert }\right)}^2\frac{\Vert {\mathit{s}}^{\mathit{*}}\Vert -\mathit{\Delta }}{\mathit{\Delta }}$;

b. Compute$\mathit{D}=\mathit{H}+{\mathit{\lambda }}_{\mathit{l}+1}{\mathit{I}}_{\mathit{n}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

c. Solve${\mathit{u}}^{\mathit{*}}$from$\mathit{Lu}=-\mathit{g}$, ${\mathit{s}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}^{\mathit{*}}$, and${\mathit{v}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{s}}^{\mathit{*}}$;

d. Compute$\Vert {\mathit{s}}^{\mathit{*}}\Vert$and$\Vert {\mathit{v}}^{\mathit{*}}\Vert$;

e. Update$\mathit{\delta }=\left|1-\frac{\Vert {\mathit{s}}^{\mathit{*}}\Vert }{\mathit{\Delta }}\right|$;

f. Set$\mathit{l}=\mathit{l}+1$.

8. Stop.

Let ${\mathit{N}}_{\mathit{DNR}}$ 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, ${\mathit{c}}_2=\left(5\mathit{n}+3{\mathit{n}}^2+\frac{{\mathit{n}}^3}{3}\right)\left({\mathit{N}}_{\mathit{DNR}}+1\right)$, including the following three operations: (1) computing $\mathit{D}=\mathit{H}+{\mathit{\lambda }}_{\mathit{k}+1}{\mathit{I}}_{\mathit{n}}$ and Cholesky decomposition $\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$ in step 2 and step 7(b) ($\mathit{n}+\frac{{\mathit{n}}^3}{3}$flops); (2) solving ${\mathit{u}}^{\text{*}}$ from $\mathit{Lu}=-\mathit{g}$, ${\mathit{s}}^{\text{*}}$ from ${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}^{\text{*}}$, and ${\mathit{v}}^{\text{*}}$ from ${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{s}}^{\text{*}}$ in step 3 and step 7(c) ($3{\mathit{n}}^2$ flops); (3) computing $\Vert {\mathit{s}}^{\text{*}}\Vert$ and $\Vert {\mathit{v}}^{\text{*}}\Vert$ in step 4 and step 7(d) ($4\mathit{n}$ 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,

${\mathit{c}}_{\mathit{DNR}}={\mathit{c}}_1+{\mathit{c}}_2=2\mathit{m}{\mathit{n}}^2+\left(7+3{\mathit{m}}^2\right)\mathit{n}+\left(5\mathit{n}+3{\mathit{n}}^2+\frac{{\mathit{n}}^3}{3}\right)\left({\mathit{N}}_{\mathit{DNR}}+1\right)+\mathbf{O}\left({\mathit{m}}^3\right)(11)$

Our numerical results indicate that the DNR TRS solver may fail when tested on the well-known Rosenbrock function, especially for problems with large $\mathit{n}$. 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 $\left|{\varphi }^{\mathit{\prime }}\left(\mathit{\lambda }\right)\right|$ 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 ${\varphi }^{\mathit{\prime }}\left(\mathit{\lambda }\right)$, Gao et al. [32] proposed a method to directly solve the TRS using inverse quadratic model interpolation (called the DIQ method), i.e., approximating $\mathit{\pi }\left(\mathit{\lambda }\right)={\left[\mathit{s}\left(\mathit{\lambda }\right)\right]}^{\mathit{T}}\mathit{s}\left(\mathit{\lambda }\right)$ by the following inverse quadratic function,

${\mathit{q}}_{\mathit{INV}}\left(\mathit{\lambda }\right)=\frac{\mathit{b}}{{\left(\mathit{\lambda }+\mathit{a}\right)}^2}(12)$

Both coefficients of $\mathit{a}$ and $\mathit{b}$ in Eq. 12 can be determined by interpolating the values of $\mathit{\pi }\left(\mathit{\lambda }\right)$ evaluated at two different points ${\mathit{\pi }}_1=\mathit{\pi }\left({\mathit{\lambda }}_1\right)$ and ${\mathit{\pi }}_2=\mathit{\pi }\left({\mathit{\lambda }}_2\right)$ with $0\le {\mathit{\lambda }}_1<{\mathit{\lambda }}_2$.

In the first iteration, we set ${\mathit{\lambda }}_1={\mathit{\lambda }}_{\mathit{min}}=0$ and accept ${\mathit{\lambda }}_{\mathit{min}}=0$ as the solution if ${\mathit{\pi }}_1\le {\mathbf{\Delta }}^2$. Otherwise, ${\mathit{\pi }}_1>{\mathbf{\Delta }}^2$ holds and we set ${\mathit{\lambda }}_2={\mathit{\lambda }}_{\mathit{max}}=\mathit{max}\left[0,\mathit{ }\frac{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }{\mathbf{\Delta }}-\mathit{\alpha }\right]$ and ${\mathit{\pi }}_2<{\mathbf{\Delta }}^2$ holds if it is not accepted as the solution. In the following iteration, either ${\mathit{\lambda }}_1$ or ${\mathit{\lambda }}_2$ will be updated accordingly such that the two conditions ${\mathit{\pi }}_1>{\mathbf{\Delta }}^2$ and ${\mathit{\pi }}_2<{\mathbf{\Delta }}^2$ always hold.

Let $\mathit{\rho }=\frac{\sqrt{{\mathit{\pi }}_1}}{\sqrt{{\mathit{\pi }}_1}-\sqrt{{\mathit{\pi }}_2}}\ge 1$ and the following solution of ${\mathit{q}}_{\mathit{INV}}\left(\mathit{\lambda }\right)={\mathbf{\Delta }}^2$ will be used as the trial search point for the next iteration,

${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_2-\mathit{\rho }\left(1-\frac{\sqrt{{\mathit{\pi }}_2}}{\mathbf{\Delta }}\right)\left({\mathit{\lambda }}_2-{\mathit{\lambda }}_1\right)(13)$

We either update ${\mathit{\lambda }}_1={\mathit{\lambda }}^{\mathit{*}}$ if $\mathit{\pi }\left({\mathit{\lambda }}^{\mathit{*}}\right)>{\mathbf{\Delta }}^2$ or update ${\mathit{\lambda }}_2={\mathit{\lambda }}^{\mathit{*}}$ if $\mathit{\pi }\left({\mathit{\lambda }}^{\mathit{*}}\right)<{\mathbf{\Delta }}^2$ iteratively until convergence. We accept ${\mathit{\lambda }}^{\mathit{*}}$ as the desired solution and terminate the iterative process either when $\left|\frac{\sqrt{\mathit{\pi }\left({\mathit{\lambda }}^{\mathit{*}}\right)}}{\mathbf{\Delta }}-1\right|<{\mathit{\delta }}_{\mathit{cr}}$, the tolerance for convergence, or when the number of iterations used to solve the TRS (${\mathit{N}}_{\mathit{DIQ}}$) reaches the user specified maximum iteration number (${\mathit{N}}_{\mathit{TRS},\mathit{max}}$), i.e.,${\mathit{N}}_{\mathit{DIQ}}\ge {\mathit{N}}_{\mathit{TRS},\mathit{max}}$.

Given the scaling factor $\mathit{\alpha }>0$, trust region size $\mathbf{\Delta }>0$, threshold of convergence ${\mathit{\delta }}_{\mathit{cr}}>0$, maximum number of iterations allowed ${\mathit{N}}_{\mathit{TRS},\mathit{max}}>0$, Hessian $\mathit{H}$ and gradient $\mathit{g}$, both the Lagrange multiplier ${\mathit{\lambda }}^{\mathit{*}}$and the trust region search step ${\mathit{s}}^{\mathit{*}}=\mathit{s}\left({\mathit{\lambda }}^{\mathit{*}}\right)$ 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

(1) Compute$\Vert \mathit{g}\Vert$;

(2) Calculate${\mathit{\lambda }}_{\mathit{max}}=\mathit{max}\left\{0,\frac{\mathit{g}}{\mathit{\Delta }}-\mathit{\alpha }\right\}$;

(3) Initialize$\mathit{k}=0$, ${\mathit{\lambda }}_1={\mathit{\lambda }}_{\mathit{min}}=0$and${\mathit{\lambda }}_2={\mathit{\lambda }}_{\mathit{max}}$;

(4) Compute$\mathit{D}=\mathit{H}+{\mathit{\lambda }}_1{\mathit{I}}_{\mathit{n}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

(5) Solve${\mathit{u}}_1$from$\mathit{Lu}=-\mathit{g}$and${\mathit{s}}_1$from${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}_1$;

(6) Compute${\mathit{\pi }}_1=\mathit{\pi }\left({\mathit{\lambda }}_1\right)={\mathit{s}}_1^{\mathit{T}}{\mathit{s}}_1$;

(7) If${\mathit{\pi }}_1\le {\mathit{\Delta }}^2$, then accept${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_1$and go to step 12;

(8) If${\mathit{\pi }}_1>{\mathit{\Delta }}^2$, then

a. Compute$\mathit{D}=\mathit{H}+{\mathit{\lambda }}_2{\mathit{I}}_{\mathit{n}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

b. Solve${\mathit{u}}_2$from$\mathit{Lu}=-\mathit{g}$and${\mathit{s}}_2$from${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}_2$;

c. Compute${\mathit{\pi }}_2=\mathit{\pi }\left({\mathit{\lambda }}_2\right)={\mathit{s}}_2^{\mathit{T}}{\mathit{s}}_2$;

(9) If$\left|1-\frac{\sqrt{{\mathit{\pi }}_1}}{\mathit{\Delta }}\right|<\left|1-\frac{\sqrt{{\mathit{\pi }}_2}}{\mathit{\Delta }}\right|$, then$\delta =\left|1-\frac{\sqrt{{\mathit{\pi }}_1}}{\mathit{\Delta }}\right|$, ${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_1$;

(10) Otherwise, $\mathit{\delta }=\left|1-\frac{\sqrt{{\mathit{\pi }}_2}}{\mathit{\Delta }}\right|$, ${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_2$;

(11) Repeat steps (a) through (i) below, until convergence ($\mathit{\delta }\le {\mathit{\delta }}_{\mathit{cr}}$ or $\mathit{k}>{\mathit{N}}_{\mathit{TRS},\mathit{max}}$):

a. Calculate$\mathit{\rho }=\frac{\sqrt{{\mathit{\pi }}_1}}{\sqrt{{\mathit{\pi }}_1}-\sqrt{{\mathit{\pi }}_2}}$;

b. Calculate${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_2-\mathit{\rho }\left(1-\frac{\sqrt{{\mathit{\pi }}_2}}{\mathit{\Delta }}\right)({\mathit{\lambda }}_2-{\mathit{\lambda }}_1)$;

c. Compute$\mathit{D}=\mathit{H}+{\mathit{\lambda }}^{\mathit{*}}{\mathit{I}}_{\mathit{n}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

d. Solve${\mathit{u}}^{\mathit{*}}$from$\mathit{Lu}=-\mathit{g}$and${\mathit{s}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{s}={\mathit{u}}^{\mathit{*}}$;

e. Compute ${\mathit{\pi }}^{\mathit{*}}=\mathit{\pi }\left({\mathit{\lambda }}^{\mathit{*}}\right)={\left[{\mathit{s}}^{\mathit{*}}\right]}^{\mathit{T}}{\mathit{s}}^{\mathit{*}}$;

f. Update$\mathit{\delta }=\left|1-\frac{\sqrt{{\mathit{\pi }}^{\mathit{*}}}}{\mathit{\Delta }}\right|$ ;

g. If$\sqrt{\mathit{ }{\mathit{\pi }}^{\mathit{*}}}>\mathit{\Delta }$, update${\mathit{\lambda }}_1={\mathit{\lambda }}^{\mathit{*}}$and${\mathit{\pi }}_1={\mathit{\pi }}^{\mathit{*}}$;

h. Otherwise, update${\mathit{\lambda }}_2={\mathit{\lambda }}^{\mathit{*}}$and${\mathit{\pi }}_2={\mathit{\pi }}^{\mathit{*}}$;

i. $\mathit{k}\leftarrow \mathit{k}+1$

(12) Stop.

Let ${\mathit{N}}_{\mathit{DIQ}}$ 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,

${\mathit{c}}_{\mathit{IQ}}=2\mathit{m}{\mathit{n}}^2+\left(7+3{\mathit{m}}^2\right)\mathit{n}+\left(5\mathit{n}+2{\mathit{n}}^2+\frac{{\mathit{n}}^3}{3}\right)\left({\mathit{N}}_{\mathit{DIQ}}+1\right)+\mathit{O}\left({\mathit{m}}^3\right)(14)$

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,

${\mathit{H}}^{(\mathit{k}+1)}=\mathit{\alpha }{\mathit{I}}_{\mathit{n}}+{\mathit{V}}^{\left(\mathit{k}\right)}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}(15)$

In Eq. 15, ${\mathit{V}}^{\left(\mathit{k}\right)}={\left[{\mathit{J}}^{(\mathit{k}+1)}\right]}^{\mathit{T}}$ is an $\mathit{n}\times \mathit{m}$ matrix, the transpose of the sensitivity matrix ${\mathit{J}}^{(\mathit{k}+1)}$ that is evaluated at the current best solution ${\mathit{x}}^{(\mathit{k}+1)}$. Here, $\mathit{m}$ denotes the number of observed data to be matched and it is not the same $\mathit{m}$ 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 ${\left[{\mathit{H}}^{(\mathit{k}+1)}+\mathit{\lambda }{\mathit{I}}_{\mathit{n}}\right]}^{-1}$ by applying the matrix inversion lemma and then solve the L-BFGS trust region search step $\mathit{s}\left(\mathit{\lambda }\right)$.

${\left[{\mathit{H}}^{(\mathit{k}+1)}+\mathit{\lambda }{\mathit{I}}_{\mathit{n}}\right]}^{-1}=\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\mathit{I}}_{\mathit{n}}-\frac{1}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^2}{\mathit{V}}^{\left(\mathit{k}\right)}{\left\{\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{V}}^{\left(\mathit{k}\right)}-{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\right\}}^{-1}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}(16)$

Let ${\mathit{P}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{V}}^{\left(\mathit{k}\right)}$ and ${\mathit{u}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$. We first solve $\mathit{v}\left(\mathit{\lambda }\right)$ from,

$\left\{\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\mathit{P}}^{\left(\mathit{k}\right)}-{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\right\}\mathit{v}\left(\mathit{\lambda }\right)={\mathit{u}}^{\left(\mathit{k}\right)}.(17)$

Then, we compute the trust region search step $\mathit{s}\left(\mathit{\lambda }\right)$ and $\mathit{\pi }\left(\mathit{\lambda }\right)=\Vert \mathit{s}\left(\mathit{\lambda }\right)\Vert { }^2$ as,

$\mathit{s}\left(\mathit{\lambda }\right)=-{\left[{\mathit{H}}^{(\mathit{k}+1)}+\mathit{\lambda }{\mathit{I}}_{\mathit{n}}\right]}^{-1}{\mathit{g}}^{\left(\mathit{k}\right)}=-\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\mathit{g}}^{\left(\mathit{k}\right)}+\frac{1}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^2}{\mathit{V}}^{\left(\mathit{k}\right)}\mathit{v}\left(\mathit{\lambda }\right),(18)$

$\mathit{\pi }\left(\mathit{\lambda }\right)=\frac{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^2}+\frac{{\left[\mathit{v}\left(\mathit{\lambda }\right)\right]}^{\mathit{T}}{\mathit{P}}^{\left(\mathit{k}\right)}\mathit{v}\left(\mathit{\lambda }\right)}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^4}-\frac{2{\left[{\mathit{u}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{v}\left(\mathit{\lambda }\right)}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^3}.(19)$

From Eq. 15, we have ${\mathit{P}}^{\left(\mathit{k}\right)}\mathit{v}\left(\mathit{\lambda }\right)=\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]\left\{{\mathit{u}}^{\left(\mathit{k}\right)}+{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\mathit{v}\left(\mathit{\lambda }\right)\right\}$ and Eq. 19 can be rewritten as,

$\mathit{\pi }\left(\mathit{\lambda }\right)=\frac{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^2}+\frac{{\left[\mathit{v}\left(\mathit{\lambda }\right)\right]}^{\mathit{T}}{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\mathit{v}\left(\mathit{\lambda }\right)-{\left[{\mathit{u}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{v}\left(\mathit{\lambda }\right)}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^3}.(20)$

We first try $\mathit{\lambda }={\mathit{\lambda }}_{\mathit{min}}=0$. If $\mathit{\pi }\left(0\right)\le {\mathit{\Delta }}^2$ holds, we accept ${\mathit{\lambda }}^{\mathit{*}}=0$ and ${\mathit{s}}^{\mathit{*}}=\mathit{s}\left(0\right)$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 $\varphi \left(\mathit{\lambda }\right)=\frac{1}{\sqrt{\mathit{\pi }\left(\mathit{\lambda }\right)}}-\frac{1}{\mathit{\Delta }}$ and we have ${\varphi }^{\mathit{\prime }}\left(\mathit{\lambda }\right)=-\frac{1}{2}\frac{{\mathit{\pi }}^{\mathit{\prime }}\left(\mathit{\lambda }\right)}{{\left[\mathit{\pi }\left(\mathit{\lambda }\right)\right]}^{3/2}}$ where ${\mathit{\pi }}^{\mathit{\prime }}\left(\mathit{\lambda }\right)$ is computed by,

${\mathit{\pi }}^{\mathit{\prime }}\left(\mathit{\lambda }\right)=\frac{{\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^3}-\frac{3\mathit{\pi }\left(\mathit{\lambda }\right)}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}+\frac{2{\left[\mathit{v}\left(\mathit{\lambda }\right)\right]}^{\mathit{T}}{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\mathit{w}\left(\mathit{\lambda }\right)-{\left[{\mathit{u}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{w}\left(\mathit{\lambda }\right)}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^3}.(21)$

In Eq. 21, $\mathit{w}\left(\mathit{\lambda }\right)={\mathit{v}}^{\mathit{\prime }}\left(\mathit{\lambda }\right)$ is solved from,

$\left\{\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\mathit{P}}^{\left(\mathit{k}\right)}-{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\right\}\mathit{w}\left(\mathit{\lambda }\right)=\frac{1}{{\left[{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }\right]}^2}{\mathit{P}}^{\left(\mathit{k}\right)}\mathit{v}\left(\mathit{\lambda }\right).(22)$

Because $\left\{\frac{1}{{\mathit{\alpha }}^{\left(\mathit{k}\right)}+\mathit{\lambda }}{\mathit{P}}^{\left(\mathit{k}\right)}-{\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}\right\}$ is an $\mathit{m}\times \mathit{m}$ symmetric and positive definite matrix, it only requires. $\frac{{\mathit{m}}^3}{3}+6{\mathit{m}}^2$ flops to solve $\mathit{v}\left(\mathit{\lambda }\right)$ and $\mathit{w}(\mathit{\lambda }$) from Eq. 17 and Eq. 22 using the Cholesky factorization.

Given the $\mathit{n}\times {\mathit{l}}_{\mathit{k}}$ variable shift matrix ${\mathit{S}}^{\left(\mathit{k}\right)}$ and gradient change matrix ${\mathit{Z}}^{\left(\mathit{k}\right)}$, we may directly compute the following three ${\mathit{l}}_{\mathit{k}}\times {\mathit{l}}_{\mathit{k}}$ matrices, ${\mathit{A}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{S}}^{\left(\mathit{k}\right)}$, ${\mathit{B}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{\left(\mathit{k}\right)}$ and ${\mathit{C}}^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{\left(\mathit{k}\right)}$, and then form the $\mathit{m}\times \mathit{m}$ matrix ${\mathit{P}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{V}}^{\left(\mathit{k}\right)}$ as,

${\mathit{P}}^{\left(\mathit{k}\right)}=\left\{\begin{array}{c}{\mathit{A}}^{\left(\mathit{k}\right)} {\mathit{B}}^{\left(\mathit{k}\right)}\\ {\left[{\mathit{B}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} {\mathit{C}}^{\left(\mathit{k}\right)}\end{array}\right\}(23)$

Directly computing the three ${\mathit{l}}_{\mathit{k}}\times {\mathit{l}}_{\mathit{k}}$ matrices, ${\mathit{A}}^{\left(\mathit{k}\right)}$, ${\mathit{B}}^{\left(\mathit{k}\right)}$ and ${\mathit{C}}^{\left(\mathit{k}\right)}$, requires $6{\mathit{l}}_{\mathit{k}}^2\mathit{n}=1.5{\mathit{m}}^2\mathit{n}$ flops. We also need to compute the vector ${\mathit{u}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$ ($2\mathit{mn}$ flops) and ${\Vert {\mathit{g}}^{\left(\mathit{k}\right)}\Vert }^2$ ($2\mathit{n}$ flops). To further reduce the computational cost, the three matrices ${\mathit{A}}^{\left(\mathit{k}\right)}$, ${\mathit{B}}^{\left(\mathit{k}\right)}$ and ${\mathit{C}}^{\left(\mathit{k}\right)}$ and the vector ${\mathit{u}}^{\left(\mathit{k}\right)}$ should be updated iteratively instead.

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

We first initialize $\mathit{k}=0$ and ${\mathit{H}}^{\left(0\right)}={\mathit{I}}_{\mathit{n}}$. If the trial search point ${\mathit{x}}_{\mathit{T}}^{\left(0\right)}={\mathit{x}}^{\left(0\right)}+{\mathit{s}}^{\mathit{*}\left(0\right)}$ does not improve the objective function, i.e., $\mathit{f}\left({\mathit{x}}_{\mathit{T}}^{\left(0\right)}\right(\ge \mathit{f}\left({\mathit{x}}^{\left(0\right)}\right)$, we set ${\mathit{x}}^{(\mathit{k}+1)}={\mathit{x}}^{\left(0\right)}$ and ${\mathit{H}}^{(\mathit{k}+1)}={\mathit{H}}^{\left(0\right)}$, recompute the sensitivity matrix ${\mathit{J}}^{(\mathit{k}+1)}=\mathit{J}\left({\mathit{x}}^{\left(0\right)}\right)$ and the gradient ${\mathit{g}}^{(\mathit{k}+1)}=\mathit{g}\left({\mathit{x}}^{\left(0\right)}\right)$ when needed (e.g., using updated training data points), and shrink the trust region size ${\mathrm{\Delta }}^{(\mathit{k}+1)}={\mathit{\gamma }}_{\mathit{d}}{\mathrm{\Delta }}^{\left(\mathit{k}\right)}$ with $0<{\mathit{\gamma }}_{\mathit{d}}<1$. We repeatedly generate trial search point ${\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}={\mathit{x}}^{\left(\mathit{k}\right)}+{\mathit{s}}^{\mathit{*}\left(\mathit{k}\right)}$ until ${\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}$ improves the objective function, i.e., $\mathit{f}\left({\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}\right)<\mathit{f}\left({\mathit{x}}^{\left(0\right)}\right)$.

If $\mathit{f}\left({\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}\right)<\mathit{f}\left({\mathit{x}}^{\left(0\right)}\right)$, we set ${\mathit{l}}_{\mathit{k}}=1$, ${\mathit{x}}^{(\mathit{k}+1)}={\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}$, evaluate ${\mathit{g}}^{(\mathit{k}+1)}=\mathit{g}\left({\mathit{x}}^{(\mathit{k}+1)}\right)$, and compute $\mathit{s}={\mathit{x}}^{(\mathit{k}+1)}-{\mathit{x}}^{\left(\mathit{k}\right)}$ and $\mathit{z}={\mathit{g}}^{(\mathit{k}+1)}-{\mathit{g}}^{\left(\mathit{k}\right)}$. In the following iterations, we update ${\mathit{l}}_{\mathit{k}}\ge 1$ and associated matrices and vectors used for solving the L-BFGS TRS for different cases accordingly. We use $\mathit{Sfla}{\mathit{g}}^{\left(\mathit{k}\right)}$ to indicate whether the new search point ${\mathit{x}}_{\mathit{T}}^{\left(\mathit{k}\right)}$ improves the objective function ($\mathit{Sfla}{\mathit{g}}^{\left(\mathit{k}\right)}="\mathit{True}"$) or not ($\mathit{Sfla}{\mathit{g}}^{\left(\mathit{k}\right)}="\mathit{False}"$).

The L-BFGS Hessian updating formulation of Eq. 4 requires ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{s}>{\mathit{c}}_3\Vert \mathit{z}\Vert \Vert \mathit{s}\Vert$ and ${\mathit{\alpha }}^{\left(\mathit{k}\right)}=\frac{{\mathit{z}}^{\mathit{T}}\mathit{z}}{{\mathit{z}}^{\mathit{T}}\mathit{s}}>{\mathit{\alpha }}_{\mathit{cr}}>0$ where ${\mathit{c}}_3$ and ${\mathit{\alpha }}_{\mathit{cr}}$ are small positive numbers. If ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}\le {\mathit{c}}_3\Vert \mathit{z}\Vert \Vert \mathit{s}\Vert$ or ${\mathit{\alpha }}^{\left(\mathit{k}\right)}\le {\mathit{\alpha }}_{\mathit{cr}}$, we simply set ${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}$, ${\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}}\right)}={\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$, ${\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}={\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$, ${\mathit{\alpha }}^{\left(\mathit{k}\right)}={\mathit{\alpha }}^{(\mathit{k}-1)}$, ${\mathit{S}}^{\left(\mathit{k}\right)}={\mathit{S}}^{\left(\mathit{k}-1\right)}$, ${\mathit{Z}}^{\left(\mathit{k}\right)}={\mathit{Z}}^{\left(\mathit{k}-1\right)}$, ${\mathit{A}}^{\left(\mathit{k}\right)}={\mathit{A}}^{(\mathit{k}-1)}$, ${\mathit{B}}^{\left(\mathit{k}\right)}={\mathit{B}}^{(\mathit{k}-1)}$, ${\mathit{C}}^{\left(\mathit{k}\right)}={\mathit{C}}^{(\mathit{k}-1)}$, and ${\mathit{P}}^{\left(\mathit{k}\right)}={\mathit{P}}^{(\mathit{k}-1)}$, with no updating. However, we need to compute ${\mathit{u}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)},\mathit{ }{\mathit{Z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$ using the new gradient ${\mathit{g}}^{\left(\mathit{k}+1\right)}$. In the following cases, we assume that both conditions ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}>{\mathit{c}}_3\Vert \mathit{z}\Vert \Vert \mathit{s}\Vert$ and ${\mathit{\alpha }}^{\left(\mathit{k}\right)}>{\mathit{\alpha }}_{\mathit{cr}}$ are satisfied.

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

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

Let ${\tilde{\mathit{Z}}}^{(\mathit{k}-1)}$ be the submatrix of ${\mathit{Z}}^{(\mathit{k}-1)}$ by removing its last column ${\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$, ${\tilde{\mathit{B}}}^{\left(\mathit{k}-1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}$ the submatrix of ${\mathit{B}}^{(\mathit{k}-1)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{(\mathit{k}-1)}$ by removing its last column, and ${\tilde{\mathit{C}}}^{\left(\mathit{k}-1\right)}={\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}$ the submatrix of ${\mathit{C}}^{(\mathit{k}-1)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{Z}}^{(\mathit{k}-1)}$ by removing its last row and last column. We first compute the following two vectors ${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$ and ${\tilde{\mathit{p}}}_4^{\left(\mathit{k}\right)}={\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{z}$ and two scalars ${\mathit{\mu }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{z}$ and ${\mathit{\tau }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}-{\left[{\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$, with computational cost $2\mathit{mn}$ flops.

We set ${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}$, ${\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}}\right)}={\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$, ${\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}=\mathit{z}$, ${\mathit{S}}^{\left(\mathit{k}\right)}={\mathit{S}}^{(\mathit{k}-1)}$, and ${\mathit{A}}^{\left(\mathit{k}\right)}={\mathit{A}}^{(\mathit{k}-1)}$, and Update ${\mathit{Z}}^{\left(\mathit{k}\right)}$ as,

${\mathit{Z}}^{\left(\mathit{k}\right)}=\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\mathit{ z}\right].(24)$

It is straightforward to derive the following equations to update ${\mathit{B}}^{\left(\mathit{k}\right)}$ and ${\mathit{C}}^{\left(\mathit{k}\right)}$ accordingly,

${\mathit{B}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\mathit{ z }\right]=\left[{\tilde{\mathit{B}}}^{\left(\mathit{k}-1\right)}\mathit{ }{\mathit{p}}_2^{\left(\mathit{k}\right)}\right],(25)$

${\mathit{C}}^{\left(\mathit{k}\right)}=\left\{\begin{array}{c}{\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\\ {\mathit{z}}^{\mathit{T}}\end{array}\right\}\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\mathit{ z }\right]=\left\{\begin{array}{c}{\tilde{\mathit{C}}}^{\left(\mathit{k}-1\right)} {\tilde{\mathit{p}}}_4^{\left(\mathit{k}\right)}\\ {\left[{\tilde{\mathit{p}}}_4^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} {\mathit{\mu }}^{\left(\mathit{k}\right)} \end{array}\right\}.(26)$

By definition,

${\mathit{u}}^{\left(\mathit{k}-1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{ }{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{ }{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)},{\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$

${\mathit{u}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{\left(\mathit{k}\right)},\mathit{ }{\mathit{Z}}^{\left(\mathit{k}\right)}\right]{\mathit{g}}^{\left(\mathit{k}+1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{ }{\tilde{\mathit{Z}}}^{\left(\mathit{k}-1\right)},\mathit{z}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$

Thus, we have,

${\mathit{u}}^{\left(\mathit{k}\right)}={\mathit{u}}^{\left(\mathit{k}-1\right)}+\left\{\begin{array}{c}{\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{ }{\tilde{\mathit{Z}}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}z\\ {\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}-{\left[{\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}\end{array}\right\}={\mathit{u}}^{\left(\mathit{k}-1\right)}+\left[\begin{array}{c}{\mathit{p}}_2^{\left(\mathit{k}\right)}\\ {\tilde{\mathit{p}}}_4^{\left(\mathit{k}\right)}\\ {\mathit{\tau }}^{\left(\mathit{k}\right)}\end{array}\right].(27)$

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

We set ${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}+1$ and ${\mathit{s}}^{\left({\mathit{l}}_{\mathit{k}}\right)}=\mathit{s}$ and ${\mathit{z}}^{\left({\mathit{l}}_{\mathit{k}}\right)}=\mathit{z}$, compute four ${\mathit{l}}_{\mathit{k}-1}$-dimensional vectors ${\mathit{p}}_1^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$ and ${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$ and ${\left[{\mathit{p}}_3^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}={\mathit{s}}^{\mathit{T}}{\mathit{Z}}^{(\mathit{k}-1)}$ and ${\mathit{p}}_4^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$ and five scalars, ${\mathit{\beta }}^{\left(\mathit{k}\right)}={\mathit{s}}^{\mathit{T}}\mathit{s}$ and ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{s}}^{\mathit{T}}\mathit{z}$ and ${\mathit{\mu }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{z}$ and ${\mathit{\gamma }}_{\mathit{s}}^{\left(\mathit{k}\right)}={\mathit{s}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$ and ${\mathit{\gamma }}_{\mathit{z}}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$, with computational cost $4\mathit{mn}$ flops. Both ${\mathit{S}}^{\left(\mathit{k}\right)}$ and ${\mathit{Z}}^{\left(\mathit{k}\right)}$ are updated by,

${\mathit{S}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\mathit{ s}\right](28)$

${\mathit{Z}}^{\left(\mathit{k}\right)}=\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\mathit{ z}\right](29)$

It is straightforward to derive the following equations to update A^(k), B^(k), and C^(k),

${\mathit{A}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\mathit{ s}\right]}^{\mathit{T}}\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\mathit{ s }\right]=\left\{\begin{array}{c}{\mathit{A}}^{\left(\mathit{k}-1\right)} {\mathit{p}}_1^{\left(\mathit{k}\right)}\\ {\left[\mathit{ }{\mathit{p}}_1^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} {\mathit{\beta }}^{\left(\mathit{k}\right)}\end{array}\right\}(30)$

${\mathit{B}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\mathit{ s}\right]}^{\mathit{T}}\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\mathit{ z }\right]=\left\{\begin{array}{c}{\mathit{B}}^{\left(\mathit{k}-1\right)} {\mathit{p}}_2^{\left(\mathit{k}\right)}\\ {\left[{\mathit{p}}_3^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} {\mathit{\epsilon }}^{\left(\mathit{k}\right)}\end{array}\right\}(31)$

${\mathit{C}}^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\mathit{ Z}\right]}^{\mathit{T}}\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\mathit{ z }\right]=\left\{\begin{array}{c}{\mathit{C}}^{\left(\mathit{k}-1\right)} {\mathit{p}}_4^{\left(\mathit{k}\right)}\\ {\left[{\mathit{p}}_4^{\left(\mathit{k}\right)}\right]}^{\mathit{T}} {\mathit{\mu }}^{\left(\mathit{k}\right)}\end{array}\right\}(32)$

We first split the $2{\mathit{l}}_{\mathit{k}-1}$-dimensional vector ${\mathit{u}}^{\left(\mathit{k}-1\right)}$ into two ${\mathit{l}}_{\mathit{k}-1}$-dimensional sub-vectors ${\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}$ and ${\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}$, i.e.,

${\mathit{u}}^{\left(\mathit{k}-1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{ }{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}=\left[\begin{array}{c}{\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}\\ {\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}\end{array}\right].$

where ${\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$ is composed of the first $l_{k-1}$ rows in ${\mathit{u}}^{\left(\mathit{k}-1\right)}$ and ${\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$ is composed of the last ${\mathit{l}}_{\mathit{k}-1}$ rows in ${\mathit{u}}^{\left(\mathit{k}-1\right)}$.

Because ${\mathit{u}}^{\left(\mathit{k}\right)}=\left[{\mathit{S}}^{\left(\mathit{k}\right)},{\mathit{Z}}^{\left(\mathit{k}\right)}\right]{\mathit{g}}^{\left(\mathit{k}+1\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)},\mathit{s}{\mathit{Z}}^{\left(\mathit{k}-1\right)},\mathit{z}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$ and ${\mathit{g}}^{\left(\mathit{k}+1\right)}={\mathit{g}}^{\left(\mathit{k}\right)}+\mathit{z}$, thus we have,

${\mathit{u}}^{\left(\mathit{k}\right)}=\left[\begin{array}{c}{\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}+{\mathit{p}}_2^{\left(\mathit{k}\right)}\\ {\mathit{\gamma }}_{\mathit{s}}^{\left(\mathit{k}\right)}\\ {\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}+{\mathit{p}}_4^{\left(\mathit{k}\right)}\\ {\mathit{\gamma }}_{\mathit{z}}^{\left(\mathit{k}\right)}\end{array}\right].(33)$

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

We set ${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}={\mathit{L}}_{\mathit{M}}$. Let ${\widehat{\mathit{S}}}^{\left(\mathit{k}-1\right)}$ and ${\widehat{\mathit{Z}}}^{\left(\mathit{k}-1\right)}$ denote the submatrices of ${\mathit{S}}^{\left(\mathit{k}-1\right)}$ and ${\mathit{Z}}^{\left(\mathit{k}-1\right)}$ by removing the first column ${\mathit{s}}^{\left(1\right)}$ and ${\mathit{z}}^{\left(1\right)}$ from ${\mathit{S}}^{(\mathit{k}-1)}$ and ${\mathit{Z}}^{(\mathit{k}-1)}$, respectively. Let ${\widehat{\mathit{A}}}^{(\mathit{k}-1)}$ and ${\widehat{\mathit{B}}}^{(\mathit{k}-1)}$ and ${\widehat{\mathit{C}}}^{(\mathit{k}-1)}$ denote the submatrices of A^(k), B^(k), and C^(k) by removing their first row and first column.

$\text{We }$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 ${\widehat{\mathit{S}}}^{\left(\mathit{k}-1\right)}$ and ${\widehat{\mathit{Z}}}^{\left(\mathit{k}-1\right)}$ and ${\widehat{\mathit{A}}}^{(\mathit{k}-1)}$ and ${\widehat{\mathit{B}}}^{(\mathit{k}-1)}$ and ${\widehat{\mathit{C}}}^{(\mathit{k}-1)}$, respectively. We can also apply Eq. 33 to update ${\mathit{u}}^{\left(\mathit{k}\right)}$ by replacing ${\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}$ with ${\widehat{\mathit{u}}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}={\left[{\widehat{\mathit{S}}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$ that is composed of the second row through the ${\mathit{l}}_{\mathit{k}-1}$-th rows in ${\mathit{u}}^{\left(\mathit{k}-1\right)}$ and replacing ${\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}$ with ${\widehat{\mathit{u}}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}={\left[{\widehat{\mathit{Z}}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$ that is composed of the last ${\mathit{l}}_{\mathit{k}-1}-1$ rows in ${\mathit{u}}^{\left(\mathit{k}-1\right)}$.

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, $\mathit{n}$-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 ${\mathit{W}}_{\mathit{I}}^{(\mathit{k}-1)}={\left[{\mathit{W}}^{(\mathit{k}-1)}\right]}^{-1}$, we can update ${\mathit{l}}_{\mathit{k}}$ and the ${\mathit{l}}_{\mathit{k}}$-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 ${\mathit{P}}^{\left(\mathit{k}\right)}$ and ${\mathit{W}}_{\mathit{I}}^{\left(\mathit{k}\right)}={\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}$, using the algorithm as summarized in Algorithm-4.

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

1. Compute${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{s}$, ${\mathit{\mu }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{z}$, ${\mathit{\beta }}^{\left(\mathit{k}\right)}={\mathit{s}}^{\mathit{T}}\mathit{s}$;

2. Compute the scaling factor${\mathit{\alpha }}^{\left(\mathit{k}\right)}={\mathit{\mu }}^{\left(\mathit{k}\right)}/{\mathit{\epsilon }}^{\left(\mathit{k}\right)}$if${\mathit{\epsilon }}^{\left(\mathit{k}\right)}>{\mathit{c}}_3\sqrt{{\mathit{\beta }}^{\left(\mathit{k}\right)}{\mathit{\mu }}^{\left(\mathit{k}\right)}}$, or set${\mathit{\alpha }}^{\left(\mathit{k}\right)}=0$otherwise;

3. If${\mathit{\alpha }}^{\left(\mathit{k}\right)}\le {\mathit{\alpha }}_{\mathit{cr}}$or${\mathit{\epsilon }}^{\left(\mathit{k}\right)}\le {\mathit{c}}_3\sqrt{{\mathit{\beta }}^{\left(\mathit{k}\right)}{\mathit{\mu }}^{\left(\mathit{k}\right)}}$, then

a. Set${\mathit{\alpha }}^{\left(\mathit{k}\right)}={\mathit{\alpha }}^{(\mathit{k}-1)}$and${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}$;

b. Update${\mathit{S}}^{\left(\mathit{k}\right)}={\mathit{S}}^{(\mathit{k}-1)}$and${\mathit{Z}}^{\left(\mathit{k}\right)}={\mathit{Z}}^{(\mathit{k}-1)}$;

c. Update${\mathit{A}}^{\left(\mathit{k}\right)}={\mathit{A}}^{(\mathit{k}-1)}$, ${\mathit{B}}^{\left(\mathit{k}\right)}={\mathit{B}}^{(\mathit{k}-1)}$, and${\mathit{C}}^{\left(\mathit{k}\right)}={\mathit{C}}^{(\mathit{k}-1)}$;

d. Update${\mathit{P}}^{\left(\mathit{k}\right)}={\mathit{P}}^{(\mathit{k}-1)}$and$\mathit{W}{\mathit{I}}^{\left(\mathit{k}\right)}=\mathit{W}{\mathit{I}}^{(\mathit{k}-1)}$;

e. Compute${\mathit{u}}^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}\right)},\mathit{ }{\mathit{Z}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$;

f. Goto step 10.

4. Else if$\mathit{Sfla}{\mathit{g}}^{\left(\mathit{k}\right)}="\mathit{False}"$(for Case-1), then

a. Set${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}$;

b. Update${\mathit{S}}^{\left(\mathit{k}\right)}={\mathit{S}}^{(\mathit{k}-1)}$;

c. Set${\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$to the last column of${\mathit{Z}}^{(\mathit{k}-1)}$;

d. Set${\tilde{\mathit{Z}}}^{\left(\mathit{k}\right)}$by removing the last column${\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}$from${\mathit{Z}}^{(\mathit{k}-1)}$;

e. Update${\mathit{Z}}^{\left(\mathit{k}\right)}=\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}\right)}\mathit{ z}\right]$;

f. Compute${\mathit{\tau }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}-{\left[{\tilde{\mathit{z}}}^{\left({\mathit{l}}_{\mathit{k}-1}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$;

g. Compute${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$and${\tilde{\mathit{p}}}_4^{\left(\mathit{k}\right)}={\left[{\tilde{\mathit{Z}}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}\mathit{z}$;

h. Compute${\mathit{u}}^{\left(\mathit{k}\right)}$usingEq. 27;

i. Update${\mathit{A}}^{\left(\mathit{k}\right)}={\mathit{A}}^{(\mathit{k}-1)}$;

j. Set${\tilde{\mathit{B}}}^{(\mathit{k}-1)}$by removing the first column from${\mathit{B}}^{(\mathit{k}-1)}$;

k. Set${\tilde{\mathit{C}}}^{(\mathit{k}-1)}$by removing the first row and the first column from${\mathit{C}}^{(\mathit{k}-1)}$;

l. Update${\mathit{B}}^{\left(\mathit{k}\right)}$and${\mathit{C}}^{\left(\mathit{k}\right)}$usingEq. 25andEq. 26;

5. Else if${\mathit{l}}_{\mathit{k}-1}<{\mathit{L}}_{\mathit{M}}$(for Case-2), then

a. Set${\mathit{l}}_{\mathit{k}}={\mathit{l}}_{\mathit{k}-1}+1$;

b. Compute${\mathit{p}}_1^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$, ${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$, ${\mathit{p}}_3^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$, and${\mathit{p}}_4^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$;

c. Compute${\mathit{\gamma }}_{\mathit{s}}={\mathit{s}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$and${\mathit{\gamma }}_{\mathit{z}}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$;

d. Update${\mathit{S}}^{\left(\mathit{k}\right)}$,${\mathit{Z}}^{\left(\mathit{k}\right)}$,${\mathit{A}}^{\left(\mathit{k}\right)}$,${\mathit{B}}^{\left(\mathit{k}\right)}$, and${\mathit{C}}^{\left(\mathit{k}\right)}$usingEq. 28throughEq. 32, respectively;

e. Set${\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}$to the first${\mathit{l}}_{\mathit{k}-1}$rows in${\mathit{u}}^{\left(\mathit{k}-1\right)}$and${\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}$to the last${\mathit{l}}_{\mathit{k}-1}$rows in${\mathit{u}}^{\left(\mathit{k}-1\right)}$;

f. Update${\mathit{u}}^{\left(\mathit{k}\right)}$usingEq. 31.

6. Else (for Case-3)

a. Set${\mathit{l}}_{\mathit{k}}={\mathit{L}}_{\mathit{M}}$;

b. Remove the first column from${\mathit{S}}^{(\mathit{k}-1)}$and${\mathit{Z}}^{(\mathit{k}-1)}$;

c. Remove the first row and first column from${\mathit{A}}^{(\mathit{k}-1)}$,${\mathit{B}}^{(\mathit{k}-1)}$, and${\mathit{C}}^{(\mathit{k}-1)}$;

d. Compute${\mathit{p}}_1^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$,${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$,${\mathit{p}}_3^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$, and${\mathit{p}}_4^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$;

e. Compute${\mathit{\gamma }}_{\mathit{s}}={\mathit{s}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$and${\mathit{\gamma }}_{\mathit{z}}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$;

f. Update${\mathit{S}}^{\left(\mathit{k}\right)}$,${\mathit{Z}}^{\left(\mathit{k}\right)}$,${\mathit{A}}^{\left(\mathit{k}\right)}$,${\mathit{B}}^{\left(\mathit{k}\right)}$, and${\mathit{C}}^{\left(\mathit{k}\right)}$usingEq. 28 throughEq. 32present, respectively;

g. Set${\mathit{u}}_{\mathit{S}}^{\left(\mathit{k}-1\right)}$to the second row through the${\mathit{l}}_{\mathit{k}-1}$-th rows in${\mathit{u}}^{\left(\mathit{k}-1\right)}$;

h. Set${\mathit{u}}_{\mathit{Z}}^{\left(\mathit{k}-1\right)}$to the last${\mathit{l}}_{\mathit{k}-1}-1$rows in${\mathit{u}}^{\left(\mathit{k}-1\right)}$;

i. Update${\mathit{u}}^{\left(\mathit{k}\right)}$usingEq. 33.

7. Form the matrix${\mathit{P}}^{\left(\mathit{k}\right)}$usingEq. 23;

8. Decompose${\mathit{B}}^{\left(\mathit{k}\right)}={\mathit{L}}^{\left(\mathit{k}\right)}+{\mathit{E}}^{\left(\mathit{k}\right)}+{\mathit{U}}^{\left(\mathit{k}\right)}$;

9. Form the matrix${\mathit{W}}_{\mathit{I}}^{\left(\mathit{k}\right)}$usingEq. 5;

10. End

The computational cost to update matrices and vectors used for solving the L-BFGS TRS as described in Algorithm-4 is at most ${\mathit{c}}_4=2\mathit{n}+4\mathit{mn}$ (taking Case-2 in step 5 as an example), including the following 4 operations: (1) computing ${\mathit{\epsilon }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{s}$, ${\mathit{\mu }}^{\left(\mathit{k}\right)}={\mathit{z}}^{\mathit{T}}\mathit{z}$, and ${\mathit{\beta }}^{\left(\mathit{k}\right)}={\mathit{s}}^{\mathit{T}}\mathit{s}$ in step 1 ($6\mathit{n}$ flops); (2) computing ${\mathit{p}}_1^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$,${\mathit{p}}_2^{\left(\mathit{k}\right)}={\left[{\mathit{S}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$,${\mathit{p}}_3^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{s}$, and ${\mathit{p}}_4^{\left(\mathit{k}\right)}={\left[{\mathit{Z}}^{\left(\mathit{k}-1\right)}\right]}^{\mathit{T}}\mathit{z}$ in step 5(b) ($8{\mathit{l}}_{\mathit{k}-1}\mathit{n}=4\mathit{mn}-8\mathit{n}$ flops),; (3) computing ${\mathit{\gamma }}_{\mathit{s}}={\mathit{s}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$ and ${\mathit{\gamma }}_{\mathit{z}}={\mathit{z}}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}+1\right)}$in step 5(c) ($4\mathit{n}$ flops).

In this paper, we apply the same idea of reducing the dimension from $\mathit{n}$ to $\mathit{m}$ 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 $\mathit{n}\times \mathit{m}$ matrix ${\mathit{V}}^{\left(\mathit{k}\right)}$ (or the transpose of the sensitivity matrix) in Eq. 15 is directly computed. Its dimensions (both the number of variables $\mathit{n}$ and the number of observed data $\mathit{m}$) are fixed, i.e., they remain the same for all iterations. Therefore, we have to directly compute the $\mathit{m}\times \mathit{m}$ matrix ${\mathit{P}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{V}}^{\left(\mathit{k}\right)}$ and the $\mathit{m}$-dimensional vector ${\mathit{u}}^{\left(\mathit{k}\right)}={\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}{\mathit{g}}^{\left(\mathit{k}\right)}$, with the computational cost $2{\mathit{m}}^2\mathit{n}+2\mathit{mn}$ flops. Using the algorithm as summarized in Algorithm-4, the computational cost can be further reduced to $2\mathit{n}+4\mathit{mn}$, by a factor of ${\mathit{l}}_{\mathit{k}}=\mathit{m}/2$, 5 to 20 roughly.

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

Given the trust region size $\mathbf{\Delta }={\mathrm{\Delta }}^{\left(\mathit{k}\right)}>0$, threshold of convergence ${\mathit{\delta }}_{\mathit{cr}}>0$, maximum number of iterations allowed ${\mathit{N}}_{\mathit{TRS},\mathit{max}}>0$, the scaling factor $\mathit{\alpha }>{\mathit{\alpha }}_{\mathit{cr}}>0$, $\mathit{m}\ge 0$, $\mathit{n}\times \mathit{m}$ matrix $\mathit{V}=\left[{\mathit{S}}^{\left(\mathit{k}\right)},{\mathit{Z}}^{\left(\mathit{k}\right)}\right]$, $\mathit{m}\times \mathit{m}$ matrices $\mathit{P}={\mathit{P}}^{\left(\mathit{k}\right)}$ and ${\mathit{W}}_{\mathit{I}}={\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}$, $\mathit{m}$-dimensional vector $\mathit{u}={\mathit{u}}^{\left(\mathit{k}\right)}$ and $\mathit{n}$-dimensional gradient vector $\mathit{g}={\mathit{g}}^{\left(\mathit{k}\right)}$ evaluated at the current best solution ${\mathit{x}}^{\left(\mathit{k}\right)}$, both the Lagrange multiplier ${\mathit{\lambda }}^{\text{*}}$and the trust region search step ${\mathit{s}}^{\text{*}}=\mathit{s}\left({\mathit{\lambda }}^{\text{*}}\right)$ 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

1. Compute $\Vert \mathit{g}\Vert$;

2. Initialize$\mathit{l}=0$, ${\mathit{\lambda }}_0=0$;

3. If$\mathit{m}=0$, then

a. Set${\mathit{\lambda }}^{\mathit{*}}=0$and${\mathit{s}}^{\ast }=-\mathit{g}/\mathit{\alpha }$if$\Vert \mathit{g}\Vert \le \mathit{\alpha }\Delta$;

b. Set${\mathit{\lambda }}^{\mathit{*}}=\frac{\Vert \mathit{g}\Vert }{\Delta }-\mathit{\alpha }$and${\mathit{s}}^{\ast }=-\mathit{g}\frac{\Delta }{\Vert \mathit{g}\Vert }$if$\Vert \mathit{g}\Vert >\mathit{\alpha }\Delta$;

c. Goto step 11.

4. Compute$\mathit{D}=\frac{1}{\mathit{\alpha }+{\mathit{\lambda }}_0}\mathit{P}-{\mathit{W}}_{\mathit{I}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

5. Solve${\mathit{p}}^{\mathit{*}}$from$\mathit{Lp}=\mathit{u}$and${\mathit{v}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{p}}^{\mathit{*}}$;

6. Compute${\mathit{u}}^{\mathit{T}}{\mathit{v}}^{\mathit{*}}$, $\mathit{p}={\mathit{W}}_{\mathit{I}}{\mathit{v}}^{\mathit{*}}$, ${\mathit{p}}^{\mathit{T}}{\mathit{v}}^{\mathit{*}}$, and${\mathit{\pi }}_0=\mathit{\pi }\left({\mathit{\lambda }}_0\right)$usingEq. 20;

7. If${\mathit{\pi }}_0\le {\mathit{\Delta }}^2$, then accept${\mathit{\lambda }}^{\mathit{*}}={\mathit{\lambda }}_0$and goto step 10;

8. Set$\mathit{\delta }=\left|1-\frac{\sqrt{{\mathit{\pi }}_0}}{\mathit{\Delta }}\right|$;

9. Repeat steps (a) through (j) below, until convergence ($\mathit{\delta }\le {\mathit{\delta }}_{\mathit{cr}}$ or $\mathit{l}>{\mathit{N}}_{\mathit{TRS},\mathit{max}}$):

a. Compute$\mathit{q}=\frac{1}{{\left(\mathit{\alpha }+{\mathit{\lambda }}_{\mathit{l}}\right)}^2}\mathit{P}{\mathit{v}}^{\mathit{*}}$ ;

b. Solve${\mathit{p}}^{\mathit{*}}$from$\mathit{Lp}=\mathit{q}$and${\mathit{w}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{w}={\mathit{p}}^{\mathit{*}}$;

c. Compute${\mathit{u}}^{\mathit{T}}{\mathit{w}}^{\mathit{*}}$,$\mathit{p}={\mathit{W}}_{\mathit{I}}{\mathit{w}}^{\mathit{*}}$,${\mathit{p}}^{\mathit{T}}{\mathit{v}}^{\mathit{*}}$, and${\mathit{\pi }}_{\mathit{l}}^{\mathit{\prime }}={\mathit{\pi }}^{\mathit{\prime }}\left({\mathit{\lambda }}_{\mathit{l}}\right)$usingEq. 21;

d. Compute${\varphi }_{\mathit{l}}=\frac{1}{\sqrt{{\mathit{\pi }}_{\mathit{l}}}}-\frac{1}{\mathit{\Delta }}$and${\varphi }_{\mathit{l}}^{\mathit{\prime }}=-\frac{1}{2}\frac{{\mathit{\pi }}_{\mathit{l}}^{\mathit{\prime }}}{{\left({\mathit{\pi }}_{\mathit{l}}\right)}^{3/2}}$;

e. Update${\mathit{\lambda }}_{\mathit{l}+1}={\mathit{\lambda }}_{\mathit{l}}-\frac{{\varphi }_{\mathit{l}}}{{\varphi }_{\mathit{l}}^{\mathit{\text{'}}}}$;

f. Set$\mathit{l}=\mathit{l}+1$;

g. Compute$\mathit{D}=\frac{1}{\mathit{\alpha }+{\mathit{\lambda }}_{\mathit{l}}}\mathit{P}-{\mathit{W}}_{\mathit{I}}$and Cholesky decomposition$\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$;

h. Solve${\mathit{p}}^{\mathit{*}}$from$\mathit{Lp}=\mathit{u}$and${\mathit{v}}^{\mathit{*}}$from${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{p}}^{\mathit{*}}$;

i. Compute${\mathit{u}}^{\mathit{T}}{\mathit{v}}^{\mathit{*}}$, $\mathit{p}={\mathit{W}}_{\mathit{I}}{\mathit{v}}^{\mathit{*}}$, ${\mathit{p}}^{\mathit{T}}{\mathit{v}}^{\mathit{*}}$, and${\mathit{\pi }}_{\mathit{l}}=\mathit{\pi }\left({\mathit{\lambda }}_{\mathit{l}}\right)$usingEq. 20;

j. Set$\mathit{\delta }=\left|1-\frac{\sqrt{{\mathit{\pi }}_{\mathit{l}}}}{\mathit{\Delta }}\right|$.

10. Compute${\mathit{p}}^{\mathit{*}}=\mathit{V}{\mathit{v}}^{\mathit{*}}$and${\mathit{s}}^{\mathit{*}}\mathit{ }=-\frac{1}{\mathit{\alpha }+{\mathit{\lambda }}^{\mathit{*}}}\mathit{g}+\frac{1}{{\left(\mathit{\alpha }+{\mathit{\lambda }}^{\mathit{*}}\right)}^2}{\mathit{p}}^{\mathit{*}}$;

11. End

Let ${\mathit{N}}_{\mathit{MIL}}$ 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, ${\mathit{c}}_5=4\mathit{n}+2\mathit{mn}+\left(9\mathit{m}+12{\mathit{m}}^2+\frac{{\mathit{m}}^3}{3}\right)\left({\mathit{N}}_{\mathit{MIL}}+1\right)$, including the following eight operations: (1) computing $\Vert \mathit{g}\Vert$ in step 1 ($2\mathit{n}$ flops); (2) computing $\mathit{D}=\frac{1}{\mathit{\alpha }+{\mathit{\lambda }}_{\mathit{l}}}\mathit{P}-{\mathit{W}}_{\mathit{I}}$ and Cholesky decomposition $\mathit{D}=\mathit{L}{\mathit{L}}^{\mathit{T}}$ in step 4 or step 9(g) ($2{\mathit{m}}^2+\frac{{\mathit{m}}^3}{3}$ flops); (3) solving ${\mathit{p}}^{\text{*}}$ from $\mathit{L}\mathit{p}=\mathit{u}$, and ${\mathit{v}}^{\text{*}}$ from ${\mathit{L}}^{\mathit{T}}\mathit{v}={\mathit{p}}^{\text{*}}$ in step 5 or step 9(h) ($2{\mathit{m}}^2$ flops); (4) computing ${\mathit{u}}^{\mathit{T}}{\mathit{v}}^{\text{*}}$, $\mathit{p}={\mathit{W}}_{\mathit{I}}{\mathit{v}}^{\text{*}}$ and ${\mathit{p}}^{\mathit{T}}{\mathit{v}}^{\text{*}}$ in step 6 or step 9(i) ($4\mathit{m}+2{\mathit{m}}^2$ flops); (5) computing $\mathit{q}=\frac{1}{{\left(\mathit{\alpha }+{\mathit{\lambda }}_{\mathit{l}}\right)}^2}\mathit{P}{\mathit{v}}^{\text{*}}$ in step 9(a) ($\mathit{m}+2{\mathit{m}}^2$ flops); (6) Solving ${\mathit{p}}^{\text{*}}$ from $\mathit{Lp}=\mathit{q}$ and ${\mathit{w}}^{\text{*}}$ from ${\mathit{L}}^{\mathit{T}}\mathit{w}={\mathit{p}}^{\text{*}}$ in step 9(b) ($2{\mathit{m}}^2$ flops); and (7) computing ${\mathit{u}}^{\mathit{T}}{\mathit{w}}^{\text{*}}$, $\mathit{p}={\mathit{W}}_{\mathit{I}}{\mathit{w}}^{\text{*}}$, and ${\mathit{p}}^{\mathit{T}}{\mathit{v}}^{\text{*}}$ in step 9(c) ($4\mathit{m}+2{\mathit{m}}^2$ flops); and (8) computing ${\mathit{p}}^{\text{*}}=\mathit{V}{\mathit{z}}^{\text{*}}$ and ${\mathit{s}}^{\text{*}}=-\frac{1}{\mathit{\alpha }+{\mathit{\lambda }}^{\text{*}}}\mathit{g}+\frac{1}{{\left(\mathit{\alpha }+{\mathit{\lambda }}^{\text{*}}\right)}^2}{\mathit{p}}^{\text{*}}$ in step 10 ($2\mathit{mn}+2\mathit{n}$).

The total computational cost (flops) used to solve the L-BFGS TRS using Algorithm-5 together with Algorithm-4 is,

${\mathit{c}}_{\mathit{MIL}}={\mathit{c}}_4+{\mathit{c}}_5=6\mathit{n}+6\mathit{mn}+\left(9\mathit{m}+12{\mathit{m}}^2+\frac{{\mathit{m}}^3}{3}\right)\left({\mathit{N}}_{\mathit{MIL}}+1\right)(34)$

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),

$\mathit{\beta }=\frac{{\mathit{c}}_{\mathit{DIQ}}}{{\mathit{c}}_{\mathit{MIL}}}=\frac{2\mathit{m}{\mathit{n}}^2+\left(7+3{\mathit{m}}^2\right)\mathit{n}+\left(5\mathit{n}+2{\mathit{n}}^2+\frac{{\mathit{n}}^3}{3}\right)\left({\mathit{N}}_{\mathit{DIQ}}+1\right)+\mathbf{O}\left({\mathit{m}}^3\right)}{6\mathit{n}+6\mathit{mn}+\left(9\mathit{m}+12{\mathit{m}}^2+\frac{{\mathit{m}}^3}{3}\right)\left({\mathit{N}}_{\mathit{MIL}}+1\right)}(35)$

Because $\mathit{n}\gg \mathit{m}$ for large-scale problems, Eq. 35 can be further simplified as,

$\mathit{\beta }=\frac{{\mathit{c}}_{\mathit{DIQ}}}{{\mathit{c}}_{\mathit{MIL}}}\approx \frac{{\mathit{N}}_{\mathit{DIQ}}+1}{18}\frac{{\mathit{n}}^2}{\mathit{m}}(36)$

For some problems the solution ${\mathit{s}}^{\mathit{*}}=\mathit{s}\left(0\right)$ is accepted with ${\mathit{N}}_{\mathit{DIQ}}={\mathit{N}}_{\mathit{MIL}}=0$. For other problems, it takes roughly ${\mathit{N}}_{\mathit{DIQ}}\approx {\mathit{N}}_{\mathit{MIL}}=5\sim 15$ iterations for a TRS solver to converge. Figure 1A,B illustrate the plots of $\mathit{\beta }$ vs. ${\mathit{n}}^2/\mathit{m}$ for different $\mathit{m}$ by fixing ${\mathit{N}}_{\mathit{DIQ}}={\mathit{N}}_{\mathit{MIL}}=0$ and ${\mathit{N}}_{\mathit{DIQ}}={\mathit{N}}_{\mathit{MIL}}=10$, respectively. As shown in both Figure 1A and Figure 1B, the MIL TRS solver may reduce the computational cost by a factor $\mathit{\beta }=10\sim {10}^5$ for $\frac{{\mathit{n}}^2}{\mathit{m}}=100\sim {10}^6$.

FIGURE 1

FIGURE 1. Plots of the ratio of computational cost β vs. n²/m, computed using Eq. 35 by fixing ${\mathit{N}}_{\mathit{DIQ}}={\mathit{N}}_{\mathit{MIL}}=0$ in (A) and ${\mathit{N}}_{\mathit{DIQ}}={\mathit{N}}_{\mathit{MIL}}=10$ in (B).

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,

$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{R}\mathbf{o}\mathbf{s}\mathbf{e}\mathbf{n}\mathbf{b}\mathbf{r}\mathbf{o}\mathbf{c}\mathbf{k}}={\displaystyle \sum }_{\mathit{i}=1}^{\mathit{n}-1}\left[100{\left({\mathit{x}}_{\mathit{i}+1}-{\mathit{x}}_{\mathit{i}}^2\right)}^2+{\left(1-{\mathit{x}}_{\mathit{i}}\right)}^2\right]\mathit{ },(37)$

$\mathit{f}{\left(\mathit{x}\right)}_{\mathbf{S}\mathbf{p}\mathbf{h}\mathbf{e}\mathbf{r}\mathbf{e}}={\displaystyle \sum }_{\mathit{i}=1}^{\mathit{n}}{\mathit{x}}_{\mathit{i}}^2\mathit{ }(38)$

The Rosenbrock function defined in Eq. 37 has one local minimum located at ${\mathit{x}}_1=-1$ and ${\mathit{x}}_{\mathit{i}}=1$ for $\mathit{i}=\mathrm{2,3},\dots \mathit{n}$ with the objective function being 4 (approximately), and one global minimum located at ${\mathit{x}}_{\mathit{i}}=1$ for $\mathit{i}=\mathrm{1,2},\dots \mathit{n}$ 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, $\mathit{n}$ as powers of two and fixed twice the maximum number of snapshots, $\mathit{M}=2{\mathit{L}}_{\mathit{M}}=10$. 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 ${\mathit{\delta }}_{\mathit{cr}}={10}^{-4}$ and the maximum number of iterations as ${\mathit{N}}_{\mathit{TRS},\mathit{max}}=16$. We set the initial trust region size at ∆ = 0.5.

TABLE 2

TABLE 2. Performance comparison of DNR vs. MIL methods testing on Rosenbrock function.

TABLE 3

TABLE 3. Performance comparison of DIQ vs. MIL methods testing on Rosenbrock function.

TABLE 4

TABLE 4. Performance comparison of DNR vs. MIL methods testing on Sphere function.

TABLE 5

TABLE 5. Performance comparison of DIQ vs. MIL methods testing on Sphere function.

We noticed that the matrix $\mathit{D}$ 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 $\mathit{m}$ 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., ${\mathit{N}}_{\mathit{DIQ}}<{\mathit{N}}_{\mathit{DNR}}$).

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

We can draw the following conclusions based on theoretical analysis and numerical tests:

1) Using the compact representation of the L-BFGS formulation in Eq. 4, the updated Hessian ${\mathit{H}}^{(\mathit{k}+1)}$ is composed of a diagonal matrix $\mathit{\lambda }{\mathit{I}}_{\mathit{n}}$ and a matrix ${\mathit{V}}^{\left(\mathit{k}\right)}{\mathit{W}}^{\left(\mathit{k}\right)}{\left[{\mathit{V}}^{\left(\mathit{k}\right)}\right]}^{\mathit{T}}$ with low rank when m < n.

2) Using the matrix inversion lemma, we are able to reduce the cost of solving the L-BFGS TRS by transforming the original equation with $\mathit{n}$ variables to a new equation with $\mathit{m}$ variables, by taking advantage of the low rank matrix updating feature.

3) Further reduction in computational cost is achieved by updating the vector u^(k) and matrices such as S^(k) and Z^(k), A^(k), B^(k), and C^(k), P^(k), and ${\left[{\mathit{W}}^{\left(\mathit{k}\right)}\right]}^{-1}$ iteratively, which effectively avoids expensive operations such as matrix-matrix and matrix-vector multiplications.

4) Through seamless integration of matrix inversion lemma with the iterative matrix and vector updating approach, the newly proposed TRS solver for the limited-memory distributed BFGS optimization method performs much more efficiently than the DNR TRS solver.

5) The MIL TRS solver can obtain the right solution with acceptable accuracy, without increasing memory usage but reducing the computational cost significantly, as confirmed by numerical results on the Rosenbrock function and Sphere function.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

GG, HF, and JV contributed to conception and algorithmic analysis of the study. HF, TW, FS, and CB contributed to design and numerical validations. GG and FH wrote the first draft of the article. All authors contributed to article revision, read, and approved the submitted version.

Conflict of Interest

Authors GG, HF, and FS were employed by Shell Global Solutions (United States) Inc. Authors JV, TW and CB were employed by Shell Global Solutions International B.V.

References

1. Ai-Mudhafar, WJ, Rao, DN, and Srinivasan, S. Robust Optimization of Cyclic CO2 Flooding through the Gas-Assisted Gravity Drainage Process under Geological Uncertainties. J Pet Sci Eng (2018). 166:490–509. doi:10.1016/j.petrol.2018.03.044

Google Scholar

2. Alpak, FO, Jin, L, and Ramirez, BA. 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 Essen, GM, Zandvliet, MJ, Den Hof, PMJ, Bosgra, OH, Jansen, JD, et al. Robust 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. Liu, Z, and Forouzanfar, F. Ensemble Clustering for Efficient Robust Optimization of Naturally Fractured Reservoirs. Comput Geosci (2018). 22:283–96. doi:10.1007/s10596-017-9689-1

CrossRef Full Text | Google Scholar

5. Liu, X, and Reynolds, AC. Augmented Lagrangian Method for Maximizing Expectation and Minimizing Risk for Optimal Well-Control Problems with Nonlinear Constraints. SPE J (2016). 21(5). doi:10.1007/s10596-017-9689-1

CrossRef Full Text | Google Scholar

6. Liu, X, and Reynolds, AC. Gradient-based Multi-Objective Optimization with Applications to Waterflooding Optimization. Comput Geosci (2016). 20:677–93. doi:10.1007/s10596-015-9523-6

CrossRef Full Text | Google Scholar

7. Oliver, DS, Reynolds, AC, and Liu, N. Inverse Theory for Petroleum Reservoir Characterization and History Matching. Cambridge, United Kingdom: Cambridge University Press (2008). doi:10.1017/cbo9780511535642

CrossRef Full Text

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

CrossRef Full Text

9. Oliver, DS. On Conditional Simulation to Inaccurate Data. Math Geol (1996). 28:811–7. doi:10.1007/bf02066348

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

11. Li, R, Reynolds, AC, and Oliver, DS. History Matching of Three-phase Flow Production Data. SPE J (2003). 8(4):328–40. doi:10.2118/87336-pa

CrossRef Full Text | Google Scholar

12. Sarma, P, Durlofsky, L, and Aziz, K. 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. Gao, G, and Reynolds, AC. An Improved Implementation of the LBFGS Algorithm for Automatic History Matching. SPE J (2006). 11(1):5–17. doi:10.2118/90058-pa

CrossRef Full Text | Google Scholar

14. Nocedal, J, and Wright, SJ. Numerical Optimization. New York City: Springer (1999).

15. Chen, C, et al. Assisted 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. Gao, G, Vink, JC, Chen, C, Alpak, FO, and Du, K. A Parallelized and Hybrid Data-Integration Algorithm for History Matching of Geologically Complex Reservoirs. SPE J (2016). 21(6):2155–74. doi:10.2118/175039-pa

CrossRef Full Text | Google Scholar

17. Gao, G, Vink, JC, Chen, C, El Khamra, Y, Tarrahi, M, et al. Distributed Gauss-Newton Optimization Method for History Matching Problems with Multiple Best Matches. Comput Geosciences (2017). 21(5-6):1325–42. doi:10.1007/s10596-017-9657-9

CrossRef Full Text | Google Scholar

18. Gao, G, Wang, Y, Vink, J, Wells, T, Saaf, F, et al. Distributed 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. Wild, SM. Derivative Free Optimization Algorithms for Computationally Expensive Functions. [PhD thesis]. Ithaca (new york): Cornell University (2009).

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

CrossRef Full Text | Google Scholar

21. Zhao, H, Li, G, Reynolds, AC, and Yao, J. Large-scale History Matching with Quadratic Interpolation Models. Comput Geosci (2012). 17:117–38. doi:10.1007/s10596-012-9320-4

CrossRef Full Text | Google Scholar

22. Audet, C, and Dennis, JE Mesh Adaptive Direct Search Algorithms for Constrained Optimization. SIAM J Optim (2006). 17:188–217. doi:10.1137/040603371

CrossRef Full Text | Google Scholar

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

CrossRef Full Text

24. Chen, C, et al. EUR 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).

CrossRef Full Text | Google Scholar

25. Chen, C, Gao, G, Li, R, Cao, R, Chen, T, Vink, JC, et al. Global-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. doi:10.2118/182639-pa

CrossRef Full Text | Google Scholar

26. Chen, C, Gao, G, Ramirez, BA, Vink, JC, and Girardi, AM. Assisted History Matching of Channelized Models by Use of Pluri-Principal-Component Analysis. SPE J (2016). 21(05):1793–812. doi:10.2118/173192-pa

CrossRef Full Text | Google Scholar

27. Armacki, A, Jakovetic, D, Krejic, N, Krklec Jerinkic, N, et al. Distributed Trust-Region Method with First Order Models. in Paper Presented at the IEEE EUROCON-2019, 18th International Conference on Smart Technologies July 1–14 2019 Novi Sad, Serbia (2019). p. 1–4.

Google Scholar

28. Moré, JJ, and Sorensen, DC. Computing a Trust Region Step. SIAM J Sci Stat Comput (1983). 4(3):553–72. doi:10.1137/0904038

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

30. Golub, GH, and Van Loan, CF. Matrix Computations. 4th ed. The Johns Hopkins University Press (2013).

31. Brust, JJ, Leyffer, S, and Petra, CG. Compact Representations of BFGS Matrices. Preprint ANL/MCS-P9279-0120. Lemont, IL: ARGONNE NATIONAL LABORATORY (2020).

32. Gao, G, Jiang, H, Vink, JC, van Hagen, PPH, Wells, TJ, et al. Performance Enhancement of Gauss-Newton Trust-Region Solver for Distributed Gauss-Newton Optimization Method. Comput Geosci (2020). 24(2):837–52. doi:10.1007/s10596-019-09830-x

CrossRef Full Text | Google Scholar

33. Burdakov, O, Gong, L, Zikrin, S, and Yuan, Y. On Efficiently Combining Limited-Memory and Trust-Region Techniques. Math Prog Comp (2017). 9:101–34. doi:10.1007/s12532-016-0109-7

CrossRef Full Text | Google Scholar

34. Burke, JV, Wiegmann, A, and Xu, L. Limited Memory BFGS Updating in a Trust-Region framework Tech. Rep. Seattle, WA: University of Washington (2008).

35. Byrd, RH, Nocedal, J, and Schnabel, RB. Representations of Quasi-Newton Matrices and Their Use in Limited Memory Methods. Math Programming (1994). 63:129–56. doi:10.1007/bf01582063

CrossRef Full Text | Google Scholar

36. Gao, G, Jiang, H, van Hagen, P, Vink, JC, and Wells, T. A Gauss-Newton Trust Region Solver for Large Scale History Matching Problems. SPE J (2017b). 22(6):1999–2011. doi:10.2118/182602-pa

CrossRef Full Text | Google Scholar

37. Gao, G, Saaf, F, Vink, J, Krymskaya, M, Wells, T, et al. Gauss-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