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


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 


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 .






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)}


 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(
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(
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(
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(
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(
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.