Implement Moore-Penrose pseudoinverse for PMMatrix #260

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