Fixed pivot in QR decomposition

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

@@ -28,39 +28,6 @@ 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 [
 
@@ -151,6 +118,51 @@ PMQRTest >> testOrthogonalize [
 	i < 10 ] whileTrue
 ]
 
+{ #category : 'tests' }
+PMQRTest >> testQRDecompositionOnRankDeficientMatrix [
+
+	| a qrDecomposition reconstruction |
+	a := PMMatrix rows: {
+			     { 1. 2. 3 }.
+			     { 4. 5. 6 }.
+			     { 7. 8. 9 } }.
+
+	qrDecomposition := PMQRDecomposition of: a.
+	qrDecomposition decompose.
+
+	self assert: qrDecomposition q rank equals: a rank.
+	self assert: qrDecomposition r rank equals: a rank.
+
+	reconstruction := qrDecomposition q * qrDecomposition r.
+	self assert: reconstruction closeTo: a
+]
+
+{ #category : 'tests' }
+PMQRTest >> testQRDecompositionWithPivotOnRankDeficientMatrix [
+
+	| a qrDecomposition expectedQR expectedPivot reconstruction |
+	a := PMMatrix rows: {
+			     { 1. 2. 3 }.
+			     { 4. 5. 6 }.
+			     { 7. 8. 9 } }.
+	expectedQR := PMMatrix rows: {
+			              { 3. 1. 2 }.
+			              { 6. 4. 5 }.
+			              { 9. 7. 8 } }.
+	expectedPivot := #( 3 1 2 ).
+
+	qrDecomposition := PMQRDecomposition of: a.
+	qrDecomposition decomposeWithPivot.
+
+	self assert: qrDecomposition q rank equals: a rank.
+	self assert: qrDecomposition r rank equals: a rank.
+	self assert: qrDecomposition q * qrDecomposition r closeTo: expectedQR.
+	self assert: qrDecomposition pivot equals: expectedPivot.
+
+	reconstruction := qrDecomposition q * qrDecomposition r * qrDecomposition permutationMatrixFromPivot inverse.
+	self assert: reconstruction closeTo: a
+]
+
 { #category : 'tests' }
 PMQRTest >> testQRFactorization [
 
@@ -225,24 +237,25 @@ PMQRTest >> testSimpleQRDecomposition [
 { #category : 'tests' }
 PMQRTest >> testSimpleQRDecompositionWithPivot [
 
-	| a qrDecomposition decomposition expected |
-	a := PMMatrix rows: { 
+	| a qrDecomposition expectedQR expectedPivot reconstruction |
+	a := PMMatrix rows: {
 			     { 12. -51. 4 }.
 			     { 6. 167. -68 }.
 			     { -4. 24. -41 } }.
+	expectedQR := PMMatrix rows: {
+			              { -51. 4. 12 }.
+			              { 167. -68. 6 }.
+			              { 24. -41. -4 } }.
+	expectedPivot := #( 2 3 1 ).
 
 	qrDecomposition := PMQRDecomposition of: a.
+	qrDecomposition decomposeWithPivot.
 
-	decomposition := qrDecomposition decomposeWithPivot.
-	decomposition first * decomposition second.
+	self assert: qrDecomposition q * qrDecomposition r closeTo: expectedQR.
+	self assert: qrDecomposition pivot equals: expectedPivot.
 
-	expected := PMMatrix rows: { 
-			            { -51. 4. 12 }.
-			            { 167. -68. 6 }.
-			            { 24. -41. -4 } }.
-	self
-		assert: decomposition first * decomposition second
-		closeTo: expected
+	reconstruction := qrDecomposition q * qrDecomposition r * qrDecomposition permutationMatrixFromPivot inverse.
+	self assert: reconstruction closeTo: a
 ]
 
 { #category : 'tests' }

src/Math-Matrix/PMQRDecomposition.class.st

@@ -9,7 +9,7 @@ Class {
 		'colSize',
 		'r',
 		'q',
-		'comparisonPrecision'
+		'pivot'
 	],
 	#category : 'Math-Matrix',
 	#package : 'Math-Matrix'
@@ -22,18 +22,6 @@ PMQRDecomposition class >> of: matrix [
 	^ self new of: matrix
 ]
 
-{ #category : 'constants' }
-PMQRDecomposition >> comparisonPrecision [
-
-	^ comparisonPrecision ifNil: [ self defaultComparisonPrecision ]
-]
-
-{ #category : 'accessing' }
-PMQRDecomposition >> comparisonPrecision: anObject [
-
-	comparisonPrecision := anObject
-]
-
 { #category : 'arithmetic' }
 PMQRDecomposition >> decompose [
 	"
@@ -73,91 +61,70 @@ https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
 
 { #category : 'arithmetic' }
 PMQRDecomposition >> decomposeWithPivot [
+	"Variant of the decompose method that introduces a pivot. At the beginning of each step it takes the largest remaining column, thus introducing a pivot. For more info, look at https://en.wikipedia.org/wiki/QR_decomposition#Column_pivoting.
+Here the pivot is stored as an array containing the new order of the columns of the input matrix. It can be used to generate the proper permutation matrix with the permutationMatrixFromPivot method"
 
-	| i vectorOfNormSquareds rank positionOfMaximum pivot matrixOfMinor |
-	vectorOfNormSquareds := matrixToDecompose columnsCollect: [ 
-		                        :columnVector | columnVector * columnVector ].
+	| i vectorOfNormSquareds rank positionOfMaximum matrixOfMinor |
+	vectorOfNormSquareds := matrixToDecompose columnsCollect: [ :columnVector | columnVector * columnVector ].
 	positionOfMaximum := vectorOfNormSquareds indexOf: vectorOfNormSquareds max.
-	pivot := Array new: vectorOfNormSquareds size.
+	pivot := (1 to: vectorOfNormSquareds size) asArray.
 	rank := 0.
-	[ 
-	| householderReflection householderMatrix householderVector columnVectorFromRMatrix |
-	rank := rank + 1.
-	pivot at: rank put: positionOfMaximum.
-	r swapColumn: rank withColumn: positionOfMaximum.
-	vectorOfNormSquareds swap: rank with: positionOfMaximum.
-	columnVectorFromRMatrix := r columnVectorAt: rank size: colSize.
-	householderReflection := self
-		                         householderReflectionOf:
-		                         columnVectorFromRMatrix
-		                         atColumnNumber: rank.
-	householderVector := householderReflection at: 1.
-	householderMatrix := householderReflection at: 2.
-	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 | 
-			| rowNumber columnNumber |
-			rowNumber := rank + index - 1.
-			columnNumber := rank + column - 1.
-			r
-				rowAt: rowNumber
-				columnAt: columnNumber
-				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: [ 
-			positionOfMaximum := (vectorOfNormSquareds
-				       copyFrom: rank + 1
-				       to: vectorOfNormSquareds size) max.
-			(positionOfMaximum closeTo: 0 precision: self comparisonPrecision) ifTrue: [ positionOfMaximum := 0 ].
-			positionOfMaximum := positionOfMaximum > 0
-				      ifTrue: [ 
-				      vectorOfNormSquareds indexOf: positionOfMaximum startingAt: rank + 1 ]
-				      ifFalse: [ 0 ] ]
-		ifFalse: [ positionOfMaximum := 0 ].
-	positionOfMaximum > 0 ] whileTrue.
+	[
+		| temp householderReflection householderMatrix householderVector columnVectorFromRMatrix |
+		rank := rank + 1.
+		temp := pivot at: rank.
+		pivot at: rank put: (pivot at: positionOfMaximum).
+		pivot at: positionOfMaximum put: temp.
+
+		r swapColumn: rank withColumn: positionOfMaximum.
+		vectorOfNormSquareds swap: rank with: positionOfMaximum.
+		columnVectorFromRMatrix := r columnVectorAt: rank size: colSize.
+		householderReflection := self householderReflectionOf: columnVectorFromRMatrix atColumnNumber: rank.
+		householderVector := householderReflection first.
+		householderMatrix := householderReflection second.
+		q := q * householderMatrix.
+		matrixOfMinor := r minor: rank - 1 and: rank - 1.
+		matrixOfMinor := matrixOfMinor
+		                 - (householderVector second tensorProduct: householderVector first * householderVector second * matrixOfMinor).
+		matrixOfMinor rowsWithIndexDo: [ :aRow :index |
+				aRow withIndexDo: [ :element :column |
+						| rowNumber columnNumber |
+						rowNumber := rank + index - 1.
+						columnNumber := rank + column - 1.
+						r rowAt: rowNumber columnAt: columnNumber 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: [
+					positionOfMaximum := (vectorOfNormSquareds copyFrom: rank + 1 to: vectorOfNormSquareds size) max.
+					(positionOfMaximum closeTo: 0) ifTrue: [ positionOfMaximum := 0 ].
+					positionOfMaximum := positionOfMaximum > 0
+						                     ifTrue: [ vectorOfNormSquareds indexOf: positionOfMaximum startingAt: rank + 1 ]
+						                     ifFalse: [ 0 ] ]
+			ifFalse: [ positionOfMaximum := 0 ].
+		positionOfMaximum > 0 ] whileTrue.
 	i := 0.
-	[ (r rowAt: colSize) isZero ] whileTrue: [ 
-		i := i + 1.
-		colSize := colSize - 1 ].
-	i > 0 ifTrue: [ 
-		r := self upperTriangularPartOf: r With: colSize.
-		i := q numberOfColumns - i.
-		pivot := pivot copyFrom: 1 to: i.
-		q := PMMatrix rows:
-			     (q rowsCollect: [ :row | row copyFrom: 1 to: i ]) ].
+	[ (r rowAt: colSize) isZero ] whileTrue: [
+			i := i + 1.
+			colSize := colSize - 1 ].
+	i > 0 ifTrue: [
+			r := self upperTriangularPartOf: r With: colSize.
+			i := q numberOfColumns - i.
+			q := PMMatrix rows: (q rowsCollect: [ :row | row copyFrom: 1 to: i ]) ].
 	^ Array with: q with: r with: pivot
 ]
 
-{ #category : 'constants' }
-PMQRDecomposition >> defaultComparisonPrecision [
-
-	^ 0.0001
-]
-
 { #category : 'private' }
 PMQRDecomposition >> householderReflectionOf: columnVector atColumnNumber: columnNumber [
 
 	| householderVector v identityMatrix householderMatrix |
 	householderVector := columnVector householder.
 	v := (PMVector zeros: columnNumber - 1) , (householderVector at: 2).
 	identityMatrix := PMSymmetricMatrix identity: colSize.
-	householderMatrix := identityMatrix
-	                     -
-	                     ((householderVector at: 1) * v tensorProduct: v).
-	^ Array with: householderVector with: householderMatrix .
+	householderMatrix := identityMatrix - (householderVector first * v tensorProduct: v).
+	^ Array with: householderVector with: householderMatrix
 ]
 
 { #category : 'private' }
@@ -183,8 +150,36 @@ PMQRDecomposition >> of: matrix [
 
 	matrixToDecompose := matrix.
 	colSize := matrixToDecompose numberOfRows.
-		r := self initialRMatrix.
-	q := self initialQMatrix.
+	r := self initialRMatrix.
+	q := self initialQMatrix
+]
+
+{ #category : 'accessing' }
+PMQRDecomposition >> permutationMatrixFromPivot [
+
+	| matrix |
+	matrix := PMMatrix zerosRows: matrixToDecompose numberOfRows cols: matrixToDecompose numberOfColumns.
+	pivot withIndexCollect: [ :column :index | matrix at: column at: index put: 1 ].
+
+	^ matrix
+]
+
+{ #category : 'accessing' }
+PMQRDecomposition >> pivot [
+
+	^ pivot
+]
+
+{ #category : 'accessing' }
+PMQRDecomposition >> q [
+
+	^ q
+]
+
+{ #category : 'accessing' }
+PMQRDecomposition >> r [
+
+	^ r
 ]
 
 { #category : 'private' }