algorithms for 'invert', was: program works with maxima 5.29.1 but freezes with maxima 5.30.0

bfloats allow unrestricted integer exponents, but you can still get the
infamous catastrophic cancellation because precision is finite.  The
numerical analysts have thought quite a bit about how to code matrix
operations (including inversion) for efficiency and accuracy using floats.
I imagine someone's thought about issues of efficiency for the rational
case (including cost of GCDs, size of intermediate results, etc.) and
efficiency and simplicity/tractability of result for the symbolic cases.

           -s


On Fri, Apr 19, 2013 at 3:45 PM, Henry Baker <hbaker1 at pipeline.com> wrote:

> I should point out that programming a really reliable, robust matrix
> inversion algorithm is a lot harder than most mathematicians would make it
> seem.
>
> For example, if the matrix is 'almost' singular (i.e., det(M) almost
> zero), then things can go haywire.  In order to avoid various kinds of
> overflow/underflow, you may want to fiddle with the floating point
> exponents in advance, etc.  BTW, don't bfloats allow unrestricted integers
> for their exponent, so that overflow/underflow is never an issue?
>
> Although the adjoint algorithm has terrible algorithmic complexity, at
> least it is perspicuous, so that the only division is the division by the
> determinant.
>
> As I pointed out in post on the topic of Matlab, many codes try to avoid
> the computation of an explicit matrix inverse, because the inverse of a
> sparse matrix is usually non-sparse.  For these reasons, Matlab & other
> packages provide the ability to compute (M^^-1).b without the explicit
> computation of (M^^-1).  I guess that Maxima could do the same thing via
> Cramer's Rule (or equivalent), so that the only division would still be
> that by det(M).
>
> At 10:46 AM 4/19/2013, Stavros Macrakis wrote:
> >I believe adjoint was chosen as the default because it works well in the
> symbolic case.
> >
> >According to Henry, LU is better for floats (but not rationals).  Perhaps
> the right policy is: if matrix is entirely float or bfloat, use LU,
> otherwise use adjoint.
> >
> >If adjoint produces simpler results, it sounds like the right algorithm
> for symbolic matrices regardless of size, as a first cut at least.  An
> unwieldy symbolic result returned fast is probably less useful to most
> users than a tractable symbolic result that takes longer.
> >
> >Perhaps there are more sophisticated policies possible (depending on
> sparseness or proportion of symbolic vs. small integer entries, etc.), but
> it seems like a good start.
> >
> >           -s
> >
> >On Fri, Apr 19, 2013 at 12:27 PM, Robert Dodier <robert.dodier at gmail.com>
> wrote:
> >On 2013-04-19, Henry Baker <hbaker1 at pipeline.com> wrote:
> >
> >> So it depends upon what you're trying to achieve, and which
> fields/rings are involved.
> >
> >> The advantages of certain 'efficient' matrix operations like LU quickly
> disappear when
> >> operating with rings other than floating point numbers.
> >
> >Can we derive from these considerations a policy for automatically
> >choosing which algorithm to use? E.g. if small then adjoint method, else
> >if numerical then LU, otherwise print a warning message ("umm, it might
> >or might not work") and use either one. Each method should be available
> >as a separate function (e.g. invert_by_adjoint) so that the user can
> >explicitly choose one or another.
> >
> >I realize that such heuristics are not foolproof, but it would be an
> >improvement on the current situation, and in any event it is
> >unreasonable to expect users to have any idea how to choose among
> >algorithms.
> >
> >best
> >
> >Robert Dodier
>
>