A Reduced-Order Gauss-Newton Method for Nonlinear Problems Based on Compressed Sensing for PDE Applications

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.

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 $x\in {\mathbb{R}}^n$ find a lower dimensional signal $\widehat{x}\in {\mathbb{R}}^r$ with $r\ll n$ to broadcast or store those lower dimensional ones $\widehat{x}.$ The most widely used techniques for data compression are based on wavelets transform. The key to using wavelets is to find a lower dimensional signal $\widehat{x}$ that relies on a known subspace properly denoted as energy compaction. It is well comprehended that the wavelets tend to accumulate energy in the low-frequency sub-band of the wavelet decomposition [3, 12, 13, 14]. The energy relates to the $L^2$-norm and is defined as $\mathrm{ϵ}=\parallel x{\parallel }^2={\sum }_{i=1}^n{x}_i^2.$ To demonstrate the energy compaction, consider Figure 1, in which one has the original image $x\in {\mathbb{R}}^{512\times 512}.$ Notice that the upper left quadrant of the wavelet decomposition is a low-dimensional approximation $\widehat{x}\in {\mathbb{R}}^{256\times 256}$ that is one-fourth the size of the original signal and resembles the low-frequency sub-band wavelet coefficients. Using the previously measured energy, the energy enclosed in $\widehat{x}$ is $95.75\%$ using just one-fourth of the coefficients. Since $\widehat{x}$ comprises most of the energy, a simple data compression scheme would execute all of the other wavelet sub-bands to zero and store only the low-frequency information as in Figure 1 (see in the bottom center). By only employing $\widehat{x}$, one can reproduce an approximation of $x\in {\mathbb{R}}^n$ by its generalized inverse of the sub-band compression [10]. Next section will provide details.

Figure 1.

Sample compression and decompression.

2.2. Reduced-order models

Let $W\in {\mathbb{R}}^{n\times n}$ describe an orthonormal wavelet. The transformation encompasses a low-pass and a high-pass submatrix that is given by

$$W_{\mathbf{n}\times \mathbf{n}}=\left[\frac{L_{\mathbf{r}\times \mathbf{n}}}{H_{\mathbf{s}\times \mathbf{n}}}\right],\mathbf{r}+\mathbf{s}=\mathbf{n}.$$E1

By orthogonality, $\mathit{rank}\left(W\right)=n,$ $\mathit{rank}\left(L\right)=r,$ and $\mathit{rank}\left(H\right)=s.$ Now, by orthogonality of W, ${\mathit{WW}}^T=I,$ then ${\mathit{LL}}^T=I_r,$ ${\mathit{HH}}^T=I_s$. The energy of a signal $x\in {\mathbb{R}}^n$ is

$$\parallel x{\parallel }^2=\parallel \mathit{Lx}{\parallel }^2+\parallel \mathit{Hx}{\parallel }^2.$$E2

Choosing $L$ that contains the majority of the energy such that the energy $\parallel \mathit{Hx}\parallel \approx 0,$ one has $\parallel x\parallel \approx \parallel \mathit{Lx}\parallel .$ This means that the energy of the original data $x\in {\mathbb{R}}^n$ is approximately equal to the energy to the compressed data $\widehat{x}=\mathit{Lx}\in {\mathbb{R}}^r.$ That is: $\parallel x\parallel \approx \parallel \widehat{x}\parallel .$ The decompressed data are $L^T\widehat{x}=\tilde{x}\in {\mathbb{R}}^{\mathbf{n}}.$ On the other hand, the compressed and decompressed energies are equal. That is $\parallel \widehat{x}\parallel =\parallel \tilde{x}\parallel$ since ${\mathit{LL}}^T=I_r.$ Therefore, the original, compressed, and decompressed data are related as follows

$$\parallel x\parallel \approx \parallel \widehat{x}\parallel =\parallel \tilde{x}\parallel ,$$E3

where $x\in {\mathbb{R}}^n$ is the original data, $\widehat{x}\in {\mathbb{R}}^r$ is the compressed data, and $\tilde{x}\in {\mathbb{R}}^n$ is the decompressed data. Thus, once an appropriate low-pass submatrix $L$ is determined, one proposes solving the corresponding optimization problem in the reduced affine subspace determined by $L^T$ and later coming back to the original size by its generalized inverse.

