Nothing attached came through. You probably need to make your attachments be .txt files, to avoid spam filters.
Thanks, -- Raul 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
