Skip to content

Implement Moore-Penrose pseudoinverse for PMMatrix #260

Description

@olekscode

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

Activity

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

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions