PolyMath icon indicating copy to clipboard operation
PolyMath copied to clipboard

Implement Moore-Penrose pseudoinverse for PMMatrix

Open olekscode opened this issue 3 years ago • 3 comments

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

olekscode avatar Apr 26 '22 18:04 olekscode

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

olekscode avatar May 05 '22 12:05 olekscode

Great job! I know this might be difficult, but do we have a way to test it?

SergeStinckwich avatar May 06 '22 01:05 SergeStinckwich

@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.

hemalvarambhia avatar Sep 05 '22 10:09 hemalvarambhia