Hi,
I hope you are all doing well. I am trying to find the inverse of this matrix, but as you can see, it is giving me problems, and these values are definitely not correct. I tried both PseudoInverse and Inverse but they gave me the same results, If any of you have an idea where the problem may come from it will help me a lot. Thank you all in advance.
I’m not seeing any error or warning message regarding the matrix computation (as distinct from its display). Is there something not shown that implies the inversion gave a bad result?
The “matrice” in your example is ill - conditioned; thus, the computation of its inverse matrix is unstable (e.g., a small change in a matrix causes huge changes in the inverse matrix). I believe you already know the followings but just a reminder. If a matrix, X, is ill - conditioned, the inverse of the matrix is a “meaningless bad approximation”. The condition number of the “matrice” in your example, (i.e., ratio of max singular value divided by min singular value), is gigantic (i.e., 2.6*10^10). Orthogonal matrix should have singular value of one.
Ill-conditioning is also known as multicollinearity in regression, confounding (or alias structure) in design of experiments, rank-deficient in matrix algebra. Note that the inverse of an orthogonal matrix is the same as its transpose (i.e., the U and V are orthogonal matrices, U’U=I, V’V=I.
X=U*S*V'.
I always compute Singular Value Decomposition (SVD) if a matrix computation is involved because SVD is the most reliable way to compute an inverse matrix. I checked Mathematica commands such as PseudoInverse, Inverse, LinearSolve, LeastSquares for ill-conditioning in the past and these commands do not provide a warning sometimes. Matlab has also the same issue.
{u, s, v} = SingularValueDecomposition[matrice]
ListPlot[Diagonal[s], PlotRange -> {All, {0, 60}}]
condNumber=Max[Diagonal[s]]/Min[Diagonal[s]] (* 2.6*10^10 *)
Yes there are no warnings or errors but the values of the inverted matrix are definitely incorrect (10^6, 10^8…)
Thank you for your informative response
Like you said, even in this case, using Least Squares and LinearSolve didn’t help solve the problem. I tried Tikhonov regularization
alpha = 10^(-6);
Ainv = Inverse[matrice + alpha IdentityMatrix[Length[matrice]]];
It seems to give better results, but it’s still not as precise as I wanted
I’m with Daniel. The following check returns zero, which is hard to improve upon: Norm[Inverse[matrice].matrice-IdentityMatrix[Length@matrice]]. If matrice has an eigenvalue with a modulus as small as 10^{-8}, the inverse will have one as large as 10^8. And the inverse will have to have some entries of a similar order of magnitude.
The following visualizes the magnitudes of the eigenvalues:
ListLogPlot[Abs[Eigenvalues[matrice]]]
If the small eigenvalues are the result of roundoff and are supposed to be zero, then the matrix is not invertible. If you want the pseudoinverse, then increase the Tolerance to a level you find acceptable.
PseudoInverse[matrice, Tolerance -> 1/1000]
You can use SingularValueList[matrice] and SingularValueList[matrice, Tolerance -> 1/1000] to examine what is discarded.