[Jprogramming] Non-scalar Newton-Raphson

[Ben Gorte - CITG]

Hello Louis, 

did you have a look at the Best Fit lab in the math section of Help -> Studio -> Labs ?
I think towards the end it is getting pretty close to what you are looking for.

Ben

________________________________________
From: Programming [programming-bounces at forums.jsoftware.com] on behalf of Louis de Forcrand [olius at bluewin.ch]
Sent: Saturday, February 13, 2016 16:50
To: programming at jsoftware.com
Subject: [Jprogramming] Non-scalar Newton-Raphson

I’ve been trying to write a conjunction that will find the zeros to a function using the
Newton-Raphson method. The simplest way to do this is probably:

English:
x_n+1 = x_n - f ( x_n ) / f ‘ ( x_n )

J:
eznewt=: 2 : ‘ ( - u % u D.1 ) ^: n ‘

This works fine for scalar -> scalar functions, but won’t work if the rank of the
result is AND the rank of the arguments are above 1.
Probably the most evident situation where this would be a problem is if one
were searching for a minimum instead of a zero, in which case the algorithm
would be applied to the derivative of the function: +/@:*: D.1 newt 30 ] 3 4
As I understand it, this would give a result shape at each iteration (with a function u) of:

( $y ) , $ u eznewt 1 y

where A is the original argument, or the result of the previous iteration. What we would
want is a result with shape $y .

First, let’s get some rules clear:
- the syntax should be: (verb) newt (# of iterations) (argument)
- the argument can be of any rank
- the shape of the argument matches the shape of the result at any iteration
- this implies that $ u y matches i. 0 or $ y

The hard part is getting u % u D.1 to have the shape of the argument.
If y is a vector, and u y is a scalar, then u D.1 y will be a vector, and
u D.1 D.1 y will be a $y by $y matrix. All that is needed then is to take
its diagonal with (<0 1)&|: .

But what if y is a matrix? Since $ u D.1 y -: ( $ y ) , $ u y , I though maybe running
y's axis together with |: might work

newt=: 2 : '(- u diag_and_divide u D.1)^:n’
diag_and_divide=: [ % ] |:~ -@#@$@[ <@}. i.@#@$@]

 but something about it doesn’t.
My head is $pinning now, and I figured I’d send this and then take a break.

Thanks in advance!
Louis
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm