One might also want to trade Householder for Givens,
at least in certain circumstances. Would be nice to
• automatically decide based on the matrix’s properties
• have a way to set flags before using %.
Am 03.05.21 um 22:25 schrieb Henry Rich:
> Can anyone make a recommendation of one (or maybe two) freely-available
> codes that I could use in place of the current QR code? LAPACK would be
> OK. I don't know the field & I don't have the time right now to look
> deeply into it.
>
> I have used Householder matrices for finding eigenvalues, & I think they
> would be cache-unfriendly for QR decomposition. Something more
> block-oriented would be preferable. But as I said I don't know the field.
>
> Henry Rich
>
> On 5/3/2021 3:56 PM, Imre Patyi 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
>
>
--
----------------------
mail written using NEO
neo-layout.org
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
