# Optimal Algorithms and the BFGS Updating Techniques for Solving Unconstrained Nonlinear Minimization Problems

Authors: Chein-Shan Liu
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

To solve an unconstrained nonlinear minimization problem, we propose an optimal algorithm (OA) as well as a globally optimal algorithm (GOA), by deflecting the gradient direction to the best descent direction at each iteration step, and with an optimal parameter being derived explicitly. An invariant manifold defined for the model problem in terms of a locally quadratic function is used to derive a purely iterative algorithm and the convergence is proven. Then, the rank-two updating techniques of BFGS are employed, which result in several novel algorithms as being faster than the steepest descent method (SDM) and the variable metric method (DFP). Six numerical examples are examined and compared with exact solutions, revealing that the new algorithms of OA, GOA, and the updated ones have superior computational efficiency and accuracy.


---

## Body

## 1. Introduction

The steepest descent method (SDM), which can be traced back to Cauchy (1847), is the simplest gradient method for solving unconstrained minimization problems. However, the SDM performs well during earlier stages and as approaching to a stationary point it converges very slowly. In this paper, we consider the following nonlinear minimization problem without considering constraint:(1)
min
⁡
f
(
x
)
,
where f
:
ℝ
n
↦
ℝ is a ℂ
2 differentiable function.In the iterative solution of (1), if x
k is the current iterative point, then we denote f
(
x
k
) by f
k, ∇
f
(
x
k
) by g
k, and ∇
2
f
(
x
k
) by A
k, which is known to be a symmetric Hessian matrix. The second-order Taylor expansion of function f
(
x
) at the point x
k is
(2)
f
(
x
)
=
f
k
+
g
k
T
Δ
x
+
1
2
(
Δ
x
)
T
A
k
Δ
x
,
where Δ
x
=
x
-
x
k. The superscript T signifies the transpose and meanwhile g
k
T
Δ
x signifies the inner product of g
k and Δ
x.Letx
=
x
k
-
λ
0
g
k, and inserting it into (2) we can obtain
(3)
f
(
x
k
-
λ
0
g
k
)
=
f
k
-
λ
0
g
k
T
g
k
+
λ
0
2
2
g
k
T
A
k
g
k
.
By requiring the minimization with respect to λ
0, we can derive
(4)
λ
0
=
∥
g
k
∥
2
g
k
T
A
k
g
k
.
Then, we have a famous steepest descent method (SDM) for solving (1).(i)
Give an initialx
0, and then compute g
0
=
∇
f
(
x
0
).
(ii)
Fork
=
0,1
,
2
,
…, we repeat the following iteration:
(5)
x
k
+
1
=
x
k
-
∥
g
k
∥
2
g
k
T
A
k
g
k
g
k
. If ∥
g
k
+
1
∥
<
ε, then stop; otherwise, go to step (ii).For the minimization problem (1) we need to solve ∇
f
=
0, and hence the residual means the value of ∥
∇
f
∥
=
∥
g
∥. The above convergence criterion ∥
g
k
+
1
∥
<
ε means that when the residual norm is smaller than a given error tolerance ε, the iterations are terminated.In the derivation of SDM for solving (1) it is easy to see that we have transformed the global minimization problem into a model problem in terms of a quadratic minimization problem of
(6)
ϕ
(
x
)
=
1
2
x
T
A
x
-
b
T
x
+
c
0
and determined the coefficient λ
0 by (4), where c
0
=
f
k
-
g
k
T
x
k
+
x
k
T
A
k
x
k
/
2
+
k
0 with k
0 being a constant to raise the level value of ϕ, and b
=
A
k
x
k
-
g
k is a constant vector within each iterative step. Here for the purpose of simple notation we omit the subscript k in (6), and we are going to modify the SDM by starting from the above locally quadratic function.Several modifications of the SDM have been addressed. These modifications have led to a new interest in the SDM that the gradient vector itself is not a bad choice but rather that the original steplengthλ
0 leads to a slow convergence behavior. Barzilai and Borwein [1] have presented a new choice of steplength through a two-point stepsize. Although their method did not guarantee the descent of the minimum function values, Barzilai and Borwein [1] were able to produce a substantial improvement of the convergence speed for a certain test of a quadratic function. The results of Barzilai and Borwein [1] have spurred many researches on the SDM, for example, Raydan [2, 3], Friedlander et al. [4], Raydan and Svaiter [5], Dai et al. [6], Dai and Liao [7], Dai and Yuan [8], Fletcher [9], and Yuan [10]. In this paper, we will approach this problem from a quite different viewpoint of invariant manifold and propose a new strategy to modify the steplength and the descent direction. Besides the SDM, there were many modifications of the conjugate gradient method for the unconstrained minimization problems, like Birgin and Martinez [11], Andrei [12–14], Zhang [15], Babaie-Kafaki et al. [16], and Shi and Guo [17].Also, there is another class method with the descent directiond in f
(
x
k
-
λ
d
) being taken to be D
∇
f
(
x
k
), where D is a positive definite matrix that approximates the inverse of the Hessian matrix A, which is usually named the quasi-Newton method. The earlier method of minimization of a nonlinear function by using this type approach is performed by Davidon [18], which was then simplified and reformulated by Fletcher and Powell [19] and was referred to as the variable metric method (DFP).The remaining portions of this paper are arranged as follows. In Section2 we describe an invariant manifold to derive the governing ordinary differential equations (ODEs). The main results are derived in Section 3, which includes the proof of convergence theorem, optimal parameter, optimal algorithm, a critical parameter, and a globally optimal algorithm. Then, in Section 4 we employ the rank-two Broyden-Fletcher-Goldfarb-Shanno (BFGS) updating techniques to update the Hessian matrix or its inversion, resulting in several novel optimal algorithms. The numerical examples are tested in Section 5 to assess the performance of the newly proposed algorithms. Finally, the conclusions are drawn in Section 6.


## 2. An Invariant Manifold

From this section on, we focus on the local minimum problem defined in terms ofϕ in (6), rather than that of f. When the novel algorithms are developed, we will return to the minimization problem (1). The present approach is different from the conventional line search method by minimizing the steplength λ
k in
(7)
f
(
x
k
-
λ
k
d
k
)
=
min
⁡
λ
>
0
f
(
x
k
-
λ
d
k
)
,
where d
k is a given search direction.At the first, we consider an iterative scheme ofx derived from the ordinary differential equations (ODEs) defined on an invariant manifold which is formed from ϕ
(
x
):
(8)
h
(
x
,
t
)
:
=
Q
(
t
)
ϕ
(
x
)
=
C
.
Here, we let x be a function of a fictitious time variable t. We do not need to specify the function Q
(
t
) a priori, of which C
/
Q
(
t
) is merely a measure of the decreasing of ϕ in time. Hence, we expect that in our algorithm if Q
(
t
)
>
0 is an increasing function of t, the iterative point x
k can tend to the minimal point. We let Q
(
0
)
=
1, and C is determined from the initial condition x
(
0
)
=
x
0 by
(9)
C
=
ϕ
(
x
0
)
>
0
.
We can suitably choose the constant k
0 and hence c
0 in (6), such that ϕ
(
x
)
>
0. Indeed, the different level of ϕ does not alter its minimal point.WhenC
>
0 and Q
>
0, the manifold defined by (8) is continuous, and thus the following differential operation being carried out on the manifold makes sense. For the requirement of consistency condition we have
(10)
Q
˙
(
t
)
ϕ
(
x
)
+
Q
(
t
)
(
A
x
-
b
)
·
x
˙
=
0
,
which is obtained by taking the time differential of (8) with respect to t, using (6), and considering x
=
x
(
t
). We suppose that x is governed by the following ODEs:
(11)
x
˙
=
-
κ
u
,
where κ is to be determined. Inserting (11) into (10) we can solve
(12)
κ
=
q
(
t
)
ϕ
g
·
u
,
where
(13)
g
:
=
A
x
-
b
,
(14)
q
(
t
)
:
=
Q
˙
(
t
)
Q
(
t
)
.We further suppose that(15)
u
=
u
1
+
α
u
2
=
:
g
+
α
B
g
,
where α is a parameter to be determined below through the solution of an optimal equation derived, and B is a descent matrix to be specified. Here, we assert that the driving vector u is an optimal linear combination of the gradient vector and a supplemental vector being the gradient vector times B.


## 3.  Numerical Methods

### 3.1. Convergence Theorem

Before the derivation of optimal algorithms we can prove the following convergence result.Theorem 1.
For an iterative scheme to solve (1) by
(16)
x
(
t
+
Δ
t
)
=
x
(
t
)
-
(
1
-
γ
)
g
·
u
u
T
A
u
u
,
which is generated from the ODEs in (11), the iterative point x on the manifold (8) has the following convergence rate:
(17)
Convergence
Rate
:
=
Q
(
t
+
Δ
t
)
Q
(
t
)
=
1
s
>
1
,
where
(18)
0
<
s
=
1
-
1
-
γ
2
2
a
0
<
1
,
(19)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
≥
1
2
,
and 0
≤
γ
<
1 is a relaxation parameter.Proof.
The proof of this theorem is quite lengthy and we divide it into three parts.
(A) Inserting (12) into (11) we can obtain an evolution equation for x:
(20)
x
˙
=
-
q
(
t
)
ϕ
g
·
u
u
.
In the algorithm if Q
(
t
) can be guaranteed to be an increasing function of t, we might have an absolutely convergent property in finding the minimum of ϕ through the following equation:
(21)
ϕ
(
t
)
=
C
Q
(
t
)
,
which is obtained from (8). Here we simplify the notation of ϕ
(
x
(
t
)
) to ϕ
(
t
).
(B) By applying the Euler method to (20) we can obtain the following algorithm:
(22)
x
(
t
+
Δ
t
)
=
x
(
t
)
-
β
ϕ
g
·
u
u
,
where
(23)
β
=
q
(
t
)
Δ
t
.
In order to keepx on the manifold defined by (21) we can insert the above x
(
t
+
Δ
t
) into
(24)
1
2
x
T
(
t
+
Δ
t
)
A
x
(
t
+
Δ
t
)
-
b
T
x
(
t
+
Δ
t
)
=
C
Q
(
t
+
Δ
t
)
-
c
0
,
and obtain
(25)
C
Q
(
t
+
Δ
t
)
-
c
0
=
1
2
x
T
(
t
)
A
x
(
t
)
-
b
T
x
(
t
)
+
β
ϕ
[
b
-
A
x
(
t
)
]
T
u
g
·
u
+
β
2
ϕ
2
u
T
A
u
2
(
g
·
u
)
2
.
Thus, by (13), (21), and (6) and through some manipulations we can derive the following scalar equation:
(26)
1
2
a
0
β
2
-
β
+
1
=
Q
(
t
)
Q
(
t
+
Δ
t
)
,
where
(27)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
≥
1
2
,
of which the inequality can be achieved by taking a suitable value of k
0 and hence c
0 in (6).
(C) Let(28)
s
:
=
Q
(
t
)
Q
(
t
+
Δ
t
)
,
and by (26) we can derive
(29)
1
2
a
0
β
2
-
β
+
1
-
s
=
0
.
From (29), we can take the solution of β to be
(30)
β
=
1
-
1
-
2
(
1
-
s
)
a
0
a
0
,
if
1
-
2
(
1
-
s
)
a
0
≥
0
.
Let
(31)
1
-
2
(
1
-
s
)
a
0
=
γ
2
≥
0
,
(32)
s
=
1
-
1
-
γ
2
2
a
0
,
and the sufficient condition 1
-
2
(
1
-
s
)
a
0
≥
0 in (30) is satisfied automatically, and thus by (30) and (31) we can obtain a preferred solution of β by
(33)
β
=
1
-
γ
a
0
.
Here 0
≤
γ
<
1 is a relaxation parameter. The inequality 0
<
s
<
1 follows from (32), 1
-
γ
2
>
0, and 2
a
0
≥
1. Inserting the above β into (22) and using (19) we can derive algorithm (16). For this algorithm we can define the local convergence rate by
(34)
Convergence
Rate
:
=
ϕ
(
t
)
ϕ
(
t
+
Δ
t
)
,
which, using (28) and (21) and 0
<
s
<
1 just proved, renders (17). This ends the proof of Theorem 1.


### 3.2. Optimization ofα

In Algorithm (22) we do not yet specify how to choose the parameter α. We can determine a suitable value of α such that s defined in (32) is minimized with respect to α, because a smaller s will lead to a faster convergence as shown by (17).Thus by inserting (19) for a
0 into (32) we can write s to be
(35)
s
=
1
-
(
1
-
γ
2
)
(
g
·
u
)
2
2
ϕ
u
·
(
A
u
)
,
where u as defined by (15) includes a parameter α. Let ∂
s
/
∂
α
=
0, and through some algebraic operations we can solve α and denote it by
(36)
α
o
=
g
·
u
1
u
1
·
(
A
u
2
)
-
g
·
u
2
u
1
·
(
A
u
1
)
g
·
u
2
u
1
·
(
A
u
2
)
-
g
·
u
1
u
2
·
(
A
u
2
)
,
where the subscript o signifies that α
o is the optimal value of α.Remark 2.
For the usual three-dimensional vectorsa, b, c
∈
ℝ
3, the following formula is famous:
(37)
a
×
(
b
×
c
)
=
(
a
·
c
)
b
-
(
a
·
b
)
c
.
Liu [20] has developed a Jordan algebra by extending the above formula to vectors in n-dimension:
(38)
[
a
,
b
,
c
]
=
(
a
·
b
)
c
-
(
c
·
b
)
a
,
a
,
b
,
c
∈
ℝ
n
.
In terms of the Jordan algebra we can write(39)
α
o
=
[
u
1
,
g
,
u
2
]
·
(
A
u
1
)
[
u
2
,
g
,
u
1
]
·
(
A
u
2
)
,
where the symmetry of A was used. It can be seen that the above equation is a more symmetric form than that in (36).


### 3.3. An Optimal Algorithm

Now we can letx
k denote the numerical value of x at the kth step and go back g to g
k, u to u
k, A to A
k, and B to B
k. Thus, by using (16) we can derive an iterative algorithm:
(40)
x
k
+
1
=
x
k
-
η
g
k
·
u
k
u
k
T
A
k
u
k
u
k
,
where
(41)
η
=
1
-
γ
.Therefore, we have the following optimal algorithm (OA).(i)
Select0
≤
γ
<
1, and give an initial x
0.
(ii)
Fork
=
0,1
,
2
,
…, we repeat the following iterations:
(42)
α
k
=
[
g
k
,
g
k
,
B
k
g
k
]
·
(
A
k
g
k
)
[
B
k
g
k
,
g
k
,
g
k
]
·
(
A
k
B
k
g
k
)
,
u
k
=
g
k
+
α
k
B
k
g
k
,
x
k
+
1
=
x
k
-
(
1
-
γ
)
g
k
·
u
k
u
k
T
A
k
u
k
u
k
.If ∥
g
k
+
1
∥
<
ε, then stop; otherwise, go to step (ii).Again we emphasize that we need to solve∇
f
=
0 for the minimization problem (1), and hence the residual means the value of ∥
∇
f
∥
=
∥
g
∥. The above convergence criterion ∥
g
k
+
1
∥
<
ε means that when the residual norm is smaller than a given error tolerance ε, the iterations are terminated.


### 3.4. A Critical Value forα

In Sections3.2 and 3.3 we have used ∂
s
/
∂
α
=
0 (or equivalently, ∂
a
0
/
∂
α
=
0) to find the optimal value of α in the descent vector u
=
g
+
α
B
g. Usually, this value of α obtained from ∂
s
/
∂
α
=
0 is not the global minimum of a
0 (or s). Here, we try another approach and attempt to derive a better value of α than α
o, such that the value of α obtained in this manner is the global minimum of a
0 (or s).In practice, we can take(43)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
=
a
s
.
When a
s is near to 0.5, the convergence speed is very fast. Inserting (15) for u into the above equation and through some elementary operations we can derive a quadratic equation to solve α:
(44)
e
1
α
2
+
e
2
α
+
e
3
=
0
,
where
(45)
e
1
:
=
ϕ
u
2
T
A
u
2
-
a
s
(
g
·
u
2
)
2
,
(46)
e
2
:
=
2
ϕ
u
1
T
A
u
2
-
2
a
s
g
·
u
1
g
·
u
2
,
(47)
e
3
:
=
ϕ
u
1
T
A
u
1
-
a
s
(
g
·
u
1
)
2
.If the following condition is satisfied:(48)
D
:
=
e
2
2
-
4
e
1
e
3
≥
0
,
then α in (44) has areal solution:
(49)
α
=
D
-
e
2
2
e
1
.Inserting (45)–(47) into thecritical equation:
(50)
D
=
e
2
2
-
4
e
1
e
3
=
0
,
we can derive an algebraic equation to determine that a
s is the lowest bound of (48).In this lowest bound
a
s
is a critical value denoted by
a
c, and for all a
s
≥
a
c it can satisfy (48) automatically. From (50) through some elementary operations, the critical value a
c can be solved as
(51)
a
c
=
ϕ
[
u
1
T
A
u
1
u
2
T
A
u
2
-
(
u
1
T
A
u
2
)
2
]
u
1
T
A
u
1
(
g
·
u
2
)
2
+
u
2
T
A
u
2
(
g
·
u
1
)
2
-
2
u
1
T
A
u
2
g
·
u
1
g
·
u
2
.
Then by inserting it for a
s into (45) and (46) we can obtain acritical value  α
c for α from (49):
(52)
α
c
=
a
c
g
·
u
1
g
·
u
2
-
ϕ
u
1
T
A
u
2
ϕ
u
2
T
A
u
2
-
a
c
(
g
·
u
2
)
2
,
where D
=
0 was used in view of (50).By inserting (51) for a
c into (52) and cancelling out the common term ϕ we can derive a final equation for α
c, of which we use the same symbols for saving notations:
(53)
a
c
=
u
1
T
A
u
1
u
2
T
A
u
2
-
(
u
1
T
A
u
2
)
2
u
1
T
A
u
1
(
g
·
u
2
)
2
+
u
2
T
A
u
2
(
g
·
u
1
)
2
-
2
u
1
T
A
u
2
g
·
u
1
g
·
u
2
,
(54)
α
c
=
a
c
g
·
u
1
g
·
u
2
-
u
1
T
A
u
2
u
2
T
A
u
2
-
a
c
(
g
·
u
2
)
2
.Here we must emphasize that, in the current descent vectoru
=
g
+
α
B
g, the above value α
c is thebest one, and the vector
(55)
u
=
g
+
α
c
B
g
(
best
descent
vector
)
is thebest descent vector. Due to its criticality, if one attempts to find a better value of the parameter
α
than
α
c
, there would be no real solution of
α. Furthermore, the best descent vector is also better than the optimal vector u
=
g
+
α
o
B
g derived in Section 3.2.


### 3.5. A Globally Optimal Algorithm

Then, we can derive the followingglobally optimal algorithm (GOA) to solve the minimization problem in (1).(i)
Select0
≤
γ
<
1 and give an initial guess of x
0.
(ii)
Fork
=
0,1
,
2
,
…, we repeat the following iterations:
(56)
u
1
k
=
g
k
,
u
2
k
=
B
k
g
k
,
a
c
k
=
(
u
1
k
·
A
k
u
1
k
u
2
k
·
A
k
u
2
k
-
(
u
1
k
·
A
k
u
2
k
)
2
)
×
(
u
1
k
·
A
k
u
1
k
(
g
k
·
u
2
k
)
2
+
u
2
k
·
A
k
u
2
k
(
g
k
·
u
1
k
)
2
-
2
u
1
k
·
A
k
u
2
k
g
k
·
u
1
k
g
k
·
u
2
k
)
-
1
,
α
k
=
a
c
k
g
k
·
u
1
k
g
k
·
u
2
k
-
u
1
k
·
A
k
u
2
k
u
2
k
·
A
k
u
2
k
-
a
c
k
(
g
k
·
u
2
k
)
2
,
(
critical
optimal
parameter
)
,
u
k
=
g
k
+
α
k
B
k
g
k
,
(
best
descent
vector
)
,
(57)
x
k
+
1
=
x
k
-
(
1
-
γ
)
g
k
·
u
k
u
k
A
k
u
k
u
k
.If x
k
+
1 converges according to a given stopping criterion ∥
g
k
+
1
∥
<
ε, then stop; otherwise, go to step (ii).Remark 3.
We have derived anovel globally optimal algorithm for solving the minimization problem in (1). In terms of the descent vector u
=
g
+
α
c
B
g, the GOA is the best one, which leads to theglobal minimum of a
0 (or s) and hence thelargest convergence rate. While the parameter γ is chosen by the user with problem dependence, the parameter α
k is exactly given by (56). Up to here we have successfully derived a novel best descent vector algorithm, with the help from (19), (53), and (54).Remark 4.
At the very beginning we have set an invariant manifoldh
(
x
,
t
)
=
Q
(
t
)
ϕ
(
x
) in (8) as our starting point to derive the iterative optimal algorithms, which includes a locally objective function ϕ in the governing equations. However, in the final stage the terms which include ϕ can be cancelled out, and thus we have obtained the optimal algorithm in (42) and the globally optimal algorithm in (57), which are both independent of ϕ.


## 3.1. Convergence Theorem

Before the derivation of optimal algorithms we can prove the following convergence result.Theorem 1.
For an iterative scheme to solve (1) by
(16)
x
(
t
+
Δ
t
)
=
x
(
t
)
-
(
1
-
γ
)
g
·
u
u
T
A
u
u
,
which is generated from the ODEs in (11), the iterative point x on the manifold (8) has the following convergence rate:
(17)
Convergence
Rate
:
=
Q
(
t
+
Δ
t
)
Q
(
t
)
=
1
s
>
1
,
where
(18)
0
<
s
=
1
-
1
-
γ
2
2
a
0
<
1
,
(19)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
≥
1
2
,
and 0
≤
γ
<
1 is a relaxation parameter.Proof.
The proof of this theorem is quite lengthy and we divide it into three parts.
(A) Inserting (12) into (11) we can obtain an evolution equation for x:
(20)
x
˙
=
-
q
(
t
)
ϕ
g
·
u
u
.
In the algorithm if Q
(
t
) can be guaranteed to be an increasing function of t, we might have an absolutely convergent property in finding the minimum of ϕ through the following equation:
(21)
ϕ
(
t
)
=
C
Q
(
t
)
,
which is obtained from (8). Here we simplify the notation of ϕ
(
x
(
t
)
) to ϕ
(
t
).
(B) By applying the Euler method to (20) we can obtain the following algorithm:
(22)
x
(
t
+
Δ
t
)
=
x
(
t
)
-
β
ϕ
g
·
u
u
,
where
(23)
β
=
q
(
t
)
Δ
t
.
In order to keepx on the manifold defined by (21) we can insert the above x
(
t
+
Δ
t
) into
(24)
1
2
x
T
(
t
+
Δ
t
)
A
x
(
t
+
Δ
t
)
-
b
T
x
(
t
+
Δ
t
)
=
C
Q
(
t
+
Δ
t
)
-
c
0
,
and obtain
(25)
C
Q
(
t
+
Δ
t
)
-
c
0
=
1
2
x
T
(
t
)
A
x
(
t
)
-
b
T
x
(
t
)
+
β
ϕ
[
b
-
A
x
(
t
)
]
T
u
g
·
u
+
β
2
ϕ
2
u
T
A
u
2
(
g
·
u
)
2
.
Thus, by (13), (21), and (6) and through some manipulations we can derive the following scalar equation:
(26)
1
2
a
0
β
2
-
β
+
1
=
Q
(
t
)
Q
(
t
+
Δ
t
)
,
where
(27)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
≥
1
2
,
of which the inequality can be achieved by taking a suitable value of k
0 and hence c
0 in (6).
(C) Let(28)
s
:
=
Q
(
t
)
Q
(
t
+
Δ
t
)
,
and by (26) we can derive
(29)
1
2
a
0
β
2
-
β
+
1
-
s
=
0
.
From (29), we can take the solution of β to be
(30)
β
=
1
-
1
-
2
(
1
-
s
)
a
0
a
0
,
if
1
-
2
(
1
-
s
)
a
0
≥
0
.
Let
(31)
1
-
2
(
1
-
s
)
a
0
=
γ
2
≥
0
,
(32)
s
=
1
-
1
-
γ
2
2
a
0
,
and the sufficient condition 1
-
2
(
1
-
s
)
a
0
≥
0 in (30) is satisfied automatically, and thus by (30) and (31) we can obtain a preferred solution of β by
(33)
β
=
1
-
γ
a
0
.
Here 0
≤
γ
<
1 is a relaxation parameter. The inequality 0
<
s
<
1 follows from (32), 1
-
γ
2
>
0, and 2
a
0
≥
1. Inserting the above β into (22) and using (19) we can derive algorithm (16). For this algorithm we can define the local convergence rate by
(34)
Convergence
Rate
:
=
ϕ
(
t
)
ϕ
(
t
+
Δ
t
)
,
which, using (28) and (21) and 0
<
s
<
1 just proved, renders (17). This ends the proof of Theorem 1.


## 3.2. Optimization ofα

In Algorithm (22) we do not yet specify how to choose the parameter α. We can determine a suitable value of α such that s defined in (32) is minimized with respect to α, because a smaller s will lead to a faster convergence as shown by (17).Thus by inserting (19) for a
0 into (32) we can write s to be
(35)
s
=
1
-
(
1
-
γ
2
)
(
g
·
u
)
2
2
ϕ
u
·
(
A
u
)
,
where u as defined by (15) includes a parameter α. Let ∂
s
/
∂
α
=
0, and through some algebraic operations we can solve α and denote it by
(36)
α
o
=
g
·
u
1
u
1
·
(
A
u
2
)
-
g
·
u
2
u
1
·
(
A
u
1
)
g
·
u
2
u
1
·
(
A
u
2
)
-
g
·
u
1
u
2
·
(
A
u
2
)
,
where the subscript o signifies that α
o is the optimal value of α.Remark 2.
For the usual three-dimensional vectorsa, b, c
∈
ℝ
3, the following formula is famous:
(37)
a
×
(
b
×
c
)
=
(
a
·
c
)
b
-
(
a
·
b
)
c
.
Liu [20] has developed a Jordan algebra by extending the above formula to vectors in n-dimension:
(38)
[
a
,
b
,
c
]
=
(
a
·
b
)
c
-
(
c
·
b
)
a
,
a
,
b
,
c
∈
ℝ
n
.
In terms of the Jordan algebra we can write(39)
α
o
=
[
u
1
,
g
,
u
2
]
·
(
A
u
1
)
[
u
2
,
g
,
u
1
]
·
(
A
u
2
)
,
where the symmetry of A was used. It can be seen that the above equation is a more symmetric form than that in (36).


## 3.3. An Optimal Algorithm

