Discussion:
derive halley's method for vector valued functions
(too old to reply)
Peter Foelsche
2017-06-18 00:58:39 UTC
Permalink
Raw Message
Assume that f(x) and x itself is a n-dimensional vector.
If one derives the newton method, we are dealing with a 2-dim Jacobian.
Can one derive this for Halley’s method for vector valued functions?
I figured how to do this only for scalar functions:

n:2;
ratsimp(subst(f0, f(x), subst(f3, diff(f(x),x,3), subst(f2, diff(f(x),x,2), subst(f1, diff(f(x),x), n*''(diff(1/f(x),x,n-1))/''(diff(1/f(x),x,n)))))));

So – how to do this for x and f(x) being a vector?

Peter
Robert Dodier
2017-06-18 07:28:40 UTC
Permalink
Raw Message
Post by Peter Foelsche
Assume that f(x) and x itself is a n-dimensional vector.
If one derives the newton method, we are dealing with a 2-dim Jacobian.
Can one derive this for Halley’s method for vector valued functions?
n:2;
ratsimp(subst(f0, f(x), subst(f3, diff(f(x),x,3), subst(f2, diff(f(x),x,2), subst(f1, diff(f(x),x), n*''(diff(1/f(x),x,n-1))/''(diff(1/f(x),x,n)))))));
The only brief description of a multidimensional Halley method which I
found is this: http://folk.uib.no/ssu029/Pdf_file/Cuyt85.pdf
Do the formulas in Section 2 match the ones you are working with?

I just skimmed through the paper but it looks like it should be
possible with Maxima. A couple of elements which could be helpful: The
second derivative f'' can be represented, I believe, as a list or maybe
a matrix of one row or one column (dunno whether a list or matrix is
preferable at this point) of matrices, each of which is the Hessian of
one of the component functions of f. The function hessian computes the
Hessian of a multivariate function.

The function "at" is a substitution function which is somewhat
specialized for representing derivatives evaluated at a point -- I guess
this is what you are trying to achive with subst above? Just guessing on
this point.

By default, operations on lists are carried out element by element --
this happens to be the definition needed for the Banach algebra on R^n
as mentioned in the paper. See the option variable listarith.

There is an outline of an approach to solve the problem given in Section
2 of the paper -- maybe that could guide a Maxima implementation.

I'll let that be enough for now -- let us know how it goes, it seems
like a good problem.

best,

Robert Dodier
Robert Dodier
2017-06-19 05:26:32 UTC
Permalink
Raw Message
I've spent a pleasant day working on this problem and got some results
which seem to be consistent with the Cuyt & Rall paper
(http://folk.uib.no/ssu029/Pdf_file/Cuyt85.pdf). I've attached a
script, halley_method.mac, and halley_method.log, the output that I
get from running batch("halley_method.mac") in Maxima.

The final results appear to be the same for the 2 example problems in
the paper, but the per-iteration values are a little bit different. I
don't know if this is single precision vs double, or inaccuracies
induced by converting floats to rationals and back again (in Maxima),
or what.

The basic method is pretty simple, but the whole thing is obscured
somewhat by a layer of dross -- stuff like needing to convert a column
matrix to a list or vice versa and forcing values to be floats. I
don't know what could be done to make it clearer.

The iteration could be more efficient -- the LU decomposition of the
Jacobian is computed repeatedly. I didn't bother trying to save and
reuse any partial results but that's an obvious thing to do.

I'd be happy to hear questions or comments.

best,

Robert Dodier
Peter Foelsche
2017-06-19 12:29:45 UTC
Permalink
Raw Message
Robert,

Thanks for your efforts!!

Peter

Sent from an universe in which Biff Tannen kept the sport almanach

From: Robert Dodier
Sent: Monday, June 19, 2017 0:26
To: Peter Foelsche
Cc: <maxima-***@lists.sourceforge.net>
Subject: Re: derive halley's method for vector valued functions

I've spent a pleasant day working on this problem and got some results
which seem to be consistent with the Cuyt & Rall paper
(http://folk.uib.no/ssu029/Pdf_file/Cuyt85.pdf). I've attached a
script, halley_method.mac, and halley_method.log, the output that I
get from running batch("halley_method.mac") in Maxima.

The final results appear to be the same for the 2 example problems in
the paper, but the per-iteration values are a little bit different. I
don't know if this is single precision vs double, or inaccuracies
induced by converting floats to rationals and back again (in Maxima),
or what.

The basic method is pretty simple, but the whole thing is obscured
somewhat by a layer of dross -- stuff like needing to convert a column
matrix to a list or vice versa and forcing values to be floats. I
don't know what could be done to make it clearer.

The iteration could be more efficient -- the LU decomposition of the
Jacobian is computed repeatedly. I didn't bother trying to save and
reuse any partial results but that's an obvious thing to do.

I'd be happy to hear questions or comments.

best,

Robert Dodier

Loading...