Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

RFC:QR decomposition in julia #5526

Merged
merged 1 commit into from
Jan 31, 2014
Merged

RFC:QR decomposition in julia #5526

merged 1 commit into from
Jan 31, 2014

Conversation

andreasnoack
Copy link
Member

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 #5178. 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.

for i = 1:mA
νAi = A[i,k]
for j = k+1:mQ
νAi += Q.factors[j,k]*A[i,j]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@ViralBShah
Copy link
Member

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
Copy link
Member

jiahao commented Jan 25, 2014

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

@andreasnoack
Copy link
Member Author

@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
Copy link
Member

mlubin commented Jan 25, 2014

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

@ViralBShah
Copy link
Member

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
Copy link
Member Author

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
Copy link
Member Author

@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.

@@ -607,7 +607,6 @@ export
eigvecs,
expm,
eye,
factorize!,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good riddance.

@andreasnoack
Copy link
Member Author

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

@jiahao
Copy link
Member

jiahao commented Jan 27, 2014

Thanks. I guess I have more bedtime reading tonight.

@andreasnoack
Copy link
Member Author

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

@andreasnoack
Copy link
Member Author

@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
Copy link
Member

jiahao commented Jan 29, 2014

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
Copy link
Member

jiahao commented Jan 30, 2014

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

@jiahao
Copy link
Member

jiahao commented Jan 31, 2014

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.

…y reflectors. Add Float16 and BigFloat to tests and test promotion. Cleaup promotion in LinAlg. Avoid promotion in mutating ! functions. Make Symmetric, Hermitian and QRs immutable. Make thin a keyword argument in svdfact. Remove cond argument in sqrtm.
andreasnoack added a commit that referenced this pull request Jan 31, 2014
RFC:QR decomposition in julia
@andreasnoack andreasnoack merged commit 8a94301 into master Jan 31, 2014
@andreasnoack andreasnoack deleted the anj/qr branch January 31, 2014 10:30
@jiahao
Copy link
Member

jiahao commented Jan 31, 2014

👯

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Stack overflows in factorization
5 participants