Gauss-Newton results for Bratu’s 3D problem presented in Section 8.2.

## Abstract

A global regularized Gauss-Newton (GN) method is proposed to obtain a zero residual for square nonlinear problems on an affine subspace built by wavelets, which allows reducing systems that arise from the discretization of nonlinear elliptic partial differential equations (PDEs) without performing a priori simulations. This chapter introduces a Petrov-Galerkin (PG) GN approach together with its standard assumptions that ensure retaining the q-quadratic rate of convergence. It also proposes a regularization strategy, which maintains the fast pace of convergence, to avoid singularities and high nonlinearities. It also includes a line-search method for achieving global convergence. The numerical results manifest the capability of the algorithm for reproducing the full-order model (FOM) essential features while decreasing the runtime by a significant magnitude. This chapter refers to a wavelet-based reduced-order model (ROM) as WROM, while PROM is the proper orthogonal decomposition (POD)-based counterpart. The authors also implemented the combination of WROM and PROM as a hybrid method referred herein as (HROM). Preliminary results with Bratu?s problem show that if the WROM could correctly reproduce the FOM behavior, then HROM can also reproduce that FOM accurately.

### Keywords

- Gauss-Newton method
- line search
- Petrov-Galerkin direction
- data compression
- wavelets

## 1. Introduction

The Newton method for solving square nonlinear problems is one of the most popular techniques used in engineering applications due to its simplicity and fast convergence rate [1, 2, 3]. However, the quality of the final numerical results is affected by the possible Jacobian’s singularity and high nonlinearity. Another drawback is that the method depends on the initial point. Therefore, it is necessary to implement a globalization strategy to get the solution independently of the initial guess. One such approach is the line-search method that relies on a suitable merit function that yields that the iterations progress toward a solution of the problem.

From the numerical point of view, the technique requires solving square linear systems several times, and it is necessary to carry out function evaluations of the order of the problem, as well as first-order information of the square law of the function to compute the Jacobians. In the case of high-dimensional nonlinear problems, the method can overcome the capacity of the computer memory or decrease speed for solving these linear systems, even in the case of a few iterations. One of the current research activities focuses on solving large-scale square nonlinear problems in real time. The purpose of this chapter is to provide an algorithm for solving large-scale square nonlinear problems, in real time, while retaining the fast convergence rate. One strategy for addressing such challenges is to characterize an affine subspace, of much lower dimension than the original one, that contains the initial solution and thus reproduces the problem’s principal features.

One procedure to characterize the affine subspace consists of solving the full-order model (FOM) in several input points whose solutions are called snapshots, then using a principal component analysis method such as singular value decomposition (SVD) to build an orthonormal basis that spans the snapshots’ majority of energy. This oblique subspace, where one seeks a solution, is projected on the original one. Good numerical results have been already reported in the literature [1, 3, 4, 5, 6, 7]. But there are still open questions about this procedure as for how and when to choose the snapshots and their number [8, 9]. It is important to emphasize that at every picture it is required solving the FOM regardless of its cost. This chapter thus promotes a new strategy that is snapshot free. The approach originated in signal processing and consists of using the notion of wavelets to compress data in a subspace of smaller dimension which retains the majority of the original energy [10, 11]. The discrete wavelet’s low-pass matrix is used as the affine subspace; then, the optimization is performed in this compressed subspace to obtain a cheaper solution that can decompress to its original size.

## 2. Reduced-order models using wavelet transformations

### 2.1. Wavelets and data compression

The rationale for using wavelet transformations for model reduction originates from the fields of image processing, image compression, and transform coding. In these fields, massive amounts of information, for example, images, are broadcasted over limited bandwidth communication lines and networks such as the internet. One needs to compress these signals for quick transmission and to diminish storage requirements [10]. In summary, data compression consists of, given a signal

### 2.2. Reduced-order models

Let

By orthogonality,

Choosing

where

## 3. Problem formulation

### 3.1. Statement of the problem

Given a nonlinear function

### 3.2. Overdetermined problem

This section formulates this problem by using the overdetermined functions

The fact that finding a solution

### 3.3. Nonlinear least-squares problem

The residual problem (4) is immediately seen to be equivalent to solving the nonlinear zero-residual least-square problem

### 3.4. First derivatives for the residual functions R and H

The Jacobian of

A direct application of the chain rule, since

### 3.5. First and second derivatives for problem (5)

The gradient of each term of the problem is

The second-order information is

## 4. Gauss-Newton method

