Merge pull request #289 from PolyMathOrg/refactor-qr-decomposition

Refactor QR Decomposition

src/Math-Matrix/PMQRDecomposition.class.st

@@ -29,64 +29,67 @@ that describes the mechanics:
 https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections
 "
 
-	| identityMatrix householderVector i matrixOfMinor |
-	identityMatrix := PMSymmetricMatrix identity: colSize.
+	| i matrixOfMinor |
 	1 to: self numberOfColumns do: [ :col | 
-		| columnVectorFromRMatrix householderMatrix v |
+		| householderReflection householderMatrix householderVector columnVectorFromRMatrix |
 		columnVectorFromRMatrix := r columnVectorAt: col size: colSize.
-		householderVector := columnVectorFromRMatrix householder.
-		v := (PMVector zeros: col - 1) , (householderVector at: 2).
-		householderMatrix := identityMatrix
-		                     -
-		                     ((householderVector at: 1) * v tensorProduct: v).
+		householderReflection := self
+			                         householderReflectionOf:
+			                         columnVectorFromRMatrix
+			                         atColumnNumber: col.
+		householderVector := householderReflection at: 1.
+		householderMatrix := householderReflection at: 2.
 		q := q * householderMatrix.
 		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 | 
+			aRow withIndexDo: [ :element :column | 
+				| rowNumber columnNumber |
+				rowNumber := col + index - 1.
+				columnNumber := col + column - 1.
 				r
-					rowAt: col + index - 1
-					columnAt: col + c - 1
-					put: ((n closeTo: 0)
+					rowAt: rowNumber
+					columnAt: columnNumber
+					put: ((element closeTo: 0)
 							 ifTrue: [ 0 ]
-							 ifFalse: [ n ]) ] ] ].
+							 ifFalse: [ element ]) ] ] ].
 	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).
+		r := self upperTriangularPartOf: r With: colSize.
 		i := q numberOfColumns - i.
 		q := PMMatrix rows:
-			     (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
+			     (q rowsCollect: [ :row | row copyFrom: 1 to: i ]) ].
 	^ Array with: q with: r
 ]
 
 { #category : #arithmetic }
 PMQRDecomposition >> decomposeWithPivot [
 
-	| identityMatrix householderVector i v vectorOfNormSquareds rank mx pivot matrixOfMinor |
+	| i vectorOfNormSquareds rank positionOfMaximum pivot matrixOfMinor |
 	vectorOfNormSquareds := matrixToDecompose columnsCollect: [ 
 		                        :columnVector | columnVector * columnVector ].
-	mx := vectorOfNormSquareds indexOf: vectorOfNormSquareds max.
+	positionOfMaximum := vectorOfNormSquareds indexOf: vectorOfNormSquareds max.
 	pivot := Array new: vectorOfNormSquareds size.
 	rank := 0.
-	identityMatrix := PMSymmetricMatrix identity: colSize.
 	[ 
-	| householderMatrix |
+	| householderReflection householderMatrix householderVector columnVectorFromRMatrix |
 	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).
+	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
@@ -95,9 +98,12 @@ PMQRDecomposition >> decomposeWithPivot [
 			                  * (householderVector at: 2) * matrixOfMinor).
 	matrixOfMinor rowsWithIndexDo: [ :aRow :index | 
 		aRow withIndexDo: [ :element :column | 
+			| rowNumber columnNumber |
+			rowNumber := rank + index - 1.
+			columnNumber := rank + column - 1.
 			r
-				rowAt: rank + index - 1
-				columnAt: rank + column - 1
+				rowAt: rowNumber
+				columnAt: columnNumber
 				put: ((element closeTo: 0)
 						 ifTrue: [ 0 ]
 						 ifFalse: [ element ]) ] ].
@@ -109,29 +115,42 @@ PMQRDecomposition >> decomposeWithPivot [
 			- (r rowAt: rank columnAt: ind) squared ].
 	rank < vectorOfNormSquareds size
 		ifTrue: [ 
-			mx := (vectorOfNormSquareds
+			positionOfMaximum := (vectorOfNormSquareds
 				       copyFrom: rank + 1
 				       to: vectorOfNormSquareds size) max.
-			(mx closeTo: 0) ifTrue: [ mx := 0 ].
-			mx := mx > 0
+			(positionOfMaximum closeTo: 0) ifTrue: [ positionOfMaximum := 0 ].
+			positionOfMaximum := positionOfMaximum > 0
 				      ifTrue: [ 
-				      vectorOfNormSquareds indexOf: mx startingAt: rank + 1 ]
+				      vectorOfNormSquareds indexOf: positionOfMaximum startingAt: rank + 1 ]
 				      ifFalse: [ 0 ] ]
-		ifFalse: [ mx := 0 ].
-	mx > 0 ] whileTrue.
+		ifFalse: [ positionOfMaximum := 0 ].
+	positionOfMaximum > 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).
+		r := self upperTriangularPartOf: r With: colSize.
 		i := q numberOfColumns - i.
 		pivot := pivot copyFrom: 1 to: i.
 		q := PMMatrix rows:
-			     (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
+			     (q rowsCollect: [ :row | row copyFrom: 1 to: i ]) ].
 	^ Array with: q with: r with: pivot
 ]
 
+{ #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 .
+]
+
 { #category : #private }
 PMQRDecomposition >> initialQMatrix [
 
@@ -158,3 +177,9 @@ PMQRDecomposition >> of: matrix [
 		r := self initialRMatrix.
 	q := self initialQMatrix.
 ]
+
+{ #category : #private }
+PMQRDecomposition >> upperTriangularPartOf: matrix With: columnSize [
+
+	^ PMMatrix rows: (r rows copyFrom: 1 to: columnSize)
+]