Now we can letx
k denote the numerical value of x at the kth step and go back g to g
k, u to u
k, A to A
k, and B to B
k. Thus, by using (16) we can derive an iterative algorithm:
(40)
x
k
+
1
=
x
k
-
η
g
k
·
u
k
u
k
T
A
k
u
k
u
k
,
where
(41)
η
=
1
-
γ
.Therefore, we have the following optimal algorithm (OA).(i)
Select0
≤
γ
<
1, and give an initial x
0.
(ii)
Fork
=
0,1
,
2
,
…, we repeat the following iterations:
(42)
α
k
=
[
g
k
,
g
k
,
B
k
g
k
]
·
(
A
k
g
k
)
[
B
k
g
k
,
g
k
,
g
k
]
·
(
A
k
B
k
g
k
)
,
u
k
=
g
k
+
α
k
B
k
g
k
,
x
k
+
1
=
x
k
-
(
1
-
γ
)
g
k
·
u
k
u
k
T
A
k
u
k
u
k
.If ∥
g
k
+
1
∥
<
ε, then stop; otherwise, go to step (ii).Again we emphasize that we need to solve∇
f
=
0 for the minimization problem (1), and hence the residual means the value of ∥
∇
f
∥
=
∥
g
∥. The above convergence criterion ∥
g
k
+
1
∥
<
ε means that when the residual norm is smaller than a given error tolerance ε, the iterations are terminated.


## 3.4. A Critical Value forα

In Sections3.2 and 3.3 we have used ∂
s
/
∂
α
=
0 (or equivalently, ∂
a
0
/
∂
α
=
0) to find the optimal value of α in the descent vector u
=
g
+
α
B
g. Usually, this value of α obtained from ∂
s
/
∂
α
=
0 is not the global minimum of a
0 (or s). Here, we try another approach and attempt to derive a better value of α than α
o, such that the value of α obtained in this manner is the global minimum of a
0 (or s).In practice, we can take(43)
a
0
:
=
ϕ
u
T
A
u
(
g
·
u
)
2
=
a
s
.
When a
s is near to 0.5, the convergence speed is very fast. Inserting (15) for u into the above equation and through some elementary operations we can derive a quadratic equation to solve α:
(44)
e
1
α
2
+
e
2
α
+
e
3
=
0
,
where
(45)
e
1
:
=
ϕ
u
2
T
A
u
2
-
a
s
(
g
·
u
2
)
2
,
(46)
e
2
:
=
2
ϕ
u
1
T
A
u
2
-
2
a
s
g
·
u
1
g
·
u
2
,
(47)
e
3
:
=
ϕ
u
1
T
A
u
1
-
a
s
(
g
·
u
1
)
2
.If the following condition is satisfied:(48)
D
:
=
e
2
2
-
4
e
1
e
3
≥
0
,
then α in (44) has areal solution:
(49)
α
=
D
-
e
2
2
e
1
.Inserting (45)–(47) into thecritical equation:
(50)
D
=
e
2
2
-
4
e
1
e
3
=
0
,
we can derive an algebraic equation to determine that a
s is the lowest bound of (48).In this lowest bound
a
s
is a critical value denoted by
a
c, and for all a
s
≥
a
c it can satisfy (48) automatically. From (50) through some elementary operations, the critical value a
c can be solved as
(51)
a
c
=
ϕ
[
u
1
T
A
u
1
u
2
T
A
u
2
-
(
u
1
T
A
u
2
)
2
]
u
1
T
A
u
1
(
g
·
u
2
)
2
+
u
2
T
A
u
2
(
g
·
u
1
)
2
-
2
u
1
T
A
u
2
g
·
u
1
g
·
u
2
.
Then by inserting it for a
s into (45) and (46) we can obtain acritical value  α
c for α from (49):
(52)
α
c
=
a
c
g
·
u
1
g
·
u
2
-
ϕ
u
1
T
A
u
2
ϕ
u
2
T
A
u
2
-
a
c
(
g
·
u
2
)
2
,
where D
=
0 was used in view of (50).By inserting (51) for a
c into (52) and cancelling out the common term ϕ we can derive a final equation for α
c, of which we use the same symbols for saving notations:
(53)
a
c
=
u
1
T
A
u
1
u
2
T
A
u
2
-
(
u
1
T
A
u
2
)
2
u
1
T
A
u
1
(
g
·
u
2
)
2
+
u
2
T
A
u
2
(
g
·
u
1
)
2
-
2
u
1
T
A
u
2
g
·
u
1
g
·
u
2
,
(54)
α
c
=
a
c
g
·
u
1
g
·
u
2
-
u
1
T
A
u
2
u
2
T
A
u
2
-
a
c
(
g
·
u
2
)
2
.Here we must emphasize that, in the current descent vectoru
=
g
+
α
B
g, the above value α
c is thebest one, and the vector
(55)
u
=
g
+
α
c
B
g
(
best
descent
vector
)
is thebest descent vector. Due to its criticality, if one attempts to find a better value of the parameter
α
than
α
c
, there would be no real solution of
α. Furthermore, the best descent vector is also better than the optimal vector u
=
g
+
α
o
B
g derived in Section 3.2.


## 3.5. A Globally Optimal Algorithm

Then, we can derive the followingglobally optimal algorithm (GOA) to solve the minimization problem in (1).(i)
Select0
≤
γ
<
1 and give an initial guess of x
0.
(ii)
Fork
=
0,1
,
2
,
…, we repeat the following iterations:
(56)
u
1
k
=
g
k
,
u
2
k
=
B
k
g
k
,
a
c
k
=
(
u
1
k
·
A
k
u
1
k
u
2
k
·
A
k
u
2
k
-
(
u
1
k
·
A
k
u
2
k
)
2
)
×
(
u
1
k
·
A
k
u
1
k
(
g
k
·
u
2
k
)
2
+
u
2
k
·
A
k
u
2
k
(
g
k
·
u
1
k
)
2
-
2
u
1
k
·
A
k
u
2
k
g
k
·
u
1
k
g
k
·
u
2
k
)
-
1
,
α
k
=
a
c
k
g
k
·
u
1
k
g
k
·
u
2
k
-
u
1
k
·
A
k
u
2
k
u
2
k
·
A
k
u
2
k
-
a
c
k
(
g
k
·
u
2
k
)
2
,
(
critical
optimal
parameter
)
,
u
k
=
g
k
+
α
k
B
k
g
k
,
(
best
descent
vector
)
,
(57)
x
k
+
1
=
x
k
-
(
1
-
γ
)
g
k
·
u
k
u
k
A
k
u
k
u
k
.If x
k
+
1 converges according to a given stopping criterion ∥
g
k
+
1
∥
<
ε, then stop; otherwise, go to step (ii).Remark 3.
We have derived anovel globally optimal algorithm for solving the minimization problem in (1). In terms of the descent vector u
=
g
+
α
c
B
g, the GOA is the best one, which leads to theglobal minimum of a
0 (or s) and hence thelargest convergence rate. While the parameter γ is chosen by the user with problem dependence, the parameter α
k is exactly given by (56). Up to here we have successfully derived a novel best descent vector algorithm, with the help from (19), (53), and (54).Remark 4.
At the very beginning we have set an invariant manifoldh
(
x
,
t
)
=
Q
(
t
)
ϕ
(
x
) in (8) as our starting point to derive the iterative optimal algorithms, which includes a locally objective function ϕ in the governing equations. However, in the final stage the terms which include ϕ can be cancelled out, and thus we have obtained the optimal algorithm in (42) and the globally optimal algorithm in (57), which are both independent of ϕ.


## 4. The Broyden-Fletcher-Goldfarb-Shanno Updating Techniques

In the above we have derived two optimal algorithms by leaving the descent matrixB to be specified by the user. First we fix B to be the exact Hessian matrix A, and then we can obtain two optimal algorithms OA and GOA.We can also apply the technique of BFGS by updating the Hessian matrixA and its inverse matrix B. To derive this updating technique let us mention the Newton iterative scheme to solve (1):
(58)
x
k
+
1
=
x
k
-
λ
k
∇
2
f
(
x
k
)
-
1
g
k
,
where λ
k is the optimal steplength along the Newton direction ∇
2
f
(
x
k
)
-
1
g
k at the kth step. In order to construct a matrix B
k to approximate ∇
2
f
(
x
k
)
-
1 we can analyze the relation between ∇
2
f
(
x
k
)
-
1 and the first order derivative g
k. We take the Taylor expansion of f
(
x
) at a point x
k
+
1, obtaining
(59)
f
(
x
)
≈
f
(
x
k
+
1
)
+
∇
f
(
x
k
+
1
)
T
(
x
-
x
k
+
1
)
+
1
2
(
x
-
x
k
+
1
)
T
∇
2
f
(
x
k
+
1
)
(
x
-
x
k
+
1
)
.
Then we have
(60)
∇
f
(
x
k
)
≈
∇
f
(
x
k
+
1
)
+
∇
2
f
(
x
k
+
1
)
(
x
k
-
x
k
+
1
)
.
Let
(61)
p
k
=
x
k
+
1
-
x
k
,
(62)
q
k
=
∇
f
(
x
k
+
1
)
-
∇
f
(
x
k
)
,
and from (60) we have
(63)
q
k
≈
∇
2
f
(
x
k
+
1
)
p
k
.
Assume that the Hessian matrix ∇
2
f
(
x
k
+
1
) is invertible and denote the inverse matrix by B
k
+
1. Then from the above equation we have the so-called quasi-Newton condition:
(64)
p
k
=
B
k
+
1
q
k
.When we take the descent matrixB to be the inverse of the Hessian matrix A, we might accelerate the convergence speed, which is however a more difficult task when the dimension of the Hessian matrix is quite large. So far, as that done in the quasi-Newton method, we can employ the following updating technique due to BFGS:
(65)
B
k
+
1
=
B
k
+
(
1
+
q
k
T
B
k
q
k
q
k
·
p
k
)
p
k
p
k
T
q
k
·
p
k
-
p
k
q
k
T
B
k
+
B
k
q
k
p
k
T
q
k
·
p
k
.
The advantage by using the above updating technique is that we can obtain an approximation of the inverse of the Hessian matrix, and we do not need to really calculate A
-
1.By the same token we can also apply the technique of BFGS to update the Hessian matrix without needing the computation of the real Hessian matrix:(66)
A
k
+
1
=
A
k
+
q
k
q
k
T
q
k
·
p
k
-
A
k
p
k
p
k
T
A
k
p
k
·
(
A
k
p
k
)
,
where at the first step we can take A
0
=
I
n. The advantage of this approach is that we do not need to calculate the Hessian matrix exactly, as developed independently by Broyden [21], Fletcher [22], Goldfarb [23], and Shanno [24]. It is easy to check that the updating technique in (66) satisfies the quasi-Newton condition:
(67)
A
k
+
1
p
k
=
q
k
,
which is an inversion relationship of (64).The present notations ofA and B are the same as those B and H used by Dai [25, 26].The above technique together with the OA by usingB
k
=
A
k generated from (66) will be named the OABFGS, and two numerical examples will be given to test its performance.When the matrixA is given exactly by the Hessian matrix in the OA and GOA, we use the above technique to compute B
k in the OA and GOA, and the resultant iterative algorithms are named OABFGS1 and GOABFGS1, respectively. Finally, if in the OA and GOA both A and B are updated by (66) and (65), the resultant iterative algorithms will be named OABFGS2 and GOABFGS2, respectively.In order to show the high performance of the optimal algorithms OABFGS1, GOABFGS1, and OABFGS2, we will also apply the DFP method [18, 19] to solve nonlinear minimization problems, of which the updating technique is
(68)
B
k
+
1
=
B
k
+
p
k
p
k
T
q
k
·
p
k
-
B
k
q
k
q
k
T
B
k
q
k
·
(
B
k
q
k
)
.
For each iterative step we search the minimization of f
(
x
k
-
λ
B
k
∇
f
(
x
k
)
) to find λ by solving the optimality equation ∇
f
(
x
k
-
λ
B
k
∇
f
(
x
k
)
)
·
(
B
k
∇
f
(
x
k
)
)
=
0 by using the half-interval method. Gill et al. [27], among several others, have shown that the modification in (65) performs more efficiently than that in (68) for most problems.


## 5. Numerical Examples

In order to evaluate the performance of the newly developed methods let us investigate the following examples. Some results are compared with those obtained from the steepest descent method (SDM) and the DFP method [18, 19]. In the above algorithms there exists a relaxation parameter γ, which is problem dependent. A good parameter value of γ can be selected easily by comparing the convergence speeds for different values of γ.Example 1.
As the first testing example of OA and GOA we use the following function given by Rosenbrock [28]:
(69)
min
⁡
{
f
=
100
(
x
2
-
x
1
2
)
2
+
(
1
-
x
1
)
2
}
.
In mathematical optimization, the Rosenbrock function is a nonconvex function used as a performance test benchmark problem for optimization algorithms. It is also known as Rosenbrock’s valley or Rosenbrock’s banana function. The minimum is zero occurring at (
x
1
,
x
2
)
=
(
1,1
). This function is difficult to be minimized because it has a steep sided valley following the parabolic curve x
1
2
=
x
2. Kuo et al. [29] have used the particle swarm method to solve this problem; however, the numerical procedures are rather complex. Liu and Atluri [30] have applied a fictitious time integration method to solve the above problem under the constraints of x
1
≥
0 and x
2
≥
0, whose accuracy can reach to the fifth order.
We apply the OA to this problem by starting from the initial point at(
3,2
) and under a convergence criterion ε
=
1
0
-
10. The SDM is run 3749 steps as shown in Figures 1(a) and 1(b) by solid lines for showing the residual and the objective function f. The SDM can reach a very accurate minimum value of f with 1.22
×
1
0
-
20. The OA with γ
=
0 converges very fast with only 6 steps, with the residual and the objective function f being shown in Figures 1(a) and 1(b) by dashed lines. The OA is faster than the SDM with about 625 times, and furthermore the minimum value of f can be reduced to 1.925
×
1
0
-
25. In Figure 1(c) we compare the solution paths generated by the SDM and OA. When the SDM is moving very slowly along the valley, the OA is moving outside the region of valley and is convergent very fast to the solution. In Figure 1(c) the red dashed line is used to represent the valley of the Rosenbrock function, which is not used to indicate the result of a numerical method. As shown in Figure 1(b), GOA is slightly better than OA, although the paths of OA and GOA seem to be identical in Figure 1(c).
The SDM usually works very well during early stages; however, as a stationary point is approached, it behaves poorly, taking small nearly orthogonal steps. On the other hand, with the help of the optimal directiong
+
α
o
A
g the OA can fast reach to the final end point with high accuracy. For this example the GOA spends the same number of steps as the OA; however, the GOA gives very accurate minimum value of f
=
1.26
×
1
0
-
29.For the Rosenbrock problem, (a) the residuals, (b) the objective functionf, and (c) the solution paths of SDM, OA, and GOA.
(a)
(b)
(c)Example 2.
Next, we consider a generalization of the Rosenbrock function [31, 32]:
(70)
min
⁡
{
f
=
∑
k
=
1
n
-
1
‍
[
100
(
x
k
+
1
-
x
k
2
)
2
+
(
1
-
x
k
)
2
]
}
.
This variant has been shown to have exactly one minimum for n
=
3 at (
x
1
,
x
2
,
x
3
)
=
(
1,1
,
1
) and exactly two minima for 4
≤
n
≤
7. The global minimum happens at (
x
1
,
x
2
,
…
,
x
n
)
=
(
1,1
,
…
,
1
) and a local minimum is near (
x
1
,
x
2
,
…
,
x
n
)
=
(
-
1,1
,
…
,
1
). This result is obtained by setting the gradient of the function equal to zero, noticing that the resulting equation is a rational function of x
i. For small n the polynomials can be determined exactly and Sturm’s theorem can be used to determine the number of real roots, while the roots can be bounded in the region of |
x
i
|
<
2.4 [32]. For larger n this method breaks down due to the size of the coefficients involved.
We apply the OA to this problem withn
=
30, starting from x
i
=
0.1
i, and under a convergence criterion ε
=
1
0
-
5. Under the above stopping criterion, the SDM is run over 13830 steps as shown in Figure 2(a) by solid lines for showing the residual and f. The SDM can reach a very accurate minimum value of f with 1.357
×
1
0
-
10. The OA with γ
=
0.2 converges with 9956 steps, with the residual and f being shown in Figure 2 by dashed lines. The OA is faster than the SDM, and similarly the minimum of f can be reduced to 1.96
×
1
0
-
11.
At the same time we show the values of the optimalα in Figure 3(a) for the OA. The value of α is much smaller than 1, which means that the term g plays a dominant role in the descent direction; however, the influence of α
A
g, although a little, is a key point to speed up the convergence.
Then we apply the GOA to this problem under the same conditions as that in OA. The GOA withγ
=
0.1 converges with 9846 steps, with the residual and f being shown in Figure 2 by dashed-dotted lines. The GOA is faster than the SDM, and the minimum value of f can be further reduced to 9.59
×
1
0
-
12. We show the values of the optimal α for the GOA in Figure 3(b). This problem is quite difficult, and many other algorithms are failure to solve this problem.For the generalized Rosenbrock problem, (a) the residuals and (b) the objective functions obtained by SDM, OA, and GOA.
(a)
(b)In the solution of the generalized Rosenbrock problem by the OA and GOA, showing the variation of optimal values ofα for (a) OA and (b) GOA.
(a)
(b)Example 3.
We consider a problem due to Powell [33]:
(71)
min
⁡
⁡
{
f
=
(
x
1
+
10
x
2
)
2
+
5
(
x
3
-
x
4
)
2
m
m
m
m
l
l
+
(
x
2
-
2
x
3
)
4
+
10
(
x
1
-
x
4
)
4
}
.
The minimum of f is zero occurring at (
x
1
,
x
2
,
x
3
,
x
4
)
=
(
0,0
,
0,0
). We apply the OA to this problem, starting from (
x
1
,
x
2
,
x
3
,
x
4
)
=
(
3
,
-
1,0
,
1
) and under a convergence criterion of ε
=
1
0
-
6. The SDM converges very slowly over 50,000 steps as shown in Figure 4 by solid lines for showing a
0, s, and the residual. The SDM can reach a very accurate value of f
=
9.938
×
1
0
-
9. At the same time, the OA with γ
=
0.15 converges with 349 steps, with a
0, s, and the residual shown in Figure 4 by dashed lines. The OA is faster than the SDM over 130 times, and furthermore we can get a more accurate f
=
8.33
×
1
0
-
10.
As shown in Figure4(c) the residual obtained by the OABFGS is rapidly growing to the order 1
0
26 and then fast decaying to the order 1
0
-
6, and through 1043 steps the OABFGS leads to a better solution with f
=
1.285
×
1
0
-
9 than that obtained by the SDM. Here the combination of the optimal algorithm with the BFGS updating technique is very useful when the exact Hessian matrix is not available or when the computation of the Hessian matrix is cumbersome.
Then we apply the GOA in Section3.5 to the Powell problem. Under the same convergence criterion of ε
=
1
0
-
6, the GOA with γ
=
0.001 converges only with 96 steps as shown in Figure 4(c), where the GOA is faster than the SDM with 500 times and than OA with 3.5 times, and furthermore we can get a very accurate minimum value f
=
1.55
×
1
0
-
9. In Figure 5 we compare the values of α obtained by the OA and the GOA, of which we can observe that both α are quite small and may be negative.
Then, we apply the OABFGS1 and the GOABGGS1 mentioned in Section4 to solve the Powell problem. Under the same convergence criterion of ε
=
1
0
-
6, the OABFGS1 with γ
=
0.1 converges only with 29 steps as shown in Figure 6(a), where we can get very accurate value of f
=
7.57
×
1
0
-
12. Similarly, the GOABFGS1 with γ
=
0.1 converges with 30 steps as shown in Figure 6(a), where we can get very accurate value of f
=
5.28
×
1
0
-
10. In Figure 6(b) we compare the values of α obtained by the OABFGS1 and GOABFGS1, of which we can observe that both α are quite large, which means that the quasi-Newton direction is a dominant factor to accelerate the convergence speed.
Then, we apply the OABFGS2 to solve the Powell problem withγ
=
0, which converges through 116 steps and with the value of f to be f
=
2.77
×
1
0
-
11, which is very accurate. In Figures 6(c) and 6(d) we show the residual and the value of α obtained by the OABFGS2.
In order to reveal the high performance of the optimal algorithms OABFGS1, GOABFGS1, and OABFGS2, we apply the DFP method [18, 19] to solve the Powell problem, which converges through 131 steps as shown in Figure 6(a). For each iteration step we search the minimization of f
(
x
k
-
λ
B
k
∇
f
(
x
k
)
) to find λ by solving the optimality equation ∇
f
(
x
k
-
λ
B
k
∇
f
(
x
k
)
)
·
(
B
k
∇
f
(
x
k
)
)
=
0 by using the half-interval method, with two end points fixed by λ
=
0 and λ
=
1, and the convergence criterion is given by ε
=
1
0
-
10. At the end of the iteration we obtain the minimum value of f to be f
=
9.79
×
1
0
-
12, which is very accurate.For the Powell problem, comparing (a)a
0 and (b) s of SDM and OA and (c) the residuals obtained by the SDM, OA, GOA, and OABFGS.
(a)
(b)
(c)Figure 5
For the Powell problem, comparingα obtained by the OA and GOA.For the Powell problem, comparing (a) the residuals and (b)α obtained by the OABFGS1 and GOABFGS1 with DFP and showing (c) the residual and (d) α obtained by OABFGS2.
(a)
(b)
(c)
(d)Example 4.
In this example we design an office block inside a structure with a curved roof given byx
=
100
-
y
2. Suppose that the number of total cuboids is n and each cuboid can have different size, and we attempt to find the dimensions of all cuboids with maximum volume which would fit inside the given roof structure; that is,
(72)
max
⁡
⁡
{
f
=
y
1
[
100
-
y
1
2
]
+
y
2
[
100
-
(
y
1
+
y
2
)
2
]
m
m
m
m
m
+
…
+
y
n
[
100
-
(
y
1
+
…
+
y
n
)
2
]
}
,
where y
i
>
0 is the height of the ith cuboid.
The maximum off is tending to 2000
/
3 when n is increasing. When n
=
95, we apply the GOA to this problem by starting from y
i
=
0.05 and under a convergence criterion of ε
=
1
0
-
5. The SDM is convergent with 6868 steps as shown in Figure 7 by solid lines for showing a
0, s, residual, and f. At the same time, the GOA with γ
=
0.35 converges with 96 steps, with a
0, s, residual, and f being shown in Figure 7 by dashed lines. Both f are tending to 661.9945. The GOA is faster than the SDM with 80 times. The heights and widths of the cuboids with respect to the number of floors are plotted in Figure 8. For this problem both OA and GOA lead to the same numerical results.
Moreover, we also apply the OABFGS and the OABFGS1 as well as the DFP to this problem under the above same initial condition and convergence criterion, where for the DFP we use the half-interval method to solve the local minimization to find the bestλ with two end points fixed by λ
=
0 and λ
=
1, and the convergence criterion is given by ε
=
1
0
-
10. The values of γ used in the OABFGS and OABFGS1 are, respectively, 0.3 and 0.05. The residuals of these three methods are compared in Figure 9(c), from which we can observe that the OABFGS converges with 35 steps, the OABFGS1 converges with 32 steps, and the DFP converges with 46 steps. They are all better than the above OA and GOA algorithms. The value of α as shown in Figure 9(a) for the OABFGS is quite small, which indicates that the descent direction is dominated by the gradient direction. The value of α as shown in Figure 9(b) for the OABFGS1 is quite large, which indicates that the descent direction is dominated by the quasi-Newton direction.For the maximum area under a given curve comparing (a)a
0, (b) s, (c) residual error, and (d) f obtained by the SDM and GOA.
(a)
(b)
(c)
(d)Figure 8
Showing the heights and widths of the floors with respect to the number of floors.For Example4, showing (a) α of OABFGS and (b) α of OABFGS1 and (c) comparing the residuals obtained by the DFP, OABFGS, and OABFGS1.
(a)
(b)
(c)Example 5.
In this example we test the minimization of the Schwefel function withn
=
100:
(73)
min
⁡
{
f
=
∑
i
=
1
n
‍
(
∑
j
=
1
i
‍
x
j
)
2
}
.
The minimum is zero at x
j
=
0
,
j
=
1
,
…
,
n.
We apply the OA to this problem by starting fromx
i
=
1 and under a convergence criterion of ε
=
1
0
-
4. The SDM does not converge within 10000 steps as shown in Figure 10 by solid lines for showing a
0, s, residual, and f. At the same time, the OA with γ
=
0.1 converges with 276 steps, with a
0, s, residual, and f being shown in Figure 10 by dashed lines. The OA is faster than the SDM over 100 times, and f is tending to 3.99
×
1
0
-
10 which is more accurate than 1.733
×
1
0
-
4 obtained by the SDM. When we apply the GOA with γ
=
0.05 to this problem, it is convergence with 299 steps and is faster than the SDM over 100 times, and f is tending to 7.38
×
1
0
-
10 which is more accurate than SDM and OA.For the minimization of Schwefel function comparing (a)a
0, (b) s, (c) residual error, and (d) f obtained by the SDM and OA.
(a)
(b)
(c)
(d)Example 6.
In this example we test the minimization of the Whitley function withn
=
8:
(74)
min
⁡
{
f
=
∑
i
=
1
n
‍
∑
j
=
1
n
‍
[
y
j
i
2
4000
-
cos
⁡
y
j
i
+
1
]
}
,
where y
j
i
=
100
(
x
i
-
x
j
2
)
2
+
(
x
j
-
1
)
2. The minimum is zero at x
j
=
1
,
j
=
1
,
…
,
n. Here, we fix n
=
8.
It is very difficult to minimize the Whitley function. The SDM diverges as shown in Figure11 by solid lines. We apply the OA to this problem, starting from x
i
=
1.12 and under a convergence criterion of ε
=
1
0
-
8. The OA with γ
=
0.06 converges with 24 steps, with residual and f being shown in Figure 11 by dashed lines and f tending to 1.47
×
1
0
-
13.For the minimization of Whitley function comparing (a) residuals and (b)f obtained by the SDM and OA, where the SDM is failure.
(a)
(b)


