From 36bc15e1e2e79871f1b46d3548ed80d470c93444 Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Tue, 7 Apr 2026 17:12:36 +0200 Subject: [PATCH 1/4] revert last pr, fixed pivot --- src/Math-Matrix/PMQRDecomposition.class.st | 127 ++++++++------------- 1 file changed, 47 insertions(+), 80 deletions(-) diff --git a/src/Math-Matrix/PMQRDecomposition.class.st b/src/Math-Matrix/PMQRDecomposition.class.st index 5fc426c..525f1fc 100644 --- a/src/Math-Matrix/PMQRDecomposition.class.st +++ b/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 [ " @@ -74,79 +62,58 @@ https://en.wikipedia.org/wiki/QR_decomposition#Using_Householder_reflections { #category : 'arithmetic' } PMQRDecomposition >> decomposeWithPivot [ - | 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 [ From cf2b92d6a3bed29093c383eeb4233e3515b43807 Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Tue, 7 Apr 2026 17:13:31 +0200 Subject: [PATCH 2/4] cleaned up methods and added accessors --- src/Math-Matrix-Tests/PMQRTest.class.st | 56 +++++----------------- src/Math-Matrix/PMQRDecomposition.class.st | 38 ++++++++++++--- 2 files changed, 44 insertions(+), 50 deletions(-) diff --git a/src/Math-Matrix-Tests/PMQRTest.class.st b/src/Math-Matrix-Tests/PMQRTest.class.st index 4f1e753..dc32ae0 100644 --- a/src/Math-Matrix-Tests/PMQRTest.class.st +++ b/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 [ @@ -225,24 +192,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' } diff --git a/src/Math-Matrix/PMQRDecomposition.class.st b/src/Math-Matrix/PMQRDecomposition.class.st index 525f1fc..5761f50 100644 --- a/src/Math-Matrix/PMQRDecomposition.class.st +++ b/src/Math-Matrix/PMQRDecomposition.class.st @@ -121,10 +121,8 @@ PMQRDecomposition >> householderReflectionOf: columnVector atColumnNumber: colum 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' } @@ -150,8 +148,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' } From 6a58d2debe178719f263cdb4eeab5ca4777d26ab Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Tue, 7 Apr 2026 17:40:19 +0200 Subject: [PATCH 3/4] added more tests --- src/Math-Matrix-Tests/PMQRTest.class.st | 45 +++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/src/Math-Matrix-Tests/PMQRTest.class.st b/src/Math-Matrix-Tests/PMQRTest.class.st index dc32ae0..d227703 100644 --- a/src/Math-Matrix-Tests/PMQRTest.class.st +++ b/src/Math-Matrix-Tests/PMQRTest.class.st @@ -118,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 [ From cdc30e2928b1b0ffebd0b30209e54776a7048c50 Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Wed, 8 Apr 2026 11:41:41 +0200 Subject: [PATCH 4/4] Added doc --- src/Math-Matrix/PMQRDecomposition.class.st | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Math-Matrix/PMQRDecomposition.class.st b/src/Math-Matrix/PMQRDecomposition.class.st index 5761f50..eb14780 100644 --- a/src/Math-Matrix/PMQRDecomposition.class.st +++ b/src/Math-Matrix/PMQRDecomposition.class.st @@ -61,6 +61,8 @@ 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 matrixOfMinor | vectorOfNormSquareds := matrixToDecompose columnsCollect: [ :columnVector | columnVector * columnVector ].