Added pseudoinverseOfDiagonal: and refactored the tests

Author: olekscode

src/Math-Numerical/PMLeastSquares.class.st

@@ -4,6 +4,28 @@ Class {
 	#category : #'Math-Numerical'
 }
 
+{ #category : #'as yet unclassified' }
+PMLeastSquares >> 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
+]
+
 { #category : #'as yet unclassified' }
 PMLeastSquares >> solveMatrixA: aMatrix matrixB: aMatrixOrVector [
 
@@ -22,7 +44,6 @@ PMLeastSquares >> solveMatrixA: aMatrix matrixB: aMatrixOrVector [
 	s := svd diagonalSingularValueMatrix.
 	v := svd rightSingularMatrix.
 
-	pseudoinverse := s inverse.
-
+	pseudoinverse := self pseudoinverseOfDiagonal: s.
 	^ v * pseudoinverse * u transpose * aMatrixOrVector
 ]

src/Math-Tests-Numerical/PMLeastSquaresFixture.class.st

@@ -17,29 +17,141 @@ PMLeastSquaresTest >> setUp [
 ]
 
 { #category : #tests }
-PMLeastSquaresTest >> testSolveIntelFortran [
+PMLeastSquaresTest >> testPseudoinverseOfDiagonalSquareMatrix [
+
+	| matrix expectedInverse inverse |
+	
+	matrix := PMMatrix rows: #(
+		(2 0 0 0)
+		(0 1 0 0)
+		(0 0 -3 0)
+		(0 0 0 -1)).
+		
+	expectedInverse := PMMatrix rows: {
+		{ 1/2 . 0 . 0 . 0 } .
+		{ 0 . 1 . 0 . 0 } .
+		{0 . 0 . -1/3 . 0} .
+		{0 . 0 . 0 . -1 }}.
+		
+	inverse := leastSquares pseudoinverseOfDiagonal: matrix.
+	self assert: inverse closeTo: expectedInverse.
+]
+
+{ #category : #tests }
+PMLeastSquaresTest >> testPseudoinverseOfDiagonalSquareMatrixWithZeros [
+
+	| matrix expectedInverse inverse |
+	
+	matrix := PMMatrix rows: #(
+		(2 0 0 0)
+		(0 1 0 0)
+		(0 0 0 0)
+		(0 0 0 0)).
+		
+	expectedInverse := PMMatrix rows: #(
+		(0.5 0 0 0)
+		(0 1 0 0)
+		(0 0 0 0)
+		(0 0 0 0)).
+		
+	inverse := leastSquares pseudoinverseOfDiagonal: matrix.
+	self assert: inverse closeTo: expectedInverse.
+]
+
+{ #category : #tests }
+PMLeastSquaresTest >> testPseudoinverseOfDiagonalTallMatrix [
 
-	| fixture solution |
+	| matrix expectedInverse inverse |
 
-	fixture := PMLeastSquaresFixtureIntelFortran new.
+	matrix := PMMatrix rows: #(
+		(2 0 0 0)
+		(0 1 0 0)
+		(0 0 -3 0)
+		(0 0 0 -1)
+		(0 0 0 0)
+		(0 0 0 0)).
+		
+	expectedInverse := PMMatrix rows: {
+		{ 1/2 . 0 . 0 . 0 . 0 . 0 } .
+		{ 0 . 1 . 0 . 0 . 0 . 0 } .
+		{0 . 0 . -1/3 . 0 . 0 . 0 } .
+		{0 . 0 . 0 . -1 . 0 . 0 }}.
+		
+	inverse := leastSquares pseudoinverseOfDiagonal: matrix.
+	self assert: inverse closeTo: expectedInverse.
+]
+
+{ #category : #tests }
+PMLeastSquaresTest >> testPseudoinverseOfDiagonalWideMatrix [
+
+	| matrix expectedInverse inverse |
+	
+	matrix := PMMatrix rows: #(
+		(2 0 0 0 0 0)
+		(0 1 0 0 0 0)
+		(0 0 -3 0 0 0)
+		(0 0 0 -1 0 0)).
+		
+	expectedInverse := PMMatrix rows: {
+		{ 1/2 . 0 . 0 . 0 } .
+		{ 0 . 1 . 0 . 0 } .
+		{0 . 0 . -1/3 . 0} .
+		{0 . 0 . 0 . -1 } .
+		{0 . 0 . 0 . 0 } .
+		{0 . 0 . 0 . 0 }}.
+		
+	inverse := leastSquares pseudoinverseOfDiagonal: matrix.
+	self assert: inverse closeTo: expectedInverse.
+]
+
+{ #category : #tests }
+PMLeastSquaresTest >> testSolveIntelFortran [
+	"An example of least squares system (AX = B) taken from Intel DGELSD Example Program in Fortran:
+https://www.intel.com/content/www/us/en/develop/documentation/onemkl-lapack-examples/top/least-squares-and-eigenvalue-problems/linear-least-squares-lls-problems/gelsd-function/dgelsd-example/dgelsd-example-fortran.html"
+	| matrixA matrixB expectedSolution solution |
+	
+	matrixA := PMMatrix rows: #(
+		( 0.12  -8.19   7.69  -2.26  -4.71)
+		(-6.91   2.22  -5.12  -9.08   9.96)
+		(-3.33  -8.94  -6.72  -4.40  -9.98)
+		( 3.97   3.33  -2.74  -7.92  -3.20)).
+		
+	matrixB := PMMatrix rows: #(
+		(7.30   0.47  -6.28)
+		(1.33   6.58  -3.42)
+		(2.68  -1.71   3.46)
+		(-9.62  -0.79   0.41)).
+		
+	expectedSolution := PMMatrix rows: #(
+		(-0.69  -0.24   0.06)
+		(-0.80  -0.08   0.21)
+		( 0.38   0.12  -0.65)
+		( 0.29  -0.24   0.42)
+		( 0.29   0.35  -0.30)).
 
 	solution := leastSquares
-		solveMatrixA: fixture matrixA
-		matrixB: fixture matrixB.
+		solveMatrixA: matrixA
+		matrixB: matrixB.
 
-	self assert: solution closeTo: fixture matrixX.
+	self assert: solution closeTo: expectedSolution.
 ]
 
 { #category : #tests }
 PMLeastSquaresTest >> testSolveSmallOneSolution [
-
-	| fixture solution |
+	"Small example of least squares system (AX = B) with one solution taken from here: https://textbooks.math.gatech.edu/ila/least-squares.html"
+	| matrixA vectorB expectedSolution solution |
 
-	fixture := PMLeastSquaresFixtureSmallOneSolution new.
+	matrixA := PMMatrix rows: #(
+		(0 1.1)
+		(1 0)
+		(0 -0.2)).
+		
+	vectorB := #(1.1 -1.1 -0.2) asPMVector.
+	expectedSolution := #(-1.1 1) asPMVector.
 
 	solution := leastSquares
-		solveMatrixA: fixture matrixA
-		matrixB: fixture matrixB.
+		solveMatrixA: matrixA
+		matrixB: vectorB.
 
-	self assert: solution closeTo: fixture matrixX.
+	self assert: solution closeTo: expectedSolution.
 ]