kaleidicassociates / lubeck

High level linear algebra library for Dlang
http://lubeck.libmir.org/
Boost Software License 1.0
66 stars 14 forks source link

Helper Function for QRResult to Calculate Q/R #17

Closed jmh530 closed 5 years ago

jmh530 commented 5 years ago

I think QRResult needs a helper function(s) that produces Q and R from matrix and tau.

I am comparing the result of one of the unittests to how I would do it in R and it doesn't match up. Compare https://run.dlang.io/is/fcMfi9 with below (in R). The lubeck version has a zero in the top left corner for matrix, but I thought it was upper trapezoidal...I know that given R, we can calculate Q as a*inv(R), but I'm not sure I would get the correct R given's lubeck's output.

> A <- matrix(c(1, 1, 0, 1, 0, 1, 0, 1, 1), 3, 3)
> X <- qr(A)
> X
$qr
           [,1]       [,2]       [,3]
[1,] -1.4142136 -0.7071068 -0.7071068
[2,]  0.7071068  1.2247449  0.4082483
[3,]  0.0000000 -0.8164966  1.1547005

$rank
[1] 3

$qraux
[1] 1.707107 1.577350 1.154701

$pivot
[1] 1 2 3

attr(,"class")
[1] "qr"
> qr.Q(X)
           [,1]       [,2]       [,3]
[1,] -0.7071068  0.4082483 -0.5773503
[2,] -0.7071068 -0.4082483  0.5773503
[3,]  0.0000000  0.8164966  0.5773503
> qr.R(X)
          [,1]       [,2]       [,3]
[1,] -1.414214 -0.7071068 -0.7071068
[2,]  0.000000  1.2247449  0.4082483
[3,]  0.000000  0.0000000  1.1547005
> 
jmh530 commented 5 years ago

It looks like orgqr/ungqr are the lapack functions you use to get Q. R can be calculated by matrix multiplication (just by transpose, rather than inverse I have above).

jmh530 commented 5 years ago

@9il I've got something working (tested at least for square matrices). The only issue right now is organization. I started off using a helper function, like extractQ (it turns out that R is actually done simply without using orgqr).

However, I am partial to putting this functionality into QRResult and then using a flag to switch between outputting the raw matrices (originally) or the actual Q/R matrices. This might reduce allocations versus having them as separate functions.

Any opinion?

9il commented 5 years ago

Sorry for the huge delay with the answer. I think additional methods that computes Q or R in the current QRResult should be fine.

jmh530 commented 5 years ago

Another delay on my part. Just submitted pull request .