Matrix Inversion ANA - Printable Version +- HP Forums ( https://www.hpmuseum.org/forum)+-- Forum: HP Software Libraries ( /forum-10.html)+--- Forum: HP Prime Software Library ( /forum-15.html)+--- Thread: Matrix Inversion ANA ( /thread-6407.html) |

Matrix Inversion ANA - churichuro - 06-12-2016 06:29 PM
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; RE: Matrix Inversion ANA - Joe Horn - 06-12-2016 09:30 PM
Is this different from the built-in inv(matrix) function? Or matrix^-1? Or 1/matrix? RE: Matrix Inversion ANA - churichuro - 06-13-2016 03:27 PM
(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 RE: Matrix Inversion ANA - rprosperi - 06-13-2016 06:30 PM
(06-13-2016 03:27 PM)churichuro Wrote: It does exactly the same 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? RE: Matrix Inversion ANA - churichuro - 06-14-2016 10:55 PM
(06-13-2016 06:30 PM)rprosperi Wrote:(06-13-2016 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. RE: Matrix Inversion ANA - rprosperi - 06-15-2016 01:32 AM
(06-14-2016 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 built-in function uses, and how can you be sure this ANA function would be faster if compiled? RE: Matrix Inversion ANA - churichuro - 06-16-2016 12:23 AM
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. RE: Matrix Inversion ANA - Paul Dale - 06-16-2016 08:12 AM
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 RE: Matrix Inversion ANA - Werner - 06-16-2016 10:27 AM
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 RE: Matrix Inversion ANA - Paul Dale - 06-16-2016 10:57 AM
(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 RE: Matrix Inversion ANA - Werner - 06-16-2016 11:06 AM
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 RE: Matrix Inversion ANA - Paul Dale - 06-16-2016 11:20 AM
(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 RE: Matrix Inversion ANA - Werner - 06-16-2016 11:21 AM
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 RE: Matrix Inversion ANA - Paul Dale - 06-16-2016 09:45 PM
There would be space in the XS for up to 8x8 matrices, so that would work. A clever place to put this information Pauli |