Been following this conversation with interest. I actually went through this exercise a long time ago (I think 7.x, could have been 6.x). I was building a framework for Kernel regression.
There were problems then, though I think the present situation is actually considerably better. The solution I came up with at the time to mitigate the issue, was lapack dgels, which I believe is what R uses. I came to the conclusion there were issues using the NIST standards for linear regression: https://www.itl.nist.gov/div898/strd/lls/lls.shtml Stuck this up if it helps with R comparisons (if it gets past grouchy Brian Ripley, it's probably fine): https://github.com/locklin/j-linreg-nist/ FWIIW I think these results are pretty reasonable. You can always produce an example that will screw up QR decomposition. If you're actually doing linear regression on data from nature, rather than coming up with pathological examples, you're presumably going to be checking for multicollinearity and using a different technique to deal with such examples. Comparing to other APLs, Kona uses/used an SVD approach as I recall which has problems with these NIST data sets. Kx uses Cholesky, which is arguably the wrong approach for linear regression; I don't remember if I checked those and my K is too stale to check. Kick me a PR if you do it; Art's pretty smart so maybe he knows more than me about linear regression. As for what technique to use: the correct answer for the general case is generally QR decomposition. J's result is decent; probably a few more iterations would give exactly the answer in lapack's dgels (there may be some exception handling thing in dgels; https://www.netlib.org/lapack/explore-html/d7/d3b/group__double_g_esolve_ga225c8efde208eaf246882df48e590eac.html if you want to dig into it). There are numerous other numeric linear algebra techniques which produce effectively the same answer. Using the same technique, or at least getting comparable results as R is pretty safe. -SL On Wed, May 5, 2021 at 5:17 AM Imre Patyi <[email protected]> wrote: > 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 > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
