Re: Nonlinear CG
"Ben FrantzDale" <[email protected]> Fri, 19 Jan 2007 14:38:00 -0500
| Newsgroups | gmane.comp.lib.mtl.devel |
|---|---|
| Message-ID | <[email protected]> |
--===============1355516514==
Content-Type: multipart/alternative;
boundary="----=_Part_144585_19807046.1169235480921"
------=_Part_144585_19807046.1169235480921
Content-Type: text/plain; charset=WINDOWS-1252; format=flowed
Content-Transfer-Encoding: quoted-printable
Content-Disposition: inline
Exactly.
Defaults could be as shown below.
The next question is a generic interface for line search. Generally line
search needs a function, f, a starting point, x, and a search direction, d.
Some line searches need f(x), which the caller often already knows, so that
should probably be passed in as well.
The CG code I have in front of me can switch between the Secant method and =
a
simpler incremental search. Both of those essentially take f, x, d, f(x),
the minimum search distance, the maximum search distance, and return the
number of function evaluations and whether or not a root was found. The
distances are in terms of the infinity norm (i.e., the maximum magnitude of
a vector component). I'm not sure if the CG algorithm could know those
values in a generic way. The code I'm looking for has them as constants and
so they could be contained in a line-search function object.
That makes it sound like a line searcher would look like this:
struct line_search {
void operator()(Fn& f,
AffinePoint& x, // x gets moved to the result position.
Direction const& d,
Scalar& f_of_x); // f_of_x is modified to be the resultin=
g
f(x).
};
In the case of Newton-Raphson, some iteration control is needed. That could
probably have a sane default that could be set by the caller of
nonlinear_cg. That is,
// Make a N-R functor with a limit of ten steps
// and a search range of 2e-4 to 0.3.
Newton_Raphson NR(10, make_pair(2e-4, 0.3));
nonlinear_cg(f, x, P, iter, NR, FletcherReeves());
That seems pretty natural and pretty generic.
=97Ben
Defaults:
// Note: x need only be an affine point, not a proper vector. (See Wikipedi=
a
for more on affine points.)
template <typename Fn, typename AffinePoint, typename Scalar>
void apply(Fn const& f, AffinePoint const& x, Scalar& result) {
result =3D f(x);
}
template <typename GradableFn, typename AffinePoint, typename Scalar,
typename GradVec>
void apply_grad(GradableFn const& f, AffinePoint const& x, Scalar& result,
GradVec& resultVec) {
apply(f, x, result);
apply_grad(f, x, resultVec);
}
template <typename HessableFn, typename AffinePoint, typename Scalar,
typename GradVec, typename HessianMat>
void apply_hessian(GradableFn const& f, AffinePoint const& x, Scalar&
result, GradVec& resultVec, HessianMat& resultHess) {
apply_grad(f, x, result, resultVec);
apply_hessian(f, x, resultHess);
}
On 1/18/07, Peter Gottschling <[email protected]> wrote:
>
>
> On 17.01.2007, at 10:24, Ben FrantzDale wrote:
>
> > On 1/16/07, Ben FrantzDale <[email protected]> wrote:
> >> Some methods require only f', others require both f and f', others
> >> require f' and f''. Any thoughts as to do that generically? One
> >> possibility would be for the operator() of f to take two, three, or
> >> four arguments. That is,
> >> f(x, d) -> d =3D f(x)
> >> f(x, d, g) -> d =3D f(x), g =3D f'(x)
> >> f(x, d, g, H) -> d =3D f(x), g =3D f'(x), H =3D f''(x).
> >>
> >> Other options include making derivative(f)(x, g) return g, but that
> >> seems a bit too clever for its own good.
> >
> >
> > Another option is to do something more like itl::mult, which has these
> > semantics:
> > itl::mult(A, x, y, z); // z =3D y + A * x
> > itl::mult(A, x, y); // y =3D A * x
> > That is something like
> > apply(f, x, scalar); // scalar =3D f(x); that could be the default
> > implementation.
> > apply_grad(f, x, scalar, grad); // scalar =3D f(x), grad =3D f'(x)
> > apply_hessian(f, x, scalar, grad, hessian); // scalar =3D f(x), grad
> > =3D f'(x), hessian =3D f''(x)
> > For some algorithms, it will be faster to compute f and its
> > derivatives in one pass, so I think the above functions may be more
> > appropriate than separate apply(f, x, scalar), apply_grad(f, x, grad),
> > and apply_hessian(f, x, hessian). On the other hand, it might make
> > sense, e.g., for apply_grad to be overloaded to be both apply_grad(f,
> > x, scalar, grad) and apply_grad(f, x, grad).
>
> I think that's a good idea. At the very least apply_grad(f, x, scalar,
> grad) can be implemented by default calling apply_grad(f, x, grad) and
> apply(f, x, scalar). In some cases (i.e. for some f) you might
> specialize it then to a more efficient computation.
>
> Peter
> >
> > I'll continue to play around with these ideas.
> >
> > =97Ben
> > _______________________________________________
> > This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/
> ------------
> Peter Gottschling
> Research Associate
> Open Systems Laboratory
> Indiana University
> 135 Lindley Hall
> Bloomington, IN 47405
> Tel.: +1-812-855-3608 Fax: +1-812-856-0853
> http://www.osl.iu.edu/~pgottsch
>
>
> _______________________________________________
> This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/
>
>
------=_Part_144585_19807046.1169235480921
Content-Type: text/html; charset=WINDOWS-1252
Content-Transfer-Encoding: quoted-printable
Content-Disposition: inline
Exactly.<br>Defaults could be as shown below. <br><br>The next question is =
a generic interface for line search. Generally line search needs a function=
, f, a starting point, x, and a search direction, d. Some line searches nee=
d f(x), which the caller often already knows, so that should probably be pa=
ssed in as well.
<br><br>The CG code I have in front of me can switch between the Secant met=
hod and a simpler incremental search. Both of those essentially take f, x, =
d, f(x), the minimum search distance, the maximum search distance, and retu=
rn the number of function evaluations and whether or not a root was found. =
The distances are in terms of the infinity norm (
i.e., the maximum magnitude of a vector component). I'm not sure if the=
CG algorithm could know those values in a generic way. The code I'm lo=
oking for has them as constants and so they could be contained in a line-se=
arch function object.
<br><br>That makes it sound like a line searcher would look like this:<br>s=
truct line_search {<br> void operator()(Fn& f,<br> &nb=
sp; =
AffinePoint& x, // x gets moved to the result position.<br=
> &n=
bsp; Direction const& d,
<br>  =
; Scalar& f_of_x); // f_of_x is modified =
to be the resulting f(x).<br>};<br><br>In the case of Newton-Raphson, some =
iteration control is needed. That could probably have a sane default that c=
ould be set by the caller of nonlinear_cg. That is,
<br> // Make a N-R functor with a limit of ten steps<br> // and=
a search range of 2e-4 to 0.3.<br>
Newton_Raphson NR(10, make_pair(2e-4, 0.3));<br> nonlinear_cg(=
f, x, P, iter, NR, FletcherReeves());<br> <br>That seems pretty natura=
l and pretty generic.<br><br>=97Ben<br><br><br>Defaults:<br><br>// Note: x =
need only be an affine point, not a proper vector. (See Wikipedia for more =
on affine points.)
<br>template <typename Fn, typename AffinePoint, typename Scalar><br>=
void apply(Fn const& f, AffinePoint const& x, Scalar& result) {=
<br> result =3D f(x);<br>}<br><br><span></span>template <typename =
GradableFn, typename AffinePoint, typename Scalar, typename GradVec>
<br>void apply_grad(GradableFn const& f, AffinePoint const& x, Scal=
ar& result, GradVec& resultVec) {<br> apply(f, x, resul=
t);<br> apply_grad(f, x, resultVec);<br>}<br><br><span></span>t=
emplate <typename HessableFn, typename AffinePoint, typename Scalar, typ=
ename GradVec, typename HessianMat>
<br>
void apply_hessian(GradableFn const& f, AffinePoint const& x, Scala=
r& result, GradVec& resultVec, HessianMat& resultHess) {<br>&nb=
sp; apply_grad(f, x, result, resultVec);<br> apply_hessia=
n(f, x, resultHess);
<br>
}<br><br><br><div><span class=3D"gmail_quote">On 1/18/07, <b class=3D"gmail=
_sendername">Peter Gottschling</b> <<a href=3D"mailto:[email protected]=
u">[email protected]</a>> wrote:</span><blockquote class=3D"gmail_quot=
e" style=3D"border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt =
0.8ex; padding-left: 1ex;">
<br>On 17.01.2007, at 10:24, Ben FrantzDale wrote:<br><br>> On 1/16/07, =
Ben FrantzDale <<a href=3D"mailto:[email protected]">bfrantz@scorec=
.rpi.edu</a>> wrote:<br>>> Some methods require only f', other=
s require both f and f', others
<br>>> require f' and f''. Any thoughts as to do that gen=
erically? One<br>>> possibility would be for the operator() of f to t=
ake two, three, or<br>>> four arguments. That is,<br>>> f(x, d)=
-> d =3D f(x)
<br>>> f(x, d, g) -> d =3D f(x), g =3D f'(x)<br>>> f(x, =
d, g, H) -> d =3D f(x), g =3D f'(x), H =3D f''(x).<br>>&g=
t;<br>>> Other options include making derivative(f)(x, g) return g, b=
ut that
<br>>> seems a bit too clever for its own good.<br>><br>><br>&g=
t; Another option is to do something more like itl::mult, which has these<b=
r>> semantics:<br>> itl::mult(A, x, y, z); // z =3D y + A * x<br>>=
itl::mult(A, x, y); // y =3D A * x
<br>> That is something like<br>> apply(f, x, scalar); // scalar =3D =
f(x); that could be the default<br>> implementation.<br>> apply_grad(=
f, x, scalar, grad); // scalar =3D f(x), grad =3D f'(x)<br>> apply_h=
essian(f, x, scalar, grad, hessian); // scalar =3D f(x), grad
<br>> =3D f'(x), hessian =3D f''(x)<br>> For some algorit=
hms, it will be faster to compute f and its<br>> derivatives in one pass=
, so I think the above functions may be more<br>> appropriate than separ=
ate apply(f, x, scalar), apply_grad(f, x, grad),
<br>> and apply_hessian(f, x, hessian). On the other hand, it might make=
<br>> sense, e.g., for apply_grad to be overloaded to be both apply_grad=
(f,<br>> x, scalar, grad) and apply_grad(f, x, grad).<br><br>I think tha=
t's a good idea. At the very least apply_grad(f, x, scalar,
<br>grad) can be implemented by default calling apply_grad(f, x, grad) and<=
br>apply(f, x, scalar). In some cases (i.e. for some f) you migh=
t<br>specialize it then to a more efficient computation.<br><br>Peter<br>&g=
t;<br>
> I'll continue to play around with these ideas.<br>><br>> =97=
Ben<br>> _______________________________________________<br>> This li=
st is archived at <a href=3D"http://www.osl.iu.edu/MailArchives/mtl-devel/"=
>
http://www.osl.iu.edu/MailArchives/mtl-devel/</a><br>------------<br>Peter =
Gottschling<br>Research Associate<br>Open Systems Laboratory<br>Indiana Uni=
versity<br>135 Lindley Hall<br>Bloomington, IN 47405<br>Tel.: +1-812-855-36=
08 Fax: +1-812-856-0853
<br><a href=3D"http://www.osl.iu.edu/~pgottsch">http://www.osl.iu.edu/~pgot=
tsch</a><br><br><br>_______________________________________________<br>This=
list is archived at <a href=3D"http://www.osl.iu.edu/MailArchives/mtl-deve=
l/">
http://www.osl.iu.edu/MailArchives/mtl-devel/</a><br><br></blockquote></div=
><br>
------=_Part_144585_19807046.1169235480921--
--===============1355516514==
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline
_______________________________________________
This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/
--===============1355516514==--