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

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