The errors I found are listed starting at
http://code.jsoftware.com/wiki/System/Interpreter/Bugs/Errors#Higher_derivatives_inaccurate_unless_symbolic
The bottom line is that you'd better test out D. on your verbs to make
sure it works.
Some years ago Roger said (as I recall it) that in retrospect d. D. etc
should not have been primitives; that a library written in J would have
been more maintainable. Since they operate on verbs to produce verbs,
any overhead in the conversion would be incurred only once, rather than
on each application of the verb; and with the table of derivatives
outside the interpreter, it could be extended easily.
That would still be a good package for someone to write. I don't think
that the bugs in D. will be, or should be, fixed.
Henry Rich
On 2/13/2016 7:36 PM, Louis de Forcrand wrote:
http://code.jsoftware.com/wiki/Vocabulary/dcapdot
<http://code.jsoftware.com/wiki/Vocabulary/dcapdot>
Towards the bottom it says that u D.n with n > 1 is less accurate
than u D.1 D.1 .
I am assuming that u D.1 D.1 isn’t very accurate either though,
as your example is correct (2nd der. of +/@:*: is 2: ).
———
Indeed, multinewt does work (better than other options in fact)
for +/@:*: .
Even though vsnewt seems to be the best so far, as you said
it won’t work if the function returns a vector with the same shape
as its argument (it doubles its shape at each iteration). Making it
work would involve running some axis together (taking diagonals)
to eliminate partial derivatives which result from differentiating
with respect to two different elements of the argument.
I can’t figure out which ones though…
Louis
On 14 Feb 2016, at 01:13, Jose Mario Quintana <[email protected]>
wrote:
Consider the following,
(+/@:*:D.1 D.1) 0 0
1 0
0 1
(+/@:*:D.1 D.1) 0 111
0 0
0 1.9861
(+/@:*:D.2 ) 0 0
2 0
0 2
(+/@:*:D.2 ) 0 111
0 _304.804
1.00391 1.99305
As far as I know, the second derivative of the sum of the squares is, in
principle, two times the identity matrix regardless of the arguments. Am I
missing something?
On Sat, Feb 13, 2016 at 4:09 PM, Marshall Lochbaum <[email protected]>
wrote:
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
----------------------------------------------------------------------
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
