1 //
2 // dfp.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Edward Seidl <seidl@janed.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #include <math.h>
29
30 #include <util/state/stateio.h>
31 #include <math/optimize/update.h>
32 #include <math/optimize/transform.h>
33 #include <util/keyval/keyval.h>
34
35 using namespace sc;
36
37 /////////////////////////////////////////////////////////////////////////
38 // DFPUpdate
39
40 static ClassDesc DFPUpdate_cd(
41 typeid(DFPUpdate),"DFPUpdate",1,"public HessianUpdate",
42 create<DFPUpdate>, create<DFPUpdate>, create<DFPUpdate>);
43
DFPUpdate()44 DFPUpdate::DFPUpdate()
45 {
46 }
47
DFPUpdate(const Ref<KeyVal> & keyval)48 DFPUpdate::DFPUpdate(const Ref<KeyVal>&keyval):
49 HessianUpdate(keyval)
50 {
51 if (keyval->exists("xprev") && keyval->exists("gprev")) {
52 Ref<SCMatrixKit> k = SCMatrixKit::default_matrixkit();
53 RefSCDimension dim = new SCDimension(keyval->count("xprev"));
54 xprev = k->vector(dim);
55 gprev = k->vector(dim);
56 for (int i=0; i<dim.n(); i++) {
57 xprev(i) = keyval->doublevalue("xprev",i);
58 gprev(i) = keyval->doublevalue("gprev",i);
59 }
60 }
61 }
62
DFPUpdate(StateIn & s)63 DFPUpdate::DFPUpdate(StateIn&s):
64 SavableState(s),
65 HessianUpdate(s)
66 {
67 Ref<SCMatrixKit> k = SCMatrixKit::default_matrixkit();
68 RefSCDimension dim;
69 dim << SavableState::restore_state(s);
70 xprev = k->vector(dim);
71 gprev = k->vector(dim);
72 xprev.restore(s);
73 gprev.restore(s);
74 }
75
~DFPUpdate()76 DFPUpdate::~DFPUpdate()
77 {
78 }
79
80 void
save_data_state(StateOut & s)81 DFPUpdate::save_data_state(StateOut&s)
82 {
83 HessianUpdate::save_data_state(s);
84 SavableState::save_state(xprev.dim().pointer(),s);
85 xprev.save(s);
86 gprev.save(s);
87 }
88
89 void
update(const RefSymmSCMatrix & ihessian,const Ref<Function> & func,const RefSCVector & xn,const RefSCVector & gn)90 DFPUpdate::update(const RefSymmSCMatrix&ihessian,const Ref<Function>&func,
91 const RefSCVector&xn,const RefSCVector&gn)
92 {
93 RefSCVector xnew, gnew;
94
95 // the update for the inverse hessian differs from the update for the
96 // hessian in that xdisp and gdisp are exchanged
97 if (inverse_hessian_) {
98 xnew = xn;
99 gnew = gn;
100 } else {
101 xnew = gn;
102 gnew = xn;
103 }
104
105 if (xprev.nonnull()) {
106 RefSCVector xdisp = xnew-xprev;
107 RefSCVector gdisp = gnew-gprev;
108 RefSCVector ihessian_gdisp = ihessian * gdisp;
109 double gdisp_ihessian_gdisp = ihessian_gdisp.scalar_product(gdisp);
110 double xdisp_gdisp = xdisp.scalar_product(gdisp);
111 ihessian.accumulate(
112 xdisp.symmetric_outer_product()*(1.0/xdisp_gdisp)
113 - ihessian_gdisp.symmetric_outer_product()*(1.0/gdisp_ihessian_gdisp)
114 );
115 xprev.assign(xnew);
116 gprev.assign(gnew);
117 } else {
118 xprev = xnew.copy();
119 gprev = gnew.copy();
120 }
121 }
122
123 void
apply_transform(const Ref<NonlinearTransform> & trans)124 DFPUpdate::apply_transform(const Ref<NonlinearTransform>& trans)
125 {
126 if (trans.null()) return;
127 HessianUpdate::apply_transform(trans);
128 trans->transform_coordinates(xprev);
129 trans->transform_gradient(gprev);
130 }
131
132 void
set_inverse(void)133 DFPUpdate::set_inverse(void)
134 {
135 HessianUpdate::set_inverse();
136 RefSCVector tmp;
137 tmp = xprev;
138 xprev = gprev;
139 gprev = tmp;
140 }
141
142 /////////////////////////////////////////////////////////////////////////
143 // BFGSUpdate
144
145 static ClassDesc BFGSUpdate_cd(
146 typeid(BFGSUpdate),"BFGSUpdate",1,"public HessianUpdate",
147 create<BFGSUpdate>, create<BFGSUpdate>, create<BFGSUpdate>);
148
BFGSUpdate()149 BFGSUpdate::BFGSUpdate()
150 {
151 }
152
BFGSUpdate(const Ref<KeyVal> & keyval)153 BFGSUpdate::BFGSUpdate(const Ref<KeyVal>&keyval):
154 DFPUpdate(keyval)
155 {
156 }
157
BFGSUpdate(StateIn & s)158 BFGSUpdate::BFGSUpdate(StateIn&s):
159 SavableState(s),
160 DFPUpdate(s)
161 {
162 }
163
~BFGSUpdate()164 BFGSUpdate::~BFGSUpdate()
165 {
166 }
167
168 void
save_data_state(StateOut & s)169 BFGSUpdate::save_data_state(StateOut&s)
170 {
171 DFPUpdate::save_data_state(s);
172 }
173
174 void
update(const RefSymmSCMatrix & ihessian,const Ref<Function> & func,const RefSCVector & xn,const RefSCVector & gn)175 BFGSUpdate::update(const RefSymmSCMatrix&ihessian,const Ref<Function>&func,
176 const RefSCVector&xn,const RefSCVector&gn)
177 {
178 RefSCVector xnew, gnew;
179
180 // the update for the inverse hessian differs from the update for the
181 // hessian in that xdisp and gdisp are exchanged
182 if (inverse_hessian_) {
183 xnew = xn;
184 gnew = gn;
185 } else {
186 xnew = gn;
187 gnew = xn;
188 }
189
190 if (xprev.nonnull()) {
191 RefSCVector xdisp = xnew-xprev;
192 RefSCVector gdisp = gnew-gprev;
193 RefSCVector ihessian_gdisp = ihessian * gdisp;
194 double gdisp_ihessian_gdisp = ihessian_gdisp.scalar_product(gdisp);
195 double xdisp_gdisp = xdisp.scalar_product(gdisp);
196 RefSCVector u = xdisp*(1.0/xdisp_gdisp)
197 - ihessian_gdisp*(1.0/gdisp_ihessian_gdisp);
198 ihessian.accumulate(
199 // DFP part
200 xdisp.symmetric_outer_product()*(1.0/xdisp_gdisp)
201 - ihessian_gdisp.symmetric_outer_product()*(1.0/gdisp_ihessian_gdisp)
202 // BFGS part
203 + u.symmetric_outer_product() * gdisp_ihessian_gdisp
204 );
205 xprev.assign(xnew);
206 gprev.assign(gnew);
207 } else {
208 xprev = xnew.copy();
209 gprev = gnew.copy();
210 }
211 }
212
213 /////////////////////////////////////////////////////////////////////////////
214
215 // Local Variables:
216 // mode: c++
217 // c-file-style: "ETS"
218 // End:
219