3. Problem formulation

3.1. Statement of the problem

Given a nonlinear function $R$ from ${\mathbb{R}}^n$ to ${\mathbb{R}}^n$, find a solution in an affine subspace determined by an initial displacement point $x_o\in {\mathbb{R}}^n$ and an orthonormal base $L_{\mathbf{n}\times \mathbf{r}}^T,$ $\mathbf{r}<\mathbf{n}.$ That is: find $x^{\ast }\in {\mathbb{R}}^n$ with $R\left(x^{\ast }\right)=0$ and $x^{\ast }\in x_o+\eta \left(L^T\right)$, where $\eta \left(L^T\right)$ is the subspace generated by the linear combination of the operator $L^T$.

3.2. Overdetermined problem

This section formulates this problem by using the overdetermined functions $H$ and $\varphi$ from ${\mathbb{R}}^r$ to ${\mathbb{R}}^n$

$$H\left(p\right)=R\left(\varphi \left(p\right)\right)=0,\varphi \left(p\right)=x_o+L^{T}p,\mathit{\text{and}}p\in {\mathbb{R}}^r.$$E4

The fact that finding a solution $p^{\ast }$ of the overdetermined problem draws attention, $H\left(p^{\ast }\right)=0$ for $p^{\ast }\in {\mathbb{R}}^r,$ is equivalent to finding the solution one is initially seeking. That is: $x^{\ast }=\varphi \left(p^{\ast }\right)$ and $R\left(x^{\ast }\right)=0$ is a solution on the affine subspace. Therefore, one studies the problem by finding a zero residual of the nonlinear least-squares problem associated to $H$. Problem (4) is called an overdetermined zero-residual problem.

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

$$\text{minimize}f\left(p\right)=\frac{1}{2}\sum _{i=1}^n{r}_i^2\left(\varphi \left(p\right)\right),p\in {\mathbb{R}}^r.$$E5

3.4. First derivatives for the residual functions $\mathit{R}$ and $\mathit{H}$

The Jacobian of $R$ at $x\in {\mathbb{R}}^n$ is given by

$$J_R\left(x\right)=J\left(x\right)={\left[{\left(\nabla r_i\left(x\right)\right)}^T\right]}_{1\le i\le n}.$$E6

A direct application of the chain rule, since $J_{\varphi }\left(p\right)=L^T$ yields:

$$J_H\left(p\right)=J\left(\varphi \left(p\right)\right)L^T.$$E7

3.5. First and second derivatives for problem (5)

The gradient of each term of the problem is $\nabla r_i^2\left(\varphi \left(p\right)\right)=2L^T{r}_i\left(\varphi \left(p\right)\right)\nabla r_i\left(\varphi \left(p\right)\right)$. Therefore the gradient of $f\left(p\right)$ is

$$\nabla f\left(p\right)=({\left(J\left(\varphi \left(p\right)\right)L^T\right)}^{T}R\left(\varphi \left(p\right)\right).$$E8

The second-order information is

$${\nabla }^{2}f\left(p\right)={\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(\mathit{J\varphi }\left(p\right)\right)L^T)+\sum _{i=1}^n{h}_i\left(p\right){\nabla }^2{h}_i\left(p\right)).$$E9

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 $q$-quadratic rate of convergence. The inconvenience to use Newton method is that the second-order information associated with the Hessian method is not easily accessible or is impractical for computational time. The latter makes the Newton method impractical for very large-scale problems.

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 $L^T\in {\mathbb{R}}^{n\times r},$ and an initial displacement $x_o\in {\mathbb{R}}^n.$

Output: Approximate solution in the affine subspace $x\in {\mathbb{R}}^n.$

1: Initial point of the problem. Given $p_o\in {\mathbb{R}}^r$.

2: Initial point in the affine subspace. $\varphi \left(p_o\right)=x_o+L^T{p}_o\in {\mathbb{R}}^n$.

