If u is (+/@:*:), then (u D. 1) returns a vector, not a matrix. From
your first email it looked like you wanted to find a zero of
(+/@:*: D. 1), that is, minimize (+/@:*:). multinewt works for that
case:
+/@:*: D. 1 multinewt (<5) 3 4
3 4
_0.142552 _0.0529989
0.00158427 0.00128811
_1.42249e_5 2.08734e_5
_2.57043e_7 _1.37868e_7
There's a more complicated formula if the function's range has fewer
dimensions than its domain. For a scalar function of a vector like
(+/@:*:) it is
v_{n+1} = v_n - f(v_n) * f'(v_n) / |f'(v_n)|^2
or
vsnewt =: 2 : '(- u * [: (% +/@:*:) u D.1)^:n'
+/@:*: vsnewt (<5) 3 4
3 4
1.5 2
0.75 1
0.375 0.5
0.1875 0.25
This converges slowly because the root has order two.
Marshall
On Sat, Feb 13, 2016 at 09:32:15PM +0100, Louis de Forcrand wrote:
> I’m sorry about that Raul, that should read “y”.
>
> ----
>
> Thanks Roger. You were right, after a quick search, I’ve found that
> indeed I would need to perform a matrix division. Several documents
> give this as a multi-variable Newton method:
>
> v_n+1 = v_n - J_f (v_n) ^-1 f(v_n)
>
> In J, this would be:
>
> multinewt=: 2 : ‘(- u %. u D.1)^:n’
>
> However, if I take a paraboloid z =: +/@:*: x y , and look for its zeros with
>
> +/@:*: multinewt (<5) 3 4
>
> I get this:
>
> 3 4
> _0.5 0.5
> _0.5 0.5
> _0.5 0.5
> _0.5 0.5
>
> which is obviously wrong. Would I need some sort of transposition? This
> could just be due to precision errors. Or because I wrote an erroneous
> conjunction?
>
> Cheers,
> Louis
>
> > On 13 Feb 2016, at 20:52, Raul Miller <[email protected]> wrote:
> >
> > Something is wrong here.
> >
> > You say "where A is the original argument, or the result of the
> > previous iteration."
> >
> > ... but I can find no use of A in the expression you are describing.
> >
> > Can you look over what you wrote and decide if it is what you meant to say?
> >
> > Thanks,
> >
> > --
> > Raul
> >
> > On Sat, Feb 13, 2016 at 10:50 AM, Louis de Forcrand <[email protected]>
> > wrote:
> >> 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
>
> ----------------------------------------------------------------------
> For information about J forums see http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
