A New Broyden rank reduction method to solve large systems of nonlinear equations

Article History: Received 04 November 2018 Accepted 03 February 2019 Available 09 April 2019 We propose a modification of limited memory Broyden methods, called dynamical Broyden rank reduction method, to solve high dimensional systems of nonlinear equations. Based on a thresholding process of singular values, the proposed method determines a priori the rank of the reduced update matrix. It significantly reduces the number of singular values decomposition calls of the update matrix during the iterations. Local superlinear convergence of the method is proved and some numerical examples are displayed.


Introduction
Let us consider the problem of finding a solution of the system of nonlinear equations F (x) = 0, F : R n → R n . ( The mapping F is assumed to fulfill the following classical assumptions (CA): • it is continuously differentiable in an open convex set D ⊂ R n , • there is an x * in D such that F (x * ) = 0, • the Jacobian F ′ is Lipschitz continuous at x * and F ′ (x * ) is nonsingular.
Newton's method (see [1][2][3]) is well suited to solve the system (1) due to its local quadratic convergence.However, this method is known to be numerically expensive.It requires the evaluation of a jacobian matrix and the solution of a linear system per iteration.An alternative to Newton's method is the Broyden's quasi-Newton method.This method uses approximations to the Jacobian matrix at each iteration by performing rank-one updates, see [4].It requires only one F -evaluation per iteration and achieves, under the classical hypotheses (CA), local superlinear convergence as shown in [5].Given an initial guess x 0 and an initial approximation B 0 of the jacobian matrix, Broyden iteration is given by where B k is updated at each iteration as *Corresponding Author with y k = F (x k+1 ) − F (x k ) and s k = x k+1 − x k .
To avoid the drawback of storing and manipulating the (n×n)-matrices of Broyden, limited memory methods put restrictions on the size of systems to solve, see references [6], [7] and [8].Equation (4) implies that if B 0 is updated p times, the resulting matrix can be written as where c j+1 = y j − B j s j s j and d j+1 = s j s j , j = 0, ..
where σ 1 ≥ . . .≥ σ p ≥ 0 are the singular values.The sets {u 1 , . . ., u p } and {v 1 , . . ., v p } are, respectively, the left and right corresponding singular vectors.By choosing B 0 = I, the matrix B p can be stored using 2p vectors of length n.Suppose the maximal rank of Q is fixed at p with p << n.Broyden rank reduction (BRR) method, see [7], is a variant of limited memory Broyden methods that approximates Q by a matrix Q with low rank q ≤ p−1 by truncating the singular value decomposition (6).
The current number of stored updates is denoted by m (0 ≤ m ≤ p).The maximal number of updates to the initial Broyden matrix is thus given by p.If p updates are stored (m = p), the singular value decomposition of Q is computed and the last (p − q) singular values are removed just before the next update is computed.In this case, the next Broyden update will proceed using the modified matrix where R = p l=q+1 σ l u l v T l is the so-called roundoff matrix.No reduction is performed as long as m < p.Thus, the liberated memory is used to store (p − q) new Broyden's updates according to the scheme ( 7) The BRR method, as presented in [7], does not give any idea how to fix a priori the rank of the matrix Q, only the smallest singular value is removed.But, in many cases there are more than one singular value that are close to zero and so they can be removed.In this case, memory will be free to store more than one Broyden's update.We propose here a new approach using a thresholding process of singular values of the update matrix by fixing a relative accuracy for the approximation of the matrix Q.In section 2 we present the new method and prove its local superlinear convergence.Section 3 is devoted to numerical results showing the efficiency of the method.

The proposed method
In many nonlinear problems, the singular values of the update matrix decay rapidly to zero, and more than one singular value can be removed.In figure 3 (see problem 3 in section 3) we present the singular values distribution of Q in case of the Spedicato function for p = 6 and p = 10.For example, when p = 6, we can see that the two last smallest singular values are zero while the third and the fourth ones are close to zero.In this example, the four last singular values can advantageously be removed as memory will be available to store four new Broyden updates and no singular values decomposition will be needed during the following four iterations.So, the question is how, in general, to choose the rank q of Q and thence the number of singular values to remove.As an answer to this question, we propose to use a tresholding process by exploiting the information about the approximation error R 2 .Given a relative accuracy ε > 0 of the approximation Q, i.e., the required rank q(ε) is given, if it exists, by The value of q(ε) is calculated each time a reduction of the update matrix is needed, and all singular values satisfying the condition σ l+1 < εσ 1 , l = 1, . . ., p − 1 are removed.If there is no k satisfying (8), only the smallest singular value is removed (we turn back to the BRR method), and in this case q = p − 1.So, the rank q will be chosen as This dynamical choice of q leads to a new method, displayed in Algorithm 1, which will be called dynamical Broyden rank reduction (DBRR) method.
In DBRR algorithm, the Sherman Morrison Woodbury formula, see [1], is used in order to compute the inverse of the Broyden matrix.If we set B 0 = I we get In this case, we have only to solve linear systems of equations with p × p matrices.Note also that the singular value decomposition of the update matrix is carried out using an economical process, see [7].Note that the error that is introduced by removing the last (p − q) singular values of Q equals Let us now prove the superlinear convergence of the proposed algorithm.
Theorem 1.Let q be defined as in (9) with where the constant α does not depend on k.Let, in addition, the classical hypotheses (CA) hold.
Then the DBRR method has local superlinear convergence.
) and e k = x k −x * for k = 0, 1, . ... To prove linear convergence of the proposed method, we need the following lemma whose proof can be found in [1], p. 77.
Equation ( 7) can be written as Using Lemma 1 we have where γ is the constant Lipschitz for F ′ .Thus we obtain, from equation ( 11), the so-called bounded deterioration property of the DBRR method since We have used the fact that is an orthogonal projection and then its norm is equal to one.Inequality (13) implies local convergence of the DBRR method.In fact, as shown in [5], any quasi-Newton method that obeys the bounded deterioration property has local linear convergence.
As consequences of the linear convergence, we have and Since DBRR method satisfies the secant equation (B k+1 s k = y k ), according to Theorem 8.2.4 in [1], a sufficient condition for {x k } k to converge superlinearly to the root x * is the so-called Dennis-Moré condition Equation (11) also implies Algorithm 1.Let x 0 ∈ R n and B 0 ∈ R n×n be given.Set the parameter p and the accuracy ε > 0.
Let m = 0. : The following lemma, where the proof can be found in [1], p. 183, will be useful in the sequel.
Lemma 2. Let s ∈ R n be nonzero and let E ∈ R n×n .Then Using inequality (14) and Lemma 2 we derive This inequality is equivalent to Using inequality (15) we show that where c > 0 is an upper bound of the sequence {E k } k .Hence, condition ( 16) is satisfied and the superlinear convergence is proved.

