A Simplified Algorithm for a Full-Rank Update Quasi-Newton Method

1. Introduction

A mathematical model $\mathsf{\varphi }\left(\omega ,\mathit{x}\right)$ with some parameters $\mathit{x}=\left[x_i\right]$ $\left(i=1,\dots ,n\right)$ is constructed for an observed system which gives observable response $\mathit{D}=\left[D_j\right]$ sampled at locations ${\omega }_j$ $\left(j=1,\dots ,m\right)$ ($m>n$) to an observable external effect. It is assumed that the simulated system response $\mathsf{\varphi }\left(\omega ,\mathit{x}\right)$ is sensitive to the perturbations $▵x_i$ of all parameters. The n adjustable parameters $\mathit{x}$ of the mathematical model are determined so that the distance

$$l_q={∥\mathit{f}\left(\mathit{x}\right)∥}_q={∥\mathsf{\varphi }\left(\omega ,\mathit{x}\right)-\mathit{D}∥}_q$$

(1.1)

between the observed and the simulated system responses ($\left[D_j\right]$ and $\mathsf{\varphi }\left(\omega ,\mathit{x}\right)$) is minimized, where $l_q:{\mathbb{R}}^n\to \mathbb{R}$ is a q norm

$$l_q={\left({\displaystyle \sum _{j=1}^m}{\left|f_j\left(\mathit{x}\right)\right|}^q\right)}^{{\displaystyle \frac{1}{q}}}$$

(1.2)

If $q=2$ is chosen, then we get the Euclidean or least-squares norm

$$l_2={\left({\displaystyle \sum _{j=1}^m}{\left|f_j\left(\mathit{x}\right)\right|}^2\right)}^{{\displaystyle \frac{1}{2}}}$$

(1.3)

that is well computable and safely usable for a great class of problems. Then minimizing the distance 1.1 leads to the problem of minimizing the least squares norm

$$l_2={∥\mathit{f}\left(\mathit{x}\right)∥}_2={∥\mathsf{\varphi }\left(\omega ,\mathit{x}\right)-\mathit{D}∥}_2$$

(1.4)

($\mathit{x}\in$${\mathbb{R}}^n$, $\mathit{f}:{\mathbb{R}}^n\to {\mathbb{R}}^m$ and $n<m$ ). The minimum value of this norm is zero, but for “real life” problems it is not reachable in most cases due to modeling inaccuracies and measurement noises. The suggested procedure gives a least-squares solution to the overdetermined system of nonlinear equations

$$\mathit{f}\left(\mathit{x}\right)=\mathbf{0}$$

(1.5)

where the solution ${\mathit{x}}^{*}$ minimizes the least squares norm 1.4 of the nonlinear residual function

$$\mathit{f}\left(\mathit{x}\right)=\mathsf{\varphi }\left(\omega ,\mathit{x}\right)-\mathit{D}$$

(1.6)

Least squares minimizations are successively made within the actual iteration steps. The basic concept of solving Equation 1.5 is that the function $\mathit{f}\left(\mathit{x}\right)$ is “linearized” by repeatedly replacing it with a linear function ${\mathit{y}}_p=\mathit{y}\left({\mathit{x}}_p\right)$ as

$$\mathit{f}\left(\mathit{x}\right)\simeq {\mathit{y}}_p={\mathit{f}}_p+{\mathit{J}}_p\left(\mathit{x}-{\mathit{x}}_p\right)$$

(1.7)

for the Newton’s method, where ${\mathit{f}}_p=\mathit{f}\left({\mathit{x}}_p\right)$ is the function value, ${\mathit{J}}_p=\mathit{J}\left({\mathit{x}}_p\right)$ is the Jacobian matrix of the function $\mathit{f}\left(\mathit{x}\right)$ (p is iteration counter) and

$$\mathit{f}\left(\mathit{x}\right)\simeq {\mathit{y}}_p={\mathit{f}}_p+{\mathit{S}}_p\left(\mathit{x}-{\mathit{x}}_p\right)$$

(1.8)

for the quasi-Newton methods methods, where ${\mathit{S}}_p=\mathit{S}\left({\mathit{x}}_p\right)$ is the finite-difference Jacobian approximate of the function $\mathit{f}\left(\mathit{x}\right)$ . Then the nonlinear problem is solved by solving a series of linear problems

$${\mathit{y}}_p=\mathbf{0}$$

(1.9)

and it follows that the new approximate ${\mathit{x}}_{p+1}$ to the solution ${\mathit{x}}^{*}$ can be given from Equations

$${\mathit{J}}_p▵{\mathit{x}}_p=-{\mathit{f}}_p$$

(1.10)

(Newton’s method) and

$${\mathit{S}}_p▵{\mathit{x}}_p=-{\mathit{f}}_p$$

(1.11)

(quasi-Newton methods) where the iteration stepsize is

$$▵{\mathit{x}}_p={\mathit{x}}_{p+1}-{\mathit{x}}_p$$

(1.12)

in both cases.

The system of simultaneous multi-variable nonlinear Equations 1.5 can be solved by the Newton’s method when the derivatives ${\mathit{J}}_p$ of $\mathit{f}\left(\mathit{x}\right)$ are available analytically and a new iterate ${\mathit{x}}_{p+1}$ can be determined. In many cases, explicit formulas for the function $\mathit{f}\left(\mathit{x}\right)$ are not available and the Jacobian ${\mathit{J}}_p$ can only be approximated. The partial derivatives of the analytic Jacobian may be replaced by suitable finite-difference quotients (discretized Newton’s iteration [12,25]) with properly chosen step sizes. For most problems, Newton’s method using analytic derivatives and Newton’s method using properly chosen divided differences are virtually indistinguishable [12]. However, the determination of the finite-difference stepsizes is not clearly defined. The suggested full-rank update procedure helps to overcome this deficit.

Let ${\mathit{S}}_p$ be a full column rank matrix (all column vectors of matrix ${\mathit{S}}_p$ are linearly independent). As $n<m$ (there are more rows than columns), ${\mathit{S}}_p$ is not invertible and the solution to the overdetermined system of Equation 1.11 doesn’t exist. However, ${\mathit{S}}_p^T{\mathit{S}}_p$ is invertible and the pseudo-inverse

$${\mathit{S}}_p^{+}={\left({\mathit{S}}_p^T{\mathit{S}}_p\right)}^{-1}{\mathit{S}}_p^T$$

(1.13)

of matrix ${\mathit{S}}_p$ is a unique left inverse which allows finding a least-squares approximate solution $\tilde{\mathit{x}}$ to Equation 1.11. The pseudo-inverse can be determined in different ways (e.g. rank factorization [14], singular value decomposition [31]). Singular value decomposition is a widely used technique, and it has an advantage, that it gives solution even if ${\mathit{S}}_p$ is singular (not a full column rank matrix). Matrix ${\mathit{S}}_p$ can be factorized as

$${\mathit{S}}_p={\mathit{U}}_p{\Sigma }_p{\mathit{V}}_p^T$$

(1.14)

where ${\mathit{U}}_p$ ($mxn$ matrix) and ${\mathit{V}}_p$ ($nxn$ matrix) are orthogonal matrices with orthonormal column eigenvectors, and

$${\Sigma }_p=diag\left({\sigma }_{i,p}\right)$$

(1.15)

$\left(i=1,\dots ,n\right)$ with ${\sigma }_{i,p}$ singular values of matrix ${\mathit{S}}_p$ and ${\sigma }_1\ge \dots \ge \sigma$. Then the least square solution to Equation 1.11 exists and the iteration stepsize is

$$▵{\mathit{x}}_p=-{\mathit{S}}_p^{+}{\mathit{f}}_p$$

(1.16)

where

$${\mathit{S}}_p^{+}={\mathit{V}}_p{\Sigma }_p^{-1}{\mathit{U}}_p^T$$

(1.17)

is the pseudoinverse of ${\mathit{S}}_p$ and

$${\Sigma }_p^{-1}=diag\left({\displaystyle \frac{1}{{\sigma }_{i,p}}}\right)$$

(1.18)

The pseudoinverse ${\left[\dots \right]}^{+}$ gives a unique least squares solution to Equation 1.11 for which the $q=2$ norm (Euclidean norm) 1.4 of the residual function 1.6 is minimum. If the rank of matrix ${\mathit{S}}_p$ is less than n (the column vectors of ${\mathit{S}}_p$ are not linearly independent and some singular values are zero), then unique solution to Equation 1.11 doesn’t exist. The spectral condition number

$$\kappa \left({\mathit{S}}_p\right)={\displaystyle \frac{{\sigma }_1}{{\sigma }_n}}\ge 1$$

(1.19)

of matrix ${\mathit{S}}_p$ measures the linear dependency of the column vectors. If they are linearly independent (non of the singular values are zero) and the condition number $\kappa$ is not much larger than one, then the matrix is well-conditioned. This is a desirable situation through the whole iteration process. Then the column vectors of matrix ${\mathit{S}}_p$ are in “general positions” and the Ortega-Rheinboldt condition [25] satisfies. If $\kappa \gg 1$ then the matrix is ill-conditioned and the solution $▵{\mathit{x}}_p$ (iteration stepsize) to the Equation 1.11 may be sensitive to small changes in ${\mathit{S}}_p$ or ${\mathit{f}}_p$.

Classic quasi-Newton linearization methods generate a sequence of improved iterates ${\mathit{x}}_p$$\left(p=0,1,2,\dots \right)$ so that the next approximate ${\mathit{x}}_{p+1}$ to the solution ${\mathit{x}}^{*}$ is determined on the basis of rank-1 update ${\mathit{S}}_{p+1}$ of the Jacobian approximate ${\mathit{S}}_p$. Such linearization update is not unique if only a single new approximate ${\mathit{x}}_{p+1}$ exists. In single variable case $\left(n=1\right)$, there are two possible updates (secant lines) from which the selection is not obvious (see details in Section 3). In multi variable case, only rank-1 update is possible if one single new approximate ${\mathit{x}}_{p+1}$ exists.

The solution to Equation 1.11 presented in [2] and applied to the identification of physically nonlinear dynamic systems [3,4,5] was based on the Wolfe-Popper procedure [29,42]. The secant equation 1.11 was splitted into two equations by introducing an auxiliary variable. The simplified solution presented in Section 4 directly solves Equation 1.11 with the pseudo-inverse 1.13 of matrix ${\mathit{S}}_p$. The simplified algorithm given in Section 5 gives an unconventional and simple strategy to improve the performance of Newton and quasi-Newton iterations for the solution of system of nonlinear equations 1.5. It can be considered both as a discretized Newton’s method or as a quasi-Newton method with full-rank update of the Jacobian approximates.

The rate of convergence $\alpha$ of the classic secant method equals to the golden section ratio $\phi$ ($\alpha =\phi \simeq 1.618$) for simple root in single variable case. The suggested method gives an additional new independent approximate ${\mathit{x}}_{p+1}^B$ to the solution x*, then rank-n linearization update can be made with the classic approximate ${\mathit{x}}_{p+1}$ (notated as ${\mathit{x}}_{p+1}^A$ in the followings) and the suggested new approximate ${\mathit{x}}_{p+1}^B$. The results of single- and multi-variable numerical examples are given in Section 6. The efficiency of the proposed method is discussed in Section 7 and it is compared with other classic rank-one update and line-search methods on the basis of available test data. Concluding remarks are summarized in Section 8.

2. Notations

Vectors and matrices are denoted by bold-face letters. Subscripts refer to components of vectors and matrices, superscripts A and B refer to interpolation base points. Notations A and B are introduced to be able to clearly distinguish between the two new approximates $x^A$ and $x^B$. Vectors and matrices may also be given by their general elements. ▵ refers to a difference of two elements. x and $\mathit{X}$ denotes unknown quantities, f and $\mathit{F}$ denotes function values and matrices. t and $\mathit{T}$ denotes multiplier scalars and scaling transformation matrices. e, $\epsilon$ and E denotes approximate error, p is iteration counter, $\alpha$ is convergence rate, ${\epsilon }^{*}$ is termination criterion. n is the number of unknowns, m is the number of function values, i, j and k and are running indexes of matrix columns and rows. Superscripts S and $TS$ refer to the traditional Secant-method and to the suggested full-rank update method (T-Secant) receptively.

3. Linearization Methods

The origin of the secant method can be traced back to ancient time to the “rule of double false position” method, described in the ${18}^{th}$ century B.C. on the Egyptian Rhind Papyrus [28] and it predates Newton’s method by more than 3000 years. It is also well-known that the local convergence rate of the classic Newton’s and the classic secant methods for simple roots are quadratic and super-linear respectively [8,11,26,38]. Weerakson and Fernando’s third order Newton variant [40] involves evaluation of a function value and two derivatives, similarly like in Traub’s suggestion [38]. Secant method variations with improved convergence have been reported in literature [13,15,16,17,23,27,33,34,40]. Mueller’s method [24] uses a quadratic approximation. A class of variations employs two function value and a derivative evaluation [37]. Third order methods are obtained by a second derivative evaluation, such as Halley’s [23] and the “super Halley” [1] method. Further improvements with higher degree polynomial approximations are proposed by Chen [9], Kanwar [17], Zhang [43] and Wang [39].

3.1. Single-Variable Case

The zero $x^{*}$ of a scalar nonlinear function $f\left(x\right)$ ($x\to f\left(x\right)$, $x\in {\mathbb{R}}^1$, $f:{\mathbb{R}}^1\to {\mathbb{R}}^1$) has to be determined, where

$$f\left(x^{*}\right)=0$$

(3.1)

Linearization means that $f\left(x\right)$ is locally replaced by a linear function $y\left(x\right)$ ($x\to y\left(x\right)$, $x\in {\mathbb{R}}^1$, $y:{\mathbb{R}}^1\to {\mathbb{R}}^1$) (secant- or tangent line) as shown on Figure 1 (Left) and operations are made on the linear function $y_p\left(x\right)$ in the $p^{th}$ iteration. $y_p\left(x\right)$ may be defined through two points $A_p\left(x_p^A{f}_p^A\right)$ and $B_p\left(x_p^B{f}_p^B\right)$ of the nonlinear function $f\left(x\right)$, where $f_p^A=f\left(x_p^A\right)$ and $f_p^B=f\left(x_p^B\right)$ are function values, or through one point $A_p\left(x_p^A{f}_p^A\right)$ and by the “slope”

$$S_p={\displaystyle \frac{▵f_p}{▵x_p}}$$

(3.2)

of the secant line $y_p\left(x\right)$, where $▵x_p=x_p^B-x_p^A$ and $▵f_p=f_p^B-f_p^A$ are differences, or by the “slope”

$$S_p=f^{\prime }\left(x_p^A\right)={\displaystyle \frac{df}{dx}}$$

(3.3)

at $x_p^A$ of the tangent line $y_p\left(x\right)$, where $f^{\prime }$ is the derivative function of $f\left(x\right)$. The successive local replacement of the nonlinear function $f\left(x\right)$ by a secant- or tangent line $y_p\left(x\right)$ gives a simple and efficient numerical root finding procedure. It follows from the condition

$$y_p\left(x_{p+1}^A\right)=0$$

(3.4)

that the zero

$$x_{p+1}^A=x_p^A-{\displaystyle \frac{f_p^A}{S_p}}$$

(3.5)

of the linear function $y_p\left(x\right)$ approximates the zero $x^{*}$ of the nonlinear function $f\left(x\right)$ and the new point $A_{p+1}\left(x_{p+1}^A{f}_{p+1}^A\right)$ can be determined for the next iteration as shown on Figure 1 (Right). With the iteration step length

