Opened 12 years ago

Last modified 7 years ago

#7392 needs_work enhancement

RDF/CDF matrices should call scipy for rank

Reported by: jason Owned by: was
Priority: major Milestone: sage-6.4
Component: linear algebra Keywords:
Cc: Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description

Right now, we rely on the horrible (for floating point) generic echelon implementation to calculate the rank of an RDF/CDF matrix.

In fact, over RR, it wouldn't hurt to call numpy/scipy either.

Change History (9)

comment:1 Changed 11 years ago by jason

  • Report Upstream set to N/A
  • Type changed from defect to enhancement

comment:2 Changed 11 years ago by dimpase

  • Status changed from new to needs_work

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

comment:3 follow-up: Changed 10 years ago by rbeezer

As a precursor to computing rank, I've written a short routine to get singular values and conveniently work with deciding the cutoff for zero values. At #10802. I'll get to this ticket next.

I think a QR decomposition can be less expensive to compute, but a bit less reliable. And I've got good advice on how to make an automated decision on "zero" entries with the SVD.

comment:4 in reply to: ↑ 3 Changed 10 years ago by jason

Replying to rbeezer:

As a precursor to computing rank, I've written a short routine to get singular values and conveniently work with deciding the cutoff for zero values. At #10802. I'll get to this ticket next.

I think a QR decomposition can be less expensive to compute, but a bit less reliable. And I've got good advice on how to make an automated decision on "zero" entries with the SVD.

That's very interesting. Should we be submitting it upstream to scipy, where it would be used by a much wider audience than Sage?

comment:5 Changed 10 years ago by leif

.

comment:6 Changed 8 years ago by jdemeyer

  • Milestone changed from sage-5.11 to sage-5.12

comment:7 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.1 to sage-6.2

comment:8 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.2 to sage-6.3

comment:9 Changed 7 years ago by vbraun_spam

  • Milestone changed from sage-6.3 to sage-6.4
Note: See TracTickets for help on using tickets.