Conjugate Gradient Methods

Richard Fears and Jason Gregersen

Purpose

To test the efficiency of different conjugate gradient methods on a variety of test problems.

Introduction

The Conjugate gradient method is an optimization technique that is very similar to a line search. We are trying to minimize some function ϕ(x) given some starting position. The general line search process says to choose a downhill direction and take some step. We then choose another direction and another step. The specifics of the Conjugate Gradient method are, how do I form a new direction and how big is my step. Written as the skeleton of an algorithm we could say that if I have some point (Subscript[x, k]) that I found using the direction (Subscript[p, k]), then I need to choose a new direction (Subscript[p, k+1]) and a step size α such that (Subscript[x, k+1])=(Subscript[x, 1])+α (Subscript[p, k+1])with the sequence of (Subscript[x, k's]) converging to the minimum value. For the method of Conjugate Gradients, to choose our next step direction, we are going to take a linear combination of the gradient at our current point and our previous direction i.e.: (Subscript[p, k+1])=-∇f((Subscript[x, 1]))+β (Subscript[p, k]) where β is chosen such that (Subscript[p, k+1]) and (Subscript[p, k]) are conjugate ((Subscript[p, k+1])A (Subscript[p, k])=0). With this approach there are multiple choices for choosing β. Deciding which of these options is best on a specific problem, or in general, is the main purpose of this paper. The options we are going to choose between are ([4] pg 2):

the Fletcher-Reeves update 		(Subscript[B, k+1]^Fr)=((∇(Subscript[f, k+1])ᵀ∇(Subscript[f, k+1]))/(∇(Subscript[f, k])ᵀ∇(Subscript[f, k]))),
the Polak-Ribiere update		(Subscript[B, k+1]^Pr)=((∇(Subscript[f, k+1])ᵀ(∇(Subscript[f, k+1])-∇(Subscript[f, k])))/(||∇(Subscript[f, k])((||)^2)))
the Hestenes-Stiefel update 		(Subscript[B, k+1]^Hs)=((∇(Subscript[f, k+1])ᵀ(∇(Subscript[f, k+1])-∇(Subscript[f, k])))/((∇(Subscript[f, k+1])-∇(Subscript[f, k]))ᵀ(Subscript[p, k])))
the Dai and Yuan update		(Subscript[B, k+1]^Dy)=((||∇(Subscript[f, k+1])((||)^2))/((∇(Subscript[f, k+1])-∇(Subscript[f, k]))ᵀ(Subscript[p, k])))
the Hager and Zhang update	 	(Subscript[B, k+1]^num)=(OverHat[((Subscript[y, k])])-2(Subscript[p, k]) ((||(OverHat[(Subscript[y, k])])((||)^2))/((OverHat[(Subscript[y, k])])ᵀ(Subscript[p, k]))))ᵀ((∇(Subscript[f, k+1]))/((OverHat[(Subscript[y, k])])ᵀ(Subscript[p, k]))) with (OverHat[(Subscript[y, k])])=∇(Subscript[f, k+1])-∇(Subscript[f, k])
the "Conjugate Descent" update 	(Subscript[B, k+1]^Cd)=((∇(Subscript[f, k+1])ᵀ∇(Subscript[f, k+1]))/((∇(Subscript[f, k+1])-∇(Subscript[f, k]))ᵀ(Subscript[p, k]))),
the Lui and Storey update 		(Subscript[B, k+1]^Ls)=((∇(Subscript[f, k+1])ᵀ(∇(Subscript[f, k+1])-∇(Subscript[f, k])))/(-(∇(Subscript[f, k]))ᵀ(Subscript[p, k])))
the "Conjugate Descent" update 	(Subscript[B, k+1]^Cd)=((∇(Subscript[f, k+1])ᵀ∇(Subscript[f, k+1]))/((∇(Subscript[f, k+1])-∇(Subscript[f, k]))ᵀ(Subscript[p, k]))),

Testing Procedures

To compare the above update methods we will be comparing the number of iterations and the wallclock times to achieve a solution . Since we will also be using a variety of problems, we will get a chance to see the robustness of the algorithms as well.

Project Outline

1. We will code the Nonlinear Conjugate Gradient method as well as the various (Subscript[β, k]) updates mentioned above.
2. We will then compare the performance of the various algorithms on the following test problems. 
	a. A high dimensional linear problem
	b. 2-D Non-linear problems (easy, Himmelblau, Siam, Rosenbrock)
	c. Higher dimensional non-linear problems.

Code

We coded a version of the Conjugate Gradient code in which we were able to input the specific (Subscript[β, k]) being used. Once the new direction in the algorithm was computed we then chose to use Mathemtica’s “Nminimize” command with an accuracy and precision goal of 4 to find the necessary step length α. In our Code we have also chosen to reset our value for (Subscript[β, k]) back to steepest descent after 20 iterations. This should account for the loss of conjugacy in directions that develops when using the algorithm for nonlinear functions.

Bk Updates

Fr[rk1_,rk2_,pk_]:=((rk2.rk2)/(rk1.rk1))

Pr[rk1_,rk2_,pk_]:=((rk2.(rk2-rk1))/(rk1.rk1))

Hs[rk1_,rk2_,pk_]:=((rk2.(rk2-rk1))/((rk2-rk1).pk))

Dy[rk1_,rk2_,pk_]:=((rk2.rk2)/((rk2-rk1).pk))

num[rk1_,rk2_,pk_]:=((rk2-rk1)-2 pk ((((rk2-rk1).(rk2-rk1))/((rk2-rk1).pk)) )).(rk2/((rk2-rk1).pk))

Cd[rk1_,rk2_,pk_]:=((rk2.rk2)/(-pk.rk1))

Ls[rk1_,rk2_,pk_]:=((rk2.(rk2-rk1))/(-(rk1).pk))

Conjugate Gradient Code

ConjGrad[y0_,f_,Tol_,Bkupd_]:=
Module[{α,pk,t,y1,Bk,pk2,rk1,rk2,df,dim=Length[y0],xvec,i=0,reset},
xvec=GenerateXVec[dim];(*Array[(Subscript[x, #])&,dim];general vector of correct dimension*)
df=D[f[xvec],{xvec}];(* gradient*)
pk=-df/.Table[xvec[[i]]->y0[[i]],{i,1,dim}];(* first pk is steepest descent, the /. is evaluating the gradient at y0*)
t=α/.FindMinimum[{f[y0+α*pk],0≤α≤1},α][[2]];(*minimizing along pk*)
y1=y0+t*pk;
Sow[y0,"points"];
Sow[y1,"points"];
While[Norm[df/.Table[xvec[[i]]->y1[[i]],{i,1,dim}]]>Tol&&i<200,
reset=0;
rk2=df/.Table[xvec[[i]]->y1[[i]],{i,1,dim}];
rk1=df/.Table[xvec[[i]]->y0[[i]],{i,1,dim}];
(*If[(rk2-rk1).pk<1*(10^(-14)),Break[]];*)
If[reset<20,Bk=Bkupd[rk1,rk2,pk],Bk=reset=0];
pk2=-rk2+Bk pk;
t2=α/.NMinimize[f[y1+α*pk2],α,PrecisionGoal->4,AccuracyGoal->4][[2]];
y2=y1+t2*pk2;
reset++;
{y1,pk,i}={y2,pk2,i+1};
Sow[y1,"points"];
];
{y1,Norm[df/.Table[xvec[[i]]->y1[[i]],{i,1,dim}]],i}
]

Miscellaneous Code

First two codes from Struthers[3]

RandomSPDMat[n_, Dist_:UniformDistribution[{1,10}]]:= Module[
{Q=RandomReal[{-1,1},{n,n}],λs=Abs[RandomReal[Dist,n]]},
Q=Eigenvectors[Q+Transpose[Q]];
Transpose[Q].DiagonalMatrix[λs].Q
]

RandomQ[n_]:=Module[
{Q},
Q=QRDecomposition[RandomReal[{-1,1},{n,n}]][[1]]
]

GenerateXVec[n_]:=Table[ToExpression[StringJoin["x",ToString[i]]],{i,1,n}]

probgen[function_,dimension_]:=Module[{B,c,dim,Q},
dim=dimension;
B=RandomSPDMat[dim-2];
c=RandomReal[5,dim-2];
big[y_]:=function[y[[1;;2]]]+(1/2)y[[3;;dim]].B.y[[3;;dim]]-c.y[[3;;dim]];
Q=RandomQ[dimension];
bigrot[y_]:=big[Transpose[Q].y];
{B,c,Q}
]

Testing Code

Testing2[f_,Schemes_]:=Module[{time,(*names*),s,sol,err,j,i},
s=Length[Schemes];
names=Table[StringJoin[ToString[Schemes[[i]]],"sol"],{i,1,s}];
j=1;
Do[
time=Timing[{sol,err,i}=ConjGrad[y0,f,Tol,Schemes[[j]]];][[1]];
names[[j]]={Schemes[[j]],time,i,err,sol};
,{j,s}
];
Grid[Join[{{"scheme","time","loops","error"}},Table[names[[i]][[1;;4]],{i,1,s}]],Frame->All]
]

cgVisGenerateOneFrame[f_,xs_,maxi_,picdimopt_:-1,borderopt_:0.1]:=Module[
{n=Length[xs[[1]]],picdim=picdimopt,border=borderopt,arrows},
(* check for stupid user *)
If[maxi<2||maxi>Length[xs],Print["Max iteration specification faulty. Must be between 2 and ",Length[xs],"."];Return[]];
(* if the user did not provide dimensions for the picture *)
If[IntegerQ[picdim],
(* create them by finding the maximum and minimum x and y values *)
picdim=Table[
{Min[xs[[1;;maxi,Table[j,{j,i,n,2}]]]],
Max[xs[[1;;maxi,Table[j,{j,i,n,2}]]]]},
{i,1,2}]
];
(* decrease min x and min y by a spacing amount *)
picdim[[1,1]]-=border;picdim[[2,1]]-=border;
(* decrease min x and min y by a spacing amount *)
picdim[[1,2]]+=border;picdim[[2,2]]+=border;
(* create the arrows for the last two points *)
arrows=Table[
Arrow[{xs[[maxi-1,(i*2-1);;(i*2)]],xs[[maxi,(i*2-1);;(i*2)]]}]
,{i,1,n/2}
];
(* finally plot it all *)
Show[
ContourPlot[f[{x,y}],{x,picdim[[1,1]],picdim[[1,2]]},{y,picdim[[2,1]],picdim[[2,2]]},
Epilog->Join[{Black,Thick},arrows]],
ListLinePlot[
Table[xs[[1;;maxi,(i*2-1);;(i*2)]],{i,1,n/2}],
PlotRange->picdim,
PlotStyle->Table[{Hue[i/(n/2)],Thick},{i,1,n/2}]
]
]
]

TestingBig[f_,dim_,schemes_]:=Module[{},
probgen[f,dim];
y0=RandomReal[{1,2},dim];
Testing2[bigrot,schemes]
]

Functions

linfun[dim_]:=Module[{A,b},
A=RandomSPDMat[dim];
b=RandomReal[10,dim];
f[xvec_]=(1/2)*xvec.A.xvec-b.xvec;
{A,b}
]

easy[{x_,y_}]:=(x^2)+(((y-1))^2)-x

H[{x1_,x2_}]:= ((((x1^2)+ x2-11))^2)+ (((x1+(x2^2)-7))^2)

Siam[{x_,y_}]:=Exp[Sin[50 x]]+Sin[60 Exp[y]]+Sin[70 Sin[x]]+Sin[Sin[80 y]]-Sin[10 (x+y)]+1/4 (x^2+y^2)

Ros[{x_,y_}]:=(1-x)^2+100(y-x^2)^2

Initializations

dimension=30;
y0=RandomReal[10,dimension];
Tol=1.*(10^(-6));
xvec=GenerateXVec[dimension];

Schemes={Pr,Fr,Hs,num,Dy,Ls,Cd};
GoodSchemes={Pr,Fr,Dy,Cd};

2-D Testing

Linear

We first chose to run our code on a 30 dimensional  linear problem. We are using a tolerance of (10^(-6)) to end our iteration and a randomly generated inital condition. The command “linfun[]” seen below creates a SPD matrix “A” via code from [4]. It also generates a random “b” vector and then defines a function “f” as the quadratic form f(x)=(1/2)xᵀ A x-xᵀ b. In the process of testing the various (Subscript[β, k]) updates we encountered overflow trouble with several of the updates. Thus the results we are going to use to compare schemes will be based on the Pr, Fr, Dy, and Cd schemes only. As shown in the table below Pr and Fr were slightly less efficient, but they all were effective.

{A,b}=linfun[30];

Testing2[f,GoodSchemes]

Since we are using a positive-definite matrix A, the minimum we find will be a global minimum. Also at the minimum the gradient will be equal to zero and since A is symmetric ∇f=A x-b=0, thus the solution to the minimization problem should be a solution to the system A x-b=0. The results below verify this.

A.names[[4]][[5]]-b (*the 4 indicates that we are using Cd to test the results*)

Non-Linear Testing

Next we want to test the schemes on some 2-D nonlinear functions, but first we need to make sure that that the CG code works on non-linear functions. Setting up the dimensions of our problem and the initial conditions we will run our code on the 2-D Himmelblau function. The results below show that Fr and Cd are the most efficient. Pr is the least efficient and we expect that this is due to it getting jammed at some point and not working properly again until the restart at the 20th iteration.

dimension=2;
y0=RandomReal[10,dimension];

Testing2[H,{Pr,Fr,Dy,Cd}]//Quiet

To test to see if the the schemes are giving us valid results we can plot the points generated during one of the schemes. The plot below uses the points generated from the Pr scheme. We can see that they are converging to a minimum value.

points=Reap[{sol,err,i}=ConjGrad[y0,H,Tol,GoodSchemes[[1]]]][[2]][[1]];
ContourPlot[H[{x,y}],{x,sol[[1]]-.01,sol[[1]]+.01},{y,sol[[2]]-.01,sol[[2]]+.01},Epilog->Point[points]]

Now confident that the Code is running properly we can compare the performance of the different (Subscript[β, k]) schemes on an easy funtion as well as the Siam and Rosenbrock functions. Here our “easy” function is easy[{x_,y_}]:=(x^2)+(((y-1))^2)-x . On this function the solution was so easy that there is no comparison between the schemes. The fact that there is zero for the number of loops indicates that the solution was found in the first step without even running the main loop of the code.

Testing2[easy,GoodSchemes]//Quiet

The results on the Siam function indicate once again that Pr seems to be less efficient (althought not the worst in this case).

Testing2[Siam,GoodSchemes]//Quiet

Te results on the Rosenbrock function indicate that this problem is difficult and all the schemes failed. Here you can also see that we have limited the number of iterations of the algorithm to 100. We believe that using a line search with more restrictions on the stepsize i.e. stronger wolfe conditions and other restrictions to ensure that the combination of (Subscript[β, k]) and α always result in a descent direction[4]

Testing2[Ros,GoodSchemes]//Quiet

Higher Dimensional Problems

Next we will test the different schemes on nonlinear functions of higher dimensions. To generate non-trivial problems we are going to put a non-linear function in the first two dimensions of a larger linear problem. Then use a rotation matrix to mix up the dimensions. These steps are all done by “probgen”. The input is the nonlinear function that goes in the first two dimensions, and the total dimension of the problem. The output is the SPD matrix and vector used to generate the linear dimensions, and the rotation matrix used to mix everything up.

dim=4;
{B,c,Q}=probgen[Siam,dim];
y0=ConstantArray[0,dim];

After generating the problem we solve it using the Dy update.

{sol,err,i}=ConjGrad[y0,bigrot,Tol,Dy];

To test the result we need to un-rotate our solution and then check to see if the solution behaves as expected. Thus after the unrotating, we calculate the residual of the linear system and then plot the solution for the non-linear system, with the results shown below. Notice that although the residual is still relatively high, the plot of the solution appears to be right on. This is a result of the large gradient values in the siam function. The conjugate gradient tends to work on the hard stuff first, thus the focus was on dealing with the nonlinear portion.

sol=Transpose[Q].sol;

B.sol[[3;;dim]]-c

ContourPlot[Siam[{x,y}],{x,sol[[1]]-.1,sol[[1]]+.1},{y,sol[[2]]-.1,sol[[2]]+.1},Epilog->{PointSize->Medium,Red,Point[sol[[1;;2]]]}]

Comfortable that the algorithm is finding solutions we then compared the different updates. The first result is for the the 4 dimensional Siam function that we just looked at. We see that Dy did the best on this function and as we discussed the althought the error might look high , the plot results were more comforting.

Testing2[bigrot,GoodSchemes]//Quiet

Given the underwhelming results of the previous example we decided to try to use an easier function but at higher dimensions. We took our “easy” function and tested it at 7, 10 and 20 dimensions. the results all suggest that Dy and Cd were much quicker and took signifigantly fewer function evaluations.

TestingBig[easy,7,GoodSchemes]

TestingBig[easy,10,GoodSchemes]

TestingBig[easy,20,GoodSchemes]

Next we tried to test the different updates on the Himmelbleau function at four dimensions. The results below showed that none of the updates yielded good results even after increasing the maximum number of iterations to 200. It does however look like while Dy and Cd appear to be divergent, the Pr and Fr schemes to be making some kind of progress.

TestingBig[H,4,GoodSchemes]

multiple copies of nonlinear

The final level of testing has us increasing the level of difficulty by taking our nonlinear function and copying it into multiple dimensions, then performing the same kind of rotation technique as was previously used. We then attempted to find the minimum using the Dy update. To get a visual of the result we created a plot of the solution for each pair of dimension. We then repeated this analysis using the Siam function. The results shown below indicate that for the Himmelbleau function when using an initial condition at the origin the solutions followed the same path in each dimension and seemed to reach a minimum. The close in results also showed thtat upon nearing the valley the points tended to bounce around for a while before converging to the solution. This was a result of the loss of conjugacy in our directions and was ultimately resolved when the reset kicked in. For the Saim however, the results were more sporadic. This was a clear result based on the difficulty of the function.

4 Copies of the Himmelbleau function

copies=4;
n=copies*2;
testpicLNLQ=RandomQ[n];
xvec=GenerateXVec[n];
testpicLNL[xv_]=Sum[H[(Transpose[testpicLNLQ].xvec)[[i*2-1;;i*2]]],{i,1,copies}]/.Table[xvec[[i]]->xv[[i]],{i,1,n}];//Quiet
testpicLNL["IC"]=ConstantArray[0.0,n];

{{sol,err,i},points}=Reap[ConjGrad[testpicLNL["IC"],testpicLNL,1.*(10^(-6)),Dy]];
points=points[[1]];
Do[points[[i]]=Transpose[testpicLNLQ].points[[i]],{i,1,Length[points]}];

cgVisGenerateOneFrame[fH,points,Length[points]]

copies=4;
n=copies*2;
testpicLNLQ=RandomQ[n];
xvec=GenerateXVec[n];
testpicLNL[xv_]=Sum[Siam[(Transpose[testpicLNLQ].xvec)[[i*2-1;;i*2]]],{i,1,copies}]/.Table[xvec[[i]]->xv[[i]],{i,1,n}];//Quiet
testpicLNL["IC"]=ConstantArray[0.0,n];

{{sol,err,i},points}=Reap[ConjGrad[testpicLNL["IC"],testpicLNL,1.*(10^(-6)),Dy]];
points=points[[1]];
Do[points[[i]]=Transpose[testpicLNLQ].points[[i]],{i,1,Length[points]}];

4 copies of Siam

cgVisGenerateOneFrame[Siam,points,Length[points]]

Conclusion

To summarize, we can say that there exist several different options for choosing the parameter (Subscript[β, k]) in the conjugate gradient method. The efficiency and stability of the different updates seems to vary on different problems. In general Dy seemed to be very efficient, but was outperformed by Fr and Pr on the 4 copy version of the Himmelbleau function. This variation of efficiency in the schemes is what has lead to the creation of hybrid schemes that switch between schemes to gain efficiency. As an example a scheme may run combinations of Pr and Fr or Dy and Hs. These choices tend to pair updates with strong convergence properties with those with may not convergen in general but if they do tend to be more efficient. For more information on Hybrid method see Hagar and Zhang[4, page6].

Reference

[1] Nocedal and Wright (2006). NumericalOptimization 2nd ed. Springer

[2] http://reference.wolfram.com/mathematica/tutorial/UnconstrainedOptimizationConjugateGradientMethods.html 

[3] Struthers (2011). LinearConjugateGradient.  http://www.mathlab.mtu.edu/~struther/Courses/5630_11/LinearConjugateGradient.nb

[4] Hagar and Zhang (2006),  A survey of Nonlinear Conjugate Gradient Methods

[5] Jonathan Richard Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain