[IPOL discuss] the problem with Eigen (and other embedded libraries)

[eml]

> CCmath (linux only)

ccmath is not linux only.  It is a set of indepentent functions
written in ANSI C.  You do not need to "install" this library.   You
copy-paste the code of the desired functions into your program.


> So I think there is not easy solution.

I think the problem is that there are too many easy solutions (for
computing the SVD, not for the general problem of linear algebra
libraries).

Eigen is a very beautiful library and probably it produces the most
efficient programs due to its clever use of expression templates.  If
one of our algorithms requires a lot of matrix computations, then
Eigen is most surely the best tool for the job.

However, if you only need it to compute the eigenvalues of a square
symmetric positive semidefinite 3x3 matrix, including the whole 2MB of
Eigen just for that is a ridiculous overkill.  The solution of this
problem has an explicit formula (using lots of square and cubic
roots), and can be also computed robustly by any naive iterative
algorithm, in few lines of portable C.  For this size of problem, the
algorithm chosen is mostly irrelevant.


I attach 3 solutions.

Solution 1 is "reallysmall_svd.c", which contains a signle function
"svd" written in ANSI C which computes the decomposition for square
symmetric positive semidefinite matrices of any size.  It is copied
almost verbatim from a book of numeric algorithms in pascal, and I
don't know about the resulting licence (if such a small piece of code
can carry a licence).

Solution 2 is "eig3.c", which contains code to compute the
eigenvectors of a symmetric 3x3 matrix.  According to its author it is
"public domain", whatever that means.

Solution 3 is "ccmath_svd.c", with the code copied from ccmath that
computes this decomposition for an arbitrary rectangular matrix
(together with an example program).  Since it is copied from ccmath,
the licence is LGPL.


These functions can be compiled without any external libraries
(besides the stantard C library, which is needed for the square
roots).  I tried their results with some random matrices and the
results for each method are exactly the same as those given by octave.



> It could be a long debate, but it's important that anyone could provide it's point of view.

My point of view is that the simplest tools should be used, even if
this results in longer or more verbose code.  For example, if I want
to solve a quadratic equation aX^2 + bX + c = 0, I think it is much
better to write down the explicit formula, than to call a full-fledged
root-finding library.  The matrix decompositions are a similar case in
my view (at least for the small dimensional case, say n=20).  It is
easier to find a trivial bug in a naive implementation, than to keep
up to date with a large library.

Now, if you need fancier or high dimensional linear algebra, the
problem is still unresolved, and difficult choices need to be made.  I
think a compromise solution is best: provide lapack, gsl, which are
stable, and allow programs that use other solutions, e.g., eigen, with
a specific version embedded.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: reallysmall_svd.c
Type: text/x-csrc
Size: 2245 bytes
Desc: not available
URL: <http://tools.ipol.im/mailman/archive/discuss/attachments/20111108/72e77016/attachment.c>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: eig3.c
Type: text/x-csrc
Size: 5954 bytes
Desc: not available
URL: <http://tools.ipol.im/mailman/archive/discuss/attachments/20111108/72e77016/attachment-0001.c>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: ccmath_svd.c
Type: text/x-csrc
Size: 9238 bytes
Desc: not available
URL: <http://tools.ipol.im/mailman/archive/discuss/attachments/20111108/72e77016/attachment-0002.c>