Skip to content

Refactor QR Decomposition #289

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 14 commits into from
Oct 24, 2022
Merged

Refactor QR Decomposition #289

Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
5e28302
Used a message that performs the same logic in its method.
hemalvarambhia Sep 7, 2022
5c09925
Performed the same refactoring in the pivot message's method.
hemalvarambhia Sep 7, 2022
3a04d72
Extracted duplicate code to an intention revealing message and method.
hemalvarambhia Sep 7, 2022
97bc32b
Clarified the name of a temporary variable.
hemalvarambhia Sep 7, 2022
df9874c
Clarified the name of a block variable.
hemalvarambhia Sep 7, 2022
0e2952f
Made similar code even more similar.
hemalvarambhia Sep 7, 2022
bc7c95d
Moved the declaration closer to the site of the computation.
hemalvarambhia Sep 11, 2022
1bd0578
Extracted a key computation step to an intention-revealing method.
hemalvarambhia Sep 11, 2022
b270613
Made similar code even more similar.
hemalvarambhia Sep 11, 2022
5fb1e1b
Reduced the scope of the identity matrix to the code block.
hemalvarambhia Sep 11, 2022
9baa7ba
We now use the new method and have reduced the scope of a temporary v…
hemalvarambhia Sep 11, 2022
c34d9cc
Removed an unused local variable.
hemalvarambhia Sep 11, 2022
009e677
Reduced the scope of a temporary variable and made similar code (vari…
hemalvarambhia Sep 11, 2022
04cc08a
Improved the name of a local variable.
hemalvarambhia Sep 25, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 63 additions & 38 deletions src/Math-Matrix/PMQRDecomposition.class.st
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I want Housholder Reflection to become a first class object.

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
Expand All @@ -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 ]) ] ].
Expand All @@ -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 [

Expand All @@ -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)
]