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: |
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
- Report Upstream set to N/A
- Type changed from defect to enhancement
comment:2 Changed 11 years ago by
- Status changed from new to needs_work
comment:3 follow-up: ↓ 4 Changed 10 years ago by
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
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
.
comment:6 Changed 8 years ago by
- Milestone changed from sage-5.11 to sage-5.12
comment:7 Changed 7 years ago by
- Milestone changed from sage-6.1 to sage-6.2
comment:8 Changed 7 years ago by
- Milestone changed from sage-6.2 to sage-6.3
comment:9 Changed 7 years ago by
- Milestone changed from sage-6.3 to sage-6.4
here is an example and some thoughts:
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:
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):
Dima