PolyMathOrg / PolyMath

Scientific Computing with Pharo
MIT License
170 stars 39 forks source link

Implement Moore-Penrose pseudoinverse for PMMatrix #260

Open olekscode opened 2 years ago

olekscode commented 2 years ago

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 commented 2 years ago

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 commented 2 years ago

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

hemalvarambhia commented 2 years ago

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