Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement Moore-Penrose pseudoinverse for PMMatrix #260

Open
olekscode opened this issue Apr 26, 2022 · 3 comments
Open

Implement Moore-Penrose pseudoinverse for PMMatrix #260

olekscode opened this issue Apr 26, 2022 · 3 comments

Comments

@olekscode
Copy link
Member

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
Copy link
Member Author

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

@SergeStinckwich
Copy link
Member

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

@hemalvarambhia
Copy link
Contributor

hemalvarambhia commented Sep 5, 2022

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants