Raul Miller wrote:
> On 8/27/07, John Randall <[EMAIL PROTECTED]> wrote:
>> QR in the Essay appears to almost work on singular matrices: it looks as
>> though a sign is wrong somewhere, making the orthogonality test fail.
>> Here is a comparison with the results from LAPACK.
>
> Here's a brute-force "fix" of the signs, to match LAPACK:
>
> NB. QR from Essays
> h =: +@|: NB. conjugate transpose
>
> QR=: 3 : 0
> n=.{:$A=.y
> if. 1>:n do.
> A (([EMAIL PROTECTED] {.@,) ; [EMAIL PROTECTED]) %:(h A) mp A
> else.
> m =.>.n%2
> A0=.m{."1 A
> A1=.m}."1 A
> 'Q0 R0'=.QR A0
> 'Q1 R1'=.QR A1 - Q0 mp T=.(h Q0) mp A1
> (-Q0,.-Q1);-(R0,.T),-(-n){."1 R1
> end.
> )
>
> But there's another difference -- the third column of Q from
> LAPACK looks very different from the third column of Q
> from QR -- LAPACK may even be noticing the singular
> state of the matrix and using a different algorithm for
> this case.
In the example we have been considering, this gives the correct value of
R, but not Q. R has to be upper triangular (which it is), but Q has to be
orthogonal: (h Q) mp Q should be the identity. In this case, since
everything is real, h is just |: .
'Q R'=:QRLA A
(h Q) mp Q
1 3.20517e_17 _4.87078e_17
3.20517e_17 1 _1.37152e_17
_4.87078e_17 _1.37152e_17 1
'Q R'=:QR A NB. modified QR from essay
(h Q) mp Q
1 _2.31678e_15 _0.0927153
_2.31678e_15 1 0.995643
_0.0927153 0.995643 1
In the QR case, there are still nonzero off-diagonal terms.
The method described in the Essay is Gram-Schmidt orthogonalization, and
this may be the problem: this requires linearly independent columns, which
obviously a singular matrix does not have. Another method is Householder
transformations, which I believe does not have this problem. For cutting
down on the number of iterations, it is common to first transform the
matrix to upper Hessenberg form, but this should not affect the outcome.
For testing, the QR factorization of a nonsingular matrix is unique if the
diagonal entries of R are positive. It is not unique for a singular
matrix.
So
QR >1 2; 3 4
+------------------+-----------------+
|0.316228 _0.948683|3.16228 4.42719|
|0.948683 0.316228| 0 _0.632456|
+------------------+-----------------+
QRLA >1 2; 3 4
+-------------------+------------------+
|_0.316228 _0.948683|_3.16228 _4.42719|
|_0.948683 0.316228| 0 _0.632456|
+-------------------+------------------+
are really giving the same result.
Best wishes,
John
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