$$▵x_p^A=x_{p+1}^A-x_p^A$$

(3.6)

from point $A_p$, Equation 3.5 can be rewritten as

$$S_p▵x_p^A=-f_p^A$$

(3.7)

that is called secant equation or we may also call it “linearization” equation as it also includes Newton’s method as a special case (tangent line). The secant based procedure has an advantage that it doesn’t need the calculation of function derivatives, it only uses function values and the order of asymptotic convergence is super-linear with convergence rate ${\alpha }^S\cong 1.618...$. However, the “slope” of the updated secant line $y_{p+1}\left(x\right)$ can be determined in two ways from the secant (quasi-Newton) conditions as

$$S_{p+1}▵x_p^A=f_{p+1}^A-f_p^A$$

(3.8)

or

$$S_{p+1}▵x_p^B=f_{p+1}^A-f_p^B$$

(3.9)

where

$$▵x_p^B=x_{p+1}^A-x_p^B$$

(3.10)

is the iteration step length from point $B_p$. The decision is far not obvious, as shown on Figure 1 (Right). The tangent line based iteration (Newton’s method) has an advantage that the “slope”

$$S_{p+1}=f^{\prime }\left(x_{p+1}^A\right)$$

(3.11)

of the updated tangent line is always well-defined. The iteration then continues with updated secant- or tangent line $y_{p+1}\left(x\right)$.

3.2. Multi-Variable Case

The scalar linearization procedure can be extended to multi-dimensions. The zero ${\mathit{x}}^{*}$ of a nonlinear vector-valued function

$$f\left(\mathit{x}\right)=\left[\begin{array}{c}{f}_1\left(\mathit{x}\right)\\ ⋮\\ f_m\left(\mathit{x}\right)\end{array}\right]$$

(3.12)

of n variables

$$\mathit{x}=\left[\begin{array}{c}{x}_1\\ ⋮\\ x_n\end{array}\right]$$

(3.13)

($\mathit{x}\to \mathit{f}\left(\mathit{x}\right)$, $\mathit{x}\in {\mathbb{R}}^n,\mathit{f}:{\mathbb{R}}^n\to {\mathbb{R}}^m,n\le m$) has to be determined, where

$$\mathit{f}\left({\mathit{x}}^{*}\right)=0$$

(3.14)

Then $\mathit{f}\left(\mathit{x}\right)$ is locally replaced by an n-dimensional hyperplane

$$\mathit{y}\left(\mathit{x}\right)(\mathit{x}\to \mathit{y}\left(\mathit{x}\right),\mathit{x}\in {\mathbb{R}}^n,\mathit{y}:{\mathbb{R}}^n\to {\mathbb{R}}^m$$

(3.15)

$\left(n\le m\right)$ and operations are made on the hyperplane ${\mathit{y}}_p\left(\mathit{x}\right)$ in the $p^{th}$ iteration. ${\mathit{y}}_p\left(\mathit{x}\right)$ may be defined through $n+1$ points $A_p\left({\mathit{x}}_p^A{\mathit{f}}_p^A\right)$ and $B_{k,p}\left({\mathit{x}}_{k,p}^B{\mathit{f}}_{k,p}^B\right)$ ($k=1,\dots ,n$) of the nonlinear function $\mathit{f}\left(\mathit{x}\right)$, where ${\mathit{f}}_p^A=\mathit{f}\left({\mathit{x}}_p^A\right)$ and ${\mathit{f}}_{k,p}^B=\mathit{f}\left({\mathit{x}}_{k,p}^B\right)$ ($k=1,\dots ,n$) are function values, or by one point $A_p\left({\mathit{x}}_p^A{\mathit{f}}_p^A\right)$ and by the “slope” (divided differences)

$${\mathit{S}}_p=▵{\mathit{F}}_p▵{\mathit{X}}_p^{-1}$$

(3.16)

$\left(i=1,\dots ,n\right)$, $\left(j=1,\dots ,m\right)$, $\left(k=1,\dots ,n\right)$ of the n-dimensional secant hyperplane ${\mathit{y}}_p\left(\mathit{x}\right)$, where $▵{\mathit{x}}_{k,p}={\mathit{x}}_{k,p}^B-{\mathit{x}}_p^A$ and $▵{\mathit{f}}_{k,p}={\mathit{f}}_{k,p}^B-{\mathit{f}}_p^A$ are difference vectors with n and m components respectively, and

$$▵{\mathit{F}}_p=\left[▵{\mathit{f}}_{k,p}\right]=\left[{\mathit{f}}_{k,p}^B-{\mathit{f}}_p^A\right]=\left[\begin{array}{ccc}▵f_{1,1,p}& \cdots & ▵f_{n,1,p}\\ ⋮& \ddots & ⋮\\ ▵f_{1,m,p}& \cdots & ▵f_{n,m,p}\end{array}\right]$$

(3.17)

$$▵{\mathit{X}}_p=\left[▵{\mathit{x}}_{k,p}\right]=\left[{\mathit{x}}_{k,p}^B-{\mathit{x}}_p^A\right]=\left[\begin{array}{ccc}▵x_{1,1,p}& \cdots & ▵x_{n,1,p}\\ ⋮& ⋮& ⋮\\ ▵x_{1,n,p}& \cdots & ▵x_{n,n,p}\end{array}\right]$$

(3.18)

The n-dimensional tangent hyperplane ${\mathit{y}}_p\left(\mathit{x}\right)$ through point $A_p\left({\mathit{x}}_p^A{\mathit{f}}_p^A\right)$ may be defined by the “slope”

$${\mathit{S}}_p=\left[{\displaystyle \frac{\partial f_{k,j,p}}{\partial x_{k,i,p}}}\right]=\left[\begin{array}{ccc}{\displaystyle \frac{\partial {\mathit{f}}_{1,1,p}}{\partial {\mathit{x}}_{1,p}}}& \dots & {\displaystyle \frac{\partial {\mathit{f}}_{\mathit{n},1,p}}{\partial {\mathit{x}}_{\mathit{n},\mathit{p}}}}\\ ⋮& ⋮& ⋮\\ {\displaystyle \frac{\partial {\mathit{f}}_{1,m,p}}{\partial {\mathit{x}}_{1,p}}}& \cdots & {\displaystyle \frac{\partial {\mathit{f}}_{\mathit{n},m,p}}{\partial {\mathit{x}}_{\mathit{n},\mathit{p}}}}\end{array}\right]$$

(3.19)

that corresponds to the Jacobian matrix ${\mathit{J}}_p$ of the nonlinear vector-valued function $\mathit{f}\left(\mathit{x}\right)$ and accordingly, Definition 3.16 corresponds to the approximation of the Jacobian matrix. It follows from the condition

$${\mathit{y}}_p\left({\mathit{x}}_{p+1}^A\right)=\mathbf{0}$$

(3.20)

that the zero

$${\mathit{x}}_{p+1}^A={\mathit{x}}_p^A-{\mathit{S}}_p^{+}{\mathit{f}}_p^A$$

(3.21)

of the n-dimensional hyperplane $\mathit{y}\left(\mathit{x}\right)$ approximates the zero ${\mathit{x}}^{*}$ of the nonlinear vector-valued function $\mathit{f}\left(\mathit{x}\right)$, where ${\left[.\right]}^{+}$ stands for the pseudo-inverse. Then the $i^{th}$element of the new approximate ${\mathit{x}}_{p+1}^A$ in the $p^{th}$ iteration will be

$$x_{i,p+1}^A=x_{i,p}^A-{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}{f}_{j,p}^A\right)$$

(3.22)

$\left(i=1,\dots ,n\right),\left(j=1,\dots ,m\right)$.With the iteration stepsize

$$▵{\mathit{x}}_p^A={\mathit{x}}_{p+1}^A-{\mathit{x}}_p^A=-{\mathit{S}}_p^{+}{\mathit{f}}_p^A=-\left[{\left({\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}{f}_{j,p}^A\right)\right)}_i\right]$$