This section presents a Gauss-Newton method to solve the nonlinear problem (5) which is equivalent to solving the overdetermined nonlinear composite function (4). It describes the standard Newton assumptions for this composite function problems that yield

### 4.1. Model order-reduction-based Gauss-Newton algorithm

This subsection presents a reduced-order Gauss-Newton algorithm for solving problem (5), which is the interest herein.

**Algorithm 1. Reduced-order Gauss-Newton (ROGN)**

**Inputs:** Given the compressed base

**Output:** Approximate solution in the affine subspace

**1:** Initial point of the problem. Given

**2:** Initial point in the affine subspace.

**3:** For

**4:** Gauss-Newton direction (compressed direction). Solve for

**5:** Compressed update:

**6:** Decompressed update:

**Remarks:**

**1.** The algorithm presents two initial points. The first

**2.** The update

**3.** Finding Gauss-Newton direction

**4.** The Gauss-Newton direction is the Petrov-Galerkin direction obtained by approximating Newton’s direction of square nonlinear problems for the following weighted problem:

### 4.2. Local convergence of reduced Gauss-Newton algorithm

It is known that the Gauss-Newton method retains

**Theorem:** Let **ROGN algorithm 1** is well defined, converges, and has

**Proof:** The residual problem

Since the Jacobian of

Second, one proves that the Jacobian of

Now, since the Jacobian of

Finally, one proves that the smallest eigenvalue of

Let

Therefore,

## 5. Regularization

Despite the advantages of the Gauss-Newton method, the algorithm will not perform well if either the problem is ill conditioned or in the presence of high nonlinearity of some components of it. The purpose of this section is to introduce two regularizations to overcome these difficulties while retaining the fast rate of convergence of the Gauss-Newton method.

### 5.1. Levenberg-Marquardt method

To prevent the Gauss-Newton algorithm to preclude in case some eigenvalues are near zero or in case of rank deficiency of the linear systems to solve, the least-squares directions are regularized by

where

Under the standard Gauss-Newton assumptions written before and choosing the regularization parameter as

### 5.2. Scaling regularization

To avoid the influence of the high order of magnitude of some components with respect to the rest of the components of the problem, one presents the following regularization:

where

This regularization prevents that large components of the problem affect the behavior of the algorithm. It is important to observe the Lipschitz constant of the problem is improved. Considering the preceding two regularizations, one has

This last regularizations prevent the smallest eigenvalue affecting the behavior of the Gauss-Newton algorithm and at the same time, through rescaling, components with small values are not considered by the influence of large components while retain its fast rate of convergence.

## 6. Globalization strategy

The good performance of the Gauss-Newton algorithm depends on a suitable initial point that must be inside its region of convergence. Rather than absorbing the computational cost associated with choosing an appropriate initial point, the chapter proposes a line-search method that provides convergence for initial points outside of the region of convergence. The goal of this approach is to obtain a sufficient decrease in the merit function. If the direction fails, then a backtracking is used until a sufficient reduction is obtained. A merit function should allow moving toward a solution of the problem.

### 6.1. Merit function

It is natural to think that the merit function for the unconstrained minimization problem (5) is itself. That is:

### 6.2. Descent direction

One proves that Gauss-Newton direction is a descent direction for the merit function

**Property:** The regularized Gauss-Newton direction

**Proof:** One proves that the directional derivative of

since

Consequently, it is possible to progress toward a solution of the problem in the

and

for fixed values

It is important to realize that these two inequalities can be reached by using a back-tracking procedure. Therefore, this work uses a line-search strategy to satisfy the inequalities. Next section proposes a line-search regularized Gauss-Newton algorithm for solving the zero-residual composite function problem.

## 7. A line-search regularized Gauss-Newton method

This section proposes the following regularized Gauss-Newton method with line search to find a solution on the affine subspace

**Algorithm 2: A reduced-order regularized Gauss-Newton (RORGN)**

**Input:** Given the compressed base

**Output:** The approximate solution in the affine subspace

**1:** Initial point of the problem. Given

**2:** Initial point in affine subspace.

**3:** For

**4:** Choose

**5:** Regularized Gauss-Newton direction. Solve for

**6:** Line search (sufficient decrease). Find

**7:** Update.

**Remarks:** The algorithm is amenable to the use of any suitable basis, not necessarily a wavelet basis. The algorithm can be tested with different initial displacement points. On the other hand, the election of the initial point of the algorithm is not limited to the origin.

## 8. Numerical examples

