Matrix Inversion ANA

06122016, 06:29 PM
Post: #1




Matrix Inversion ANA
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; 

06122016, 09:30 PM
Post: #2




RE: Matrix Inversion ANA
Is this different from the builtin inv(matrix) function? Or matrix^1? Or 1/matrix?
<0ɸ0> Joe 

06132016, 03:27 PM
Post: #3




RE: Matrix Inversion ANA  
06132016, 06:30 PM
Post: #4




RE: Matrix Inversion ANA
(06132016 03:27 PM)churichuro Wrote: It does exactly the same Is there some advantage using your method over the builtin function? In other words, why would one want to enter a function that duplicates something builtin? Bob Prosperi 

06142016, 10:55 PM
Post: #5




RE: Matrix Inversion ANA
(06132016 06:30 PM)rprosperi Wrote:(06132016 03:27 PM)churichuro Wrote: It does exactly the same 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. 

06152016, 01:32 AM
Post: #6




RE: Matrix Inversion ANA
(06142016 10:55 PM)churichuro Wrote: That is a good question, 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 builtin function uses, and how can you be sure this ANA function would be faster if compiled? Bob Prosperi 

06162016, 12:23 AM
Post: #7




RE: Matrix Inversion ANA
I made a program for this method for the TI89 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 TINspire 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. 

06162016, 08:12 AM
Post: #8




RE: Matrix Inversion ANA
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 

06162016, 10:27 AM
Post: #9




RE: Matrix Inversion ANA
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 inplace. The 15C has only 64 registers and would otherwise be unable to invert an 8x8 matrix. Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE 

06162016, 10:57 AM
Post: #10




RE: Matrix Inversion ANA
(06162016 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 inplace, the equations are then solved and it is converted back (possibly later). Pauli 

06162016, 11:06 AM
(This post was last modified: 06162016 11:06 AM by Werner.)
Post: #11




RE: Matrix Inversion ANA
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE 

06162016, 11:20 AM
Post: #12




RE: Matrix Inversion ANA
(06162016 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 

06162016, 11:21 AM
(This post was last modified: 06162016 11:23 AM by Werner.)
Post: #13




RE: Matrix Inversion ANA
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 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE 

06162016, 09:45 PM
Post: #14




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

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)