Inverting a Matrix with MatrixOp and Matrixinverse
sevastionger
Mon, 03/16/2015 - 02:35 am
My problem seems quite simple but I just don't manage to get the solution.
My main goal is to get the inverse of a Matrix which is at least 1000x1000 big. Each row represents the x-value going from -1 to 1, columns represent y-values, also going from -1 to 1.
I fill the matrix with the derivative of the Fermi-Distribution which so far works fine.
The problems appear, when I try to invert the Matrix and check, if the Product of the Inverse Matrix with the Initial Matrix gives the Unity Matrix.
In general it shouldnt matter how I multiply so A*A^-1 = A^-1 * A = 1
When using the Matrixop Inv() and check the product for unity i get absolute rubbish, no matter if I multiply from left or right.
Using MatrixInverse A, however, turns out to give me a reasonable inverted Matrix BUT it matters how I multiply them. Only multiplying the inverse from the right works out, so 1 = A* A^-1 switching the two matrices gives a complete different result.
Why do I get no reasonable result at all for the Matrixop operation and why does the way I multiply them matter?
Thank you for your help.
This is how I filled the square Matrix, coef is just a coeficient wave, range defines the x range
function fermimatrix(sizematrix, range, matrixname)
Variable sizematrix, range, coef
String matrixname
//Wave squarematrix = $matrixname
coef = 0.01
Make/O/N = (sizematrix,sizematrix) $matrixname // creates empty matrix
Wave squarematrix = $matrixname
setscale x,-range, range, squarematrix // scales x
setscale y,-range, range, squarematrix // scales y
//x = energy
//y = voltage
squarematrix[][] = (1/(exp((x+y)/coef)+1)^2*exp((x+y)/coef)/coef)+0.000001 // fills matrix with derivative of fermi dist; adding 1e-06 allows inverting
end
fermimatrix(coef, 1000, 1, square)
MatrixInverse square
Matrixop/o checkunity1 = square x M_inverse // WORKS
Matrixop/o checkunity2 = M_inverse x square // DOESNT WORK
Matrixop/o inverse = Inv(square)
Matrixop/o checkunity3 = inverse x square // Both don't work
Matrixop/o checkunity4 = square x inverse
Variable sizematrix, range, coef
String matrixname
//Wave squarematrix = $matrixname
coef = 0.01
Make/O/N = (sizematrix,sizematrix) $matrixname // creates empty matrix
Wave squarematrix = $matrixname
setscale x,-range, range, squarematrix // scales x
setscale y,-range, range, squarematrix // scales y
//x = energy
//y = voltage
squarematrix[][] = (1/(exp((x+y)/coef)+1)^2*exp((x+y)/coef)/coef)+0.000001 // fills matrix with derivative of fermi dist; adding 1e-06 allows inverting
end
fermimatrix(coef, 1000, 1, square)
MatrixInverse square
Matrixop/o checkunity1 = square x M_inverse // WORKS
Matrixop/o checkunity2 = M_inverse x square // DOESNT WORK
Matrixop/o inverse = Inv(square)
Matrixop/o checkunity3 = inverse x square // Both don't work
Matrixop/o checkunity4 = square x inverse
March 16, 2015 at 02:36 am - Permalink
First, a curious look at your data shows the 1e-6 offset on an SP wave. That is usually not a good idea because the general LAPACK accuracy in SP would not be very far from that value so I recommend working with DP waves only.
As a second step I'd look at the SVD of the matrix. You can use MatrixSVD or use the /D flag of the MatrixInverse operation. In either case you will find that your eigenvalues fall off very rapidly so that even if you did not know the full details of the construction of your matrix you would be able to convince yourself that the inversion of this matrix is not numerically stable.
I would go back to the theoretical equations that you are trying to solve and rewrite them in a way that is more conducive for numerical computation. In practice you will find that matrix inverse are not used as often as they appear in theoretical manipulations.
A.G.
WaveMetrics, Inc.
March 16, 2015 at 10:16 am - Permalink