Numerical results
We present now numerical tests by applying the proposed method to some classical test functions from the literature and we present a comparison of this method with the classical BRR method (q = p − 1).The numerical experiments were carried out using the scientific computing software MATLAB.We use the following stopping criterion for our computer programs where ε a = 10 −15 and ε r = 10 −15 (respectively, absolute and relative tolerances).For both BRR and DBRR methods most of the computational time is spent in evaluations of the function F and computation of the singular values decomposition of the update matrix.

Problem 1
Let us consider the trigonometric system The size of this problem is n = 1000000, and the initial guess is given by x 0 = (1.2, . . ., 1.2) T .We plot in figure 1 the distribution of singular values of the update matrix for p = 5, 8, 10 and p = 15.For these values of p, the update matrix has rank two in all nonlinear iterations.Hence, the DBRR method requires less singular values decomposition calls.The size of this problem is n = 1000000, and the initial guess is given by x 0 = (1.2, . . ., 1.2) T .We plot in figure 1 the distribution of singular values of the update matrix for p = 5, 8, 10 and p = 15.For these values of p, the update matrix has rank two in all nonlinear iterations.Hence, the DBRR method requires less singular values decomposition calls.
For this example both BRR and DBRR do not converge for p ≤ 4. Performances of DBBR and BRR methods, for different values of p, are presented in tables 1 and 2, respectively.For all p values, the choice of ε does not significantly affect the convergence rate and the computational time.

