1 //------------------------------------------------------------------------
2 // Copyright (C) 1993:
3 // J.C. Meza
4 // Sandia National Laboratories
5 // meza@california.sandia.gov
6 //------------------------------------------------------------------------
7
8 #ifdef HAVE_CONFIG_H
9 #include "OPT++_config.h"
10 #endif
11
12 #ifdef HAVE_STD
13 #include <cstring>
14 #include <cmath>
15 #else
16 #include <string.h>
17 #include <math.h>
18 #endif
19
20 #include "OptConstrQNewton.h"
21 #include "cblas.h"
22 #include "ioformat.h"
23 #include <float.h>
24
25
26 using Teuchos::SerialDenseVector;
27 using Teuchos::SerialDenseMatrix;
28 using Teuchos::SerialSymDenseMatrix;
29
30
31 namespace OPTPP {
32
33 //------------------------------------------------------------------------
34 //
35 // Quasi-Newton Method member functions
36 // checkDeriv()
37 // updateH()
38 //------------------------------------------------------------------------
39
40 // static char* class_name = "OptConstrQNewton";
41
42 // Check the analytic gradient with FD gradient
checkDeriv()43 int OptConstrQNewton::checkDeriv()
44 { return checkAnalyticFDGrad(); }
45
46
47 //----------------------------------------------------------------------------
48 //
49 // Update Hessian using a Quasi-Newton update
50 //
51 //----------------------------------------------------------------------------
updateH(SerialSymDenseMatrix<int,double> & Hk,int k)52 SerialSymDenseMatrix<int,double> OptConstrQNewton::updateH(SerialSymDenseMatrix<int,double>& Hk, int k)
53 {
54
55 double mcheps = DBL_EPSILON;
56 double sqrteps = sqrt(mcheps);
57
58 int i;
59
60 NLP1* nlp = nlprob();
61 int nr = nlp->getDim();
62 SerialDenseVector<int,double> grad(nr), xc(nlp->getXc().length());
63 xc = nlp->getXc();
64 grad = nlp->getGrad();
65
66 SerialDenseVector<int,double> D(nr);
67 // BFGS formula
68
69 if (k == 0) { // Initial Hessian is set equal to the Identity Matrix
70 Hessian = 0.0;
71 // D = sx.AsDiagonal()*sx.AsDiagonal();
72 // D = sfx.AsDiagonal();
73 double typx, xmax, gnorm;
74 // double gamma;
75 //gnorm = Norm2(grad);
76 gnorm = sqrt(grad.dot(grad));
77
78 // Initialize xmax, typx and D to default values
79 xmax = -1.e30; typx = 1.0; D = 1.0;
80
81 for (i=0; i < nr; i++) xmax = max(xmax,fabs(xc(i)));
82 if( xmax != 0.0) typx = xmax;
83 if( gnorm!= 0.0) D = gnorm/typx;
84 if (debug_) {
85 *optout << "UpdateH: gnorm0 = " << gnorm
86 << "typx = " << typx << "\n";
87 }
88 for (i=0; i < nr; i++) Hessian(i,i) = D(i);
89 return Hessian;
90 }
91
92 SerialDenseVector<int,double> yk(nr), sk(nr), Bsk(nr);
93 SerialSymDenseMatrix<int,double> Htmp(nr), Htmp2(nr);
94
95 yk = grad;
96 yk -= gprev;
97 sk = xc;
98 sk -= xprev;
99
100 if (debug_) {
101 Print(yk);
102 Print(sk);
103 }
104
105 double gts =gprev.dot(sk);
106 double yts = yk.dot(sk);
107
108 double snorm = sqrt(sk.dot(sk));
109 double ynorm = sqrt(yk.dot(yk));
110
111 if (debug_) {
112 *optout << "UpdateH: gts = " << gts
113 << " yts = " << yts << "\n";
114 *optout << "UpdateH: snorm = " << snorm
115 << " ynorm = " << ynorm << "\n";
116 }
117
118 if (yts <= sqrteps*snorm*ynorm) {
119 if (debug_) {
120 *optout << "UpdateH: <y,s> = " << e(yts,12,4) << " is too small\n";
121 *optout << "UpdateH: The BFGS update is skipped\n";
122 }
123 Hessian = Hk; return Hk;
124 }
125
126 SerialDenseVector<int,double> res(nr);
127 //res = yk - Hk*sk;
128 SerialDenseVector<int,double> restmp(Hk.numRows());
129 restmp.multiply(Teuchos::LEFT_SIDE,1.0,Hk,sk,0.0);
130 res = yk;
131 res -= restmp;
132 if (res.normInf() <= sqrteps) {
133 if (debug_) {
134 *optout << "UpdateH: <y,s> = " << e(yts,12,4) << " is too small\n";
135 *optout << "UpdateH: The BFGS update is skipped\n";
136 }
137 Hessian = Hk; return Hk;
138 }
139
140 // Bsk = Hk*sk;
141 Bsk.multiply(Teuchos::LEFT_SIDE,1.0,Hk,sk,0.0);
142 double sBs = sk.dot(Bsk);
143 double etol = 1.e-8;
144
145 if (sBs <= etol*snorm*snorm) {
146 if (debug_) {
147 *optout << "UpdateH: <y,s> = " << e(yts,12,4) << " is too small\n";
148 *optout << "UpdateH: The BFGS update is skipped\n";
149 }
150 Hk = 0;
151 for (i=0; i < nr; i++) Hk(i,i) = sx(i)*sx(i);
152 Hessian = Hk; return Hk;
153 }
154
155 // Otherwise update the Hessian approximation
156 if (debug_) {
157 // *optout << "\nUpdateH: before update, k = " << k << "\n";
158 // FPrint(optout, Hk);
159 }
160
161 // //Htmp = - (Bsk * Bsk.t()) / sBs;
162 // for(int i=0; i<nr; i++)
163 // for(int j=0; j<=i; j++)
164 // {Htmp(i,j) = -Bsk(i)*Bsk(j)/sBs;}
165
166 // //Htmp = Htmp + (yk * yk.t()) / yts;
167 // for(int i=0; i<nr; i++)
168 // for(int j=0; j<nr; j++)
169 // {Htmp2(i,j) = yk(i)*yk(j)/yts;}
170 // Htmp += Htmp2;
171
172 // Htmp = Hk;
173 // Htmp += Htmp;
174 // Hk = Htmp;
175
176 SerialDenseMatrix<int,double> Htmp3(Htmp.numRows(),Htmp.numCols());
177 Htmp3.multiply(Teuchos::NO_TRANS,Teuchos::TRANS,-1/sBs,Bsk,Bsk,0.0);
178
179 for(int i=0;i<Htmp3.numRows();i++)
180 for(int j=0; j<=i;j++)
181 {Htmp(i,j) = Htmp3(i,j);}
182
183 // Htmp = Htmp + (yk * yk.t()) / yts;
184 Htmp3.multiply(Teuchos::NO_TRANS,Teuchos::TRANS,1/yts,yk,yk,0.0);
185 for(int i=0; i<nr; i++)
186 for(int j=0; j<=i; j++)
187 {Htmp2(i,j) = Htmp3(i,j);}
188
189 Htmp += Htmp2;
190
191
192 Htmp += Hk;
193
194
195 Hk = Htmp;
196 // Htmp.Release();
197 SerialDenseVector<int,double> Bgk(nr);
198 //Bgk = Hk*grad;
199 Bgk.multiply(Teuchos::LEFT_SIDE, 1.0, Hk, grad, 0.0);
200 double gBg = grad.dot(Bgk);
201 double gg = grad.dot(grad);
202 double ckp1= gBg/gg;
203 if (debug_) {
204 // *optout << "\nUpdateH: after update, k = " << k << "\n";
205 // FPrint(optout, Hk);
206 *optout << "UpdateH: sBs = " << sBs << "\n";
207 *optout << "UpdateH: ckp1 = " << ckp1 << "\n";
208 }
209 Hessian = Hk;
210 return Hk;
211 }
212
213 } // namespace OPTPP
214