how to get the inverse of a symmetric matrix and ensure accurancy

22 views (last 30 days)
Pengpeng HE
Pengpeng HE on 26 Feb 2019
Commented: Walter Roberson on 1 Mar 2019
Hi,
I have a flexibility matrix (20*20), F, which is symmetric and positively defined. I need to reverse it to get the stiffness matrix, K=F_inv, and then to obtain eigenvalues using K. Theoretically speaking, K should also be symmetric and positively defined, which also results in positive real eigenvalues. However, in Matlab, the inverse of F does not yield to a symmetric matrix anymore by using K=inv(F) probability due to numerical errors. Moreover, this phenomenon leads to a few complex eigenvalues due to the non-symmetric K.
I tested norm(K*F) very close to 1. So how can I deal with this situation? I also tried to use Cholesky decomposition to get the inverse matrix instead of build-in inv. This approach can definitely provides symmetric inverse matrix of F, however, the accurancy is reduced as well. norm(F_inv*F) using Cholesky is around 1.2, and F_inv*F is close to the identity matrix, but not accurate enough. So how can get the inverse of F to be symmetrix and accurate as well? Anyone has suggestion on this? Thank you very much!
  2 Comments
Walter Roberson
Walter Roberson on 1 Mar 2019
Note that in cases of very slight asymmetry due to round off error, it is common to use
F = (F + F.')/2;
which is guaranteed to be symmetric.

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Stephan
Stephan on 26 Feb 2019
Edited: Stephan on 26 Feb 2019
Hi,
you could consider to do the calculation symbolic. Here is a small example:
A = [9 3 -6 12; 3 26 -7 -11; -6 -7 9 7; 12 -11 7 65]
B = inv(A)
C = A*B
% do the same symbolic
D = sym(A)
E = inv(D)
F = D*E
Note that maybe things are getting ugly, when using this approach on a 20x20 matrix.
Best regards
Stephan
  9 Comments
Show 7 older comments
Pengpeng HE
Pengpeng HE on 27 Feb 2019
Yes. It can work, but when I convert F_inv_sym to double by
F_inv = double(F_inv_sym);
F_inv is still assymetric. You method really works. The reason is this F matrix... Sorry about that.
Stephan
Stephan on 27 Feb 2019
Edited: Stephan on 27 Feb 2019
No problem - your question teached me something about the Symbolic Toolbox ;-)
Even if it does not solve your problem, if you liked the answer, please accept it.

Sign in to comment.

More Answers (2)

John D'Errico
John D'Errico on 27 Feb 2019
Edited: John D'Errico on 27 Feb 2019
You want to compute the inverse of a symmetric matrix F. So let me first check your claim...
norm(F,inf)
ans =
9.7249e-06
norm(F-F',inf)
ans =
2.5801e-09
F is not in fact symmetric in the first place! In fact, given that norm(F) is pretty small, the norm of F-F' is embarassingly large in comparison. So F is not even remotely close to being a symmetric matrix.
That the inverse of F is not symmetric is not even a remote surprise, since F is not in fact symmetric! I think you are trying to fix the problem from the wrong end. Start by learning why it is that F is NOT in fact symmetric. Once you fix F, then you can worry about what you do with it.
  4 Comments
Show 2 older comments
Pengpeng HE
Pengpeng HE on 27 Feb 2019
Yes. I found it. Elements in F matrix are computed using integral. The default absolute error tolerance of 'integral' is 1e-10, however, the order of magnitudes of F is around 1e-6 ~ 1e-10, which is very close to the AbsTol. I changed AbsTol and RelTol to be much smaller and it works now!
John D'Errico
John D'Errico on 1 Mar 2019
Very good!
Easy to not think about the tolerances there, not realizing the tolerance used can have an impact on other computations down the line. In this case, the tolerance was huge in context.
I've seen similar problems arise, usually whenever two different numerical methods are used, one within the other. So anything that employs a tolerance.
A typical example is where someone wants to optimize a result that includes a numerical integration inside. Since the optimization uses a gradient computation internally, small errors in the objective function can become highly important. Now the objective function no longer appears to be a continuous, differentiable function of the inputs.

Sign in to comment.

Bjorn Gustavsson
Bjorn Gustavsson on 27 Feb 2019
If you want an inverse that is nice in that sense you could always simply do an eigenvalue decomposition of F. Then the eigenvalues of the inverse of F is simply the inverse of the eigenvalues of F, and you'll also be able to produce an invF that is symmetric and positive definite (within finite precision numerical accuracy):
[V,D] = eig(F);
Finv = V*inv(D)*V.';
D_of_inv_F = inv(D);
HTH
  1 Comment
Pengpeng HE
Pengpeng HE on 27 Feb 2019
Thank you, Bjorn. I tested the F matrix and found that the issue should be caused by the very slight asymmetry of F although it should be symmetric. Now I replace the lower triangular part of F by its upper one to make it completely symmetric. This can work now. I also attach the F and M matrix in the original post for you to test if you are interested.
Regards,
Pengpeng

Sign in to comment.

Categories

Find more on Matrix Indexing in Help Center and File Exchange

Tags

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!