Cholseski update of rank one

Open discussion regarding features, bugs, issues, vendors, etc.

Cholseski update of rank one

Postby nice44 » Wed Aug 17, 2011 3:31 pm

Since they are still not in Lapack (or I'm blind?), and since I had a need for it, I have a made port of Choleski update/downdate from the original LINPACK dchud and dchdd routines.
http://www.netlib.org/linpack/dchud.f
http://www.netlib.org/linpack/dchdd.f
I followed the suggestions and guidelines of Matthias Seeger from:
Low Rank Updates for the Cholesky Decomposition - Matthias Seeger
http://infoscience.epfl.ch/record/16146 ... update.pdf
(I noted that he previouslyposted on the subject in 2007)
The main changes are :
- uses of BLAS drot and drotg for Given's rotations
- forces positive elements on the diagonal of R
- accepts upper or lower Cholesky factor R
I also changed sqrt(1-norm**2) into a sin(acos(norm)) or cos(asin(norm) so as to be more stable if norm is close to 1.
This might be questionable because there is no use of inverse trigonometric functions inthe whole Lapack, but if you have a better way...

I named the routines dchur1 and dchdr1 for choleski update/downdate of rank 1, but of course, this can be changed at will...
You can find the code attached, it is released under MIT license.
Tell me if it is of any interest.

Cheers

Nicolas Cellier
Attachments
dchdr1.txt
Choleski downdate of rank one
(9.63 KiB) Downloaded 68 times
dchur1.txt
Choleski update of rank one
(7.92 KiB) Downloaded 155 times
nice44
 
Posts: 3
Joined: Wed Aug 17, 2011 3:00 pm

Re: Cholseski update of rank one

Postby admin » Wed Aug 24, 2011 5:50 pm

Thank you Nicolas,
As it is, it is not possible to integrate this code directly into LAPACK.
- The routine needs a better description, name ,...
- We need the double complex version of the code (we have generator to get the single precision and the complex one)
- We need to have some TESTING that goes with it
Also, if we decide to integrate the routine, we may want to look for a code that use a "blocking" technic. But it would be fine to start with yours.

Nicolas, if you feel like finishing the polishing, we will be more than happy to count you as one of the several LAPACK contributors.
If you want more information about LAPACK Contribution, please take a look at http://www.netlib.org/lapack/#_strong_c ... ors_strong and especially the LAPACK Program Style document (Sorry, it has not been updated for a while).
Let us know if you have any questions
Julie
admin
Site Admin
 
Posts: 486
Joined: Wed Dec 08, 2004 7:07 pm

Re: Cholseski update of rank one

Postby sven » Thu Aug 25, 2011 5:13 am

Dear Nicolas,

Just a small point, but sqrt(1-norm**2) can be computed stably as sqrt((1+norm)*(1-norm)).

I agree with Julie; updating by rank k instead of rank one would allow the use of Householder transformations and block partitioning and hence the use of Level 3 BLAS. But of course we appreciate that this would be much more work. Daniel Kressner has software in this spirit for updating the QR factorization at

http://www.math.ethz.ch/~kressner/qrupdate.php

Best wishes,

Sven Hammarling.
sven
 
Posts: 144
Joined: Wed Dec 22, 2004 4:28 am

Re: Cholseski update of rank one

Postby nice44 » Thu Aug 25, 2011 1:17 pm

Julie, yes the names would have to be changed.
I have some tests for the real/double cases but haven't any complex version available.
Also, I didn't start thinking about blocking.
I needed a rank-1 update in the real case, the version from Seeger was not available, so I decided to share this little piece.
My added value was minimal. I might have a deeper look if there is any interest, but I don't promise anything.

Sven, thanks for the sqrt(1-n**2) tips, that's dead simple!
And yes, of course, updates of rank n can be performed with QR (or its LQ transpose in case of lower triangular factors I think).
Do you know how to handle a rank n "downdate"?

Nicolas
nice44
 
Posts: 3
Joined: Wed Aug 17, 2011 3:00 pm

Re: Cholseski update of rank one

Postby sven » Fri Aug 26, 2011 7:29 am

Dear Nicolas,

The downdating problem needs care. Some ideas, in the context of QR are given in:

[url]eprints.ma.man.ac.uk/1192/[/url]

The best reference for much of the background is probably

A Bjorck. Numerical Methods for Least Squares Problems. SIAM, 1996.

Best wishes,

Sven.
sven
 
Posts: 144
Joined: Wed Dec 22, 2004 4:28 am

Re: Cholseski update of rank one

Postby nice44 » Fri Aug 26, 2011 3:43 pm

Sven, thanks, this is a much more advanced analysis!
What are you plans about developed code?
nice44
 
Posts: 3
Joined: Wed Aug 17, 2011 3:00 pm


Return to User Discussion

Who is online

Users browsing this forum: Yahoo [Bot] and 1 guest

cron