Hi Jason Grout (and others),
we started something on IRC last night which I really wanted to finish,
since you mentioned that you might pick up work again on dense RDF/CDF
matrices at some point.
Since I'm working on this from the sparse end (but can't upload my
changes for some time) here's an attempt to get us coordinated -- I'd
love your opinion on the following when you have time for it.
Background:
a) The Sage-implemented solvers (and determinant finder!) are
inherited by RDF/CDF matrices are totally unsuitable for numerics
b) I believe Sage should by default give correct, dependable answers
or an exception, also for numerics
c) Which solver to use for a matrix can be seen as a property of the
matrix.
So my thoughts are:
Do not fix, but remove, solve_left_LU (#4932). solve_left/solve_right
should be used instead.
Have solve_X take an "algorithm" argument, which can be 'cholesky',
'lu', 'qr' and so on. This:
sage: A.solve_right(B, algorithm=['cholesky', 'lu'])
would first attempt Cholesky (which would fail early if the matrix
didn't have all-positive (or all-negative) diagonal, a little later if
it wasn't Hermitian, and only after attempting if the matrix was
indefinite).
By default, solvers would use iterative refinement:
sage: A.solve_right(B, algorithm=['cholesky', 'lu'], rtol=1e-5, atol=1e-5)
This would use iterative refinement with a given decomposition to get
below the specificed relative and absolute errors, see numpy.allclose
[1], and an exception could be raised if everything was horribly
unstable. One could pass iterate=False or similar to skip this.
solve_X would accept any number of additional keyword arguments, any
unknown arguments are simply ignored.
Since which solver is suitable is in some sense a property of the
matrix, there would be a set_algorithms method which would set the
default order of algorithms to try and their options. Also all the same
options would apply to other methods (log_determinant, inverse, and so on).
sage: A.set_algorithms(['cholesky', 'cholesky', 'qr', 'lu'], rtol=1e-3,
cholmod_permutation_finder='metis')
sage: A.log_determinant() # say this ends up with QR
sage: A.solve_right(B) # uses cached QR
The defaults would depend on the type of matrix as discussed yesterday.
Many have direct solutions -- the explicit Hermitian matrices would
default to attempting Cholesky first as a simple heuristic.
[1]
http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html
--
Dag Sverre
--
To post to this group, send an email to sage-devel@googlegroups.com
To unsubscribe from this group, send an email to
sage-devel+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sage-devel
URL: http://www.sagemath.org
