Matrix Inversion
Posted: Sat Jun 16, 2018 11:54 am
Needed to invert some large matrices, so rolled me some PB to do so, based on the method presented in this paper:
https://arxiv.org/ftp/arxiv/papers/1304/1304.6893.pdf
I've only verified the results for <= 16x16 matrices so far (and it yielded the right answer in each case).
In the code, I use the same 4x4 source matrix example as in the paper, and although the inverse result is correct here, it is actually NOT correct in the paper
Anyway, there's not a shred of 'advanced PB' in the code, in fact it is all pretty much 'entry level' - so hopefully very easy to understand, but it sacrifices efficiency and there's plenty of opportunity for optimization...
https://arxiv.org/ftp/arxiv/papers/1304/1304.6893.pdf
I've only verified the results for <= 16x16 matrices so far (and it yielded the right answer in each case).
In the code, I use the same 4x4 source matrix example as in the paper, and although the inverse result is correct here, it is actually NOT correct in the paper
Anyway, there's not a shred of 'advanced PB' in the code, in fact it is all pretty much 'entry level' - so hopefully very easy to understand, but it sacrifices efficiency and there's plenty of opportunity for optimization...
Code: Select all
Declare UpdateNumbersArray(pVal.l, kVal.l)
Declare DeterminePK()
Global NewList N.s() ;Dictionary List of columns
Global NewList B.s() ;Dictionary List of rows
Global Dim A.d(4,4) ;The 'Numbers Array' ie Source Matrix to be inverted. Change dimensions for different size matrix
Global Dim InvArr.d(4,4) ;Holds inverse of the Source Matrix. Change dimensions to that of Numbers Array, if different size matrix
;Change number of elements in B and N according to size of Numbers Array:
AddElement(B()) : B() = "y;1"
AddElement(B()) : B() = "y;2"
AddElement(B()) : B() = "y;3"
AddElement(B()) : B() = "y;4"
AddElement(N()) : N() = "x;1"
AddElement(N()) : N() = "x;2"
AddElement(N()) : N() = "x;3"
AddElement(N()) : N() = "x;4"
A(1,1) = 2
A(1,2) = 5
A(1,3) = 3
A(1,4) = 2
A(2,1) = 4
A(2,2) = 10
A(2,3) = 1
A(2,4) = 7
A(3,1) = 1
A(3,2) = 5
A(3,3) = 2
A(3,4) = 2
A(4,1) = 2
A(4,2) = 1
A(4,3) = 2
A(4,4) = 1
;Method for inverting a matrix per paper "A space efficient flexible pivot selection approach to evaluate determinant and inverse of a matrix":
;First determine p (row val) and k (column val):
;1. iterate through the B list until you find the first letter a 'y'. Make the position of p the value after the colon
;2. look in matrix of numbers at position (p,p). If it is non zero, then look at the entry in N at position p, and make k the number after the colon.
;If the value is zero, then keep iterating to the right: (p,p+n) until a non zero is found. then, look at the entry in N at position p+n, and make k the number after the colon
;Now that the p and k values are known:
;3. Refresh the values in lists B and N according to following formula: B = B + [xk] - [yp] and N = N - [xk] + [yp]. Add the entry at the same position as the entry removed
;4. For all values in the Numbers matrix other than those in row p and column k, calculate new values as follows:
;mi = a(i,k)/a(p,k) and then: a(i,j) = a(i,j) + a(p,j)*mi
Procedure DeterminePK()
p.l ;row position of pivot point
k.l ;column position of pivot point
ResetList(B())
ResetList(N())
;Find position of p:
Pos = 0
While NextElement(B())
Pos = Pos + 1
If StringField(B(),1,";") <> "x"
Break
EndIf
Wend
p = Pos
;Find position of k by finding first non-zero value from Numbers Array:
For t = Pos To ArraySize(A(),2)
If A(Pos,t) <> 0
Break
EndIf
Next
SelectElement(N(), t-1) ;we do t-1 because for 'SelectElement', element numbering begins at 0. So Pos=1 is actually 0 for SelectElement
k = Val(StringField(N(),2,";"))
;Now we refresh the values in lists B and N, as per rule 3:
ResetList(B())
While NextElement(B())
If StringField(B(),1,";") = "y" And Val(StringField(B(),2,";")) = p
B() = "x" + ";" + Str(k)
Break
EndIf
Wend
ResetList(N())
While NextElement(N())
If StringField(N(),1,";") = "x" And Val(StringField(N(),2,";")) = k
N() = "y" + ";" + Str(p)
Break
EndIf
Wend
UpdateNumbersArray(p,k)
EndProcedure
Procedure UpdateNumbersArray(pVal.l, kVal.l)
;The pivot point is (p,k) ie (pVal,kVal)
;Divide pivot row by value at pivot point, and pivot column by -ve of value at pivot point. Set the value at pivot point as 1 divided by the value at pivot point:
For t = 1 To ArraySize(A(),2)
If t<>kVal
A(pVal,t) = A(pVal,t)/A(pVal,kVal)
EndIf
Next
For t = 1 To ArraySize(A(),1)
If t<>pVal
A(t,kVal) = -1*A(t,kVal)/A(pVal,kVal)
EndIf
Next
A(pVal,kVal) = 1/A(pVal,kVal)
;Now calculate new values for the Numbers Array as per rule 4:
For t = 1 To ArraySize(A(),1)
For v = 1 To ArraySize(A(),2)
If t<>pVal And v<>kVal
m.d = A(t,kVal)/A(pVal,kVal)
A(t,v) = A(t,v) + (A(pVal,v)*m)
EndIf
Next
Next
EndProcedure
;If the matrix has an inverse, the number of pivot elements would be equal to the order of the matrix. So, we call DeterminePK() as many times as the order of the matrix.
;The order of the matrix is the number of rows or columns of the matrix (doesn't matter which, as it has to be square to have an inverse):
If ArraySize(A(),1) = ArraySize(A(),2)
For t = 1 To ArraySize(A(),1)
DeterminePK()
Next
EndIf
;Re-arrange the values of A() into InvArr according to the row and column sequence indicated by the dictionary lists B() and N():
For t = 1 To ArraySize(A(),1)
SelectElement(B(), t-1)
RowVal = Val(StringField(B(),2,";"))
For v = 1 To ArraySize(A(),2)
SelectElement(N(), v-1)
ColVal = Val(StringField(N(),2,";"))
InvArr(RowVal,ColVal) = A(t,v)
Next
Next
;Print out the values in InvArr. This is the inverse of Matrix A():
For t = 1 To ArraySize(InvArr(),1)
For v = 1 To ArraySize(InvArr(),2)
Debug StrD(InvArr(t,v),4) ;Change 4 to however many decimal places you require in the result
Next
Next