## 6. Conclusions

By formulating the minimization problem into a continuous manifold, we can derive a governing system of ODEs for deriving the iterative algorithms which being proven convergence to find the unknown minimum pointx of a given nonlinear minimization function. The novel algorithm is named “an optimal algorithm (OA)," because in the local frame we have derived the optimal parameter of α in the descent direction, which is a linear combination of the gradient vector and a supplemental vector. We have demonstrated a critical descent vector to derive aglobally optimal algorithm (GOA), which can substantially accelerate the convergence speed in the numerical solution of nonlinear minimization problem. It was proven that the critical value α
c in the critical descent vector leads to thelargest convergence rate among all the descent vectors specified by u
=
g
+
α
B
g. Due to its criticality, if one attempts to find a better descent vector than u
=
g
+
α
c
B
g, there would be no real descent vector of u. Through several numerical tests we found that the both the OA and the GOA outperformed very well. Then we have proposed novel algorithms based on OA and GOA by replacing the Hessian matrix A and the descent matrix B with the updating techniques of BFGS for one or for both of these two matrices. Two numerical examples were given to test the performances of these algorithms, which are faster than the original OA and GOA algorithms, due to the enhancement by using the quasi-Newton conditions on the updated matrices.



---

*Source: 324181-2014-03-12.xml*