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

Thanks, Barton.

Yes, I was just learning about condition numbers.

 From what Prof. Strang said on his MIT YouTube video, for a symmetric positive definite matrix, the condition number is the ratio of the largest to the smallest (both necessarily real&positive) eigenvalues.  According to Prof. Strang, the loss of precision is related to log(condition#).

For other matrices, the condition number is sqrt(condition(transpose(M).M)).  (I assume that this still works even when M is symmetric positive definite.)

What I didn't understand from Prof. Strang is what happens when transpose(M).M is still singular.

BTW, it looks like the ill-conditioned matrices you're showing are precisely those generated by 'GCD-with-multipliers' applied to large relatively prime numbers; is this correct?

At 04:46 AM 4/20/2013, Barton Willis wrote:
>> For example, if the matrix is 'almost' singular (i.e., det(M) almost zero), then things can go haywire.
>
>A matrix can be ill conditioned with a determinant of one; famous example:
>   (%i1)  m : matrix([63245986, 102334155],   [102334155,   165580141])$
>
>The determinant of m is 1, but its condition number is about 10^16
>
>   (%i2) [determinant(m), mat_norm(m,'inf) * mat_norm(m^^-1,'inf)];
>
>   (%o2) [1,71778070001175616] Let m be a coefficient matrix
>
>(%i4) eq : m . matrix([x],[y])- matrix([1],[1])$
>
>(%i5) eq  : xreduce('append,args(eq))$
>
>Solve using rational numbers
>
>(%i6) linsolve(eq,[x,y]);
>
>(%o6) [x=63245986,y=-39088169]
>
>Solve using binary64
>
>(%i7) linsolve(float(eq),[x,y]),keepfloat : true;
>
>(%o7) [x=3.1622993000000007*10^7,y=-1.95440845*10^7]
>
>Partial pivoting doesn't save the calculation--at least the estimated condition number will warn the alert user:
>
>(%i8) linsolve_by_lu(m,[1,1],'floatfield);
>
>(%o8) [matrix([4.1475558898394227*10^7],[-2.5633305101605773*10^7]),4.7070743271309312*10^16]
>
>--Barton