(3.23)

Equation 3.21 can be rewritten as

$${\mathit{S}}_p▵{\mathit{x}}_p^A=-{\mathit{f}}_p^A$$

(3.24)

that is called “secant equation” or “linearization equation” in multi-variable case. With the new approximate ${\mathit{x}}_{p+1}^A$ and with the corresponding function value ${\mathit{f}}_{p+1}^A$, an updated Jacobian or Jacobian approximate ${\mathit{S}}_{p+1}$ has to be determined. Table 1 summarizes the basic equations of the above detailed linearization methods.

If all partial derivatives of the function $\mathit{f}\left(x\right)$ is known, then the Jacobian update can easily be calculated. Newton’s method is one of the most widely used algorithm with very attractive theoretical and practical properties and with some limitations. The computational costs of Newton’s method is high, since the Jacobian ${\mathit{J}}_p$ and the solution to the linear system 3.24 must be computed at each iteration. It is well-known that the local convergence of Newton’s method is q-quadratic if the initial trial approximate ${\mathit{x}}_0$ is close enough to the solution ${\mathit{x}}^{*}$, if $\mathit{J}\left({\mathit{x}}^{*}\right)$ is non-singular and if $\mathit{J}\left(\mathit{x}\right)$ satisfies the Lipschitz condition

$$∥\mathit{J}\left(\mathit{x}\right)-\mathit{J}\left({\mathit{x}}^{*}\right)∥\le L∥\mathit{x}-{\mathit{x}}^{*}∥$$

(3.25)

for all $\mathit{x}$ close enough to ${\mathit{x}}^{*}$.

However, in many cases, the function $\mathit{f}\left(\mathit{x}\right)$ is not an analytical function, the partial derivatives are not known or difficult to evaluate and Newton’s method cannot be applied. Quasi-Newton methods are widely used for solving systems of nonlinear equations when the Jacobian is not known or difficult to determine. The Jacobian is approximated by divided differences 3.16, and the system of nonlinear equations 1.5 is solved by repeatedly solving systems of linear equations 3.24 providing the new approximates ${\mathit{x}}_{p+1}^A$.

The Jacobian approximate 3.16 should be updated according to the fundamental equation of the quasi-Newton methods (“quasi-Newton condition” or “secant condition”)

$${\mathit{S}}_{p+1}▵{\mathit{x}}_p^A={\mathit{f}}_{p+1}^A-{\mathit{f}}_p^A$$

(3.26)

for all $p=0,1,2,\dots$ . However, the condition 3.26 doesn’t uniquely specify the update ${\mathit{S}}_{p+1}$, so the iterative procedure is not well-defined and further constraints are needed. Different methods offer their own specific solution, but new quasi-Newton approximate ${\mathit{x}}_{p+1}$ will never allow full-rank update of the Jacobian approximate ${\mathit{S}}_{p+1}$ (Equation 3.26 is an underdetermined system of m linear equations with $mn$ unknowns). A large family of methods are available with different additional conditions or assumptions. Martinez [21] has been made a thorough survey on the family of practical quasi-Newton methods.

The partial derivatives of the Jacobian matrix may be replaced by suitable difference quotients (discretized Newton iteration, see [12,25])

$$\left[{\displaystyle \frac{▵\mathit{f}}{▵\mathit{x}}}\right]=\left[{\displaystyle \frac{f_j\left(\mathit{x}+▵x_k{\mathit{d}}^k\right)-f_j\left(\mathit{x}\right)}{▵x_k}}\right]=\left[{\displaystyle \frac{▵f_j\left(\mathit{x}\right)}{▵x_k}}\right]$$

(3.27)

$\left(k=1,\dots ,n\right)$, $\left(j=1,\dots ,m\right)$ with n additional function value evaluations, where ${\mathit{d}}^k$ is the $k^{th}$ Cartesian unit vector. However, it is difficult to choose the stepsize $▵\mathit{x}$. If any $▵x_k$ is too large, then Expression 3.27 can be a bad approximation to the Jacobian so the iteration converges much more slowly, if it converges at all. On the other hand, if any $▵x_k$ is too small, then $▵f_j\left(\mathit{x}\right)\simeq 0$, and cancellations can occur which reduces the accuracy of the difference quotients 3.27 (see [35]). Another modification is the inexact-Newton approach, when the nonlinear equation is solved by an iterative linear solver (see [6,10,22]). Wolfe [42] and Popper [29] suggested column update of matrix 3.17.

One of the most widely used Broyden-formula [7] makes a rank-one update for the Jacobian approximate as

$${\mathit{S}}_{p+1}={\mathit{S}}_p+{\mathit{Q}}_p={\mathit{S}}_p+{\displaystyle \frac{▵{\mathit{f}}_p^A-{\mathit{S}}_p▵{\mathit{x}}_p^A}{{∥▵{\mathit{x}}_p^A∥}^2}}{\left[▵{\mathit{x}}_p^A\right]}^T$$

(3.28)

where ${\mathit{Q}}_p$ is a rank-one matrix and by using the Sherman-Morrison formula, the inverse Jacobian approximate update is given as

$${\mathit{S}}_{p+1}^{+}={\mathit{S}}_p^{+}+{\displaystyle \frac{▵{\mathit{x}}_p^A-{\mathit{S}}_p^{+}▵{\mathit{f}}_p^A}{{\left[▵{\mathit{x}}_p^A\right]}^T{\mathit{S}}_p^{+}▵{\mathit{f}}_p^A}}{\left[▵{\mathit{x}}_p^A\right]}^T{\mathit{S}}_p^{+}$$

(3.29)

Broyden’s secant condition 3.26 can be re-written as

$${\mathit{S}}_p▵{\mathit{x}}_p^A+{\mathit{Q}}_p▵{\mathit{x}}_p^A={\mathit{f}}_{p+1}^A-{\mathit{f}}_p^A$$

(3.30)

and with the secant Equation 3.24 we have

$${\mathit{Q}}_p▵{\mathit{x}}_p^A={\mathit{f}}_{p+1}^A$$

(3.31)

As ${\mathit{Q}}_p$ is a rank-one matrix, then this equation has infinitely large number of solutions (the left nullspace of matrix ${\mathit{Q}}_p$ is an $n-1$ dimensional vector space [36]) and the Ortega-Rheinboldt [25] condition will not be satisfied (the column vectors of matrix ${\mathit{S}}_p$ should be linearly independent and have to be “in general position” through the whole iteration process) and it may result an ill-conditioned Broyden update ${\mathit{S}}_{p+1}$.

4. T-Secant Method

A numerical procedure has been developed for solving an overdetermined system of nonlinear equations 1.5. It can be considered both as a discretized Newton’s method or as a quasi-Newton method with full-rank update of the Jacobian approximates. The suggested full-rank update procedure (“T-Secant” method [2]) provides an unconventional and simple strategy to improve the performance of quasi-Newton iterations. It is based on the classic secant linearization, that was completed with a new independent approximate ${\mathit{x}}_{p+1}^B$. The secant equation was modified by a suitably chosen non-uniform scaling transformation and a new independent approximate ${\mathit{x}}_{p+1}^B$ was determined from the modified secant equation. The solution to the secant equation presented in [2] was based on the Wolfe-Popper procedure [29,42] so that the secant equation was splitted into two equations by introducing an auxiliary variable. The simplified algorithm, presented in the paper, directly solves the secant equation by using the pseudo-inverse of the transformed Jacobian approximate.

4.1. Single-Variable Case

The suggested procedure uses the information on the improvement of the classic secant procedure and gives a new approximate $x_{p+1}^B$ in the vicinity of the classic secant approximate $x_{p+1}^A$ as follows. Given two initial approximates $x_p^A$ and $x_p^B$ with function values $f_p^A$ and $f_p^B$ and the new approximate $x_{p+1}^A$ with function value $f_{p+1}^A$ is known from the solution of the classic secant Equation 3.7. An independent new approximate