3: For $k=0:$ until convergence $(\parallel R\left(\varphi \left(p_k\right)\parallel \le \mathrm{ϵ}\right)$.

4: Gauss-Newton direction (compressed direction). Solve for $\mathrm{\Delta }{p}_{k+1}$

$${\left(J\left(\varphi \left(p_k\right)\right)L^T\right)}^T\left(\mathit{J\varphi }\left(p_k\right)\right)L^T)\Delta p_{k+1}=-{\left(J\left(\varphi \left(p_k\right)\right)L^T\right)}^{T}R\left(\varphi \left(p_k\right)\right).$$E10

5: Compressed update: $p_{k+1}=p_k+\Delta p_{k+1}$.

6: Decompressed update: $\varphi \left(p_{k+1}\right)=x_o+L^T{p}_{k+1}$.

Remarks:

1. The algorithm presents two initial points. The first $x_o$ is the displacement to characterize the affine subspace, and the second one $p_o$ is the initial point for the algorithm.

2. The update $\varphi \left(p_{k+1}\right)$ is the approximation to the solution one is looking for which one denotes by $x_{k+1}$.

3. Finding Gauss-Newton direction $\Delta p_{k+1}$ is equivalent to solving the following linear least-squares problem:

$$\underset{\Delta p_{k+1}}{\text{min}}\left\{\frac{1}{2}\parallel \left(\mathit{J\varphi }\left(p_k\right),L^T\right)\mathit{\Delta p}+R\left(\varphi \left(p_k\right)\right){\parallel }^2\right\}.$$E11

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:

$$\underset{\Delta p_{k+1}}{\text{min}}\left\{\frac{1}{2}\parallel (L^T\mathit{\Delta p}+R\left(x\right){\parallel }_{Q^{-1}}^2,Q=J{\left(x\right)}^{T}J\left(x\right)>0\right\}.$$E12

4.2. Local convergence of reduced Gauss-Newton algorithm

It is known that the Gauss-Newton method retains $q$-quadratic rate of convergence under standard assumptions for zero-residual single-function problems [15]. The natural question is: What are the standard assumptions that guarantee the Gauss-Newton conditions for the composite function one is working with, that conserve $q$-quadratic rate of convergence? The next theorem establishes these assumptions.

Theorem: Let $H$ from ${\mathbb{R}}^r$ to ${\mathbb{R}}^n$ be defined by $H\left(p\right)=R\left(x_o+L^{T}p\right),$ $x_o\in {\mathbb{R}}^n,$ $p\in {\mathbb{R}}^r,$ and $L\in {\mathbb{R}}^{r\times n}$ are orthonormal operators with $r<n.$ Assume there exists a solution $p^{\ast }\in \tilde{D}\subset {\mathbb{R}}^r$, with $\tilde{D}$ convex and open. Define $D=\left\{x_o\right\}+L^T\left(\tilde{D}\right)$, where $L^T\left(\tilde{D}\right)$ is the image of $\tilde{D}$ under $L^T\in {\mathbb{R}}^{n\times r}$. Assume that $J_R\in L_{\gamma }\left(D\right),$ $J_R$ is bounded on $D,$ and the minimum eigenvalue of $J{\left(x^{\ast }\right)}^{T}J\left(x^{\ast }\right)$ is positive. Then, the sequence $\left\{p_{k+1}\right\}$ given ROGN algorithm 1 is well defined, converges, and has $q$-quadratic rate of convergence. That is:

$$\parallel p_{k+1}-p^{\ast }\parallel \le \frac{1}{2}\parallel p_k-p^{\ast }\parallel$$E13

$$\parallel p_{k+1}-p^{\ast }\parallel \le \frac{\tilde{c}\tilde{\gamma }}{2\tilde{\sigma }}\parallel p_k-p^{\ast }{\parallel }^2\mathit{\text{with}}\tilde{c},\tilde{\gamma },\tilde{\sigma }\in {\mathbb{R}}_{+}.$$E14

Proof: The residual problem $R$ has solution $x^{\ast }$ on $D$. First, one proves that the Jacobian of $H$ is Lipchitz on $\tilde{D}.$ Since $J_H\left(p\right)=J\left(\varphi \left(p\right)\right)L^T,$ then

$$\begin{array}{l}\parallel J_H\left(p_1\right)-J_H\left(p_2\right)\parallel =\parallel \left(J\left(\varphi \left(p_1\right)\right)-J\left(\varphi \left(p_2\right)\right)\right)L^T\parallel ,\mathit{\text{for}}{p}_1,p_2\in \tilde{D}.\end{array}$$E15

Since the Jacobian of $R$ is Lipschitz on $D,$ one concludes

$$\begin{array}{l}\parallel J_H\left(p_1\right)-J_H\left(p_2\right)\parallel \le \tilde{\gamma }\parallel p_1-p_2\parallel ,\tilde{\gamma }={\gamma }^{\ast }\parallel L^T{\parallel }^2,\mathit{\text{for}}{p}_1,p_2\in \tilde{D}.\end{array}$$E16

Second, one proves that the Jacobian of $H$ is bounded on $\tilde{D}.$ Since the Jacobian of $H$ at $p$ is $J(\varphi \left(p\right)L^T$, then

$$\parallel J_H\left(p\right)\parallel =\parallel J(\varphi \left(p\right)L^T\parallel \mathit{\text{for}}p\in \tilde{D}.$$E17

Now, since the Jacobian of $R$ is bounded on $D,$ one concludes

$$\parallel J_H\left(p\right)\parallel \le \tilde{c},\tilde{c}=c^{\ast }\parallel L^T\parallel \mathit{\text{for}}p\in \tilde{D}.$$E18

Finally, one proves that the smallest eigenvalue of $J_H{\left(p^{\ast }\right)}^T{J}_H\left(p^{\ast }\right)$ is greater than zero.

$$J_H{\left(p^{\ast }\right)}^T{J}_H\left(p^{\ast }\right)={\left(J\left(x^{\ast }\right)L^T\right)}^T\left(J\left(x^{\ast }\right)L^T\right)\mathit{\text{with}}{x}^{\ast }=x_o+L^T{p}^{\ast }.$$E19

Let $p\ne 0\in {\mathbb{R}}^k$ and $\sigma \in \mathbb{R}$ be an eigenvector and eigenvalue associated with the last symmetric matrix. Then

$$\parallel L^{T}p{\parallel }_Q^2=\sigma \parallel p{\parallel }^2,Q=J{\left(x^{\ast }\right)}^{T}J\left(x^{\ast }\right)>0.$$E20

Therefore, $\sigma >0$ since $L^T$ is a full rank and $p\ne 0$. The convergence and its fast rate of convergence given by the last two inequalities follow from the Theorem 10.2.1 in the Dennis and Schnabel book [15].

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

$$\begin{array}{l}\underset{\mathit{\Delta p}}{\text{min}}\left\{\frac{1}{2}\parallel \left(J\left(\varphi \left(p\right)\right)L^T\right)\mathit{\Delta p}+R\left(\varphi \left(p\right)\right){\parallel }^2+\frac{\mu }{2}\parallel \mathit{\Delta p}{\parallel }^2\right\}\end{array}$$E21

where $\mu >0$. The solution is given by

$$\begin{array}{l}\left({\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(J\left(\varphi \left(x\right)\right)L^T\right)+\mathit{\mu I}\right)\mathit{\Delta p}=-{\left(J\left(\varphi \left(p\right)\right)L^T\right)}^{T}R\left(\varphi \left(p\right)\right).\end{array}$$E22

Under the standard Gauss-Newton assumptions written before and choosing the regularization parameter as $\mu =O(\parallel {\left(J\left(\varphi \left(p\right)\right)L^T\right)}^{T}R\left(\varphi \left(p\right)\parallel \right)$, the regularized Gauss-Newton algorithm converges and the $q$-quadratic rate of convergence is retained; see Theorem 10.2.6 [15].

$$\parallel p_{k+1}-p^{\ast }\parallel \le \frac{1}{2}\frac{\tilde{\alpha }\tilde{\gamma }}{\tilde{\lambda }}\parallel p_k-p^{\ast }{\parallel }^2.$$E23

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:

$$\begin{array}{l}\underset{\mathit{\Delta p}}{\text{min}}\left\{\frac{1}{2}\parallel \left(J\left(\varphi \left(p\right)\right)L^T\right)\mathit{\Delta p}+R\left(\varphi \left(p\right)\right){\parallel }^2+\frac{\sigma }{2}\parallel L^T\mathit{\Delta p}{\parallel }_Q^2\right\},\end{array}$$E24

where $Q={\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(J\left(\varphi \left(p\right)\right)L^T\right)$. The solution is given by

$$\begin{array}{l}\left({\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(J\left(\varphi \left(p\right)\right)L^T\right)\right)\mathit{\Delta p}=-{\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\frac{R\left(\varphi \left(p\right)\right)}{1+\sigma \parallel R\left(\varphi \left(p\right)\right)\parallel }.\end{array}$$E25

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

$$\begin{array}{l}\left({\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(J\left(\varphi \left(p\right)\right)L^T\right)+\mathit{\mu I}\right)\mathit{\Delta p}=-{\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\frac{R\left(\varphi \left(p\right)\right)}{1+\sigma \parallel R\left(\varphi \left(p\right)\right)\parallel }.\end{array}$$E26

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: $M\left(p\right)=f\left(p\right)$.

6.2. Descent direction

One proves that Gauss-Newton direction is a descent direction for the merit function $M\left(p\right)=f\left(p\right)$.

Property: The regularized Gauss-Newton direction $\mathit{\Delta p}$ given by (26) is a descent direction for the merit function $M\left(p\right)=f\left(p\right)$.

Proof: One proves that the directional derivative of $f$ at the direction $\mathit{\Delta p}$ is less than zero. The gradient of $f\left(p\right)$ is given by $\nabla f\left(x\right)={\left(J\left(\varphi \left(p\right)\right)L^T\right)}^{T}R\left(\varphi \left(p\right)\right)$. Therefore

$$\nabla f{\left(p\right)}^T\mathit{\Delta p}=-\parallel {\left(J\left(\varphi \left(p\right)\right)L^T\right)}^{T}R\left(\varphi \left(p\right)\right){\parallel }_{Q^{-}1}^2<0$$E27

since $Q=\left({\left(J\left(\varphi \left(p\right)\right)L^T\right)}^T\left(J\left(\varphi \left(p\right)\right)L^T\right)+\mu \right)$ is positive definite.

Consequently, it is possible to progress toward a solution of the problem in the $\mathit{\Delta p}$ direction. The purpose is to find a step length $\alpha \in \left(0,1\right]$ that yields a sufficient decrease. To that effect one follows the Armijo-Goldstein conditions given by

$$f\left(p+\mathit{\alpha \Delta p}\right)\le f\left(p\right)+\alpha \left({\lambda }^{\ast }\nabla f{\left(p\right)}^T\mathit{\Delta p}\right)$$E28

and

$$\nabla f{\left(p+\mathit{\alpha p}\right)}^T\mathit{\Delta p}\ge \beta \nabla f{\left(p\right)}^T\mathit{\Delta p},$$E29

for fixed values $\lambda ,\beta \in \left(0,1\right).$ The first inequality allows sufficient decrease of the merit function, and the second one avoids step lengths that are very small. It is important to observe that if $\beta$ is chosen, $\beta \in \left[\lambda ,1\right],$ then the two inequalities can be satisfied simultaneously. Wolfe proved that if $f$ is continuously differentiable on ${\mathbb{R}}^r,$$\mathit{\Delta p}$ is a descent direction, and assuming the set $\{f\left(p+\mathit{\alpha \Delta p}\right);\alpha \in (0,1\left]\right\}$ is bounded below, then there exists an ${\alpha }^{\ast }\in \left(0,1\right]$ such that the two inequalities be satisfied simultaneously [16].

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 $x_o+\eta \left(L^T\right)$ for problem (4).

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

Input: Given the compressed base $L^T\in {\mathbb{R}}^{n\times r},$ and a displacement $x_o\in {\mathbb{R}}^n$.

Output: The approximate solution in the affine subspace $x\in {\mathbb{R}}^n.$

1: Initial point of the problem. Given $p_o\in {\mathbb{R}}^r.$

2: Initial point in affine subspace. $x_1=x_o+L^T{p}_o\in {\mathbb{R}}^n.$

3: For $k=1:$ until convergence ($\parallel R\left(x_k\right)\parallel \le \epsilon .$).

4: Choose ${\mu }_k={\sigma }_k\parallel R\left(x_k\right)\parallel$, and ${\sigma }_k\in \left(0,1\right].$

5: Regularized Gauss-Newton direction. Solve for $\Delta p_k$

$$(\left({\left(J\left(x_k\right)L^T\right)}^T\left(J\left(x_k\right)L^T\right)+{\mu }_{k}I\right)\Delta p_k=-\frac{{\left(J\left(x_k\right)L^T\right)}^{T}R\left(x_k\right)}{1+{\sigma }_k\parallel R\left(x_k\right)\parallel }.$$E30

6: Line search (sufficient decrease). Find ${\alpha }_k\in \left(0,1\right]$ such that

$$\parallel R\left(x_k+{\alpha }_k{L}^T\Delta p_k\right){\parallel }^2<\parallel R\left(x_k\right){\parallel }^2+2\cdot {10}^{-4}{\alpha }_k\nabla f{\left(p_k\right)}^T\Delta p_k.$$E31

7: Update. $x_{k+1}=x_k+{\alpha }_k{L}^T\Delta p_k$.

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 ${\mu }_k={\sigma }_k\Vert \mathbf{R}\left(x_k\right)\Vert ,$ with ${\sigma }_k\in \left(0,1\right)$. One employs as stopping criteria for the algorithm; either the norm of the residual, ${\mathrm{ϵ}}_k=\parallel \mathbf{R}\left(x_k\right)\parallel$, is less than some small positive real value given, $\mathrm{ϵ}$, or a maximum number of iterations reached, $k_{\mathit{max}}$.

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 $\mathrm{\Omega }={\left[0,1\right]}^n,n=2,3$ be a unitary square or cube, where $x_i\in \left[0,1\right],i=1,\dots ,n$ are the spatial variables while $n$ is the space dimension.

$$\begin{array}{c}\mathrm{\Delta }u+\zeta \cdot e^u=0\mathrm{on}\mathrm{\Omega };u=u\left(\underset{¯}{x}\right),\\ u=0\mathrm{on}\mathrm{\partial \Omega },\end{array}$$E32

and $\zeta \in \mathbb{R}$ is a coefficient. The Laplacian is defined by

$$\mathrm{\Delta }\left(\cdot \right)=\sum _{i=1}^n\frac{{\partial }^2\left(\cdot \right)}{\partial x_i^2}.$$E33

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 $10\times 10\times 10$ and $\zeta =1.5$ is fixed. In the figure, one has from left to right the FOM, the WROM, and the absolute error. Neither of these WROM models could reproduce the FOM behavior beyond $\zeta >9$. They could not get into neither the second branch nor close to the bifurcation point, where the system becomes highly nonlinear. On the other hand, Figure 3 compares FOM versus WROM at 90% in order to show that these WROM could properly reproduce the FOM behavior in the whole cube. Table 1 summarizes the performance of the family of Gauss-Newton algorithms applied to FOM Bratu’s 3D problem. One employs these numerical values, ${\mathrm{ϵ}}_{\mathit{tol}}={10}^{-3}$ and $k_{\mathit{max}}=32$. Once again, one gets close to the bifurcation point by choosing, $\zeta =9.9$, to pose a challenging nonlinear system while $\sigma$ was tuned to achieve performance for a given rank.

Figure 2.

FOM vs. WROM parameter continuation solutions of (32), for $n=3$, are shown.

Figure 3.

FOM and ROM are compared: FOM (left), WROM at 90% (center), and error (right).

Newton method	${\mathrm{ϵ}}_k$	${\delta }_k$	#Iter	Success
N = 512, $\sigma =0.03$
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, $\sigma =0.025$
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, $\sigma =0.020$
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, $\sigma =0.015$
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

Table 1.

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

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:

$$\begin{array}{l}\mathbf{R}:{\mathbb{R}}^n\to {\mathbb{R}}^n,x\in {\mathbb{R}}^n,R_i\left(x\right)=0,1\le i\le n,\\ R_i\left(x\right)=2.5x_i^3-10.5x_i^2+11.8x_i-i+\sum _{i=1}^n{x}_i=0.\end{array}$$E34

where $n$ is the size of the nonlinear system. One implemented the algorithms with ${\mathrm{ϵ}}_{\mathit{tol}}={10}^{-6}$, $k_{\mathit{max}}=128,$ and $\sigma =0.025$.

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 $n=32$. Conversely, the regularized method converged for all realizations. This latter method outperformed all others. On the other hand, the regularized and line-search method could consistently converge for all realizations but $n=2048$. For larger ranks, that is, 512 and 1024, the latter method is the more efficient bottom line; the regularized method performed well in this example.

Newton method	${\mathrm{ϵ}}_k$	${\delta }_k$	#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

Table 2.

Gauss-Newton results for Yamamura problem.

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.

Figure 4.

Bratu?s 1D, 10% compression.

Figure 5.

Bratu?s 1D, 5% compression.

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.

Figure 6.

Bratu?s 2D, 10% compression.

Figure 7.

Bratu?s 2D, 5% compression.

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, $\zeta$, for Bratu’s 1D problem. When one approaches the bifurcation point, $\zeta =3.5$, one should gradually bump up energy as shown. Notice that the distribution tends to concentrate more points toward the bifurcation point. With this energy distribution into account, one obtains the AWROM results in Figure 9. This ROM accurately reproduces the FOM behavior as noted. Let now conduct the following experiment. One must run the FOM, and at every snapshot, one needs to store the number of Newton iterations and the resulting error tolerance. One then computes the energy distribution that Figure 10 depicts. The following formula was employed to do so:

$$E_i=\left(1-\theta \right)+\theta \left({\mathit{\text{nIters}}}_i/\mathit{\text{nMaxIters}}\right),$$E35

where ${\mathit{\text{nIters}}}_i$ is the number of iterations at the current location, and $\mathit{\text{nMaxIters}}$ is the maximum number of iterations reported by the FOM and $\theta \in \left[0,1\right]$. Figure 11 depicts excellent accordance between FOM and AWROM $\left(\theta =0.2\right)$. One should expect that if this last AWROM model is inputted for an HROM simulation, HROM should reproduce the original FOM correctly.

Figure 8.

Linear energy distribution.

Figure 9.

10–20% variable compression.

Figure 10.

Variable energy distribution.

Figure 11.

Variable compression.

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 $\left(\theta =0.2\right)$. These ROM reproduce the FOM behavior accordingly, which proves that there is potential to study the performance of HROM further. Another important question that arises from further research is how to improve the compression ratio of the AWROM scheme.

Figure 12.

Linear AWROM.

Figure 13.

Variable AWROM.

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 $q$-quadratic rate of convergence. The findings include that the Petrov-Galerkin projection directions for the Newton method are no other than the Gauss-Newton ones for a composite function. The technique uses two initial points, one that determines the affine subspace and the other is the starting guess for solving the composition mentioned earlier. The notion of compressed sensing with wavelets produces the characterization of an affine subspace that comprises the majority of the energy of the problem. The chapter showed some numerical experimentations that back up the proposed globalization methodology for solving highly nonlinear dynamic systems in real time. These last ones reproduce the principal features of the FOM. Results underline the fact that one does not need to employ information at any particular point. The Bratu’s 3D FOM results prove that the proposed RORGN algorithm outperforms the standard GN method while retaining its $q$-quadratic rate of convergence. This chapter concludes that the regularized and line-search-enabled scheme is the most robust and efficient algorithm for the problems presented herein. The numerical results imply that this approach performs well and it does not significantly increase the CPU time.

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.