Thanks for the info. I agree with you about the necessity of having
methods beyond %. available for handling data.
If I can find a code that I can drop in to replace %., I will use it.
LAPACK seems to have a couple of candidates.
Henry Rich
On 5/6/2021 8:15 AM, Scott Locklin wrote:
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
--
This email has been checked for viruses by AVG.
https://www.avg.com
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
