I'm working on some dynamic problems, and often we need to determine the inverse of a matrix of order 50x50 and larger. I need to speed up the process.
Do just some parts of your matrix elements change from one step to another? If so you could look up the woodburry identity. This allows you to keep the information of your inverse matrix and use it again when implementing slight changes in the elements, without making a hole new inverse matrix. Otherwise if you have a sparse matrix, there are plenty of preprocessing steps to speed up for a follow up gauss elimination algorithm.
Well, how do you do it now? Iterative solvers can find approximate solutions very fast, so one important information lacking in your problem description is exactly how accurate each inverse needs to be.
I currently use Wolfram Mathematica, and it is slow. No, my matrix isn't sparse, all elements are nonzero, but it is a symmetrical. Thank you for the idea of parallel computation..
If you have a symmetric matrix then you can easily apply a skyline method for storing your matrix and the Gauss method for the inverse and solution (I assume that you solve the dynamic behavior of a structure). If you manage to write this code in any low-level programming language like C, C++ or fortran (you can find also open source subroutines that will do it for you), then you will have the inverse in less than a second (personally I prefer fortran 90/95 with the latest Intel Visual Fortran Compiler). 50x50 or even 1000x1000 is not even close to what we call a large matrix. According to the current numerical methods for inverting and solving a linear system of equations this is considered to be small, especially when you are handling numerically the matrix as an array during the analysis procedure.
Now on the other hand, if you decide to continue with Mathematica then you need to go parallel as suggested by Deepraj. Nevertheless, parallel processing will not improve dramatically your computational cost given that you will base the optimization on your hardware and not the code itself. For example having a quad processor, it will allow you to speed up the procedure by 3 to 3.5 times approximately given the nature of the implemented parallel code. If you have to inverse several times during a single analysis (or perform numerous analyses) then the duration of the computational time will remain relatively large. By incorporating my first proposal, you will spend some time programming but eventually perform your analysis in terms of seconds.
Conclusively, you have to consider how many analyses you are planing to perform in the future and decide which path to take.
Deepraj Paul and George Markou are definitely right. I had similar problem with matlab and ODE solver with large nonhomogenous differential equations... After C++ programing, calculations are reduced from aprox 45 mins to few mins..
In C or FORTRAN I would just use a good implementation of LAPACK (Intel Math Kernel Library, ATLAS a.o.) and choose the algorithm that is best suited for the matrix at hand. There are stable and fast methods for the inversion of general, symmetric, positive-definite dense and sparse matrices and parallelized code is usually available (but parallel execution might only pay of for large matrices, for very large ones you could also consider GPU-computing).
You could also examine how your matrix changes between the inverse computations: if it can be described as a (series of) rank updates, it might not be necessary to invert the matrix each time. Instead e.g. Sherman–Morrison formula might be applicable. You can also find some methods for rank-updates in LAPACK.
Check out PETSc (Portable Extensible Toolkit for Scientific Computing). You have options to use preconditioners and a variety of Krylov Subspace solvers here.
I doubt you need the inverse explicitly. Probably you need to apply it to another objexct, e.g. a vector. If so, the inversion should be avoided completely! If you really need the inverse explicitly, a fast method exploiting modern computer achitecture as available in current noteboks and desktops, read "Matrix Inversion on CPU-GPU Platforms with Applications in Control Theory" by Peter Benner, Pablo Ezzatti, Enrique S. Quintana-Ortí, Alfredo Remón, see http://onlinelibrary.wiley.com/doi/10.1002/cpe.2933/abstract
As Peter Benner and Gang Mei suggest, I'd advise the use of GPU (CUDA or OpenCL) for maximum performance on a single computer. As they suggested CUDA frameworks, I'd suggest you the ViennaCL package for OpenCL: http://viennacl.sourceforge.net/
Also, are you solving dense or sparse matrix inversion?
(1) As Peter Benner notes, it is rare that one actually requires a matrix inverse. Solving A.x=b for matrix A and vector B is typically 2-3 times faster than inverting A.
(2) If your matrix is positive definite as well as symmetric, then you should get further speed gain over basic LinearSolve by using Method -> "Cholesky".
(3) For a matrix of machine real or complex numbers, Mathematica uses Lapack with whatever BLAS is on your particular platform (and typically this is multithreaded so you get some benefit from parallelism automatically). So this should be fast. If you have an exact or symbolic matrix then that would be a different matter. Assuming that you are working with machine precision numbers, if you are finding the linear algebra to be excessively slow then that is something worth showing to the Tech Support people ([email protected]). Also feel free to send a smallish example my way. (If RG makes it difficult to locate do this, it should not be too hard to track down ny contact info on the Internet. Every spammer seems to find it...)
The first question. Asked by others already: why do you need the inverse explicitly?
Maybe you need the diagonal only (gives _some_ speedup in the positive definite case)?
The second:: what is ''large" for you.
Below some 1000 where is little sense to bother about specialities: take lapack
with the optimized blas, and if you don't have an optimized blas make one yourself using ''atlas''. If dimension is larger, then _warning_ : the inverse might be full even if your matrix ia very sparse, take a (sparse) lu decomposition (umfpack for example) and use forward/backward substitution and store the results where they fit in (maybe you need a diskfile for this)
illcondioning has nothing to deal with largeness, this is a second trap.
Use the 'method of successive bordering' to invert matrices of large order. Starting from the North-West corner, you can invert large matrices step by step without any problem. This method is included in standard books on Computational Mathematics'.
I recommened you to use fortran/C which is the fastest. but if you use If u are using mathlab, then using fortran/c code it can speed-up the inverse computation.
Depending upon your type of matrix to be inverted, if its a sparse matrix, I think there will be specific codes available which can do sparse matrix inversion.
If your matrix isn't a sparse one, then you can go for Parallel Computiation using mathlab/python.
You did not state your system specification,namely the motherboard type(Intel etc), RAM size (512 MB, 1GB, 2 GB etc), the processor speed (512 MHz, 1GHz, 2 GHz). All these factors influence the speed of mathematical computations. Secondly, you have not mentioned how long it takes your system to invert a 50x50, 100x100, 500x500 matrices. This will enable the classification of computational speed.
The best and fastest means of computing matrix inverse is C or C++ based program because they have virtual memory capability whereby part of hard disk space is accessed as memory, and hence they can handle very large matrices. However, you need to increase the page file size by. going to control panel, system, advanced, performance, advanced, and then virtual memory.
There does not seem to be an efficient means of computing the inverse of a non symmetric matrix because you need to form an augmented matrix of the given matrix and a unit matrix and proceed to reduce the given matrix to a unit matrix by row reduction operations.
Many researchers have questioned the need for a matrix inverse. In the field of soil-structure interaction (Civil Engineering) matrix inversion is needed in order to convert the foundation flexibility matrix to a foundation stiffness matrix.
They were not saying you do not need the inverse - they were saying usually you want the inverse to do something else like solving a set of equations, and often the something else is faster than computing the inverse explicitly - of course, if you are using it to solve many systems of equations for different sets of values, then computing the inverse explicitly is faster
it all depends on which problem you are solving and how often
In the soil-structure problem I am talking about (known as a surface element method) the inverse of the flexibility matrix is not being computed as part of the solution of a set of simultaneous equations. The inverse of the foundation flexibility matrix, K_ff, becomes the foundation stiffness matrix K_fs. The foundation stiffness matrix is then added to the superstructure, (raft perhaps) stiffness matrix to form the global stiffness K_g. Thereafter the simultaneous equations are set up namely K_g*D = P.
This can be solved using Gaussian elimination. You do not need to obtain the inverse of K_g. However, you must compute the inverse of K_f.
In addition in the solution of dynamic problems it may be necessary to perform the so called Guyan reduction in order to reduce the size of the stiffness matrix by eliminating the slave degrees of freedom. It is then compulsorily necessary to obtain the inverse of a matrix K_sm.
It seems, someone did not like my answer! There is this method of inverting a matrix by partitioning it. This method of inverting using partitions was modified to invert the matrix step by step, which is known as the method of successive bordering. When the matrix is large, this method works very well indeed.
The formula of the inverse can't be applied to any matrix : you can check this with the example given by Yang in "finite element structurall analysis" where there is an Error in calculating the inverse.
I have faced the problem to invert Successfully a 3x3 matrix for the mixture of materials (concrete+steel) which is like a composite in order to obtain the material properties matrix [D] given by {sigma}=[D]{epsilon}. This matrix is not well conditioned. So I have solved three systems of equations from [Dinv][D]=[Identity] with the Efficient Partial Pivoting Gauss LU Algorithm such that [Dinv]{One column of D}={One column of Identity}.
The fastest language to operate with vectors/matrix is, by far, CUDA. If you operate intensively with matrix and you have access to any modern NVIDIA card try this language.
"Pros ":
- There exist a wide variety of languages that support CUDA via wrappers (C, C++, Fortran, IDL...)
- Speed, speed and speed: a multiplication of two 4000x4000 matrices took me obout 20 seconds using IDL, 0.30 seconds (yes 30/100 seconds) using a CUDA based library (gpuLib) from IDL (CPU=Xeon E5620 @2.4Ghz and 8 threads working in parallel , GPU=Tesla M2050)
"Cons":
- Long learning curve: It is a new concept in programming languages. When using CUDA you are programming a GPU not a CPU and the programming procedure is very different
- It is necessary to completely rewrite the old Fortran and C code (but only the functions involving CPU expensive matrices operations)
Of course, this is not a "short-term solution" ...
Dear Aleksandar. Perhaps your question is better posted in the form: how to speed up inverse of large matrix with a given error. You could calculate a relative large matrix very fast with large errors, but the time may be dramatically extended if a more accurate results you pretend. I have developed a two-stage method in which a "coarse" inverse is firstly obtained and then a refination method by sucessive aproximations reduce the errors to a required value. This method may be more valuable for bad-conditioned matrix. If you are interested I could include this work in my personal publications and send to you a copy of the algorithm and computation codes.
That is why it is sometimes important to use Monte Carlo methods to help make a large problem into a smaller problem but I don't know how to incorporate those techniques into matrix algebra but there is likely a way
Carlos R. Gonzalez is quite right! It is not the speed of computation but the accuracy of the result that A.V. Nikolic should be worried about. It is really amazing how many inverses of sparse matrices and ill conditioned matrices are computed wrongly. Twelve years ago, in the process of developing a finite element program for plate bending, I got stuck on the way because my program could not compute the inverse of a sparse 12 by 12 matrix. My investigations revealed that along the line the program was dividing with a very small number around exp(-9). What I did was to implement partial row pivoting whereby whenever a small pivot is encountered (less than 0.001) the program interchanges the current row with another row with a larger pivot. The best result is obtained if full pivoting (bot row and column) is implemented. However, there is much secondary book keeping that accompanies column pivoting though it be handled easily.
Please, find my modest contribution to the matrix inversion methods in the list of my publications and give me your opinions. The algorithm of refination is very simple and I give a Delphi code, but a MATLAB code could be implemented easily. Sorry, it is in Spanish. It you consider it useful please enjoy it and later give me your impressions.
inverting a matrix is always expensive to do and should be avoided. I would do solve instead as in matlab: inv(A) = A \ I, where I is the identity matrix
Fist of all we need to study the nature of the matrix, and accordingly we implement the formula. For example: if the given matrix is a triangle or banded matrix, then we can find their inverses quickly.
Dear colleagues, I just insert into my publication "Rutina para el refinamiento de la inversión de matrices de cualquier rango con la precisión deseada" an addendum with a MATLAB code and some experiments carried out with it inverting matrices of 100x100 and 200x 200 size. The calculations were implemented calculating random inverse matrices with condition numbers between 4701 and 121730 in two steps: a Gauss-Jordan inversion without pivoting and a later refining with my algorithm. The whole inversion time for 100x100 was about 0.05 s and for 200x200 was about 0.3- 0.5 s. The errors of the inverses were about 2.10(-13) and 2.10(-12) respectively. Just prove it and tell me your own results.
This is with Mathematica but I would expect similar results from Matlab since it iand Mathematica both use Lapack for this. I show timings and max error in computing the residuals for several random examples at dimension 200.
Error in individual entries of the inverse matrices will be different from the above, though I would not expect them to be orders of magnitude different from what the residuals show.
Carlos Gonzalez-Gonzalez matrix inversion routine is quite interesting. However, it has been verified for matrices as large as 200 by 200. Would it work for larger matrices? In my work on soil-structure interaction using a C-program I succeeded in inverting a 289 by 289 matrix but failed when the matrix size increased to 1089 by 1089 that is obtained when the grid size changes from 16 by 16 to 32 by 32. Could Prof. Gonzalez be kind enough to try his routine on a random 1089 by 1089 matrix?
Dear Amaechi, I tried to rise the size to 500 but the MATLAB version I have is only able to manage old EXCEL editions (xls) with a maximum of 256 columns and 65000 lines. The use of EXCEL is not obligatory, is only a question of facility for data entering. It can be used an alternative file format less restrictive like .dat or .txt for larger matrices, but in these cases a convertion rutine from EXCEL to such other data formats is to be needed. In the next days I will try to do this. Also, I will do something in the sense of the answer of Daniel, because he reported times and errors better as mines; perhaps we are near of the limits that could be reached with up-to-day PC's for this problem.
There is an our article in (http://www.sciencedirect.com/science/article/pii/S0898122107000922 ). In that study, we have given an iterative decreasing dimension method (IDDM) which decreases by one dimension at every step for solution of a system of linear algebraic equations without any pre-processing and an iterative decreasing dimension algorithm (IDDA) which is based on this method.
In our other article (http://edergi.sdu.edu.tr/index.php/sdufd/article/viewFile/579/671, but Turkish), we have developed a method to calculate inverse of given regular matrix with Iterative Decreasing Dimension Method (IDDM). We have given a algorithm which is based on this method and Maple program.
Dears Amaechi and Daniel. Please, see the attachment in which is demonstrated how is possible to obtain the inverse of almost all matrix (up to 1089x1089) in 2 s with an error of ~ 10exp(-11) The programm to do this is also attached. Dear Kemal. your algorithm and method may be very good, but numerical results you have to show in order to compare.
For that size of matrices I would recommend you to write your own code using C/C++ or Fortran and then to speed-up using BLAS functions especially for the rank-one updates (dger function) and the rows interchange (swap function). An easier alternative could be to use LAPACK only for square matrices because their functions to solve rectangular systems are expensive.
Obviously I was not saying to compute directly the inverse but factorize by means of LU method with pivoting technique. That is to say, if your system is Ax=B, by pivoting you will have PAx=PB. Applying the LU factorization you will have LUx=PB. Then, using MATLAB inverse bar (for notation only), your solution will be x=U\(L\PB).
Bueno Andrés, yo estoy un poco liado con la docencia aquí en Angola, pero si te interesa este tema le podemos meter mano entre los dos y quizás también introduzca a algún profesor joven angolano de Mat. Numérica para que se vaya vinculando a investigaciones. Qué te parece? Yo tengo buenos amigos allá en la Autónoma.
Absent any further information on the class of matrices, nothing really interesting can be said. A 50 x 50 matrix isn't ``large'' by today's standards: in double precision, i.e. 8 bytes per entry, it represents 20 kbytes-for real entries, or 40 kbytes for complex entries.
So size isn't a problem. Of course it's necessary to have an idea of how much larger the size can become.
The structure of the matrix is of relevance, as has been mentioned-but, also, its spectral properties: do we expect these matrices to have eigenvalues of very small absolute value (almost zero)? Here's where the physics of the problem comes in and plays a vital role.
Se puede pensar en usar la tarjeta gráfica (GPU) para el procesamiento del cálculo de la inversa (solución del sistema lineal), y así acelerar el tiempo de cómputo del resto del proceso. Tambien usar estructuras paralelas como MPI o openMP, ya que se sabe que este problema posee una baja escalabilidad...
There several published algorithms on ACM Digital Library for large sparse matrices. But i have never used and tested them. You should take a look if you have programming skills.
There are generally two classes of optimization possible: based on knowledge and based on tech.
Based on knowledge:
- try exploiting the structure or properties of the matrix (if any) (e.g. symmetry, sparse structure, etc.)
- try exploiting the knowledge about rhe underlying problem (e.g. multi grid merhods in FEM exploit a smoothing property of the problem, FFT exploits devide and conquer)
- exploit other knowledge (e.g. does the matrix need to be inverted in each step or is some property constant?)
Technical
- try using optimized libraries (but I fear BLAS would not help with the given size of the problem)
- try parallelization on CPU or "hyper parallelization" through CUDA or OpenCL (GPU)
- try using or interfacing a more perfomant language (C or FORTRAN)
You didn't say whether you're doing numerical or symbolic work. If it's numerical, there are many tools to try, e.g. below I show how to do this in R. The computation for a 50x50 matrix takes just 0.2 milliseconds on a simple laptop. I also show below how to check a few values, and if 15 digits is enough for you, then R or matlab will be sufficient. (PS: sorry the formatting is odd. I don't know how to set monospace font in research gate comments.)
Yes, 50x50 is not large, it is actually very small. Moreover, except for possible rare occurrences, you do not compute inverses of matrices! In 99% of all cases where an inverse appears it is just meant as "solving a linear system AX=B for X". Although theoretically this still involves the inverse, numerically this is something different. There is a vast variety of powerful direct methods available for this task.
A good reference might be http://www.nr.com/. The efficiency of calculating inverse of large matrix quite depends on the exact properties of the matrix. You may provide more detailed information about your matrix.
Abdellah Rabbah: In some cases it is not possible to avoid matrix inverse computation. Consider the matrix equation [K_s + Inv(K_f) ]d = P
Where Inv(K_f) is the inverse of matrix K_f and it is required to compute the vector d . I can solve this equation for d using Gaussian elimination without finding the inverse of the compounded matrix K_s +Inv(K_f).
Please could you tell me how to determine d without first computing Inv(K_f)?
Amaechi Anyaegbunam could solve that is (K_f K_s+I) d = K_f P where I is the identity matrix. Whether the matrix multiplication is an improvement over inverting K_s might depend on specifics of these matrices.
ATTN! Daniel Lichtblau: Your answer is appreciated. It is theoretically correct and is certain to give the correct result. Your equation implies that one (unity) is to be added to all the diagonal entries of the matrix product K_f K_s. However, is it not possible for the solution of (K_f K_s+I) d = K_f P by Gaussian elimination to suffer instability if the diagonal entries of the matrix product K_f K_s were to be of the order of unity?
Amaechi Anyaegbunam: Certainly what I mentioned could have numeric stability problems. Also the original version, which requires an inversion, might have stability problems e.g. if K_f is ill-conditioned.
One advantage to the modification I showed is that it is likely to be better if one is using an indirect method such as conjugate gradient. In that case multiplications would only be of the matrix-times-vector type. Also depending on the particulars of the matrices involves, preconditioning might be an option.
Your solution can be found with the Kidder's Method by using the expansion of the inverse of the matrix : [G]=[ [ Ks*Kf ] + [ I ] ] when multiplying your system by [Kf] where {d}=[Ginv]*[Kf]{ P} : References : Kidder / Golub / Moussaoui. I have used in the development and tested in my last research work a second order term in the cited expansion. If [ A ]= [ Ks*Kf ] you will get [Ginv]=[ I ]-[ A ]+[ A ]^2 + .......
You need to exploit the properties of your matrix. See for example, if he matrix is sparce, band or symmetric definite positive. See also if it is not ill-conditionned. Check also if you really need to inverse matrices or just solve some algebraic equations. For this you have iterative and direct methods. Also a matrix of order 50 is not big unless your are looking for a technique that can be applied also to solve the same problem for matrices of order of 1000 when for example you reduce some problem parameters like dicretization step-seize.
If you have a positive-definite matrix, you can use the Cholesky decomposition, so you get a lower triangular matrix and upper triangular matrix. If you have to solve a linear equation, using this decompossition you don´t need the inverse.
Is the some iterative method for searching the inverse of matrix available? For example, do we can try to extend the Newton's iterative method for computing the inverse of a given matrix, refer to Hou-Biao Li et al. Chebyshev-type methods and preconditioning techniques, Appl. Math. Comput. 2011.
Also can use algorithms of kind fixed point, where can iteratively approximate de inverse matrix using elemental operations over de original matrix without modifications.
It is better to redesign the algorithm to reduce the rows, as example in a 3d simulation with one axis of symetry by translate to a 2D model.
Here is a opencl method: https://github.com/sschaetz/nvidia-opencl-examples/blob/master/OpenCL/src/oclTranspose/transpose.cl
Unfortunately opencl/cuda uses usually floats that generates rounding errors if the algorithm adds large arrays elements one to last instead using a tree system
computing the inverse of a matrix (when it exists) is difficult, because the underlying mathematics is not very nice (small matrix changes can imply VERY LARGE inverse changes) and not conformant with our usual intuitions about continuity
that said, it is essential for the important computations for Kalman filters (and many variants thereof), so it cannot simply be avoided in many cases - there are versions of the Kalman filter equations that take this into account and use different formulas
it means that each instance of a matrix that you want to invert should be very carefully analyzed for its proximity to non-invertibility, and if the equations allow, the inverse should be avoided
but there is no fast and easy inversion algorithm (in some applications, a few seconds is WAY too long, since there are hundreds of millions of these computations)
I do agree with Christopher Landauer that computing the inverse of a matrix (when it exists) is not an easy task, because the underlying mathematics is not very nice.
There are several methods of finding the inverse of a matrix. However, none of them is as easy as we expect. Mathematics researchers have not yet been able to discover/invent a method, which is as easy as we expect, of finding the inverse of a matrix
So far my knowledge is concerned, I assess that determination of inverse of a matrix by the application of elementary transformations is the easiest among the available difficult methods of finding the inverse of matrix provided the inverse exists.
Gauss Jordan is one of the algorithms which can be used to inverse of large and complex matrices. In MATLAB and FORTRAN we can get inverse of large matrix.