diff --git a/src/Math-Matrix-Tests/PMSingularValueDecompositionExample.class.st b/src/Math-Matrix-Tests/PMSingularValueDecompositionExample.class.st new file mode 100644 index 0000000..533fc9a --- /dev/null +++ b/src/Math-Matrix-Tests/PMSingularValueDecompositionExample.class.st @@ -0,0 +1,71 @@ +Class { + #name : 'PMSingularValueDecompositionExample', + #superclass : 'Object', + #instVars : [ + 'matrix', + 'actualU', + 'actualV', + 'actualS' + ], + #category : 'Math-Matrix-Tests', + #package : 'Math-Matrix-Tests' +} + +{ #category : 'instance creation' } +PMSingularValueDecompositionExample class >> newMatrix: aMatrix u: uMatrix v: vMatrix s: sMatrix [ + + ^ self new + matrix: aMatrix; + actualU: uMatrix; + actualV: vMatrix; + actualS: sMatrix; + yourself +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualS [ + + ^ actualS +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualS: anObject [ + + actualS := anObject +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualU [ + + ^ actualU +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualU: anObject [ + + actualU := anObject +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualV [ + + ^ actualV +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> actualV: anObject [ + + actualV := anObject +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> matrix [ + + ^ matrix +] + +{ #category : 'accessing' } +PMSingularValueDecompositionExample >> matrix: anObject [ + + matrix := anObject +] diff --git a/src/Math-Matrix-Tests/PMSingularValueDecompositionTest.class.st b/src/Math-Matrix-Tests/PMSingularValueDecompositionTest.class.st index a1bafc1..bbfc783 100644 --- a/src/Math-Matrix-Tests/PMSingularValueDecompositionTest.class.st +++ b/src/Math-Matrix-Tests/PMSingularValueDecompositionTest.class.st @@ -28,131 +28,91 @@ Internal Representation and Key Implementation Points. " Class { #name : 'PMSingularValueDecompositionTest', - #superclass : 'TestCase', + #superclass : 'ParametrizedTestCase', #instVars : [ - 'matrix', - 'actualU', - 'actualV', - 'actualS' + 'example' + ], + #classInstVars : [ + 'example1', + 'example2', + 'example3', + 'example4' ], #category : 'Math-Matrix-Tests', #package : 'Math-Matrix-Tests' } -{ #category : 'as yet unclassified' } -PMSingularValueDecompositionTest >> loadExample1 [ - "Simple example. Square matrix with real values" - - matrix := PMMatrix rows: #( - (0 4) - (0 0)). - - actualU := PMMatrix rows: #( - (1 0) - (0 1)). - - actualV := PMMatrix rows: #( - (0 1) - (1 0)). - - actualS := PMMatrix rows: #( - (4 0) - (0 0)). +{ #category : 'initi' } +PMSingularValueDecompositionTest class >> loadExamples [ + + example1 := PMSingularValueDecompositionExample + newMatrix: (PMMatrix rows: #( #( 0 4 ) #( 0 0 ) )) + u: (PMMatrix rows: #( #( 1 0 ) #( 0 1 ) )) + v: (PMMatrix rows: #( #( 0 1 ) #( 1 0 ) )) + s: (PMMatrix rows: #( #( 4 0 ) #( 0 0 ) )). + example2 := PMSingularValueDecompositionExample + newMatrix: (PMMatrix rows: #( #( 0 1 0 0 ) #( 0 0 2 0 ) #( 0 0 0 3 ) #( 0 0 0 0 ) )) + u: (PMMatrix rows: #( #( 0 0 1 0 ) #( 0 1 0 0 ) #( 1 0 0 0 ) #( 0 0 0 1 ) )) + v: (PMMatrix rows: #( #( 0 0 0 1 ) #( 0 0 1 0 ) #( 0 1 0 0 ) #( 1 0 0 0 ) )) + s: (PMMatrix rows: #( #( 3 0 0 0 ) #( 0 2 0 0 ) #( 0 0 1 0 ) #( 0 0 0 0 ) )). + example3 := PMSingularValueDecompositionExample + newMatrix: (PMMatrix rows: #( #( -1 1 0 ) #( 0 -1 1 ) )) + u: (PMMatrix rows: { { (1 / 2 sqrt). (1 / 2 sqrt) }. { (-1 / 2 sqrt). (1 / 2 sqrt) } }) + v: (PMMatrix rows: { { (-1 / 6 sqrt). (-1 / 2 sqrt). (1 / 3 sqrt) }. { (2 / 6 sqrt). 0. (1 / 3 sqrt) }. { (-1 / 6 sqrt). (1 / 2 sqrt). (1 / 3 sqrt) } }) + s: (PMMatrix rows: { { 3 sqrt. 0. 0 }. { 0. 1. 0 } }). + example4 := PMSingularValueDecompositionExample + newMatrix: (PMMatrix rows: #( #( 1 0 0 0 2 ) #( 0 0 3 0 0 ) #( 0 0 0 0 0 ) #( 0 2 0 0 0 ) )) + u: (PMMatrix rows: #( #( 0 1 0 0 ) #( 1 0 0 0 ) #( 0 0 0 1 ) #( 0 0 1 0 ) )) + v: (PMMatrix rows: { { 0. 0.2 sqrt. 0. 0.8 sqrt. 0 }. { 0. 0. 1. 0. 0 }. { 1. 0. 0. 0. 0 }. { 0. 0. 0. 0. 1 }. { 0. 0.8 sqrt. 0. 0.2 sqrt negated. 0 } }) + s: (PMMatrix rows: { { 3. 0. 0. 0. 0 }. { 0. 5 sqrt. 0. 0. 0 }. { 0. 0. 2. 0. 0 }. { 0. 0. 0. 0. 0 } }). ] -{ #category : 'as yet unclassified' } -PMSingularValueDecompositionTest >> loadExample2 [ - "Simple example. Square matrix with real values" - - matrix := PMMatrix rows: #( - (0 1 0 0) - (0 0 2 0) - (0 0 0 3) - (0 0 0 0)). - - actualU := PMMatrix rows: #( - (0 0 1 0) - (0 1 0 0) - (1 0 0 0) - (0 0 0 1)). - - actualV := PMMatrix rows: #( - (0 0 0 1) - (0 0 1 0) - (0 1 0 0) - (1 0 0 0)). - - actualS := PMMatrix rows: #( - (3 0 0 0) - (0 2 0 0) - (0 0 1 0) - (0 0 0 0)). +{ #category : 'building suites' } +PMSingularValueDecompositionTest class >> testParameters [ + + self loadExamples. + + ^ ParametrizedTestMatrix new forSelector: #example addOptions: { + example1. + example2. + example3. + example4 } ] -{ #category : 'as yet unclassified' } -PMSingularValueDecompositionTest >> loadExample3 [ - "Rectangular matrix with real values" - - matrix := PMMatrix rows: #( - (-1 1 0) - (0 -1 1)). - - actualU := PMMatrix rows: ( - Array - with: (Array with: -1 / (2 sqrt) with: 1 / (2 sqrt)) - with: (Array with: 1 / (2 sqrt) with: 1 / (2 sqrt)) ). - - actualV := PMMatrix rows: ( - Array - with: (Array with: 1 / (6 sqrt) with: -1 / (2 sqrt) with: 1 / (3 sqrt)) - with: (Array with: -2 / (6 sqrt) with: 0 with: 1 / (3 sqrt)) - with: (Array with: 1 / (6 sqrt) with: 1 / (2 sqrt) with: 1 / (3 sqrt)) ). - - actualS := PMMatrix rows: ( - Array - with: (Array with: 3 sqrt with: 0 with: 0) - with: (Array with: 0 with: 1 with: 0) ). +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> actualS [ + + ^ self example actualS ] -{ #category : 'as yet unclassified' } -PMSingularValueDecompositionTest >> loadExample4 [ - "Rectangular matrix with real values" - - matrix := PMMatrix rows: #( - (1 0 0 0 2) - (0 0 3 0 0) - (0 0 0 0 0) - (0 2 0 0 0)). - - actualU := PMMatrix rows: #( - (0 0 1 0) - (0 1 0 0) - (0 0 0 -1) - (1 0 0 0)). - - actualV := PMMatrix rows: (Array - with: ( Array with: 0 with: 0 with: 0.2 sqrt with: 0 with: -1 * (0.8 sqrt)) - with: ( Array with: 1 with: 0 with: 0 with: 0 with: 0) - with: ( Array with: 0 with: 1 with: 0 with: 0 with: 0) - with: ( Array with: 0 with: 0 with: 0 with: 1 with: 0) - with: ( Array with: 0 with: 0 with: 0.8 sqrt with: 0 with: 0.2 sqrt) - ). - - actualS := PMMatrix rows: (Array - with: (Array with: 2 with: 0 with: 0 with: 0 with: 0) - with: (Array with: 0 with: 3 with: 0 with: 0 with: 0) - with: (Array with: 0 with: 0 with: 5 sqrt with: 0 with: 0) - with: (Array with: 0 with: 0 with: 0 with: 0 with: 0) - ) +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> actualU [ + + ^ self example actualU +] + +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> actualV [ + + ^ self example actualV ] -{ #category : 'initialization' } -PMSingularValueDecompositionTest >> setUp [ - super setUp. - self loadExample1 - "self loadExample2." - "self loadExample3." - "self loadExample4." +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> example [ + + ^ example +] + +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> example: anObject [ + + example := anObject +] + +{ #category : 'accessing' } +PMSingularValueDecompositionTest >> matrix [ + + ^ self example matrix ] { #category : 'tests' } @@ -160,8 +120,8 @@ PMSingularValueDecompositionTest >> testIsExampleCorrect [ "Verifying that the example SVD is correct" | reconstructed | - reconstructed := actualU * actualS * actualV transpose. - self assert: reconstructed closeTo: matrix. + reconstructed := self actualU * self actualS * self actualV transpose. + self assert: reconstructed closeTo: self matrix. ] { #category : 'tests' } @@ -170,10 +130,10 @@ PMSingularValueDecompositionTest >> testIsExampleOrthonormalU [ | m identity | - m := actualU numberOfRows. + m := self actualU numberOfRows. identity := PMSymmetricMatrix identity: m. - self assert: (actualU transpose * actualU) closeTo: identity + self assert: (self actualU transpose * self actualU) closeTo: identity ] { #category : 'tests' } @@ -182,10 +142,10 @@ PMSingularValueDecompositionTest >> testIsExampleOrthonormalV [ | n identity | - n := actualV numberOfRows. + n := self actualV numberOfRows. identity := PMSymmetricMatrix identity: n. - self assert: (actualV transpose * actualV) closeTo: identity + self assert: (self actualV transpose * self actualV) closeTo: identity ] { #category : 'tests' } @@ -194,7 +154,7 @@ PMSingularValueDecompositionTest >> testIsOrthonormalU [ "Verifying that matrix U is orthonormal. That is, all columns of U are unit vectors and orthogonal to each other." | u m identity | - u := matrix decomposeSV leftSingularMatrix. + u := self matrix decomposeSV leftSingularMatrix. m := u numberOfRows. identity := PMSymmetricMatrix identity: m. self assert: u transpose * u closeTo: identity @@ -206,7 +166,7 @@ PMSingularValueDecompositionTest >> testIsOrthonormalV [ "Verifying that matrix U is orthonormal. That is, all columns of U are unit vectors and orthogonal to each other." | v n identity | - v := matrix decomposeSV rightSingularMatrix. + v := self matrix decomposeSV rightSingularMatrix. n := v numberOfRows. identity := PMSymmetricMatrix identity: n. self assert: v transpose * v closeTo: identity @@ -218,7 +178,7 @@ PMSingularValueDecompositionTest >> testIsSquareU [ "U should be a square matrix" | u | - u := matrix decomposeSV leftSingularMatrix. + u := self matrix decomposeSV leftSingularMatrix. self assert: u numberOfRows equals: u numberOfColumns ] @@ -228,7 +188,7 @@ PMSingularValueDecompositionTest >> testIsSquareV [ "V should be a square matrix" | v | - v := matrix decomposeSV rightSingularMatrix. + v := self matrix decomposeSV rightSingularMatrix. self assert: v numberOfRows equals: v numberOfColumns ] @@ -236,12 +196,12 @@ PMSingularValueDecompositionTest >> testIsSquareV [ PMSingularValueDecompositionTest >> testReconstruction [ | svd u v s reconstructed | - svd := matrix decomposeSV. + svd := self matrix decomposeSV. u := svd leftSingularMatrix. v := svd rightSingularMatrix. s := svd diagonalSingularValueMatrix. reconstructed := u * s * v transpose. - self assert: reconstructed closeTo: matrix + self assert: reconstructed closeTo: self matrix ] { #category : 'tests' } @@ -250,9 +210,9 @@ PMSingularValueDecompositionTest >> testShapeS [ "If A is an (m x n) matrix, then S should be an (m x n) matrix" | s m n | - s := matrix decomposeSV diagonalSingularValueMatrix. - m := matrix numberOfRows. - n := matrix numberOfColumns. + s := self matrix decomposeSV diagonalSingularValueMatrix. + m := self matrix numberOfRows. + n := self matrix numberOfColumns. self assert: s numberOfRows equals: m. self assert: s numberOfColumns equals: n ] @@ -263,8 +223,8 @@ PMSingularValueDecompositionTest >> testShapeU [ "If A is an (m x n) matrix, then U should be an (m x m) matrix" | u m | - u := matrix decomposeSV leftSingularMatrix. - m := matrix numberOfRows. + u := self matrix decomposeSV leftSingularMatrix. + m := self matrix numberOfRows. self assert: u numberOfRows equals: m. self assert: u numberOfColumns equals: m ] @@ -275,8 +235,8 @@ PMSingularValueDecompositionTest >> testShapeV [ "If A is an (m x n) matrix, then V should be an (n x n) matrix" | v n | - v := matrix decomposeSV rightSingularMatrix. - n := matrix numberOfColumns. + v := self matrix decomposeSV rightSingularMatrix. + n := self matrix numberOfColumns. self assert: v numberOfRows equals: n. self assert: v numberOfColumns equals: n ] @@ -287,8 +247,8 @@ PMSingularValueDecompositionTest >> testValueS [ "Comparing S to its known value" | s | - s := matrix decomposeSV diagonalSingularValueMatrix. - self assert: s closeTo: actualS + s := self matrix decomposeSV diagonalSingularValueMatrix. + self assert: s closeTo: self actualS ] { #category : 'tests' } @@ -297,8 +257,8 @@ PMSingularValueDecompositionTest >> testValueU [ "Comparing U to its known value" | u | - u := matrix decomposeSV leftSingularMatrix. - self assert: u closeTo: actualU + u := self matrix decomposeSV leftSingularMatrix. + self assert: u closeTo: self actualU ] { #category : 'tests' } @@ -307,6 +267,6 @@ PMSingularValueDecompositionTest >> testValueV [ "Comparing V to its known value" | v | - v := matrix decomposeSV rightSingularMatrix. - self assert: v closeTo: actualV + v := self matrix decomposeSV rightSingularMatrix. + self assert: v closeTo: self actualV ] 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]. ] 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 ]