operatorsjnewtons-methodtacit-programming

J: Tacit adverb of Newton's method


I've found in 'addons/math/misc/brent.ijs' implementation of Brent's method as an adverb. I would like to build a Newton's method as an adverb too but it's much harder than building tacit verbs.

Here is a explicit version of Newton's iteration:

   newton_i =: 1 : '] - u % u d.1'

With such usage:

   2&o. newton_i^:_ (1) NB. (-: 1p1) must be found
1.5708
   2 o. 1.5708 NB. after substitution we get almost 0
_3.67321e_6

And of course, for convenience:

    newton =: 1 : 'u newton_i^:_'

What's a tacit equivalent?


Solution

  • TL;DR

    Per the comments, a short answer; the tacit equivalent to the original, explicit newton_i and newton are, respectively:

    n_i =: d.0 1 (%/@:) (]`-`) (`:6) 
    newton =: n_i (^:_)
    

    Some techniques for how such translations are obtained, in general, can be found on the J Forums.

    Construction

    The key insights here are that (a) that a function is identical to its own "zeroeth derivative", and that (b) we can calculate the "zeroeth" and first derivative of a function in J simultaneously, thanks to the language's array-oriented nature. The rest is mere stamp-collecting.

    In an ideal world, given a function f, we'd like to produce a verb train like (] - f % f d. 1). The problem is that tacit adverbial programming in J constrains us to produce a verb which mentions the input function (f) once and only once.

    So, instead, we use a sneaky trick: we calculate two derivatives of f at the same time: the "zeroth" derivative (which is an identity function) and the first derivative.

       load 'trig'
       sin              NB. Sine function (special case of the "circle functions", o.)
    1&o.
    
       sin d. 1 f.      NB. First derivative of sine, sin'.
    2&o.
    
       sin d. 0 f.      NB. "Zeroeth" derivative of sine, i.e. sine.
    1&o."0
    
       sin d. 0 1 f.    NB.  Both, resulting in two outputs.
    (1&o. , 2&o.)"0
    
       znfd =: d. 0 1   NB. Packaged up as a re-usable name.
       sin znfd f.
    (1&o. , 2&o.)"0
    

    Then we simply insert a division between them:

       dh =: znfd (%/@) NB. Quotient of first-derivative over 0th-derivattive
    
       sin dh
    %/@(sin d.0 1)
    
       sin dh f.
    %/@((1&o. , 2&o.)"0)
    
       sin dh 1p1  NB. 0
    _1.22465e_16
    
       sin 1p1     NB. sin(pi) = 0
    1.22465e_16
       sin d. 1 ] 1p1  NB. sin'(pi) = -1
    _1
       sin dh 1p1  NB. sin(pi)/sin'(pi) = 0/-1 = 0
    _1.22465e_16
    

    The (%/@) comes to the right of the znfd because tacit adverbial programming in J is LIFO (i.e. left-to-right, where as "normal" J is right-to-left).

    Stamp collecting

    As I said, the remaining code is mere stamp collecting, using the standard tools to construct a verb-train which subtracts this quotient from the original input:

       ssub  =: (]`-`) (`:6)     NB. x - f(x)
    
       +: ssub                   NB. x - double(x)
    ] - +:
       -: ssub                   NB. x - halve(x)
    ] - -:
    
       -: ssub 16                NB. 16 - halve(16)
    8
       +: ssub 16                NB. 16 - double(16)
    _16
       *: ssub 16                NB. 16 - square(16)
    _240
       %: ssub 16                NB. 16 - sqrt(16)
    12
    

    Thus:

        n_i =: znfd ssub         NB. x - f'(x)/f(x)
    

    And, finally, using "apply until fixed point" feature of ^:_, we have:

        newton =: n_i (^:_)
    

    Voila.