Merge pull request #286 from jordanmontt/least-squares

Added Least squares implementation

Author: SergeStinckwich

src/Math-Numerical/PMLeastSquares.class.st

@@ -0,0 +1,49 @@
+"
+SVD-based implementation of a minimum-norm solution to the overdetermined least squares problem.
+"
+Class {
+	#name : #PMLeastSquares,
+	#superclass : #Object,
+	#category : #'Math-Numerical-Math-LinearAlgebra'
+}
+
+{ #category : #accessing }
+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 : #api }
+PMLeastSquares >> solveMatrixA: aMatrix matrixB: aMatrixOrVector [
+
+	"If b is a vector:
+	x' = minimize || b - Ax ||
+	
+	If B is a matrix:
+	X' = minimize || B - AX ||"
+
+	| svd u s v pseudoinverse |
+	svd := PMSingularValueDecomposition decompose: aMatrix.
+
+	u := svd leftSingularMatrix.
+	s := svd diagonalSingularValueMatrix.
+	v := svd rightSingularMatrix.
+
+	pseudoinverse := self pseudoinverseOfDiagonal: s.
+	^ v * pseudoinverse * u transpose * aMatrixOrVector
+]

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

@@ -0,0 +1,126 @@
+"
+A PMLeastSquaresTest is a test class for testing the behavior of PMLeastSquares
+"
+Class {
+	#name : #PMLeastSquaresTest,
+	#superclass : #TestCase,
+	#instVars : [
+		'leastSquares'
+	],
+	#category : #'Math-Tests-Numerical'
+}
+
+{ #category : #running }
+PMLeastSquaresTest >> setUp [
+
+	super setUp.
+	leastSquares := PMLeastSquares new
+]
+
+{ #category : #tests }
+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 [
+
+	| matrix expectedInverse inverse |
+	
+	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 >> testSolveSmallOneSolution [
+	"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 |
+	
+	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: matrixA
+		matrixB: vectorB.
+	
+	self assert: solution closeTo: expectedSolution.
+]