From 331ee84fe40730238519614c6c859de995226d54 Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Fri, 13 Mar 2026 17:11:39 +0100 Subject: [PATCH 1/2] add different version of the jacobi transformation algorithm --- .../PMJacobiTransformationBis.class.st | 136 ++++++++++++++++++ .../PMJacobiTransformationHelper.class.st | 13 +- 2 files changed, 147 insertions(+), 2 deletions(-) create mode 100644 src/Math-Matrix/PMJacobiTransformationBis.class.st diff --git a/src/Math-Matrix/PMJacobiTransformationBis.class.st b/src/Math-Matrix/PMJacobiTransformationBis.class.st new file mode 100644 index 0000000..a9504b9 --- /dev/null +++ b/src/Math-Matrix/PMJacobiTransformationBis.class.st @@ -0,0 +1,136 @@ +Class { + #name : 'PMJacobiTransformationBis', + #superclass : 'PMJacobiTransformation', + #instVars : [ + 'n', + 'd', + 'v', + 'nrot', + 'nMax' + ], + #category : 'Math-Matrix', + #package : 'Math-Matrix' +} + +{ #category : 'as yet unclassified' } +PMJacobiTransformationBis >> eigsrt [ + + | k p | + 1 to: n - 1 do: [ :i | + k := i. + p := d at: i. + i + 1 to: n do: [ :j | + (d at: j) >= p ifTrue: [ + k := j. + p := d at: j ] ]. + k ~= i ifTrue: [ + d at: k put: (d at: i). + d at: i put: p. + 1 to: n do: [ :j | + p := v at: j at: i. + v at: j at: i put: (v at: j at: k). + v at: j at: k put: p ] ] ] +] + +{ #category : 'operation' } +PMJacobiTransformationBis >> evaluate [ + + ^ d +] + +{ #category : 'initialization' } +PMJacobiTransformationBis >> initialize: aMatrix [ + + | a b z sm tresh c s g h t theta tau | + n := aMatrix numberOfColumns. + a := PMMatrix new: n. + v := PMMatrix zerosRows: n cols: n. + 1 to: n do: [ :i | + 1 to: n do: [ :j | + a at: i at: j put: (aMatrix at: i at: j). + i = j ifTrue: [ v at: i at: j put: 1 ] ] ]. + nrot := 0. + nMax := 500. + + b := PMVector new: n. + d := PMVector new: n. + + 1 to: n do: [ :ip | + b at: ip put: (a at: ip at: ip). + d at: ip put: (b at: ip) ]. + z := PMVector zeros: n. + + 1 to: 50 do: [ :i | + sm := self sumUpperTriangle: a. + sm = 0 ifTrue: [ + self eigsrt. + ^ self ]. + tresh := i < 4 + ifTrue: [ 0.2 * sm / n squared ] + ifFalse: [ 0 ]. + 1 to: n - 1 do: [ :ip | + ip + 1 to: n do: [ :iq | + g := 100 * (a at: ip at: iq) abs. + i > 4 & ((d at: ip) abs + g = (d at: ip) abs) & ((d at: iq) abs + g = (d at: iq) abs) + ifTrue: [ a at: ip at: iq put: 0 ] + ifFalse: [ + (a at: ip at: iq) abs > tresh ifTrue: [ + h := (d at: iq) - (d at: ip). + h abs + g = h abs + ifTrue: [ t := (a at: ip at: iq) / h ] + ifFalse: [ + theta := 0.5 * h / (a at: ip at: iq). + t := 1 / (theta abs + (1.0 + theta squared) sqrt). + theta < 0 ifTrue: [ t := t negated ] ]. + c := 1 / (1 + t squared) sqrt. + s := t * c. + tau := s / (1 + c). + h := t * (a at: ip at: iq). + z at: ip put: (z at: ip) - h. + z at: iq put: (z at: iq) + h. + d at: ip put: (d at: ip) - h. + d at: iq put: (d at: iq) + h. + a at: ip at: iq put: 0. + 1 to: ip - 1 do: [ :j | + g := a at: j at: ip. + h := a at: j at: iq. + a at: j at: ip put: g - (s * (h + (g * tau))). + a at: j at: iq put: h + (s * (g - (h * tau))) ]. + ip + 1 to: iq - 1 do: [ :j | + g := a at: ip at: j. + h := a at: j at: iq. + a at: ip at: j put: g - (s * (h + (g * tau))). + a at: j at: iq put: h + (s * (g - (h * tau))) ]. + iq + 1 to: n do: [ :j | + g := a at: ip at: j. + h := a at: iq at: j. + a at: ip at: j put: g - (s * (h + (g * tau))). + a at: iq at: j put: h + (s * (g - (h * tau))) ]. + 1 to: n do: [ :j | + g := v at: j at: ip. + h := v at: j at: iq. + v at: j at: ip put: g - (s * (h + (g * tau))). + v at: j at: iq put: h + (s * (g - (h * tau))) ]. + nrot := nrot + 1 ] ] ] ]. + 1 to: n do: [ :ip | + b at: ip put: (b at: ip) + (z at: ip). + d at: ip put: (b at: ip). + z at: ip put: 0 ] ]. + self eigsrt. + ^ self +] + +{ #category : 'as yet unclassified' } +PMJacobiTransformationBis >> sumUpperTriangle: a [ + + | sm | + sm := 0. + 1 to: n - 1 do: [ :ip | ip + 1 to: n do: [ :iq | sm := sm + (a at: ip at: iq) abs ] ]. + ^ sm +] + +{ #category : 'private - transforming' } +PMJacobiTransformationBis >> transform [ + + ^ v +] diff --git a/src/Math-Matrix/PMJacobiTransformationHelper.class.st b/src/Math-Matrix/PMJacobiTransformationHelper.class.st index 28daefa..e88054a 100644 --- a/src/Math-Matrix/PMJacobiTransformationHelper.class.st +++ b/src/Math-Matrix/PMJacobiTransformationHelper.class.st @@ -24,12 +24,21 @@ PMJacobiTransformationHelper class >> new [ ^ self error: 'Illegal creation message for this class' ] +{ #category : 'as yet unclassified' } +PMJacobiTransformationHelper >> cleanValues: aCollection [ + + ^ aCollection collect: [ :each | + (each closeTo: 0) + ifTrue: [ 0 ] + ifFalse: [ each ] ] +] + { #category : 'initialization' } PMJacobiTransformationHelper >> initialize: aSymmetricMatrix [ | jacobi | - jacobi := PMJacobiTransformation matrix: aSymmetricMatrix. - eigenvalues := jacobi evaluate. + jacobi := PMJacobiTransformationBis matrix: aSymmetricMatrix. + eigenvalues := self cleanValues: jacobi evaluate. eigenvectors := jacobi transform columnsCollect: [ :each | each]. ] From 50be58c843bf959a9b04a3511fd0578362d9333e Mon Sep 17 00:00:00 2001 From: DurieuxPol Date: Fri, 13 Mar 2026 17:11:58 +0100 Subject: [PATCH 2/2] change how svd is computed --- .../PMSingularValueDecomposition.class.st | 50 +++++++++++++++---- 1 file changed, 41 insertions(+), 9 deletions(-) diff --git a/src/Math-Matrix/PMSingularValueDecomposition.class.st b/src/Math-Matrix/PMSingularValueDecomposition.class.st index 50a740a..c01729a 100644 --- a/src/Math-Matrix/PMSingularValueDecomposition.class.st +++ b/src/Math-Matrix/PMSingularValueDecomposition.class.st @@ -45,6 +45,31 @@ PMSingularValueDecomposition class >> decompose: aMatrix [ ^ self new initialize: aMatrix ] +{ #category : 'adding' } +PMSingularValueDecomposition >> addMissingLeftVectorsFor: aMatrix [ + + | m k | + m := leftSingularMatrix numberOfRows. + k := aMatrix rank. + + 1 to: m do: [ :i | + | w norm | + k = m ifTrue: [ ^ self ]. + + w := PMVector zeros: m. + w at: i put: 1.0. + 1 to: k do: [ :j | + | uj proj | + uj := leftSingularMatrix columnAt: j. + proj := uj * w. + w := w - (uj * proj) ]. + + norm := w norm. + norm > 0 ifTrue: [ + k := k + 1. + leftSingularMatrix atColumn: k put: w / norm ] ]. +] + { #category : 'accessing' } PMSingularValueDecomposition >> diagonalSingularValueMatrix [ "Diagonal matrix S in A = USV decomposition. All elements are 0 except those on the main diagonal. The values on the diagonal are singular values of matrix A. By convention, they are sorted in descending order" @@ -55,20 +80,27 @@ PMSingularValueDecomposition >> diagonalSingularValueMatrix [ { #category : 'initialization' } PMSingularValueDecomposition >> initialize: aMatrix [ - | m n symU symV eigenU eigenV diag | + | m n symV eigenV diag singularValues | m := aMatrix numberOfRows. n := aMatrix numberOfColumns. - symU := aMatrix * aMatrix transpose. + symV := aMatrix transpose * aMatrix. + eigenV := symV eigen. "Sort eigenpairs by descending eigenvalue" "Clamp tiny negative eigenvalues to 0" - "Expensive computation" - eigenU := symU eigen. - eigenV := symV eigen. - leftSingularMatrix := (PMMatrix rows: eigenU vectors) transpose. + singularValues := eigenV values collect: [ :lam | (lam max: 0) sqrt ]. rightSingularMatrix := (PMMatrix rows: eigenV vectors) transpose. - diag := m < n - ifTrue: [ eigenU values ] - ifFalse: [ eigenV values ]. + leftSingularMatrix := PMMatrix rows: m columns: m random: 0. + 1 to: singularValues size do: [ :i | + | sigma vi ui | + sigma := singularValues at: i. + sigma > 0 ifTrue: [ + vi := rightSingularMatrix columnAt: i. + ui := aMatrix * vi / sigma. + leftSingularMatrix atColumn: i put: ui ] ]. + + self addMissingLeftVectorsFor: aMatrix. + + diag := eigenV values. diagonalSingularValueMatrix := PMMatrix rows: m columns: n random: 0. diagonalSingularValueMatrix setDiagonal: diag sqrt ]