Merge pull request #287 from PolyMathOrg/qr-decomposition-extract-class

Refactor: Extract QR Decomposition Class

Author: SergeStinckwich

src/Math-Matrix/PMMatrix.class.st

@@ -409,6 +409,12 @@ PMMatrix >> columnAt: anInteger [
 	^ rows collect: [ :each | each at: anInteger ]
 ]
 
+{ #category : #'cell accessing' }
+PMMatrix >> columnVectorAt: col size: dimension [
+
+	^ (self columnAt: col) copyFrom: col to: dimension 
+]
+
 { #category : #iterators }
 PMMatrix >> columnsCollect: aBlock [
 	"Perform the collect: operation on the rows of the receiver."
@@ -751,41 +757,7 @@ PMMatrix >> productWithVector: aVector [
 { #category : #'as yet unclassified' }
 PMMatrix >> qrFactorization [
 
-	| identMat q r hh colSize i |
-	self numberOfRows < self numberOfColumns ifTrue: [ 
-		self error: 'numberOfRows<numberOfColumns' ].
-	r := PMMatrix rows: rows deepCopy.
-	colSize := self numberOfRows.
-	q := PMSymmetricMatrix identity: colSize.
-	identMat := q deepCopy.
-	1 to: self numberOfColumns do: [ :col | 
-		hh := ((r columnAt: col) copyFrom: col to: colSize) householder.
-		i := (PMVector new: col - 1 withAll: 0) , (hh at: 2).
-		q := q * (identMat - ((hh at: 1) * i tensorProduct: i)). "not really necessary, should be simplified"
-		i := PMMatrix rows:
-			     ((r rows allButFirst: col - 1) collect: [ :aRow | 
-				      aRow allButFirst: col - 1 ]).
-		i := i - ((hh at: 2) tensorProduct: (hh at: 1) * (hh at: 2) * i).
-		i rows withIndexDo: [ :aRow :index | 
-			aRow withIndexDo: [ :n :c | 
-				r
-					rowAt: col + index - 1
-					columnAt: col + c - 1
-					put: ((n closeTo: 0)
-							 ifTrue: [ 0 ]
-							 ifFalse: [ n ]) ] ]
-		"col <colSize ifTrue: [i :=(hh at: 2) copyFrom: 2 to: colSize -col +1.  i withIndexDo: [:n :index| r rowAt: col columnAt: index put: n  ]  ]""and this part is not correct, dont uncomment before the bug is corrected! useful if q is not explicitely necessary" ].
-	"r rows allButFirst withIndexDo: [:aRow :ri|1 to: (ri min: self numberOfColumns ) do: [:ci|aRow at: ci put:0  ]   ]  ""not necessary with equalsTo:0"
-	i := 0.
-	[ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ] whileTrue: [ 
-		i := i + 1.
-		colSize := colSize - 1 ].
-	i > 0 ifTrue: [ 
-		r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
-		i := q numberOfColumns - i.
-		q := PMMatrix rows:
-			     (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
-	^ Array with: q with: r
+	^ (PMQRDecomposition of: self) decompose
 ]
 
 { #category : #'as yet unclassified' }

src/Math-Matrix/PMQRDecomposition.class.st

@@ -0,0 +1,91 @@
+"
+I am responsible for the decomposition of a matrix, A say, into a product A = QR of an orthogonal matrix Q and an upper triangular matrix R
+"
+Class {
+	#name : #PMQRDecomposition,
+	#superclass : #Object,
+	#instVars : [
+		'matrixToDecompose',
+		'colSize'
+	],
+	#category : #'Math-Matrix'
+}
+
+{ #category : #'instance creation' }
+PMQRDecomposition class >> of: matrix [
+	matrix numberOfRows < matrix numberOfColumns ifTrue: [ 
+		self error: 'numberOfRows<numberOfColumns' ].
+	^ self new of: matrix
+]
+
+{ #category : #private }
+PMQRDecomposition >> decompose [
+
+	"
+The method appears to be using Householder reflection. There is a wiki page
+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 := PMSymmetricMatrix identity: colSize.
+	1 to: self numberOfColumns do: [ :col | 
+		| columnVectorFromRMatrix householderMatrix |
+		columnVectorFromRMatrix := r columnVectorAt: col size: colSize.
+		householderVector := columnVectorFromRMatrix householder.
+		i := (PMVector zeros: col - 1) , (householderVector at: 2).
+		householderMatrix := identityMatrix
+		                     -
+		                     ((householderVector at: 1) * i tensorProduct: i).
+		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 | 
+			aRow withIndexDo: [ :n :c | 
+				r
+					rowAt: col + index - 1
+					columnAt: col + c - 1
+					put: ((n closeTo: 0)
+							 ifTrue: [ 0 ]
+							 ifFalse: [ n ]) ] ] ].
+	i := 0.
+	[ (r rowAt: colSize) allSatisfy: [ :n | n = 0 ] ] whileTrue: [ 
+		i := i + 1.
+		colSize := colSize - 1 ].
+	i > 0 ifTrue: [ 
+		r := PMMatrix rows: (r rows copyFrom: 1 to: colSize).
+		i := q numberOfColumns - i.
+		q := PMMatrix rows:
+			     (q rows collect: [ :row | row copyFrom: 1 to: i ]) ].
+	^ Array with: q with: r
+]
+
+{ #category : #private }
+PMQRDecomposition >> initialQMatrix [
+
+	^ PMSymmetricMatrix identity: colSize
+]
+
+{ #category : #private }
+PMQRDecomposition >> initialRMatrix [
+
+	^ PMMatrix rows: matrixToDecompose rows deepCopy
+]
+
+{ #category : #private }
+PMQRDecomposition >> numberOfColumns [
+
+	^ matrixToDecompose numberOfColumns
+]
+
+{ #category : #'instance creation' }
+PMQRDecomposition >> of: matrix [
+
+	matrixToDecompose := matrix.
+	colSize := matrixToDecompose numberOfRows.
+]

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

@@ -27,6 +27,20 @@ PMQRTest >> assert: inverse isMoorePenroseInverseOf: aMatrix [
 		closeTo: aMatrix mpInverse transpose.
 ]
 
+{ #category : #tests }
+PMQRTest >> testHorizontalRectangularMatrixCannotBeDecomposed [
+
+	| horizontalRectangularMatrix |
+	horizontalRectangularMatrix := PMMatrix rows: { 
+			                               { 1. 2. 4 }.
+			                               { 5. 6. 7 } }.
+
+	self
+		should: [ 
+			(PMQRDecomposition of: horizontalRectangularMatrix) ]
+		raise: Error
+]
+
 { #category : #tests }
 PMQRTest >> testMoorePenroseInverseOfLargeNonRandomMatrixAndItsTranspose [
 	| a inverse transposeOfA |