$$x_{p+1}^B=x_{p+1}^A+▵x_{p+1}$$

(4.1)

is planned to be determined in the vicinity of the classic secant approximate $x_{p+1}^A$. Let the ratio

$$t_p^f={\displaystyle \frac{f_{p+1}^A}{f_p^A}}$$

(4.2)

$\left(f_p^A\ne 0\right)$ of the function value improvement and the ratio

$$t_p^x={\displaystyle \frac{x_{p+1}^B-x_{p+1}^A}{x_{p+1}^A-x_p^A}}={\displaystyle \frac{▵x_{p+1}}{▵x_p^A}}$$

(4.3)

$\left(▵x_p^A\ne 0\right)$ of the desired iteration stepsize change be defined. The single variable secant equation 3.7 is then modified as

$$S_{Tp}▵x_p^A=-f_p^A$$

(4.4)

where

$$S_{Tp}={\displaystyle \frac{t_p^f}{t_p^x}}{S}_p$$

(4.5)

Then the new approximate $x_{p+1}^B$ can be expressed from Equation 4.4

$$t_p^f{\displaystyle \frac{x_{p+1}^A-x_p^A}{x_{p+1}^B-x_{p+1}^A}}{S}_p▵x_p^A=-f_p^A$$

(4.6)

After re-arrangement

$$x_{p+1}^A-x_{p+1}^B={\displaystyle \frac{t_p^f}{f_p^A}}{S}_p{\left(▵x_p^A\right)}^2$$

(4.7)

and with the secant Equation 3.7

$$x_{p+1}^B=x_{p+1}^A-{\displaystyle \frac{t_p^f}{f_p^A}}{S}_p{\left(▵x_p^A\right)}^2=x_{p+1}^A+t_p^f▵x_p^A$$

(4.8)

the desired new iteration stepsize

$$▵x_{p+1}=t_p^f▵x_p^A$$

(4.9)

can be given. The suggested procedure can also be applied to the Newton’s method (T-Newton method). Then the “slope” $S_p$ corresponds to the derivative $f^{\prime }$ of function $f\left(x\right)$.

The geometrical representation of the suggested method can be derived as follows ([2]). Let the new approximate $x_{p+1}^B$ be the zero of a function $z_p\left(x\right)$ :

$$z_p\left(x_{p+1}^B\right)=0$$

(4.10)

Then by replacing $x_{p+1}^B$ with x in Equation 4.8 gives the function

$$z_p\left(x\right)={\displaystyle \frac{t_p^f{S}_p{\left(▵x_p^A\right)}^2}{x-x_{p+1}^A}}+f_p^A=0$$

(4.11)

with zero at

$$x=x_{p+1}^B$$

(4.12)

Equation 4.11 is a hyperbolic function with vertical and horizontal asymptotes $x_{p+1}^A$ and $f_p^A$, and its root $x_{p+1}^B$ will be in the vicinity of $x_{p+1}^A$ in “appropriate distance” that is regulated by the function value $f_{p+1}^A$ (see Figure 2). This virtue of the suggested procedure ensures an automatic mechanism for having the actual approximates $x_p^A$ and $x_p^B$ in general positions through the whole iteration process providing stable and efficient numerical performance. The suggested method’s performance is demonstrated with the results of numerical tests on test functions in Section 6.

4.2. Multi-Variable Case

Let two independent approximates ${\mathit{x}}_p^A=\left[x_{k,p}^A\right]$ and ${\mathit{x}}_p^B=\left[x_{k,p}^B\right]$ to the zero ${\mathit{x}}^{*}$ of the nonlinear vector-valued function $\mathit{f}\left(\mathit{x}\right)$ be given in the $p^{th}$ iteration $\left(k=1,\dots ,n\right)$. Let the series of n approximates ${\mathit{x}}_{k,p}^B$ be constructed by individually increment the elements $x_{k,p}^A$ of the approximate ${\mathit{x}}_p^A$ by an increment

$$▵x_{k,p}=x_{k,p}^B-x_{k,p}^A$$

(4.13)

as

$${\mathit{x}}_{k,p}^B={\mathit{x}}_p^A+▵x_{k,p}{\mathit{d}}^k$$

(4.14)

where ${\mathit{d}}^k$ is the $k^{th}$ Cartesian unit vector. It follows from this special construction of the approximates ${\mathit{x}}_{k,p}^B$ , that $x_{k,i,p}^B-x_{k,p}^A=0$ for $i\ne k$ and $x_{k,i,p}^B-x_{k,p}^A=▵x_{k,p}$ for $i=k$ and matrix 3.18 will be a diagonal matrix :

$$▵{\mathit{X}}_p=\left[▵{\mathit{x}}_{k,p}\right]=\left[{\mathit{x}}_{k,p}^B-{\mathit{x}}_p^A\right]=\left[\begin{array}{ccc}▵x_{1,p}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & ▵x_{n,p}\end{array}\right]=diag\left(▵x_{k,p}\right)$$

(4.15)

Let the ratios

$${\mathit{T}}_p^F=diag\left(t_{j,p}^F\right)=diag\left({\displaystyle \frac{f_{j,p+1}^A}{f_{j,p}^A}}\right)$$

(4.16)

$\left(f_{j,p}^A\ne 0\right)$ of the function value improvements $\left(j=1,\dots ,m\right)$ and the ratios

$${\mathit{T}}_p^X=diag\left(t_{i,p}^X\right)=\mathit{diag}\left({\displaystyle \frac{{\mathit{x}}_{\mathit{i},\mathit{p}+\mathbf{1}}^{\mathit{B}}-{\mathit{x}}_{\mathit{i},\mathit{p}+\mathbf{1}}^{\mathit{A}}}{{\mathit{x}}_{\mathit{i},\mathit{p}+\mathbf{1}}^A-{\mathit{x}}_{\mathit{i},\mathit{p}}^{\mathit{A}}}}\right)=\mathit{diag}\left({\displaystyle \frac{▵{\mathit{x}}_{\mathit{i},\mathit{p}+\mathbf{1}}}{▵{\mathit{x}}_{\mathit{i},\mathit{p}}^{\mathit{A}}}}\right)$$

(4.17)

$\left(▵x_{i,p}^A\ne 0\right)$ of the desired iteration stepsize changes $\left(i=1,\dots ,n\right)$ be defined similarly like in the single variable case. Let the secant equation 3.24 be modified as

$${\mathit{S}}_{Tp}▵{\mathit{x}}_p^A=-{\mathit{f}}_p^A$$

(4.18)

where

$${\mathit{S}}_{Tp}={\mathit{T}}_p^F{\mathit{S}}_p{\left({\mathit{T}}_p^X\right)}^{-1}=diag\left(t_{j,p}^F\right)\left[{\displaystyle \frac{▵f_{k,j,p}}{▵x_{i,p}}}\right]diag\left({\displaystyle \frac{1}{t_{i,p}^X}}\right)=\left[{\displaystyle \frac{t_{j,p}^F}{t_{i,p}^X}}{\displaystyle \frac{▵f_{k,j,p}}{▵x_{i,p}}}\right]$$

(4.19)

$\left(i=1,\dots ,n\right)$, $\left(j=1,\dots ,m\right)$, $\left(k=1,\dots ,n\right)$ or in explicit form :

