Quasi-Newton Methods

Introduction to Data Science Unconstrained Optimization Quasi-Newton Methods

Quasi-Newton methods is a framework of algorithims that is built on the Newton-Rhason method and the line search in this method is based on a $n*n$ matrix S, which serves as the same purpose as the inverse Hessian matrix in the Newton-Rhapson Method. The Quasi-Newton method is yet another framework of algorithims that do not rquire the Hessian. The expression in Quasi-Newton method can be given as

$x_{k+1}=x_k-\alpha_kS_kg_k$

When $S_K=I$ , it is the algorithim for the steepest-descent method. When $S_k=H^{-1}$ , the algorithim is for the Newton-Rhapson method. The above equation of $x_{k+1}$ converges fastest when $S_k=H^{-1}$ . Quasi-Newton methods combines the advantages of the steepest-descent and the conjugate-direction methods. They are among the most efficient methods and are used extensively in many applications.

Assume that the difference in successive $x_k$ is $\delta_k=x_{k+1}-x_k$ and the difference in successive gradient is $g_k$ is $\gamma_k=g_{k+1}-g_k$ , then we have $\gamma_k=H\delta_k$ . The Newton-Rhapson method determines the direction vector as $d_k=-S_kg_k$ . The magnitude $\alpha_k$ can be determined by a line search minimizing the function $f(x_k+\alpha d_k)$ . Now here we want to avoid inverting a matrix such as computing $S_k=H^{-1}$ and checking that $S_k$ is positive definite in each iteration, rather we simply want an updating rule for $S_{k+1}$ such as

$S_{k+1}S_k+C_k$ .

Here $C_k$ is a $n*n$ correction matrix that can be computed form available data. In any definition of $C_k$ , $S_{k+1}$ must satisfy these properties:

1. $S_{k+1}\gamma_k=\delta_k$ .

2. Vectors $\{\delta_0,\delta_1,...,\delta_{n-1}\}$ must be conjugate direction vectors.

3. A positive definite $S_k$ must produce a positive definite $S_{k+1}$ .

The second property states that the advantages of tyhe conjugat-direction method method applies to Quasi-Newton method. The third method states that the direction vector is valid in each iteration. There are many ways to define $C_k$ and they give rise to different variants of the Quasi-Newton method.

Rank-One Method

The Rank-Method method is featured by the correction matrix $C_k$ has a unity rank. Now we assume that

$S_{k+1}\gamma_k=\delta_k$ .

Let

$S_{k+1}=S_k\gamma_k+\beta\xi_k\xi_k^T\gamma_k$

And

$\gamma^T_k(\delta_k-S_k\gamma_k)=\beta_k\gamma^T_k\xi_k\xi_k^T\gamma_k=\beta_k(\xi^T_k\gamma_k)^2$

$(\delta_k-S_k\gamma_k)(\delta_k-S_k\gamma_k)^T=\beta_k(\xi_k^T\gamma_k)^2\beta_k\xi_k\xi_k^T$

Then we have

$\displaystyle \beta_k\xi_k\xi_k^T=\frac{(\delta_k-S_K\gamma_k)(\delta_k-S_k\gamma_k)^T}{\gamma_k^T(\delta_k-S_k\gamma_k)}$

Therefore,

$S_{k+1}=S_k+\beta_k\xi_k\xi_k^T$ .

There are two problems with this method. First, the positive definite $S_k$ may not give a positive definite $S_{k+1}$ . In such a case, the next direction is not a good direction. Second, the denominator in the correction formula may approach or become equal to zero. If that happens, the method breaks down as $S_{k+1}$ is undefined.

Davidon-Fletcher-Powell Method

The Davidon-Fletcher-Powell (DPF) method is similar to the rank-one method but it has one important advantage, that if the initial matrix $S_0$ s positive definite, then the subsequent matrices are always positive definite. Unlike the Raank-One method, every new direction is a descent direction.

The DPF formula is

$\displaystyle S_{k+1}=S_k+\frac{\delta_k\delta_k^T}{\delta_k^T\gamma_k}-\frac{S_k\gamma_k\gamma_k^Ts_k}{\gamma_k^TS_k\gamma_k}$

The correction is an $n*n$ symmetric matrix of rank two.

Broyden-Fletcher-Goldfarb-Shanno Method

The BFGS method formula is

$\displaystyle S_{k+1}=S_k+(1+\frac{\gamma^T_kS_k\gamma_k}{\gamma_k^T\delta_k})\frac{\delta_k\delta_k^T}{\gamma^T_k\delta_k}-\frac{\delta_k\gamma_k^TS_k+S_k\gamma_k\delta_k^T}{\gamma_k^T\delta_k}$

The BFGS method has the following properies:

1. $S_{k+1}$ becomes identical to $H^{-1}$ for $=n-1$ .

2. Direction vectors $\{\delta_i\}_{i1,2,...,n-1}$ form a conjugate set.

3. $S_{k+1}$ is the positive definite if $S_k$ is positive definite.

4. $\delta_k^T\gamma_k=\delta^T_kg_{k+1}-\delta^T_kg_k>0$ .

BFGS method is the best method among all.

Huang Family Method

Huang method is a general formula that includes the rank-one, DFP, BFGS, Pearson and McCormick ethods. It has the following form

