Thanks. Perhaps, this is the kind of QR routine that is implemented inside J in C.
It shows some numerical instability. Here is an example. All the best, Imre Patyi NB. From Raul Miller, "[email protected]". QR1=: 3 : 0 n=.{:$A=.y if. 1>:n do. A ((% {.@,) ; ]) %:(+@|: A) +/ .* A else. m =.>.n%2 A0=.m{."1 A A1=.m}."1 A 'Q0 R0'=.QR1 A0 'Q1 R1'=.QR1 A1 - Q0 +/ .* T=.(+@|: Q0) +/ .* A1 (Q0,.Q1);(R0,.T),(-n){."1 R1 end. ) NB. We can use the above QR routine to demonstrate a loss of NB. orthogonality, which is a form of numerical instability. NB. Here is a matrix of primes. A=:p:i. 32 32 'q r'=:QR1 A NB. The matrix q is not as orthogonal as we might wish in double precision. >./|,((|:q)(+/ .*)q)-(=/~i.32) NB. 4.12269e_9 >./|,(q(+/ .*)|:q)-(=/~i.32) NB. 1.48817e_9 NB. Here is the same thing in Octave with the built-in command qr(). NB. A=reshape(list_primes(32*32),32,32) ; NB. [q,r]=qr(A) ; NB. norm(q'*q-eye(size(q)),"inf") NB. ans = 3.5686e-15 NB. norm(q*q'-eye(size(q)),"inf") NB. ans = 4.0289e-15 On Sun, May 2, 2021 at 12:54 PM Raul Miller <[email protected]> wrote: > Here's something interesting. > > When I compare the 128!:0 implementation of QR with a straightforward > J implementation of QR, I get slightly different results. > > QR0=: 128!:0 > > QR1=: 3 : 0 > n=.{:$A=.y > if. 1>:n do. > A ((% {.@,) ; ]) %:(+@|: A) +/ .* A > else. > m =.>.n%2 > A0=.m{."1 A > A1=.m}."1 A > 'Q0 R0'=.QR1 A0 > 'Q1 R1'=.QR1 A1 - Q0 +/ .* T=.(+@|: Q0) +/ .* A1 > (Q0,.Q1);(R0,.T),(-n){."1 R1 > end. > ) > > > In various J90x version (I tested J901, J902 and J903 beta): > > (QR0 -&> QR1) p:i.3 3 > 0 _2.22045e_16 2.33147e_15 > 5.55112e_17 _8.88178e_16 9.74221e_15 > 1.11022e_16 _2.16493e_15 _6.96318e_15 > > 0 3.55271e_15 7.10543e_15 > 0 1.33227e_15 _1.5099e_14 > 0 0 _2.22045e_16 > > x=:_0.5+(i.1000)%999 > c=:(^x) %. X=:(^ 0j1 * x */ i:10) > (+/%#)*:|(X (+/ .*) c)-(^x) > 6.57613e_5 > > Compare J807: > > (QR0 -&> QR1) p:i.3 3 > 0 0 0 > 0 0 0 > 0 0 0 > > 0 0 0 > 0 0 0 > 0 0 0 > > x=:_0.5+(i.1000)%999 > c=:(^x) %. X=:(^ 0j1 * x */ i:10) > (+/%#)*:|(X (+/ .*) c)-(^x) > 0.0279006 > > I imagine that this has to do with floating point precision > limitations, though I have not looked deeply enough to characterize > the issues here. That said, when I look at the current J source, I see > that it is using 'long double' (which has a 64 bit mantissa): > > jsource$ find . -iname '*.h' | xargs grep -l 'long double' | xargs > grep -c 'long double' > ./jsrc/jtype.h:1 > ./sleef/include/sleef.h:4 > ./sleef/src/arch/helpervecext.h:15 > ./sleef/src/common/misc.h:4 > > So that's probably part of it. > > FYI, > > -- > Raul > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