$${\mathit{S}}_{Tp}=\left[\begin{array}{ccc}{\displaystyle \frac{{\mathit{f}}_{\mathbf{1},\mathit{p}+\mathbf{1}}^{\mathit{A}}}{{\mathit{f}}_{\mathbf{1},\mathit{p}}^{\mathit{A}}}}& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\displaystyle \frac{{\mathit{f}}_{\mathit{m},\mathit{p}+\mathbf{1}}^{\mathit{A}}}{{\mathit{f}}_{\mathit{m},\mathit{p}}^{\mathit{A}}}}\end{array}\right]\left[\begin{array}{ccc}{\displaystyle \frac{▵{\mathit{f}}_{\mathbf{1},\mathbf{1},\mathit{p}}}{▵{\mathit{x}}_{\mathbf{1},\mathit{p}}}}& \dots & {\displaystyle \frac{▵{\mathit{f}}_{\mathbf{1},\mathit{n},\mathit{p}}}{▵{\mathit{x}}_{\mathit{n},\mathit{p}}}}\\ ⋮& ⋮& ⋮\\ {\displaystyle \frac{▵{\mathit{f}}_{1,m,p}}{▵{\mathit{x}}_{\mathbf{1},\mathit{p}}}}& \cdots & {\displaystyle \frac{▵{\mathit{f}}_{\mathit{n},\mathit{m},\mathit{p}}}{▵{\mathit{x}}_{\mathit{n},\mathit{p}}}}\end{array}\right]\left[\begin{array}{ccc}{\displaystyle \frac{▵x_{1,p}^A}{▵x_{1,p+1}}}& 0& 0\\ 0& \ddots & 0\\ 0& 0& {\displaystyle \frac{▵x_{n,p}^A}{▵x_{n,p+1}}}\end{array}\right]$$

(4.20)

and

$${\mathit{S}}_{Tp}^{+}={\mathit{T}}_p^X{\mathit{S}}_p^{+}{\left({\mathit{T}}_p^F\right)}^{-1}$$

(4.21)

is the pseudoinverse of ${\mathit{S}}_p$. Then Equation 4.18 can be re-written as

$$▵{\mathit{x}}_p^A=-{\mathit{S}}_{Tp}^{+}{\mathit{f}}_p^A$$

(4.22)

Then the $i^{th}$ element $x_{i,p+1}^B$ of the new approximate ${\mathit{x}}_{p+1}^B$ in the $p^{th}$ iteration will be

$$x_{i,p+1}^B=x_{i,p+1}^A+▵x_{i,p+1}=x_{i,p+1}^A-{\displaystyle \frac{{\left(▵x_{i,p}^A\right)}^2}{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}\frac{f_{j,p}^A}{t_{j,p}^F}\right)}}$$

(4.23)

where $t_{p,j}^F\ne 0$, $\left(j=1,\dots ,m\right)$ and $\left(i=1,\dots ,n\right)$. Let the ratios

$${\mu }_{i,p}={\displaystyle \frac{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}{f}_{j,p}^A\right)}{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}\frac{f_{j,p}^A}{t_{j,p}^F}\right)}}={\displaystyle \frac{-▵x_{i,p}^A}{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}\frac{f_{j,p}^A}{t_{j,p}^F}\right)}}$$

(4.24)

be introduced . By using Equation 3.23, the new iteration stepsize can be expressed from Equation 4.23 as

$$▵x_{i,p+1}=x_{i,p+1}^B-x_{i,p+1}^A=-{\displaystyle \frac{{\left(▵x_{i,p}^A\right)}^2}{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}\frac{f_{j,p}^A}{t_{j,p}^F}\right)}}={\mu }_{i,p}▵x_{i,p}^A$$

(4.25)

Table 2 summarizes the basic equations of the above detailed multi-variable update method. The suggested procedure can also be applied to the Newton’s method (T-Newton method). Then matrix ${\mathit{S}}_p$ corresponds to the Jacobian matrix ${\mathit{J}}_p$ of function $\mathit{f}\left(\mathit{x}\right)$. The function value improvement parameter ${\mathsf{t}}_{\mathsf{j},\mathsf{p}}^{\mathsf{F}}$ $\left(j=1,\dots ,m\right)$ in Equation 4.24, defined as 4.16 has a key role in the suggested iteration process and the absolute value of the denominator $\left|f_{j,p}^A\right|$ has to be low bounded with a pre-given value $T_{min}$ to avoid division by zero. When the values $\left|f_{j,p}^A\right|$ tend to zero with increasing iteration counter p, then

$${\mu }_{i,p}\to T_{min}$$

(4.26)

Figure 3 shows the effect of iteration parameter $T_{min}$ on the iteration process ($N=1000$). It is clear that too high value of $T_{min}$ has negative effect on the iteration efficiency. It can be seen on the example that a value $0.001$ causes the iteration process failure, and there is almost no effect with lower values $0.0001\leqq T_{min}$.

The suggested procedure may also be repeated without using the classic secant approximates and without updating the previously known Jacobian approximate ${\mathit{S}}_p$ as

$$▵x_{i,p+1+h}=x_{i,p+1+h}^B-x_{i,p+1+h}^A={\displaystyle \frac{-{\left({\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}{f}_{j,p+h}^A\right)\right)}^2}{{\displaystyle \sum _{j=1}^m}\left(S_{i,j,p}^{+}\frac{f_{j,p+h}^A}{t_{j,p+h}^F}\right)}}$$

(4.27)

where

$$t_{j,p+h}^F={\displaystyle \frac{f_{j,p+1+h}^A}{f_{j,p+h}^A}}$$

(4.28)

$\left(h=1,\dots \right)$ until the new approximates

$$x_{i,p+1+h}^B=x_{i,p+1+h}^A+▵x_{i,p+1+h}$$

(4.29)

$\left(i=1,\dots ,,n\right)$ are sufficiently improved. A numerical example is shown in Table 3 with test function

$$f\left(x\right)=x^3-2x-5$$

(4.30)

for demonstration purpose. The results indicate linear convergence with convergence rate $\alpha =1$. $\alpha$ is the computed convergence rate, $N_f$ is the number of function value evaluations and L is the mean convergence rate, suggested by Broyden [7] in Table 3.

5. Algorithm

Let ${\epsilon }^{*}$ be the error bound for termination criterion and if ${\mathit{x}}^{*}$ is known then

$${\mathit{e}}_p^A={\mathit{x}}_p^A-{\mathit{x}}^{*}$$

(5.1)

is the error vector of approximate ${\mathit{x}}_p^A$ in the $p^{th}$ iteration with elements $e_{i,p}^A$ $\left(i=1,\dots ,n\right)$. Let the error norm

$${\epsilon }_p={\displaystyle \frac{{∥{\mathit{e}}_p^A∥}_2}{n}}={\displaystyle \frac{\sqrt{{\displaystyle \sum _{i=1}^n}{\left(e_{i,p}^A\right)}^2}}{n}}$$

(5.2)

be defined, where ${∥.∥}_2$ is Euclidean norm and let the iteration be terminated when

$${\epsilon }_p<{\epsilon }^{*}$$

(5.3)

holds. Choose $T_{min}$ as lower bound for $\left|t_{j,p}^F\right|\left(j=1\dots ,m\right)$ and let $f_{min}$ be a lower bound for $\left|f_{j,p}^A\right|$ $\left(j=1,\dots ,m\right)$. Let $p=0$ and let the approximate ${\mathit{x}}_p^A$ and the difference vector $▵{\mathit{x}}_p$ be given. Calculate the corresponding function values ${\mathit{f}}_p^A$ and assure that $f_{min}<\left|f_{j,0}^A\right|$ $\left(j=1,\dots ,m\right)$. Iteration constants $f_{min}$ and $T_{min}$ are necessary to avoid division by zero and to avoid computed values be near the numerical precision.

