wilmott.com

Two questions:1.) I'm looking for an implementation (or description) of an algorithm for calculating the inverse of a symmetric matrix. The routine that I have at the moment makes no distinction between symmetric and non-symmetric matrices, which means that the resulting inverse matrix is "not" always symmetric! (when examined element by element to the full available number of decimal places - of course it is symmetric really - this is just numerical error!). It would seem to me that there therefore must exist a more efficient algorithm. What is it?!! 2.) My matrix library is implemented in C#. This is fine most of the time, but when I get a 500 x500 matrix which I have to multiply against a 500 x 1 vector 64 x 130 times things start to get rather slower! What are my options? I'm aware that I could use my library in it's original C++ form, but are there other options - like special fast matrix libraries coded in an even lower language, or a library provided by a chip manufacturer, or perhaps there are approximation routines for such calculations? ! Any advice/ help appreciated!Thanks Charlie

I'm not sure about porting libraries to C#, as I am a C++.NET developer, but you might want to look up the MATLAB compiler and use functions from their libraries.

Charlie,if you have an ill-conditioned matrix it is not surprising yourinverses are not symmetrical due to accuracy problems. Have you had a look at Sec. 2 of The Numerical Recipes in C:http://www.library.cornell.edu/nr/cbookcpdf.htmlI am not sure if this contains an optimized version for symmetricmatrices. I doubt you will be able to get below O(N^3) in computingtime.I do not know how C# exactly compares to C/C++ for matrix operations butwould suggest you either use the Numerical recipes code (availableon the Web) or have a look at NAG etc. libraries.http://www.nag.co.uk/numeric/fn/manual/ ... ,Hans-Marc

Most of the time you dont need to compute the inverse of a matrix explicitky. It is much more efficient to work with the Cholesky factors of the matrix. There are algorithms that can be used to work with symmetric and asymmetric matrices. For symmetric, we usually use the LLt and for asymmetric we use LU factorisation.

Thanks for all your replies. I will look into them all. A bit of googling and I also found some other speed enhancements:for A^r where A is an NxN matrix the number of matrix multiplications can be reduces from r to O(log2(r)) by observing that (for example) A^2 = A x A A^4 = (A^2) x (A^2)A^8 = (A^4) x (A^4)A^15 = A^8 x A^4 x A^2 x A=6 matrix multiplications Likewsie matrix multiplication (in general) can be reduced from an O(n^3) operation to an O(n^log2(7)) operation (<a target=new class=ftalternatingbarlinklarge href="http://www.library.cornell.edu/nr/bookc ... >Numerical Recipes</a>. A little further research and I found that if you are prepared to do it, different memory storage schemes for the matrix can also deliver further improvements for large matrices.The easiest way to get a speed increase will be to get and use a graphics which is dedicated to linear algebra (as many of them are). And/or depoly the code on a MacAnother (to me very surprising) route that I found is to use reflection in C# to write out the matrix multiplication semi-explictly at run time. Apparantly this will deliver around 50% speed increases for medium sized matrices (and I presume with tweaking could deliver similar improvements for larger (i.e. N>100) matrices).I looked at our current matrix library (based on LAPACK) and found that it was (for instance) creating a new matrix for the transpose. By overloading the indexing I got a more efficient transpose property/ method. For large matrices some kind of matrix memory pooling may deliver some further improvements?!The LLt algorithm? I don't recognise the name - but will check my copy of nr when I get home and see what is there! If I'm really lucky it may already be somewhere in my current library!!!regarding C++/ C# - one would expect that C++, especially unmanaged, would be faster, but I then loose the depth of debugging that I get with C#, since everything else is C#! I can step into the matrix operations and see what has gone wrong. (No doubt I will now be informed that this is rubbish and that I can use the two together in the same Visual Studio Solution ??!!!)

LLt is specialised form of LU for symmetric matrices. the Lt is L transposed - sometimes its also called LDLt. check tnetlibhttp://www.netlib.no/netlib/c/

The reason you can't find it is because "LLt" isn't recognised as the name of an algorithm. The method daveangel is trying to talk about is Cholesky factorization. And it doesn't apply to symmetric matricies, it applies to symmetric positive definite matricies! If A is s.p.d it will give you unique L such that A = LL^T and L is a lower triangular matrix.