#7392: RDF/CDF matrices should call scipy for rank
------------------------------+---------------------------------------------
Reporter: jason | Owner: was
Type: enhancement | Status: needs_work
Priority: major | Milestone: sage-4.5.3
Component: linear algebra | Keywords:
Author: | Upstream: N/A
Reviewer: | Merged:
Work_issues: |
------------------------------+---------------------------------------------
Changes (by dimpase):
* status: new => needs_work
Comment:
here is an example and some thoughts:
{{{
sage: m=matrix(2,[1.5,1.75,1.5,1.75])
sage: m
[1.50000000000000 1.75000000000000]
[1.50000000000000 1.75000000000000]
sage: m.rank()
2
sage:
}}}
This can be traced to _echelon_in_place_classical in
$SAGE_ROOT/devel/sage/sage/matrix2.pyx, which attempts to compute
row-echelon form of this matrix (basically, Gaussian elimination), and
fails due to a precision loss.
Numerics people recommend using the QR-decomposition, (i.e. into product
of an orthogonal and an upper-triangular matrices q and r) followed by
inspection of r. This is apparently more stable.
For instance:
{{{
sage: from scipy.linalg import qr
sage: q,r=qr(m)
sage: q
array([[-0.70710678, 0.70710678],
[ 0.70710678, 0.70710678]])
sage: r
array([[ -2.12132034e+00, -2.47487373e+00],
[ -0.00000000e+00, 3.96275894e-16]])
}}}
One then has to decide how to round entries of r to 0.0.
This would take care of RDF and the like.
Interestingly, the computation of the charpoly is more robust, and can
perhaps provide an alternative (just count the powers of x it is divisible
by):
{{{
sage: m.charpoly()
x^2 - 3.25000000000000*x
sage:
}}}
Dima
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/7392#comment:2>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