• Step 1 : Generate a set of n additional approximates ${\mathit{x}}_{k,p}^B$ (Equation 4.14) and evaluate function values ${\mathit{f}}_{k,p}^B$ $\left(k=1,\dots ,n\right)$. Assure that $f_{min}<\left|f_{k,j,p}^B\right|$ $\left(j=1,\dots ,m\right)$.
• Step 2 (Classic secant update method) : Construct the Jacobian approximate matrix ${\mathit{S}}_p$, determine its pseudo-inverse ${\mathit{S}}_p^{+}$ [31], calculate ${\mathit{x}}_{p+1}^A$ from Equation 3.21 and ${\epsilon }_p$ from Equation 5.2.
• Step 3 : If ${\epsilon }_p<{\epsilon }^{*}$ then terminate iteration, else continue with Step 4.
• Step 4 (suggested update method) : Calculate ${\mathit{f}}_{p+1}^A$ (assure that $f_{min}<\left|f_{j,p}^A\right|$), ${\mathit{T}}_p^F$ from Equation 4.16 and ${\mu }_i$ from Equation 4.24. Let $T_{\min }<\left|t_{j,p}^F\right|$ and determine ${\mathit{x}}_{p+1}^B$ from Equation 4.25.
• Step 5 : Continue iteration from Step 1 with $p=p+1$, ${\mathit{x}}_{p+1}^A$ and ${\mathit{x}}_{p+1}^B$.

If $p_{max}$ is the number of necessary iterations for satisfying the termination criterion ${\epsilon }_p<{\epsilon }^{*}$ and n is the number of unknowns to be determined, then the suggested update method needs $n+1$ function evaluations in each iterations and altogether

$$N_f=p_{max}\left(n+1\right)$$

(5.4)

function evaluations to reach the desired termination criterion. $p_{max}$ is depending on many circumstances such as the nature of the function $f\left(x\right)$, termination criteria (${\epsilon }^{*}$ or others), the distance of the initial approximate ${\mathit{x}}^A$ from the solution ${\mathit{x}}^{*}$ and from the iteration constants $f_{min}$ and $T_{min}$.

6. Numerical Tests Results

6.1. Single Variable Test Function

A numerical example is given with a single-variable test function

$$f\left(x\right)=cos\left(x\right)-x$$

(6.1)

with root $x^{*}\cong 0.73908513...$ The results of the classic secant iteration process are shown on Table 4 and Figure 4. Iterations were made with initial approximates $x_0^A=-2.0$ and $x_0^B=2.0$ $\left(p=0\right)$ providing $f_0^A=1.584$. The first secant approximate $x_1^A=-0.416...$ is found as the zero of the first secant line $y_0\left(x\right)$ providing $f_1^A=1.331...$ $\left(p=1\right)$. The next iteration ($p=2$) will continue with approximate $x_2^A=0.442......$and with $f_2^A=0.462...$ .

The results of the suggested iteration process are shown on Table 5 and Figure 5. Iterations were made with the same initial approximates as in case of the classic secant iteration and the first secant approximate $x_1^A=-0.416...$ is found. The first T-secant approximate $x_1^B=0.915...$ is given as the zero of the first hyperbola function $z_0\left(x\right)$ (Figure 5, left). Iteration then goes on with approximates $x_1^A=-0.416...$ and $x_1^B=0.915...$ providing $f_1^A=1.331...$ $\left(p=1\right)$, and new approximates $x_2^A=0.667...$ and $x_2^B=0.764...$ as the zeros of the second secant and hyperbola functions $y_1\left(x\right)$ and $z_1\left(x\right)$ receptively (Figure 5, left). The next iteration ($p=2$) will then continue with approximates $x_2^A=0.667......$ and $x_2^B=0.764...$ and with $f_2^A=0.119...$ and gives $x_3^A=0.7387......$ and $x_3^B=0.7391...$ (Figure 5, right).

6.2. Solution of an Inverse Problem

An example is given for an imaginary modeling problem. Modeling of systems on the basis of observed or specified system responses plays a key role in scientific research. A central problem of modeling is to fit a structured mathematical model to available data so that the values of the unknown model parameters are determined providing that the simulated data is as near to the available data as possible. A system response can generally be represented by a curve in a two (or three) dimensional space and can be simulated by a computer program. Parameter identification corresponds to minimize the distance between observed and simulated system responses. Parameter identification problems can generally be formulated as non-linear least squares problems so that some unknown model parameters are determined providing minimal deviation between observed and simulated system responses. Solving non-linear least squares problems is often a difficult task especially if the number of unknowns is high.

A rational system under investigation produces measurable response to a known measurable external effect. A mathematical model with n unknown parameters $\left(p_1,\dots ,p_n\right)$ simulates the behavior of the system providing response to the same external effect as acted on the system during observation. Let, the observed and the simulated system responses be represented by two dimensional curves and sampled in k discrete points with coordinates $x_i$ and $y_i$$\left(i=1,\dots k\right)$. The coordinates of a synthetic observed response were generated as

$$x\left(z\right)=\left(a+b\right)cos\left(z+\alpha \right)-ccos\left({\displaystyle \frac{a+b}{b}}z+\alpha \right)+x_0$$

(6.2)

$$y\left(z\right)=\left(a+b\right)sin\left(z+\alpha \right)-csin\left({\displaystyle \frac{a+b}{b}}z+\alpha \right)+y_0$$

(6.3)

(an epicycloid) with parameters $x_0=10$, $y=8$, $a=4$, $b=2$, $c=3.5$ and $\alpha =1$ and with $0\leqq z\leqq 2\pi$. The parameters

$$p=\left[\begin{array}{ccccc}{x}_0& y_0& a& b& c\end{array}\right]$$

(6.4)

were were considered to be unknown, and the system response was simulated with initial approximate

$$p_0=\left[\begin{array}{ccccc}{x}_0& y_0& a& b& c\end{array}\right]=\left[\begin{array}{ccccc}8& 11& 3.5& 2.5& 3\end{array}\right]$$

(6.5)

Distance between observed and simulated system responses was defined for arbitrary two-dimensional curves. This definition can be extended to n dimensions making wide range of parameter identification problems possible to formulate. The value of the defined distance is dimensionless and it expresses the ratio of area between system response graphs in normalized coordinate system to the area of a rectangle with unit area.

The distance between the observed and the simulated system responses was quantified by area $D\left(p\right)$ between the graphs of system responses in a normalized coordinate system. Partial area $r_j\left(p\right)$ was defined as the area of a quadrangle among the ${\left(j-1\right)}^{th}$ and $j^{th}$ division points on the system response graphs $\left(j=1,\dots ,m\right)$ as shown on Figure 6. These division points were selected equidistantly along system response graphs.

Thus, parameter identification problems can be formulated as solving a system of nonlinear equations

$$\mathit{r}\left(\mathit{p}\right)=0$$

(6.6)

where $\mathit{p}\to \mathit{r}\left(\mathit{p}\right)$, $\mathit{p}\in {\mathbb{R}}^5$, $\mathit{r}:{\mathbb{R}}^5\to {\mathbb{R}}^m$ and $\mathit{r}$ is residual vector with components $r_j$$\left(j=1,\dots ,m\right)$. The unknown parameters $\mathit{p}$ of an epicycloid 6.2 and 6.3 were determined according to the suggested update algorithm. The variation of parameters $\mathit{p}$ through iterations and the observed and the simulated system responses for the initial approximate are shown on Figure and Figure 6. The observed and simulated system response pairs are shown on Figure 7 7, 13, 19, 25, 31 and 43 simulations (S).

The results of numerical tests with different initial approximates showed that the optimal solutions were not reached in all cases. The convergence was especially sensitive for the initial values of parameters a, b and c. Several local optimal solutions were detected as shown on Figure 8.

7. Efficiency

The efficiency of an algorithm for the solution of nonlinear equations is thoroughly discussed by Traub [38] as follows. Let $\beta$ be the order of the iteration sequence such that for the approximate errors $e_i=x_i-x^{*}$, there exists a nonzero constant C (asymptotic error constant) for which

$${\displaystyle \frac{\left|e_{i+1}\right|}{{\left|e_i\right|}^{\beta }}}\to C$$

(7.1)

A natural measure of the information used by an algorithm is the “informational usage” d, which is defined as the number of new pieces of information (values of the function and its derivatives) required per iteration (called “horner” by Ostrowski [26]). Then the efficiency of the algorithm within one iteration can be measured by the “informational efficiency”

$$EFF={\displaystyle \frac{\beta }{d}}$$

