[Jprogramming] Non-scalar Newton-Raphson

[Louis de Forcrand]

Thanks everyone!

The lab was pretty much exactly what I was looking for, and indeed
matrix division was the way to go. For satisfactory results, ~1000
iterations are needed in some cases, but they are spectacular:



The code is probably minimal, but it’s surprisingly versatile; the code used to generate the curves is:

newt=: 2 : '(- (0{::n) * u %. u D.1)^:(1{::n)’ NB. fetch is used here to be able to box the 2nd parameter
								  NB. for a list of iterations
chi2=: 1 : '[: +/@:*: {:@[ - ] u {.@[‘
points=: (i: 5) ,: ? 10 $ 10
green=: points&(p. chi2) D.1 newt 0.5 1000 ] 1 1 1 1 1
red=: points&(p. chi2) D.1 newt 0.5 1000 ] 1 1 1 1 1 1 1 1 1 1

The optimised parameters were then stored in green or red, and I simply
plotted the two functions on top of the points with the plot package.

Theoretically (ignoring possible precision problems) any function which takes parameters on its left
and a value to its right can be used as a placeholder to minimise the square of the differences with
the points to fit; with only two operators!
I really am glad I decided to look into APL and then J; you won’t see me doing this any time soon in C.

Thanks again for your generous help,
Louis


> On 14 Feb 2016, at 15:14, Ben Gorte - CITG <B.G.H.Gorte at tudelft.nl> wrote:
> 
> 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
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm