Hi all,
I have to say that I am all in for implementations of common algebra algorithms in uBlas for the following reasons:
* Newcomers really need to get their projects prototyped fast, having implementations in place gives them all the heads up they need
* There are projects (think for example a small opengl 3d visualization library) that relying on algebra algorithms based on bindings from the classical algebra libraries is not really worth the logistics trouble. In small projects or libraries I also find myself preferring making custom implementations of a few algorithms (inversions, linear solvers), even if they are not optimal, or even if they have some development cost.
* The numeric bindings is a very good library and whoever needs that extra something, or is in the research sector can always take advantage of it. I use it mainly in isolated research projects, rather than deploying them in full programs that I have to share with others (I don't suggest that someone doesn't want to do this though, it is a personal taste relating to the logistics and maintenance issues relating with such an approach, especially in cross-platform projects). The bindings way is better for larger projects/teams.
Relating to the nice comment by Matwey:
0. uBlas style is good for algorithms that are efficiently implemented with expression templates. Others (take for example matrix multiplication) do not benefit and I would suggest to cut the bait and go for a classical functional approach.
2. We may want to implement the algorithms in an stl kind of way. i.e. defining them through iterators that work pretty much like the bindings. This has the benefit of (maybe) giving a common interface and also differentiate algorithms from containers (a nice vision of the c++ leaders and I think of some of the ppl in here). This approach though may create a very nasty algebraic interface and may not be appropriate after all.
3. It looks like we will soon have the ublas::tie implementations (I am not sure if the name is ideal but we can start suggesting somethings). This would give some extremely nice interfaces, that can be chosen by the compiler for us. For example:
tie(U, S, V) = svd(A);
//and
tie(S) = svd(A); //or just S=svd(A);
I am hopeful that not only this is very convenient but can be implemented in such a way that the returned objects actually define which algorithm is chosen (up to an ambiguity level)! Furthermore this would be an one of its kind extension of C++ and a great showcase. Thanx Jesse for thinking this, and I hope you see the possibilities it opens!
To recapitulate my view is that we should provide custom implementations, open the development as possible so that people can contribute (maybe make an easy to access branch for all), and most importantly first define the interface it should have. The interface choices that come into my mind are the following (I will take SVD as an example):
1. Procedural: svd(A,U,S,V);
2. Functional - procedural: S = svd(A,U,V);
3. tieable: tie(U, S, V) = svd(A); tie(S) = svd(A); //or just S=svd(A);
4. Classical/or stl like: SVD(A.begin<1>(), A.end<1>(), S.begin<1>(), U.begin<1>(), V.begin<1>()) (ugly uh? fairly generic though).
Overall I would vote for (3), but we need first to implement the "tie" interface.
Thanx for reading all of this!
Nasos
> To: ublas@lists.boost.org
> From: matwey.kornilov@gmail.com
> Date: Sat, 24 Jul 2010 16:39:40 +0400
> Subject: Re: [ublas] todo list
>
>
> As far as I understand when we talk about implementation of such algorithms
> in uBLAS we mean that API should be 'ublas-style'. For instance, when we
> write inner_prod( v, u ) we get an object of this operation which looks like
> a scalar for us. The question is what is the 'ublas-style' for SVD,
> Cholesky, LU, QR, etc. algorithms? It is not clear for me.
>
> Most of the algorithms are in-place. For instance, we have to alter our
> matrix in order to get its SVD decomposition. The next example. Cholesky
> decomposition is defined only for positive-defined matrices. When matrix is
> not positive-defined there is a square root of negative number in the
> algorithm. The good C++ style is to throw exception here and to leave the
> matrix in initial state (i.e not altered). But known Cholesky implementation
> assumes undefined matrix at output in such cases.
>
> For instance the result of SVD decomposition is a triple of matrices: one
> singular diagonal S matrix and two unitary matrices U and V. Imagine we want
> to know only singular values to calculate conditioning number of matrix or
> it's rank. What happens when we don't want to know U and V matrices? As
> uBLAS user I expect that decomposition saves calculations and memory for
> placing unitary matrices. What happens when we want to use SVD for matrix
> inversion? In this case we have to know U and V and thus they should be
> calculated and allocated.
>
> The third question is the estimations of rounding errors.
>
> I suppose there are a several dozens of questions left.
>
>
> I tried to implement SVD and Cholesky for my purposes. The code is
> available:
> http://curl.sai.msu.ru/~matwey/lsp
> Subversion:
> svn://curl.sai.msu.ru/lsp
> But IMO the API definitely is not the good one.
>
>
> Marco Guazzone wrote:
>
> > # SVD, Cholesky, LU, QR, Shur
> > Do you think this should be written from scratch or maybe using
> > existent solvers like the one provided by LAPACK (possibly using the
> > boost::numeric::bindings lib)?
>
>
>
> _______________________________________________
> ublas mailing list
> ublas@lists.boost.org
> http://lists.boost.org/mailman/listinfo.cgi/ublas
> Sent to: nasos_i@hotmail.com
The New Busy is not the too busy. Combine all your e-mail accounts with Hotmail. Get busy.