(7.2)

An alternative definition of efficiency is

$$*EFF={\beta }^{{\displaystyle \frac{1}{d}}}$$

(7.3)

called “efficiency index” by Ostrowski [26]. Another measure of efficiency, called “computational efficiency” takes into account the “cost” of calculating different derivatives. The concept of informational efficiency ($EFF$) and efficiency index ($*EFF$) doesn’t take into account the cost of evaluating f and its derivatives, nor does they take into account the total number of pieces of information needed to achieve a certain accuracy in the root of the function. If f is composed of elementary functions, then the derivatives are also composed of elementary functions, thus the cost of evaluating the derivatives is merely the cost of combining the elementary functions.

Broyden [7] suggested the mean convergence rate

$$L={\displaystyle \frac{1}{N_f}}ln{\displaystyle \frac{R\left({\mathit{x}}_0^A\right)}{R\left({\mathit{x}}_{p_{\max }}^A\right)}}$$

(7.4)

as a measure of efficiency of an algorithm for solving a particular problem, where $N_f$ is the total number of function evaluations, ${\mathit{x}}_0^A$ is the initial approximate, $x_{p_{\max }}^A$ is the last approximate to the solution $x^{*}$ when the termination criteria is satisfied after $p_{max}$ iterations. $R\left(x\right)$ is the Euclidean norm of $f\left(x\right)$.

Zang’s two-step method [43] is similar to King-Werner method [18,41] with convergence rate $2.414\dots$ [32]. Chen’s [9] and Wang’s [39] three-step methods show asymptotic convergence order $2.732...$. Chen’s [9] method require two function and one derivative evaluations by using a second-order polynomial, proposed in [39] and shows $1.174\dots$ Ostrowski’s [26] efficiency index. Wang’s [39] method require three function and two derivative evaluations shows $1.101\dots$ Ostrowski-index. Zhang [43] published a finite difference based secant method with $2.618\dots$ asymptotic convergence order. Later Ren [32] showed that Zhang’s [43] method has only $2.414\dots$ convergence order due to some mistakes in the derivation. Table 6 compares the efficiencies of classic (Rows 1-2 : Secant, Newton) and improved algorithms (Rows 3, 5 : T-Secant, T-Newton with the suggested full-rank update and from references [9,39]), the T-Secant method with constant Jacobian approximate ${\mathit{S}}_p$ (row 4 : TS-const. ${\mathit{S}}_p$, see data in Table 3), Chen’s [9] and Wang’s [39] methods.

Very limited data are available to compare the performance of the suggested update method (T-Secant) with other classic methods, especially for large number of unknowns. Efficiency results were given by Broyden [7] for the Rosenbrock function for $N=2$. The calculated convergence rates for the two Broyden method variants [7], for the Powell’s method [30], for the adaptive coordinate descent method [19] and for the Nelder-Mead simplex method [20] were compared with the calculated values for the T-secant method in Table 7 [2]. Rows 1-5 are data from referenced papers, rows 6-8 are T-secant results with the referenced initial approximates and rows 9-15 are calculated data for $N>2$. If the value of $R\left(\mathit{x}\right)=0$ (for $p_{max}$) is zero, then the mean convergence rates (L and $L_N$) are not countable (zero in the denominator). A substitute value ${10}^{-25}$ was used when $R\left(\mathit{x}\right)=0$ (for $p_{max}$) in rows 6, 7, 10 and 13.

Results show that the mean convergence rate L (Equation 7.4) for $N=2$ is much higher for the T-secant method ($\simeq 5.5-6.9$) than for the other listed methods ( $\simeq 0.1-0.6$). However it is obvious that the mean convergence rate values decrease rapidly with increasing N values (more unknowns need more function evaluations). A modified convergence rate

$$L_N=N*L={\displaystyle \frac{N}{N_f}}ln{\displaystyle \frac{R\left({\mathit{x}}_0^A\right)}{R\left({\mathit{x}}_{p_{\max }}^A\right)}}$$

(7.5)

is suggested to use as an $"N"$ independent measure of efficiency (see Table 7). The values of L and $L_N$ are at least 10 times larger for the T-secant method than for the referenced classic methods for $N=2$ (see Table 7).

The efficiency measures (L and $L_N$) are also depending on the initial conditions (distance of the initial approximate from the optimal solution, termination criterion). Results from large number of numerical tests indicate an average $L_N\simeq 7.4$ with standard deviation around $3.7$ for the T-secant method even for large N values.

8. Conclusions

A numerical procedure has been developed for solving an overdetermined system of nonlinear equations 1.5. It can be considered both as a discretized Newton’s method or as a quasi-Newton method with full-rank update of the Jacobian approximates. Quasi-Newton methods are widely used for solving systems of nonlinear equations when the function derivatives (Jacobian) are not known or difficult to determine. The derivatives (Jacobian) are approximated by divided differences 3.16, and the system of nonlinear equations 1.5 is solved by repeatedly solving systems of linear equations 3.24. The new approximate ${\mathit{x}}_{p+1}^A$ is determined from the classic secant equation and the divided differences are updated according to the secant condition 3.26. However, the secant condition doesn’t uniquely specify the Jacobian approximate, so the update procedure is not well-defined and further constraints are needed. Different methods offer specific update solutions.

The suggested numerical procedure (“T-Secant” method [2]) provides an unconventional and simple strategy to improve the performance of quasi-Newton iterations by full-rank update of the Jacobian approximates. It is based on the classic secant linearization and allows full-rank update of the Jacobian approximates. The classic secant iteration was completed with a new independent approximate ${\mathit{x}}_{p+1}^B$, that was determined from a modified secant equation 4.18. Modification was made by a suitably chosen non-uniform scaling transformation. Scaling was made by the difference quotients of the classic secant function value improvements 4.16 and the difference quotients of the desired new iteration stepsizes and the classic secant stepsizes 4.17. A new independent approximate ${\mathit{x}}_{p+1}^B$ was determined from the modified secant equation 4.18. The solution to the secant equation presented in [2] was based on the Wolfe-Popper procedure [29,42] so that the secant equation was splitted into two equations by introducing an auxiliary variable. The simplified algorithm, presented in the paper, directly solves the secant equation by using the pseudo-inverse 4.21 of the transformed Jacobian approximate 4.19.

It has been shown that the new T-secant approximate ${\mathit{x}}_{p+1}^B$ will be in the vicinity of the classic secant approximate ${\mathit{x}}_{p+1}^A$ if the classic secant iterates converge to the root of the nonlinear function [2]. The Jacobian approximate was then full-rank updated by constructing new divided differences (Jacobian approximates) from the classic secant and from the T-secant approximates.

It was shown that the iterative procedure possesses super-quadratic asymptotic convergence property with convergence rate $\alpha =\phi +1={\phi }^2\simeq \mathbf{2}.\mathbf{618}$ for simple root in single variable case ([2]), where $\phi$ is the golden section ratio. The suggested procedure can also be applied to the Newton’s method (matrix ${\mathit{S}}_p$ corresponds to the Jacobian matrix ${\mathit{J}}_p$). Numerical test results indicate that the quadratic convergence ($\alpha =\mathbf{2}$) of the Newton’s method increases to cubic ($\alpha =\mathbf{3}$).

The efficiency has been studied in multi-variable case and compared with other classic rank-one update and line-search methods on the basis of available data. Results show that its efficiency is considerably better than the efficiency of other classic low-rank update methods. The suggested method’s performance was demonstrated by the results of numerical tests with widely used single- and multi-variable benchmark test functions. A Rosenbrock test function was used with up to 1000 variables [2]. The method has also been successfully applied for the solution of different “real life” inverse problems [3,4,5] on physically nonlinear dynamic systems. Further studies can be made on convergence properties in case of multiple roots and on the efficiency characteristics in cases of more benchmark test functions. The suggested procedure may be used for other applications especially with large number of unknowns and when low number of function evaluations is crucial.