Dear Eric Iverson, Thank you for trying to put improving %. on the agenda. I am only a complex analyst not a numerical analyst. As far as I know, Givens rotations make a cache friendlier approach than Householder reflections. I will try to find some good quality C code for QR via Givens on the web or in a book, but it may take me some time. I am not sure how J works on the inside, and I am also not a professional programmer. I will try to find at least a starting point for QR.
Yours sincerely, Imre Patyi On Mon, May 3, 2021 at 4:25 PM Eric Iverson <[email protected]> wrote: > 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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
