Michael Orlitzky wrote:
> 
> You're going to have a rather hard time convincing me that, upon 
> encountering a decimal point, any other math software is going to flip 
> out and silently return incorrect results for basic matrix operations.
> 
> Recall:
> 
> sage: m = matrix([ [-0.3, 0.2, 0.1],
>                     [0.2, -0.4, 0.4],
>                     [0.1, 0.2, -0.5] ])
> 
> sage: m.echelon_form()
> 
>        [ 1.00000000000000 0.000000000000000 0.000000000000000]
>        [0.000000000000000  1.00000000000000 0.000000000000000]
>        [0.000000000000000 0.000000000000000  1.00000000000000]
> 


The problem here is that no one has really implemented a numerically 
stable version of echelon_form for RR.  I thought we called scipy for 
rank calculations over RDF, but apparently we don't even do that!  Scipy 
answers correctly:

sage: import scipy
sage: scipy.rank(m.numpy())
2


Very little attention has been paid to numerically stable linear algebra 
for non-exact rings.  We do get some numerically stable things for 
matrices over RDF and CDF because those are based on numpy and scipy 
(which call blas and lapack routines).  However, apparently we don't 
call numpy/scipy to calculate the rank, but instead rely on our unstable 
echelon form computation!

Scipy and numpy don't have facilities for calculating the echelon form 
because it does not make sense to calculate the echelon form in 
numerical settings (paraphrasing their words to the best of my memory).

That said, there is plenty of work for someone who would like to 
implement numerically stable, arbitrary precision linear algebra in 
Sage.  I've talked with at least one person who has worked on such a C 
library (arbitrary precision, numerically stable linear algebra) about 
contributing to Sage in the last year, but politics have held them up 
from contributing any code.  We don't even need huge experts in the 
area; even someone implementing basic textbook algorithms would be a 
huge improvement, I think.

If someone wanted to take make a good library of numerically stable 
linear algebra routines based on mpfr, I think it would be absolutely 
awesome.

Thanks,

Jason



-- 
Jason Grout


--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to 
sage-support-unsubscr...@googlegroups.com
For more options, visit this group at 
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
-~----------~----~----~----~------~----~------~--~---

Reply via email to