The authors run on a MacBook Pro laptop equipped with an Intel(R) Quad-Core(TM) i7-2720QM CPU @ 2.20GHz and 8 GB of RAM. Section 8.2 presents Bratu’s 3D problem. This problem is more challenging since the nonlinear systems become ill conditioned where one approaches the bifurcation point. Therefore, the RORGN algorithm was tested to solve them efficiently. All these Bratu’s problems utilized wavelets-based ROM. One takes the regularization by

### 8.1. Bratu’s 2D problem

The Bratu’s 1D equation can be generalized by replacing the second derivative by a Laplacian [1, 17]. This section numerically studies the nonlinear diffusion equation with exponential source term in two and three dimensions. Let

and

One can discretize (32) by means of central finite differences on regular tensor product meshes. Homogenous Dirichlet boundaries are enforced conditions in all square or cube faces, see [17] for details.

### 8.2. Bratu’s 3D problem

Figures 2 and 3 present results from Bratu’s 3D problem. Figure 2 shows the parameter continuation problem while Figure 3 compares FOM and WROM results in the whole domain. For visualization purposes, one cuts away the front half of the cube to see inside it. Figure 2 shows the FOM in continuous line while dashed blue and magenta lines correspond to WROM at 85 and 90%, respectively. The mesh size is

Newton method | #Iter | Success | ||
---|---|---|---|---|

N = 512, | ||||

Standard | 1.637075E−04 | 4.497266E−08 | 7 | True |

Combined | 5.089393E−04 | 8.391213E−07 | 5 | True |

Regularized | 5.089393E−04 | 8.391213E−07 | 5 | True |

Scaled | 1.637075E−04 | 4.497266E−08 | 7 | True |

Line search | 1.637075E−04 | 4.497266E−08 | 7 | True |

Com. and Line search | 5.089393E−04 | 8.391213E−07 | 5 | True |

N = 1000, | ||||

Standard | 3.303599E−04 | 1.682251E−07 | 7 | True |

Combined | 1.424983E−05 | 7.450719E−10 | 6 | True |

Regularized | 1.424983E−05 | 7.450719E−10 | 6 | True |

Scaled | 3.303599E−04 | 1.682251E−07 | 7 | True |

Line search | 3.303599E−04 | 1.682251E−07 | 7 | True |

Com. and Line search | 1.424983E−05 | 7.450719E−10 | 6 | True |

N = 1728, | ||||

Standard | 5.792054E−04 | 4.806531E−07 | 7 | True |

Combined | 1.197055E−04 | 6.096723E−08 | 5 | True |

Regularized | 1.197055E−04 | 6.096723E−08 | 5 | True |

Scaled | 5.792054E−04 | 4.806531E−07 | 7 | True |

Line search | 5.792054E−04 | 4.806531E−07 | 7 | True |

Com. and Line search | 1.197055E−04 | 6.096723E−08 | 5 | True |

N = 2744, | ||||

Standard | 1.143068E−05 | 1.756738E−10 | 8 | True |

Combined | 1.678400E−04 | 1.052920E−07 | 6 | True |

Regularized | 1.678400E−04 | 1.052920E−07 | 6 | True |

Scaled | 1.143068E−05 | 1.756738E−10 | 8 | True |

Line search | 1.143068E−05 | 1.756738E−10 | 8 | True |

Com. and Line search | 1.678400E−04 | 1.052920E−07 | 6 | True |

One observes for all ranks reported herein that the regularized method provides convergence tolerances likewise but it usually spent two iterations less than standard Newton and hence CPU time reduces as well. One notices for this particular problem that line search is equivalent to standard Newton as well as combined matches the performance of the combined plus line-search method.

### 8.3. Nonlinear benchmark problems

One also considers a benchmark nonlinear problem from the literature in order to challenge the proposed algorithms. The Yamamura [18] problem is a nonlinear system of equations defined by:

where

The objective is to challenge the algorithms presented in this research, and the results are reported in Table 2. One can infer from this table that the standard Newton method could not converge in any of these realizations except

Newton method | #Iter | Success | ||
---|---|---|---|---|

N = 32 | ||||

Standard | 6.643960E−06 | 1.282327E−07 | 47 | True |

Regularized | 3.547835E−05 | 8.310947E−07 | 32 | True |

Line search | 8.835078E−09 | 2.415845E−13 | 78 | True |

Reg. and Line search | 4.271577E−06 | 8.288108E−09 | 35 | True |

N = 256 | ||||

Standard | 2.992259E−01 | 3.018468E+07 | 128 | False |

Regularized | 4.260045E−06 | 5.643326E−08 | 21 | True |

Line search | 2.569192E−02 | 2.618397E+03 | 128 | False |

