Open
Description
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