PolyMath
PolyMath copied to clipboard
Implement Moore-Penrose pseudoinverse for PMMatrix
https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse
A simple (but slightly time-consuming) implementation can be based on SVD:
PMMatrix >> pseudoinverse
| svd u s v sPseudoinverse |
svd := PMSingularValueDecomposition decompose: matrix.
u := svd leftSingularMatrix.
s := svd diagonalSingularValueMatrix.
v := svd rightSingularMatrix.
sPseudoinverse := self pseudoinverseOfDiagonal: s.
^ v * sPseudoinverse * u transpose
PMMatrix >> pseudoinverseOfDiagonal: aMatrix [
"To get pseudoinverse of a diagonal rectangular matrix, we take reciprocal of any no-zero element of the main diagonal, leaving all zeros in place. Then we transpose the matrix."
| pseudoinverse diagonalSize |
"Rows become columns and columns become rows because we transpose"
pseudoinverse := PMMatrix
zerosRows: aMatrix numberOfColumns
cols: aMatrix numberOfRows.
"The size of the main diagonal of a rectangular matrix is its smallest dimension"
diagonalSize := aMatrix numberOfRows min: aMatrix numberOfColumns.
"Inverting the elements on the main diaginal"
1 to: diagonalSize do: [ :i |
pseudoinverse at: i at: i put: ((aMatrix at: i at: i) = 0
ifTrue: [ 0 ] ifFalse: [ 1 / (aMatrix at: i at: i) ]) ].
^ pseudoinverse
There is this but it has to be improved:
PMMatrix >> mpInverse
"Moore Penrose Inverse. "
|f g|
self numberOfRows < self numberOfColumns
ifTrue:[ f := self transpose qrFactorizationWithPivoting.
g := f first.
f := f second inversePivotColumns: (f at:3) ]
ifFalse: [ f := self qrFactorizationWithPivoting.
g := (f second inversePivotColumns: (f at:3)) transpose.
f := f first transpose ].
^ g * ((f *self *g) inverse) *f
Great job! I know this might be difficult, but do we have a way to test it?
@olekscode @SergeStinckwich the tests for these (MP Inverse) fail because of QR Decomposition with pivoting. I am attempting to fix it with #288. Specifically, the Pivot has nils in it for some matrices and I focusing on why very closely.