I need to multiply the inverse of a matrix A to another matrix B. If B were to be a vector, I would simply solve a linear system Ax = B to get the product x=inv(A)B. With B being a matrix, I don't know the most efficient way|method to compute inv(A)B.
Kindly share your experience.
You can use MATLAB : https://stackoverflow.com/questions/15374441/matlab-inverse-of-large-matrixhttps://www.reddit.com/r/matlab/comments/cb6mz4/inverse_of_a_large_sparse_matrix_in_matlab/
As your system matrix most probably is sparse, you might be better off with numerical linear algebra routines especially designed for large, sparse matrices. I did stuff on that in the late 80s, through the Cray-2 machine, but there must be more streamlined ways to do it nowadays. (In those days one had to "cut" the lengths of the vectors in order to match the smartest size of sub-vectors for the linear algebra routines, but nowadays that will most surely be automatic, to a degree, nowadays.) Talk with numerical linear algebra experts, and they will help you.
In any case, you will definitely not create the inverse of the system matrix - it's a waste of time.
Thank you, Michael and Peter. Yes, the matrices are sparse, large, tridiagonal, and arising from finite difference of an unsteady PDE.
Now, the full scheme which I want to implement is of the form:
( E- kA^{-1}B ) u^{n+1} = u^n + ( E - k A^{-1}B ) G( w),
with
w = ( E - (t^{n+1}-s) A^{-1}B ) A^{-1}b,
where
u is a vector of solution, n is the time level, k = time step.
A^{-1} is the inverse of matrix A, and B is also a matrix.
E is the identity matrix with the same dimensions as A and B,
b is a vector and G() is a vector function of the indicated argument, and
s depends on t^n and t^{n+1}.
If possible, kindly suggest a way to avoid computing the inverse of A.
I code in C++ and only interested in/use C++ libraries, not Matlab or any of those interpreted languages.
Thanks.
Hi: i think that the best way is using the "Product Form of Inverse - PFI", showed in Advanced Linear-Programming computing Techniques, by William Orchard Hays.
As your matrix is sparse, you could use an iterative method like the generalized minimal residual method (GMRES) which is implemented in many library of linear algebra.
Seeing your goal was to solve a PDE, which is discretized into a system where
matrix representation are applied. Since the usual numerical PDEs can be reformulated via large-scale sparse matrices along with vectorization (discrete truncation) of the nonlinear or inhomogeneous terms when available, C plus compiled with Matlab are typically a good efficient way to implement the "iteration algorithm". (this will often be subject to the specific solution you try to construct.) I agree with the iterative numerics applying "asymptotic convergence (adaptive Schulz algorithm) for solving it with the generalized inverse" which seems efficient and flexible algorithm.
For more advanced development in the immensely broad area, I would recommend the following article and the references therein
https://www.pnas.org/content/110/17/6634
Hi Chinedu,
Since you mentioned that your matrices are tridiagonal, you can use the Thomas algorithm. You will get your solution in O(n) operations.
Best regards,
Timur
Chinedu Nwaigwe It definitely seems that you do not need to compute the inverse of matrix A explicitly. You can solve for AX=B and Ax=b using the numerical linear algebra packages in C++ (SuperLU, MUMPS, etc are available packages) and obtain X and x. Then you can rewrite your equations as,
( E- k X ) u^{n+1} = u^n + ( E - k X ) G( w)
with, w = ( E - (t^{n+1}-s) X ) x
@Anoop, yes, I already sorted it that way using the eigen library.
Thanks.
Chinedu thanks. The world has advanced a great deal with the advent of computers and software development. With the availability of computerized ways of solving this problem it is difficult to find the most efficient technique. The parameters involved are computer time involved and user friendliness.
You can do the following :
A(m by n )X(n by p) = B(m by p)
Applying vec operation on both sides we get :
kron(Identity,A)(mp by np)*vecX(np by 1) = vecB (mp by 1) which is a matrix- vector linear system and can be solved by gaussian elimination and related algorithms.
Kron(U,V) stands for kronecker product U × V, the vec operation converts a matrix Z into a vector vec Z by column wise stacking one on top of the other, the topmost vector block is the first column of matrix Z, 2nd block is the 2nd column and the last column of Z is the lowest vector block in vec Z.
solving for vec X, you can reconstruct X by restructuring the solution vector blocks as respective columns in matrix X.
https://www.researchgate.net/post/How_to_find_the_inverse_of_a_large-scale_sparse_matrix_in_an_efficient_way
- Badges
- Science method