Problem 2
We consider now the so-called extended system of Byeong The size of this problem is n = 1000000, and the initial guess is given by x 0 = (0.0087, . . ., 0.0087) T .For this example, the update matrix has rank one during all nonlinear iterations, see figure 2.
So, whenever a reduction of the update matrix is needed, (p − 1) singular values are removed freeing the memory to store new updates.For a given value of p, the choice of the parameter ε does not affect the performance of DBRR method as shown in table 3. A comparison of tables 3 and 4 shows the efficiency of the singular values thresholding process.

Problem 3
In this example, we compute the root of the socalled Spedicato function The size of this problem is n = 1000000, and the initial guess is given by x 0 = (−1.2, . . ., −1.2) T .Both BRR and DBRR methods do not converge for p ≤ 4. The rank of Q increases with the nonlinear iterations as shown in figure 3. Performances of DBBR and BRR methods are presented in tables 5 and 6, respectively.

Problem 4
We consider the nonlinear convection-diffusion partial differential equation    We set C = 20 and u 0 = 0. Equation ( 18) is discretized on a 500 × 500 grid using centered differences.Evolution of the sigular values of Q is displayed in figure 4 for different values of p.
For p ≤ 5 the rank of Q is p, and only the p th singular value is removed.In this case, BRR and DBRR methods behave similarly, see tables 7 and 8.For p = 10 the rank of Q is ten but condition ( 8) is satisfied in the first iterations for ε = 10 −1  since, in this case, q = 6 at iteration 11 and q = 8 at both iterations 15 and 17.For p = 20 the rank of Q increases with nonlinear iterations, but condition ( 8) is satisfied for ε = 10 −1 , 10 −3 , 10 −5 .

Conclusion
We have introduced a new Broyden rank reduction method to solve systems of nonlinear equations.The method is based on a thresholding

Figure 1 .
Figure 1.Singular values of the update matrix Q for p = 5, 8, 10 and p = 15 for Problem 1.

Figure 2 .
Figure 2. Singular values of the update matrix Q for p = 3, 5, 10 and p = 15 for Problem 2.

Figure 3 .
Figure 3. Singular values of the update matrix Q for p = 6, 7, 10 and p = 15 for Problem 3.

Figure 4 .
Figure 4. Singular values of the update matrix Q for p = 5, 10, 15 and p = 20 for Problem 4.

Mohammed
Ziani received his PhD in applied Mathematics at the University of Rennes (France) and Mohammed V University in Rabat (Morocco).He is professor of applied Mathematics at Faculty of sciences in Rabat.His research interests include issues related to solving large nonlinear systems of equations and partial differential equations based image processing.An International Journal of Optimization and Control: Theories & Applications (http://ijocta.balikesir.edu.tr)This work is licensed under a Creative Commons Attribution 4.0 International License.The authors retain ownership of the copyright for their article, but they allow anyone to download, reuse, reprint, modify, distribute, and/or copy articles in IJOCTA, so long as the original authors and source are credited.To see the complete license contents, please visit http://creativecommons.org/licenses/by/4.0/.

Table 1 .
Performance of the DBRR method for Problem 1.

Table 2 .
Performance of the BRR method for Problem 1.

Table 3 .
Performance of the DBRR method for Problem 2.

Table 4 .
Performance of the BRR method for Problem 2.

Table 5 .
Performance of the DBRR method for Problem 3.

Table 6 .
Performance of the BRR method for Problem 3.

Table 7 .
Performance of the DBRR method for Problem 4.

Table 8 .
[8]formance of the BRR method for Problem 4.process of the Broyden matrix singular values.All singular values of the update matrix that are smallest than a given threshold are removed.We have proved the local superlinear convergence and numerically tested the method on a variety of problems.As compared to classical Broyden rank reduction, our method induces a significant reduction of the execution time (CPU time) by reducing the number of calls of the singular values decomposition.As a perspective of this work is to combine the proposed method with that presented in[8].Souissi is Professor at the faculty of sciences, Mohammed V University, Rabat, Morocco.His main areas of interest are Scientific Computing, Optimal shape Design and Numerical Analysis of Partial Differential Equations.He received the PhD in Applied Mathematics from Laval University, Canada.He has publications in well-known journals and supervised many PhD theses.