On Jan 7, 2010, at 9:56 PM, Dag Sverre Seljebotn wrote:
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).
+1
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.
I'm not usually a fan of silently ignoring arguments, but it may make
sense to do so if you can't always know what algorithms will be used.
rtol and atol seem a bit cryptic to me, but perhaps that's because I'm
not in the numerical field?
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')
I'm not as sure about this idea (though not completely opposed
either). It seems more natural have a way of asserting, or even
requesting, a property of a matrix (maybe with a check flag) and
having it deduce what algorithms to use based on that. In terms of
naming the method could be called set_hints, and could take an
algorithms parameter.
sage: A.log_determinant() # say this ends up with QR
sage: A.solve_right(B) # uses cached QR
Are you implying that it would use the QR if it is already cached, but
otherwise may use a different algorithm (which could be handy)?
- Robert
--
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
