Matrix Inversion ANA

[churichuro]

This program will invert Matrixes

EXPORT ANA(a)
BEGIN
local i,j,k,r,p,b,n;
r:=1;
n:=rowDim(a);
b:=MAKEMAT(0,n,n);
FOR k FROM 1 TO n DO
p:=a[k,k];
FOR j FROM 1 TO n DO
b[k,j]:= -a[k,j];
b[j,k]:= a[j,k];
END;
b[k,k]:=r;
FOR i FROM 1 TO n DO
IF i=k THEN CONTINUE; END;
FOR j FROM 1 TO n DO
IF j=k THEN CONTINUE; END;
b[i,j]:=(p*a[i,j]-a[i,k]*a[k,j])/r;
END;
END;
a:=b;
r:=p;
END;
return(a/p);
END;

[Joe Horn]

Is this different from the built-in inv(matrix) function? Or matrix^-1? Or 1/matrix?

[churichuro]

(06-12-2016 09:30 PM)Joe Horn Wrote:  Is this different from the built-in inv(matrix) function? Or matrix^-1? Or 1/matrix?

It does exactly the same
only that by a different method

[rprosperi]

(06-13-2016 03:27 PM)churichuro Wrote:  It does exactly the same
only that by a different method

Is there some advantage using your method over the built-in function? In other words, why would one want to enter a function that duplicates something built-in?

[churichuro]

(06-13-2016 06:30 PM)rprosperi Wrote:  

(06-13-2016 03:27 PM)churichuro Wrote:  It does exactly the same
only that by a different method

Is there some advantage using your method over the built-in function? In other words, why would one want to enter a function that duplicates something built-in?

That is a good question,
It is like asking why use a calculator
when you can use a computer?

on the other hand, this method is more resource efficient and faster
Unfortunately, as basic is interpreted and the primitive in the calculator
is a compiled version results in the primitive is faster.

if the method ANA was compiled, it would be at least twice as fast than
the original function.

for the moment there is no special reason for you to use this little program.

[rprosperi]

(06-14-2016 10:55 PM)churichuro Wrote:  That is a good question,
It is like asking why use a calculator
when you can use a computer?

on the other hand, this method is more resource efficient and faster
Unfortunately, as basic is interpreted and the primitive in the calculator
is a compiled version results in the primitive is faster.

if the method ANA was compiled, it would be at least twice as fast than
the original function.

for the moment there is no special reason for you to use this little program.

I agree that having the source code for a function can be useful, as it could help you understand how and why unexpected answers are produced, and other insights.

How do you know the algorithm the built-in function uses, and how can you be sure this ANA function would be faster if compiled?

[churichuro]

I made a program for this method for the TI-89 calculator in C language, which is compiled, and the result is that I got a function that is 6 to 10 times faster than the original function.

in the TI-Nspire try to invert a matrix of 50x50 elements but resources were exhausted and did not give me a result, but using the same method ANA, calculator achievement give a correct result.

[Paul Dale]

This algorithm might be faster but I have doubts it will be as numerically stable as the builtin function. There is no pivoting involved which is a good sign it will be unstable for some inputs. I'm also unsure it will necessarily be more resource efficient, you need a second matrix of the same size as the input and it is possible to do a matrix inversion insitu with a temporary vector of length n.

The 34S firmware does both of these -- it uses the Doolittle LU decomposition with partial pivoting as the basis for its matrix inversion. This is what HP used in 15C and I suspect it is what they still use -- writing numeric code is difficult and theirs was designed by one of, if not, the best.

Writing programs like this is always a bit of fun


- Pauli

[Werner]

Even without a temporary vector of length n (I suppose you're not referring to the permutation vector)

F=inv(L), G=inv(U) and GxF can all be computed in-place.

The 15C has only 64 registers and would otherwise be unable to invert an 8x8 matrix.

Cheers, Werner

[Paul Dale]

(06-16-2016 10:27 AM)Werner Wrote:  Even without a temporary vector of length n (I suppose you're not referring to the permutation vector)

The extra vector of length n is the permutation vector contain indices that define the partial pivots.

I believe, but am not certain, that these values are stored in the matrix descriptor itself rather than occupying other registers. Each permutation index fits in a nibble and there is plenty of space in the matrix description. Look at Synthetic Matrices. This means that the matrix can be quickly converted from the LU decomposition back into its original form. E.g. solving a system of equations doesn't destroy the matrix but there isn't space to store a second copy containing the decomposition -- I suspect it is converted in-place, the equations are then solved and it is converted back (possibly later).


Pauli

[Werner]

In the 15C they are encoded in the pivots (A(i,i), I mean) - in the exponent sign, I seem to remember. Read it once, somewhere.
Werner

[Paul Dale]

(06-16-2016 11:06 AM)Werner Wrote:  In the 15C they are encoded in the pivots (A(i,i), I mean) - in the exponent sign, I seem to remember. Read it once, somewhere.

The pivots can't go in the diagonal of the decomposed matrix. The U part of the LU decomposition uses the diagonal elements. The L part has 1's down its diagonal so they don't need to be stored.

Probably, the easiest way to confirm if the descriptor contains the pivots would be to modify an older version on nonpareil and inspect the registers.


Pauli

[Werner]

I meant that the permutation info was recorded in the diagonal elements, sorry for not expressing myself more clearly!

Found it:
"The row interchanges
are recorded in the otherwise underused XS
format fields of the result matrix's diagonal elements"

on pg. 32 of

HP Journal May 1983

Cheers, Werner

[Paul Dale]

There would be space in the XS for up to 8x8 matrices, so that would work.
A clever place to put this information


Pauli