Managed SVD fails / results in incorrect vt when applied to non-symmetric matrices [BUG]

[ZimmerD]

Describe the bug

Given the non-symmetric matrix A :

let A = Matrix.ofJaggedArray [| [|2.;-1.;0.|]; [|4.;3.;-2.|]|]

the result of a SVD is expected to yield:

val vt : Matrix<float> = matrix [[0.7877263614; 0.5012804118; -0.358057437]
                                 [0.5883484054; -0.7844645406; 0.1961161351]
                                 [0.1825741858; 0.3651483717; 0.9128709292]]
val u : Matrix<float> = matrix [[0.1961161351; 0.9805806757]
                                [0.9805806757; -0.1961161351]]
val s : Matrix<float> = matrix [[5.477225575; 0.0; 0.0]
                                [0.0; 2.0; 0.0]]

but instead a call to FSharp.Stats.Algebra.LinearAlgebraManaged.SVD yields:

val vtM : matrix = matrix [[-0.7877263614; -0.5883484054; -0.1825741858]
                           [-0.4803209521; 0.5739984443; 0.2226514461]
                           [0.3507052555; -0.1222886149; -1.119059972]]
val uM : matrix = matrix [[-0.1961161351; -0.9805806757]
                          [-0.9805806757; 0.1961161351]]
val sM : Matrix<float> = matrix [[5.477225575; 0.0; 0.0]
                                 [0.0; 2.0; 0.0]]

Comparing the matrices we can see:

1. The values of sigma (s vs. sM) are equal
2. The results of u (u vs UM) differ in sign (which should be ok, since the direction of the unit vectors is rather arbitrary and can vary between svd implementations)
3. The result of vt (vt, vs vtM) show larger discrepancies.

While the first col of v is equal to the first row of vtM, this is not the case for the other rows/columns.

I suspect that vtM is actually not transposed and other than that something else seems to be wrong when calculating the SVD.

This is backed by the following observation: calculating u s vt should yield A. It does for the reference example but it doesn't for vtM, uM and sM. However, if you transpose vtM, the first column of the matrix product is equal to the first column of A.

Additionally, I tried to compute the svd of let B' = Matrix.ofJaggedArray [| [|2.;-1.;2.;-1.|]; [|4.;3.;4.;3.|]; |] which resulted in :

System.IndexOutOfRangeException: Index was outside the bounds of the array.
   at FSharp.Stats.Algebra.SVD.computeInPlace(Double[,] a) in C:\Users\bvenn\source\repos\FSharp.Stats\src\FSharp.Stats\Algebra\SVD.fs:line 22
   at FSharp.Stats.Algebra.LinearAlgebraManaged.SVD(Matrix`1 a) in C:\Users\bvenn\source\repos\FSharp.Stats\src\FSharp.Stats\Algebra\LinearAlgebraServiceManaged.fs:line 22
   at FSI_0075.svdManaged(Matrix`1 A)
   at <StartupCode$FSI_0091>.$FSI_0091.main@()

The behaviour could not be reproduced using symmetric matrices.

OS and framework information:

• OS: Windows 10
• FSharp.Stats, 0.4.2 used in script via #r "nuget: FSharp.Stats, 0.4.2"

Known workaround Locate FSharp.Stats lib folder using: FSharp.Stats.ServiceLocator.setEnvironmentPathVariable Register LinearAlgebra service calling: LinearAlgebra.Service() Perform SVD calling: FSharp.Stats.Algebra.LinearAlgebra.SVD