$\displaystyle S_{k+1}=S_k+\frac{\delta_k(\theta\delta_k+\phi S^T_k\gamma_k)^T}{(\theta\delta_k+\phi S^T_k\gamma_k)^T\gamma_k}-\frac{S_k\gamma_k(\psi\delta_k+\omega S^T_k\gamma_k)^T}{(\psi\delta_k+\omega S^T_k\gamma_k)^T\gamma_k}$

$\theta,\, \phi,\, \psi$ and $\omega$ are independent parameters.

By letting $theta=1,\, \phi=-1,\, \psi=1,\, \omega=-1,$ we have the Rank-One formula.

By letting $theta=1,\, \phi=0,\, \psi=0,\, \omega=1,$ we have the DFP formula.

By letting $theta=1,\, \phi=0,\, \psi=1,\, \omega=0,$ we have the McCormick formula:

$\displaystyle S_{k+1}=S_k+\frac{(\delta_k-S_k\gamma_k)\delta^T_k}{\delta_k^T\gamma_k}$

By letting $theta=0,\, \phi=1,\, \psi=0,\, \omega=1,$ we have the Pearson Formula:

$\displaystyle S_{k+1}=S_k+\frac{(\delta_k-S_k\gamma_k)\gamma_k^TS_k}{\gamma_k^TS_k\gamma_k}$

By letting $\displaystyle \frac{\phi}{\theta}=\frac{-\gamma_k\delta^T_k}{\gamma_k\delta_k^T+\gamma^T_kS_k\gamma_k},\,\psi=1,\, \omega=0,$ we have the BFGS formula.

The following code minimizes the Himmelblau function $f(x_1,x_2)=(x_1^2+x_2-11)^2+(x_1+x_2^2-7)^2$ using different Quasi-Newton methods.

				
					//The Himmelblau function: 
// f(x) = (x1^2 + x2 - 11)^2 + (x1 + x2^2 - 7)^2
 
val f: RealScalarFunction f= object : new RealScalarFunction() {
override fun Double evaluate(Vector x) {
double x1 = x.get(1),
double x2 = x.get(2),

double result = pow(x1 * x1 + x2 - 11, 2),
result += 12 * pow(x1 + x2 * x2 - 7, 2),
return result
}
override public int dimensionofDomain() {
return 2
}
override public int dimensionofRange() {
return 1
}
}

val g: RealVectorFunction g = object : new RealVectorFunction() {
override fun evaluate(Vector x): {
double x1 = x.get(1),
double x2 = x.get(2),

double w1 = x1 * x1 + x2 - 11,
double w2 = x1 + x2 * x2 -7,
    
double[] result = new double[2],
result[0] = 4 * w1 * x1 + 2 * w2,
result[1] = 2 * w1 + 4 * w2 * x2,
return new DenseVector (result)
}  
override public int dimensionofDomain() {
return 2
}
override public int dimensionofRange() {
return 2
}
}

val problem: C2OptimProblem = new C2OptimProblemImpl(f, g)
val solver : QuasiNewtonMinimizer RankOneMinimizer = new RankOneMinimizer(
1e-16,
15)
val soln1: RankOneMinimizer.solve(problem)
val xmin1: soln1.search(new DenseVector(new double[]{6, 6}))
val fmin1= f.evaluate(xmin1)
System.out.println(String.format("f(%s) = %.16f", xmin1.toString(), fmin1))

val solver : QuasiNewtonMinimizer DFPMinimizer = new DFPMinimizer(
1e-16,
15)
val soln2: DFPMinimizer.solve(problem)
val xmin2: soln2.search(new DenseVector(new double[]{6, 6}))
val fmin2= f.evaluate(xmin2)
System.out.println(String.format("f(%s) = %.16f", xmin2.toString(), fmin2))

val solver : QuasiNewtonMinimizer BFGSMinimizer = new BFGSMinimizer(
false,
1e-16,
15)
val soln3: BFGSMinimizer.solve(problem)
val xmin3: soln3.search(new DenseVector(new double[]{6, 6}))
val fmin3= f.evaluate(xmin3)
System.out.println(String.format("f(%s) = %.16f", xmin3.toString(), fmin3))

val solver : QuasiNewtonMinimizer HuangMinimizer = new HuangMinimizer(
0,
1,
0,
1,
1e-16,
15)
val soln4: HuangMinimizer.solve(problem)
val xmin4: soln4.search(new DenseVector(new double[]{6, 6}))
val fmin4= f.evaluate(xmin4)
System.out.println(String.format("f(%s) = %.16f", xmin4.toString(), fmin4))

val solver : QuasiNewtonMinimizer PearsonMinimizer = new PearsonMinimizer(
1e-16,
15)
val soln5: PearsonMinimizer.solve(problem)
val xmin5: soln5.search(new DenseVector(new double[]{6, 6}))
val fmin5= f.evaluate(xmin5)
System.out.println(String.format("f(%s) = %.16f", xmin5.toString(), fmin5))

The output is:

				
					f([3.000000, 2.000000]) = 0.000000000000000
f([3.000000, 2.000000]) = 0.000000000000000
f([3.000000, 2.000000]) = 0.000000000000000
f([3.000000, 2.000000]) = 0.000000000000000
f([3.000000, 2.000000]) = 0.000000000000000

The Quasi-Newton methods are much more efficient than any of the Conjugate-Drection methods. They take less iterations than the Conjugate-Direction methods and give etter accuracy and precision.

Previous Topic

Back to Lesson

Next Topic