[racket] math/matrix

2014-05-14 15:31 GMT+02:00 Konrad Hinsen <konrad.hinsen at fastmail.net>:
> Neil Toronto writes:
>
>  > One thing we should really do is get your LAPACK FFI into the math
>  > library and have `flmatrix-solve` use that, but fail over to Racket code
>  > systems that don't have LAPACK. If I remember right, it would have to
>  > transpose the data because LAPACK is column-major.
>
> Sounds good. Is that LAPACK FFI already available somewhere?

See my post in this thread 3 days ago. It is attached.
The code works (to my knowledge), but improvements are
certainly possible. For example, it ought to be reasonable
straightforward to support more than Flonum matrices.

That said, I hope someone takes the code and turns it into
something more.

I have tested it on OS X, on other systems you might
need add names/paths of the libraries in question.

>  > Some thoughts, in no particular order:
>  >
>  >   1. Because of transposition and FFI overhead, there's a matrix size
>  > threshold under which we ideally should use the code below, even on
>  > systems with LAPACK installed.
>
> Transposition can be avoided in many practically relevant situations
> (e.g. symmetric matrices), so such decisions should be taken on a
> case-by-case basis.

The problem is that (Matrix Flonum) and LAPACK has different
representation for matrices. In order to take a (Matrix Flonum)
and turn it into a LAPACK one, the entries in the underlying vector
must be transposed.

If we(I) had realized it sooner, the representation in (Matrix Flonum)
would have been the same.

> As for the FFI overhead, even after reading the introduction to the FFI
> documentation I have no clear idea of how important it is.  It seems
> (but I am not at all certain) that there are different vector-like
> data structures, some of which are optimized for access from Racket and
> others for access via C pointers. If that's true, it may be of interest
> to have a constructor for arrays and matrices that live in "C-optimized"
> space. Or even in "Fortran-optimized" space, using column-major storage.

I think it is negligible for all but very small matrices.

>  >   2. Because of small differences in how it finds pivots, LAPACK's
>  > solver can return slightly different results. Should we worry about that
>  > at all?
>
> I'd say documenting the issue is sufficient.

Agree.

>  >   3. A design decision: if a matrix contains just one flonum, should we
>  > convert it to (Matrix Flonum) and solve it quickly with
>  > `flmatrix-solve`, or use the current `matrix-solve` to preserve some of
>  > its exactness?
>  >
>  > I lean toward regarding a matrix with one flonum as a flonum matrix.
>  > It's definitely easier to write library code for, and would make it
>  > easier for users to predict when a result entry will be exact.
>
> I agree. I'd like to see a clear use case for any other approach.
> Predictability is  important.

Following the inexactness-is-contagious principle, I agree.

--
Jens Axel Søgaard