Roger Hui wrote:
> 128!:0 underlies %. and %. rejects singular matrices.
>
> http://www.jsoftware.com/jwiki/Essays/QR_Decomposition
> presents a close model of 128!:0 . What modifications to
> the verb QR are needed so that it works on singular matrices?
I would not mess with %.
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.
9!:7 '+++++++++|-'
mp =:+/ .*
det=:-/ .*
NB. LAPACK version of QR reduction
require '~addons/math/lapack/lapack.ijs'
require '~addons/math/lapack/dgeqrf.ijs'
QRLA=:[: clean_jlapack_ dgeqrfQ_jlapack_
NB. QR from Essays
h =: +@|: NB. conjugate transpose
QR=: 3 : 0
n=.{:$A=.y
if. 1>:n do.
A ((% {.@,) ; ]) %:(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.
)
QRcheck=: 3 : 0
A=. y
assert. >:/$A
'Q R'=.z=. QR A
'm n'=. $A
assert. A -: Q mp R
assert. (($Q) -: m,n) *. 5e_14 > >./, | (=i.n) - (h Q) mp Q
assert. (($R) -: n,n) *. (0=R) >: >/~i.n
z
)
test=:3 : 0 NB. test assertion
Q=.y
n=.{. $ y
(=i.n) - (h Q) mp Q
)
A=:>1 2 3;4 5 6;5 7 9 NB. singular matrix
QRLA A
+----------------------------+--------------------------+
|_0.154303 0.801784 _0.57735|_6.48074 _8.79529 _11.1098|
|_0.617213 _0.534522 _0.57735| 0 0.801784 1.60357|
|_0.771517 0.267261 0.57735| 0 0 0|
+----------------------------+--------------------------+
QR A
+----------------------------+---------------------------+
|0.154303 0.801784 _0.806874|6.48074 8.79529 11.1098|
|0.617213 _0.534522 0.480691| 0 0.801784 1.60357|
|0.771517 0.267261 _0.343351| 0 0 2.5868e_14|
+----------------------------+---------------------------+
QRcheck A
|assertion failure: QRcheck
| (($Q)-:m,n)*.5e_14>>./,|(=i.n)-(h Q)mp Q
test >{. QRLA A
0 _3.20517e_17 4.87078e_17
_3.20517e_17 _4.44089e_16 1.37152e_17
4.87078e_17 1.37152e_17 0
test >{. QR A
1.11022e_16 _2.31678e_15 0.0927153
_2.31678e_15 0 0.995643
0.0927153 0.995643 _2.22045e_16
Best wishes,
John
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