Reg. and Line search | 2.893999E−06 | 2.380562E−08 | 35 | True |

N = 512 | ||||

Standard | 2.280299E−02 | 3.448757E+05 | 128 | False |

Regularized | 2.384802E−07 | 2.579578E−10 | 46 | True |

Line search | 1.127484E−01 | 1.845466E+07 | 128 | False |

Reg. and Line search | 3.551655E−08 | 2.368659E−11 | 42 | True |

N = 1024 | ||||

Standard | 9.540053E−05 | 4.501800E+01 | 128 | False |

Regularized | 1.730306E−06 | 8.414453E−09 | 41 | True |

Line search | 7.199610E−05 | 5.095672E+00 | 128 | False |

Reg. and Line search | 8.857816E−06 | 4.250476E−07 | 35 | True |

N = 2048 | ||||

Standard | 1.769950E+00 | 8.504131E+14 | 128 | False |

Regularized | 7.270952E−08 | 1.578029E−10 | 68 | True |

Line search | 6.195869E−02 | 4.180347E+09 | 128 | False |

Reg. and Line search | 1.069042E−02 | 3.229689E−01 | 128 | False |

## 9. Hybrid method: HROM

The idea is simple since the wavelet subspace is not a function of a priori known snapshots, it can be determined without executing the so-called, computationally expensive, off-line stage, in which one thoroughly studies the FOM and can sample the input space to record a representative set of snapshots. HROM considers as input snapshots those outputted by a WROM procedure. As shown below, if this WROM happens to reproduce the FOM behavior correctly, then one should expect the resulting PROM, that is, HROM, to replicate the original FOM behavior accurately. At first, glance, depending upon the WROM compression ratio, this procedure reduces the runtime of the off-line stage.

This section conducts a series of preliminary numerical experiments on the well-known Bratu’s nonlinear benchmark problem, in particular in one and two dimensions [5] to sell this case. Figures 4 and 5 depict results for the 1D continuation problem. They utilize the following WROM compression ratios: 10 and 5%, where the compression ratio is constant during the continuation problem. All plots display two distinct HROM compression ratios. The WROM at 20% that is not depicted completely misses the bifurcation zone and the second branch thus the HROM is also way off. However, as the WROM starts to catch up with the FOM then HROM too does. Indeed, it is observed that HROM yields comparable results when comparing it to the version that takes the original snapshots, that is, PROM.

One can repeat a similar experiment with Bratu’s 2D continuation problem, which produces the same trend as before. Indeed, if the input WROM is way off targeting then, HROM is off as well. For instance, 27% compression implies that all HROM miss the second branch, but still, the 21% model could slightly reproduce the proper FOM trend. Things significantly improve on models with 21 and 20% as shown in Figures 6 and 7. However, these models still miss the bifurcation zone. They just render a flat profile there. These insights suggest that if the WROM can correctly reproduce the FOM behavior, then HROM can do so. One is probably able to improve the accuracy of the WROM by changing the compression ratio during the online stage. For Bratu’s problem, one can assume a graded energy distribution which provides more of it while approaching the bifurcation point. This approach is referred as “adaptive WROM” or AWROM as shorthand.

This section presents an alternative strategy that can rely on insights from the FOM, such as Newton tolerances and number of iterations, which are in turn indirect error estimates. Figure 8 plots the energy distribution that was utilized as a function of the continuation parameter,

where

Figures 12 and 13 depict preliminary results of HROM applied over a couple of AWROM models whose energy distribution was described in Figures 8 and 10, respectively

## 10. Concluding remarks

This chapter introduced a global regularized Gauss-Newton method for resolving square nonlinear algebraic systems in an affine subspace that enable real-time solutions without a priori simulations. Solving a nonlinear least-squares composite function poses the problem where the answer is derived from the inside argument. The authors thus presented the standard Newton assumptions that guarantee a

From the numerical results presented in Section 9 for Bratu’s 1D and 2D problems, the data fusion procedure (HROM) can be used as an alternative procedure when the simulation time for a problem can be limited. This method uses two different sceneries. In the first, one implies that the data come through out of the model, while in the second, the data are obtained independently of the latter. That is, in the first case the model governs the data, and the other the input information rules the model.

## Acknowledgments

The authors recognize the financial support of the project “Reduced-Order Parameter Estimation for Underbody Blasts” financed by the Army Research Laboratory, through the Army High-Performance Computing Research Center under Cooperative Agreement W911NF-07-2-0027, and also acknowledge Dr. Martine Ceberio at UTEP for proofreading the manuscript.