QR Decomposition Method Object Refactoring (#288)

* Added an example discovered in wikipedia that we can use to get a deeper understanding.

* Added a test to see how QR Decomposition works and then began migrating logic to the new method object.

* Clarified the classification of a method.

* Copied the logic to a more sensible place.

* Replaced logic with a method that does the same thing.

* Made the code more intention-revealing.

* Extracted a computation to an intention revealing variable.

* Broke up the computation of the indices. The computation i appears to be a matrix of minor.

* Extracted a computation to a more sensible place in an intention-revealing method.

* Added a test that computes the minor of a square matrix.

* Added another assertion to make logic of method clear.

* Added another assertion to raise confidence.

* Another another interesting test where we compute the minor of a non-square matrix.

* Improved code formatting, replaced duplicate code with a message send.

* Inlined some local variables, made the decomposeWithPivot method look more similar to decompose.

* Used a more sensible name, derived from the wikipedia entry.

* Improved the name of the temporary variable.

* Added a test that gets to the route of the problem.

* Promoted the Q and R to state variables as they are duplicated in the decomposition methods.

* Extracted a pivot temporary variable.

* Inproved the names of some temporary variables.

* Improved the name of the temporary variables.

* Corrected the name of the temporary variable.

* Use an already existing intention-revealing message.

* Moved code to a more sensible place.

* Inproved the names of some variables.

* Improved the name of a variable.

* Improved the name of a temporary variable.

* Extracted logic to a more sensible place, enhancing encapsulation.

* Made similar code even more similar.

* Made similar code even more similar.

* Made similar code even more similar.

* Made the name of the local variable the same.

Author: hemalvarambhia

src/Math-Core/PMVector.class.st

@@ -220,6 +220,11 @@ PMVector >> isReal [
 	^ self allSatisfy: [ :each | each isRealNumber ].
 ]
 
+{ #category : #testing }
+PMVector >> isZero [
+ ^ self allSatisfy: [ :element | element = 0 ]
+]
+
 { #category : #operation }
 PMVector >> log [
 	"Apply log function to every element of a vector"

src/Math-Matrix/PMMatrix.class.st

@@ -660,6 +660,14 @@ PMMatrix >> lupInverse [
 		ifNotNil: [ :i | ^ self class rows: i ]
 ]
 
+{ #category : #operation }
+PMMatrix >> minor: rowIndex and: columnIndex [
+
+	^ PMMatrix rows:
+		  ((self rows allButFirst: columnIndex) collect: [ :aRow | 
+			   aRow allButFirst: rowIndex ])
+]
+
 { #category : #'as yet unclassified' }
 PMMatrix >> mpInverse [
 	"Moore Penrose Inverse. "
@@ -883,6 +891,12 @@ PMMatrix >> rowsDo: aBlock [
 	^ rows do: aBlock
 ]
 
+{ #category : #iterators }
+PMMatrix >> rowsWithIndexDo: aBlock [
+
+	^ rows withIndexDo: aBlock
+]
+
 { #category : #transformation }
 PMMatrix >> scaleBy: aNumber [
 

src/Math-Matrix/PMQRDecomposition.class.st

@@ -6,7 +6,9 @@ Class {
 	#superclass : #Object,
 	#instVars : [
 		'matrixToDecompose',
-		'colSize'
+		'colSize',
+		'r',
+		'q'
 	],
 	#category : #'Math-Matrix'
 }
@@ -18,7 +20,7 @@ PMQRDecomposition class >> of: matrix [
 	^ self new of: matrix
 ]
 
-{ #category : #private }
+{ #category : #arithmetic }
 PMQRDecomposition >> decompose [
 
 	"
@@ -27,25 +29,23 @@ that describes the mechanics:
 https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
 "
 
-	| identityMatrix q r householderVector i |
-	r := self initialRMatrix.
-	q := self initialQMatrix.
+	| identityMatrix householderVector i matrixOfMinor |
 	identityMatrix := PMSymmetricMatrix identity: colSize.
 	1 to: self numberOfColumns do: [ :col | 
-		| columnVectorFromRMatrix householderMatrix |
+		| columnVectorFromRMatrix householderMatrix v |
 		columnVectorFromRMatrix := r columnVectorAt: col size: colSize.
 		householderVector := columnVectorFromRMatrix householder.
-		i := (PMVector zeros: col - 1) , (householderVector at: 2).
+		v := (PMVector zeros: col - 1) , (householderVector at: 2).
 		householderMatrix := identityMatrix
 		                     -
-		                     ((householderVector at: 1) * i tensorProduct: i).
+		                     ((householderVector at: 1) * v tensorProduct: v).
 		q := q * householderMatrix.
-		i := PMMatrix rows:
-			     ((r rows allButFirst: col - 1) collect: [ :aRow | 
-				      aRow allButFirst: col - 1 ]).
-		i := i - ((householderVector at: 2) tensorProduct:
-			      (householderVector at: 1) * (householderVector at: 2) * i).
-		i rows withIndexDo: [ :aRow :index | 
+		matrixOfMinor := r minor: col - 1 and: col - 1.
+		matrixOfMinor := matrixOfMinor
+		                 - ((householderVector at: 2) tensorProduct:
+				                  (householderVector at: 1)
+				                  * (householderVector at: 2) * matrixOfMinor).
+		matrixOfMinor rowsWithIndexDo: [ :aRow :index | 
 			aRow withIndexDo: [ :n :c | 
 				r
 					rowAt: col + index - 1
@@ -54,7 +54,7 @@ https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
 							 ifTrue: [ 0 ]
 							 ifFalse: [ n ]) ] ] ].
 	i := 0.
-	[ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ] whileTrue: [ 
+	[ (r rowAt: colSize) isZero ] whileTrue: [ 
 		i := i + 1.
 		colSize := colSize - 1 ].
 	i > 0 ifTrue: [ 
@@ -65,6 +65,73 @@ https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
 	^ Array with: q with: r
 ]
 
+{ #category : #arithmetic }
+PMQRDecomposition >> decomposeWithPivot [
+
+	| identityMatrix householderVector i v vectorOfNormSquareds rank mx pivot matrixOfMinor |
+	vectorOfNormSquareds := matrixToDecompose columnsCollect: [ 
+		                        :columnVector | columnVector * columnVector ].
+	mx := vectorOfNormSquareds indexOf: vectorOfNormSquareds max.
+	pivot := Array new: vectorOfNormSquareds size.
+	rank := 0.
+	identityMatrix := PMSymmetricMatrix identity: colSize.
+	[ 
+	| householderMatrix |
+	rank := rank + 1.
+	pivot at: rank put: mx.
+	r swapColumn: rank withColumn: mx.
+	vectorOfNormSquareds swap: rank with: mx.
+	householderVector := (r columnVectorAt: rank size: colSize)
+		                     householder.
+	v := (PMVector zeros: rank - 1) , (householderVector at: 2).
+	householderMatrix := identityMatrix
+	                     -
+	                     ((householderVector at: 1) * v tensorProduct: v).
+	q := q * householderMatrix.
+	matrixOfMinor := r minor: rank - 1 and: rank - 1.
+	matrixOfMinor := matrixOfMinor
+	                 - ((householderVector at: 2) tensorProduct:
+			                  (householderVector at: 1)
+			                  * (householderVector at: 2) * matrixOfMinor).
+	matrixOfMinor rowsWithIndexDo: [ :aRow :index | 
+		aRow withIndexDo: [ :element :column | 
+			r
+				rowAt: rank + index - 1
+				columnAt: rank + column - 1
+				put: ((element closeTo: 0)
+						 ifTrue: [ 0 ]
+						 ifFalse: [ element ]) ] ].
+	rank + 1 to: vectorOfNormSquareds size do: [ :ind | 
+		vectorOfNormSquareds
+			at: ind
+			put:
+			(vectorOfNormSquareds at: ind)
+			- (r rowAt: rank columnAt: ind) squared ].
+	rank < vectorOfNormSquareds size
+		ifTrue: [ 
+			mx := (vectorOfNormSquareds
+				       copyFrom: rank + 1
+				       to: vectorOfNormSquareds size) max.
+			(mx closeTo: 0) ifTrue: [ mx := 0 ].
+			mx := mx > 0
+				      ifTrue: [ 
+				      vectorOfNormSquareds indexOf: mx startingAt: rank + 1 ]
+				      ifFalse: [ 0 ] ]
+		ifFalse: [ mx := 0 ].
+	mx > 0 ] whileTrue.
+	i := 0.
+	[ (r rowAt: colSize) isZero ] whileTrue: [ 
+		i := i + 1.
+		colSize := colSize - 1 ].
+	i > 0 ifTrue: [ 
+		r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
+		i := q numberOfColumns - i.
+		pivot := pivot copyFrom: 1 to: i.
+		q := PMMatrix rows:
+			     (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
+	^ Array with: q with: r with: pivot
+]
+
 { #category : #private }
 PMQRDecomposition >> initialQMatrix [
 
@@ -88,4 +155,6 @@ PMQRDecomposition >> of: matrix [
 
 	matrixToDecompose := matrix.
 	colSize := matrixToDecompose numberOfRows.
+		r := self initialRMatrix.
+	q := self initialQMatrix.
 ]

src/Math-Tests-Matrix/PMMatrixTest.class.st

@@ -652,6 +652,41 @@ PMMatrixTest >> testMatrixTrace [
 	self assert: a tr equals: 15
 ]
 
+{ #category : #'linear algebra' }
+PMMatrixTest >> testMinorsOfNonSquareMatrix [
+
+	| matrix expected |
+	matrix := PMMatrix rows: { 
+			          { 1. 2 }.
+			          { 6. 7 }.
+			          { 5. 4 } }.
+
+	expected := PMMatrix rows: { 
+			            { 7 }.
+			            { 4 } }.
+	self assert: (matrix minor: 0 and: 0) equals: matrix.
+	self assert: (matrix minor: 1 and: 1) equals: expected
+]
+
+{ #category : #'linear algebra' }
+PMMatrixTest >> testMinorsOfSquareMatrix [
+
+	| matrix expected |
+	matrix := PMMatrix rows: { 
+			          { 1. 2. 3 }.
+			          { 6. 7. 8 }.
+			          { 5. 4. 1 } }.
+
+	expected := PMMatrix rows: { 
+			            { 7. 8 }.
+			            { 4. 1 } }.
+	self assert: (matrix minor: 0 and: 0) equals: matrix.
+	self assert: (matrix minor: 1 and: 1) equals: expected.
+	self
+		assert: (matrix minor: 2 and: 2)
+		equals: (PMMatrix rows: { { 1 } })
+]
+
 { #category : #tests }
 PMMatrixTest >> testOnesMatrix [
 	| a |

src/Math-Tests-Matrix/PMQRTest.class.st

@@ -27,6 +27,39 @@ PMQRTest >> assert: inverse isMoorePenroseInverseOf: aMatrix [
 		closeTo: aMatrix mpInverse transpose.
 ]
 
+{ #category : #tests }
+PMQRTest >> testDecompositionOfMatrixCausingErraticFailure [
+
+	| a qrDecomposition matricesAndPivot q r expectedMatrix pivot |
+	a := PMSymmetricMatrix rows:
+		     #( #( 0.41929313699681925 0.05975350554089691
+		           0.2771676258543356 0.35628773381760703 )
+		        #( 0.05975350554089691 0.12794227252152854
+		           0.3257742693302102 0.28814463284245906 )
+		        #( 0.2771676258543356 0.3257742693302102 0.8468441832097453
+		           0.9101872061892353 )
+		        #( 0.35628773381760703 0.28814463284245906
+		           0.9101872061892353 0.5163744224777326 ) ).
+
+	qrDecomposition := PMQRDecomposition of: a.
+	matricesAndPivot := qrDecomposition decomposeWithPivot.
+
+	expectedMatrix := PMMatrix rows:
+		                  #( #( 0.2771676258543356 0.35628773381760703
+		                        0.41929313699681925 0.05975350554089691 )
+		                     #( 0.3257742693302102 0.28814463284245906
+		                        0.05975350554089691 0.12794227252152854 )
+		                     #( 0.8468441832097453 0.9101872061892353
+		                        0.2771676258543356 0.3257742693302102 )
+		                     #( 0.9101872061892353 0.5163744224777326
+		                        0.35628773381760703 0.28814463284245906 ) ).
+	q := matricesAndPivot at: 1.
+	r := matricesAndPivot at: 2.
+	pivot := matricesAndPivot at: 3.
+	self assert: q * r closeTo: expectedMatrix.
+	self assert: pivot equals: #( 3 4 3 nil )
+]
+
 { #category : #tests }
 PMQRTest >> testHorizontalRectangularMatrixCannotBeDecomposed [
 
@@ -169,6 +202,45 @@ PMQRTest >> testRank [
 			self assert: matrix transpose rank equals: 3 ] ].
 ]
 
+{ #category : #tests }
+PMQRTest >> testSimpleQRDecomposition [
+
+	| a qrDecomposition decomposition |
+	a := PMMatrix rows: { 
+			     { 12. -51. 4 }.
+			     { 6. 167. -68 }.
+			     { -4. 24. -41 } }.
+
+	qrDecomposition := PMQRDecomposition of: a.
+
+	decomposition := qrDecomposition decompose.
+	decomposition first * decomposition second.
+	self assert: decomposition first * decomposition second equals: a
+]
+
+{ #category : #tests }
+PMQRTest >> testSimpleQRDecompositionWithPivot [
+
+	| a qrDecomposition decomposition expected |
+	a := PMMatrix rows: { 
+			     { 12. -51. 4 }.
+			     { 6. 167. -68 }.
+			     { -4. 24. -41 } }.
+
+	qrDecomposition := PMQRDecomposition of: a.
+
+	decomposition := qrDecomposition decomposeWithPivot.
+	decomposition first * decomposition second.
+
+	expected := PMMatrix rows: { 
+			            { -51. 4. 12 }.
+			            { 167. -68. 6 }.
+			            { 24. -41. -4 } }.
+	self
+		assert: decomposition first * decomposition second
+		closeTo: expected
+]
+
 { #category : #tests }
 PMQRTest >> testVectorHouseholder [