Re: Matrix decompositions.
Nasos Iliopoulos <[email protected]>
| Newsgroups | gmane.comp.lib.boost.ublas |
|---|---|
| Message-ID | <[email protected]> |
Yes, additionally with the boost transition to github we are planning to make contributing to uBlas much easier. We would like though to stress quality and find a way to define performance requirements specifications. When the github repo (a fork of boostorg/ublas) is ready (I hope I can have it ready this weekend) we will include instructions on how to work with pull requests (to submit an algorithm or a small bugfix for example), or in cases where larger involvement is required maybe give access to the ublas development repo to ppl that want to contribute. -Nasos On 05/28/2013 10:59 AM, Salman Javaid wrote: > Are there any plans to integrate Karl's QR Decomposition > implementation into Boost.uBLAS? I will be grateful for a response. > > > > Best Regards, > Salman Javaid > > > On Thu, May 16, 2013 at 7:06 AM, Karl Rupp <[email protected] > <mailto:[email protected]>> wrote: > > Hi guys, > > I have a working implementation of QR for uBLAS in ViennaCL for > about a year already: > > https://github.com/viennacl/viennacl-dev/blob/master/viennacl/linalg/qr.hpp > > Feel free to copy&paste and relicense as needed, I agree to > whatever is necessary to integrate into uBLAS if of interest. > > There is a bit of duplication for the ViennaCL types (requiring > host<->GPU transfers) and the uBLAS types (operating in main RAM), > yet it gives an idea how things can be implemented. It can > certainly be improved here and there (more details on request), > yet it addresses most of the points raised by Oswin. And it's > faster than a standard LAPACK for sizes above ~1k times 1k. > > I recommend extracting the Householder reflections into a nice > separate interface, since this functionality will also be needed > for other algorithms like SVD or GMRES. As a nice side effect, it > makes the implementation for QR more compact. > > Generally, in order to get *any* reasonable performance, one > really needs to use matrix-matrix multiplications where possible > in order to avoid the memory bandwidth bottleneck. > > Best regards, > Karli > > > > On 05/16/2013 01:06 AM, oswin krause wrote: > > Hi, > > These are further good points! > > I also came up with a few new ones(and tips!): > > - QR needs pivoting. It's main usages are SVD and > pseudo-inverses. In > both cases the input does not necessary have full rank. Also > pivoting > helps for matrices with high condition numbers. > > - For the same reasons H was not formed explicitly, Q should > not be > formed. Instead there should be a version of the algorithm > which does > only return the reflection vectors forming Q. > > - For dense matrices at least, it is possible to do the QR > in-place by > storing the R part as lower triangular and the householder > transformation vectors in the upper triangular. (That is very > similar to > how LU is implemented). > > - The best sources for algorithmic information are the LAPACK > working notes. > http://www.netlib.org/lapack/lawns/ > > In your case lawn114 sems to be the most relevant, even though it > assumes a fast BLAS3. > > Greetings, > Oswin > > > On 16.05.2013 05:32, Nasos Iliopoulos wrote: > > That's not a bad start. > > I think Oswin covered a whole lot of items here, but a > more complete > algorithm needs to satisfy some or all of the following: > > - The algortihm should have a dispatch mechanism so that > optimized > versions for various matrix types can be provided. > (sparse, banded, > etc.). You don't necessarily need to provide them all to > start with. > - The algorithm should operate on matrix expressions > rather than > matrices (so it can be applied to for example subranges). > Static > dispatch or overload if for some reason this seems to > reduce performance. > - Const correctness is important. Try using const > reference on > immutable types. > - Instead of 0.00025 provide a value based on a user > choice.If it is > hard coded by the user, the compiler will probably > convert it into a > const value. > - Don't use ints for indexing, use either std::size_t, or > container::size_type. If you need a signed type (i.e. to > count for > differences on unsigned types) use ptrdiff_t. uBlas > containers provide > a difference_type typedef for that purpose (i.e. > matrix<double>::difference_type). > - use noalias(..) = in every assignment that the lhs is > not a part of > rhs, or when the algebraic operation is mapped 1-1. (i.e. > A=2*A+B can > be written as noalias(A)=2*A+B, but A=prod(B,A)+D cannot > atm). This > provides more benefits than just avoiding temporaries. > > > The QR decomposition of a 100x100 matrix should take no > more than a > few miliseconds (or even less than a milisecond) to run. > A 1000x1000 should take around 1/3 to 1/10 of a sec. > > Compile with: > g++ -O3 -DNDEBUG Main.cpp -o qrtest > > Then you'll see that your code runs pretty fast, but it > doesn't scale > well as Oswin noted. > > Best regards, > -Nasos > > > > > On 05/15/2013 10:12 PM, Salman Javaid wrote: > > Thank you, Oswin for the detailed response. I am going > to update the > code. > > David, Nasos, any advise on coding conventions? Or > anything else that > you can possible suggest? I will stand grateful. > > > > > > Best Regards, > Salman Javaid > > > On Tue, May 14, 2013 at 10:53 PM, oswin krause > <[email protected] > <mailto:[email protected]> > <mailto:[email protected] > <mailto:[email protected]>>> wrote: > > Hi, > > in the order I stumbled over the things: > > main.cpp > line 44-54: you don't need a copy, instead you > should use a > combination of row/subrange. > line 58-60: you should take a look at inner_prod > line 63: 0.00025 is too big. > line 66: You should never create H explicitly. > line 67: because you formed H, this step is O(n^3) > which makes > the whole algorithm O(n^4). This can be done in O(n^2) > line 73-79: same applies here. > > Greetings, > Oswin > > On 14.05.2013 22:12, Salman Javaid wrote: > > Hello uBLAS Contributors: > > I have > applied to GSoC 2013 > and pitched implementation of SVD > factorization for uBLAS. In > order to better prepare myself and to get my > hands dirty at > uBLAS, I ended up implementing QR > Factorization employing > Householder Reflections using uBLAS. This is > only the first > draft and will be needing significant > improvement, e.g., > computation of QR decomposition of 100 * 100 > matrix takes around > 30 seconds. But I guess just to get familiar > with code base, it > was a good exercise. Over the next week or two > I will be trying > to optimize the code. > > > I will be > absolutely > grateful if contributors can have a quick > glance at the code, > and point me to any improvements they can > suggest. Particularly > in speeding up matrix multiplication. > > I used Visual Studio 2010 to compile the code. > I will try to get > the code running on my Ubuntu machine in a > couple of days hopefully. > > Here the header file: > https://github.com/salmanjavaid/QR_Decomposition/blob/master/QR_Header.hpp > > > The main file: > https://github.com/salmanjavaid/QR_Decomposition/blob/master/Main.cpp > > > Best Regards, > Salman Javaid > > > _______________________________________________ > ublas mailing list > [email protected] > <mailto:[email protected]> > <mailto:[email protected] > <mailto:[email protected]>> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to:[email protected] > <mailto:to%[email protected]> > <mailto:[email protected] > <mailto:[email protected]>> > > > > _______________________________________________ > ublas mailing list > [email protected] <mailto:[email protected]> > <mailto:[email protected] > <mailto:[email protected]>> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: [email protected] > <mailto:[email protected]> > <mailto:[email protected] > <mailto:[email protected]>> > > > > > > _______________________________________________ > ublas mailing list > [email protected] <mailto:[email protected]> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to:[email protected] > <mailto:to%[email protected]> > > > > > _______________________________________________ > ublas mailing list > [email protected] <mailto:[email protected]> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to:[email protected] > <mailto:to%[email protected]> > > > > > _______________________________________________ > ublas mailing list > [email protected] <mailto:[email protected]> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: [email protected] <mailto:[email protected]> > > > _______________________________________________ > ublas mailing list > [email protected] <mailto:[email protected]> > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: [email protected] <mailto:[email protected]> > > > > > _______________________________________________ > ublas mailing list > [email protected] > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: [email protected] _______________________________________________ ublas mailing list [email protected] http://lists.boost.org/mailman/listinfo.cgi/ublas Sent to: [email protected]