Robert Bradshaw wrote:
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.
That's what I was thinking.
rtol and atol seem a bit cryptic to me, but perhaps that's because I'm
not in the numerical field?
I don't care much about that, just copied what NumPy uses.
relative_tolerance is fine.
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.
Hmm. This gets tougher if one wants to introduce iterative methods as
well. Also the typical usecase is "the set of methods I want to call". I
guess the decompositions are roughly covered in this way:
- Cholesky: Hermitian
- QR: Make sure to call least_squares or QR before solve_right
- SVD: Make sure to singular_values or SVD before solve_right
Does anybody know if one would like to use QR instead of LU for
numerical stability in any setting? Anyway, I think we want to avoid
sage: A.check_condition_number_less_than(2) # force use of frobnicate
iteration method for solving
sage: my_code_using_solve(A)
Anyway, perhaps one can wait with introducing set_hints until there's a
proven need for it. The particular case of Cholesky is covered by the
"subset of matrix space" idea. I.e.:
sage: A.set_immutable()
sage: A = A.hermitian_matrix() # checks, and makes element of subset
Making it immutable first ensures sharing of data.
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)?
I was just implying the matrix wasn't Hermitian, really. But I envision
that the default algorithm is to search for and use any cached
decomposition first.
So just calling QR first is a way to make QR the default. Perhaps that
is enough, though it doesn't deal with the case of
sage: code_which_only_may_call_solve(A)
--
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
