Your request for an improvement to %. is well stated and makes sense. We will take a serious look at this and will try to get it in as part of the 903 beta and for inclusion in the final 903 release.
As you have already done much of the leg work, are you ready to state a strong preference for which code or example J should adopt. On Mon, May 3, 2021 at 3:56 PM Imre Patyi <[email protected]> wrote: > Thank you. > > I have played with implementing a QR algorithm in J, > called qrhh for QR via Householder reflections. > Please find it attached along with some tests. > > The given QR algorithm qrhh solves all the 5 test cases attached > while the built in %. solves none of them with comparable or even > acceptable accuracy. > > I am more and more convinced that the built-in %. is numerically very > unstable. It would be better to find a reputable QR implementation > that is numerically stable. I am sure there are nice implementations > of QR in C freely available that are numerically stable. Any > mathematically correct > QR algorithm can solve small examples, but to solve examples of > useful size, we need to have a numerically stable algorithm. Most > plain vanilla math algorithms have very little chance of being numerically > stable (i.e., they magnify the roundoff errors rather than keeping them > more or less constant), it is not for nothing that numerical analysts > exist. > > The rosetta code site has a C implementation with Householder reflections. > I have not tried testing it, but it is likely to beat J's built-in %. by > several > orders of magnitude of accuracy if the implementation is correct. > The recursive algorithm given also on the same page of rosetta code > for J (and mentioned in this thread by Raul Miller) is based upon a clever > rewriting of the classical Gram--Schmidt orthogonalization process, > which is well-known to be numerically unstable and every book on > numerical analysis recommends against using it. While a numerically > unstable algorithm (of which 128!:0 is probably an example) can be accurate > on some examples, it often fails on similar examples for no clear reason. > See Example 3 in the attachment, where > > NB. 'q r'=: 128!:0 X > NB. >./|,((|:q)(+/ .*)q)-(=/~i.1{$q) > NB. 0.999904 > > 128!:0 produces an "orthogonal" matrix Q with Q'Q-I having an entry > nearly 1 (all should be nearly zero). Here X is a 1000 by 21 real matrix, > with not too large or too small numbers. When %. solves a linear > least squares problem Ax=b, or QRx=b, it relies on Q'Q=I in order to solve > Rx=Q'b > by backward substitution. If Q'Q is not nearly I, it ( %. ) is likely to > fail badly. > > https://rosettacode.org/wiki/QR_decomposition#C > > https://code.jsoftware.com/wiki/Essays/QR_Decomposition > > https://en.wikipedia.org/wiki/QR_decomposition > > Book, "Numerical recipes in C". > > I would like to request that the J maintainers replace the current > %. with a numerically stable one. QR decomposition is a well > entrenched small program, and it would likely be easy enough to > find a suitable numerically stable replacement. It would cost you > little to fix it. > > Thank you very much. > > Yours sincerely, > Imre Patyi > > > On Sun, May 2, 2021 at 11:37 PM Henry Rich <[email protected]> wrote: > > > J does not use long double. > > > > It is possible that 8.07 used 8087 instructions with 80-bit > > floating-point; I don't remember. > > > > The QR algorithm in JE is a clever recursive algorithm that Roger > > invented. It uses +/ . * to do almost all the calculations; reciprocals > > are performed only at the bottom of the recursion, on 2x2 or 1x1 > matrices. > > > > At the lower levels of recursion it works with short, wide matrices that > > the blocking matrix-multiply code can't grip efficiently, so I added a > > version of +/ . * that does matrix multiplication via inner products two > > at a time: just the thing for such cases. The overall algorithm is > > unchanged. > > > > The differences you report look pretty small to me, and could be > > artifacts of differing orders of summation. > > > > I have seen references to what Python does, but I don't know what its > > qr() code is defined to do. If the matrix is ill-conditioned, perhaps > > some extra processing is performed. > > > > Henry Rich > > > > > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
