RFC:QR decomposition in julia

[andreasnoack]

This pull request started as a QR decomposition in julia such that e.g. Matrix{BigFloat} can be decomposed. However, while adding the tests I also decided to go through more combinations of types which revealed some problems and I have therefore changed a lot of the promotion rules. This request should therefore also fix JuliaLang/LinearAlgebra.jl#54. I decided that the ! versions shouldn't promote but assume that it has called with a reasonable type.

The option for returning the condition number in sqrtm was broken and undocumented so I decided to remove it. There was also a few bugs in leading dimension calculations in lapack.jl and more the factorization and special matrix types have been changed from type to immutable.

I am aware that this pr includes the balance option to eigfact, but I expect that the balance pull request will be merged very soon.

[ViralBShah]

I am not sure if this gets hoisted by the compiler, but you may want to have Qfactors = Q.factors outside the loops.

[andreasnoack]

@Viral It made a significant difference in speed to do this.

[ViralBShah]

This is amazing stuff, but also quite a big PR, making it a bit difficult to review.

@jiahao should take a look as well.

[jiahao]

This looks correct upon a first reading. However, I'm AFC until later in the afternoon.

[andreasnoack]

@ViralBShah Yes it got too big. I didn't anticipate the consequences of the changes to the tests.

I forgot that I have changed the names of the QR types. There are three now: the usual unpivoted (QR), pivoted (QRPivoted) and the WY version (without pivoting) (QRWY). I renamed the old QR to QRWY and added a traditional QR for the julia implementation. All versions return a least squares solution when m>n and the minimum norm solution when m<n, but only QRPivoted estimates the rank.

I should also add that there is a lot of rank one updating written out in loops in the QR code. When we have fast array views, it should be possible to remove many code lines.

[mlubin]

Not directly related, but will this make it easy to implement functionality like in qrupdate?

[ViralBShah]

For fast qrupdate, we need to implement a block version of the algorithm for good performance, which is not what you necessarily need for BigFloat.

[andreasnoack]

I have rebased. @jiahao During the rebase I dedollared de6c1f6. I think the line count is about the same and the code maybe a little easier to read.

[andreasnoack]

@mlubin I don't know what it will take to implement updating. I tried to look at your link, but I think the only documentation is the source code.

[jiahao]

Good riddance.

[andreasnoack]

@jiahao I have updated the pull request with your suggestions so far.

[jiahao]

Thanks. I guess I have more bedtime reading tonight.

[andreasnoack]

One more update of the code. The direct test against LAPACK was useful and caught a conjugation error.

[andreasnoack]

@jiahao In order to get the same result for the julia LAPACK LU I had to switch to the BLAS style square root free absolute value abs(real(z))+abs(imag(z)), but after thinking about it I'd prefer the more generic abs(z). The two methods give different ordering and therefore permutations. What do you think?

[jiahao]

Just want to point out that your observation can be rephrased (with slight mathematical abuse) as LAPACK using the l1-norm over (real(z), imag(z)) and you would prefer the l2-norm. I don't really have a strong feeling about one over the other. In fact, you could even choose to specify the norm as a parameter.

[jiahao]

Since factorize! is no longer exported, that probably deserves a note in NEWS.md

[jiahao]

I've read through the core QR algorithms twice now and they seem correct. So if no one else has any objections, I think you can declare victory as you see fit.

[jiahao]

👯