``PMPrincipalComponentAnalyser`` - sqrt undefined for number less than zero

[nikhilpinnaparaju]

Code to Reproduce Error

|a pca |
a := PMMatrix rows: #(#(-1 -1 1) #(-2 -1 2) #(-3 -2 3)).
pca := PMPrincipalComponentAnalyserSVD new componentsNumber: 2.
pca fit: a.
pca transformMatrix.

sqrt attempted on negative eigenvalues.

[khinsen]

A shorter code snippet with the same effect:

|a |
a := PMMatrix rows: #(#(-1 -1 1) #(-2 -1 2) #(-3 -2 3)).
a decompose.

The SVD algorithm implemented in PMSingularValueDecomposition is far from optimal, and in particular is subject to loss of precision. This is what happens here: the input matrix has a zero singular value which ends up being computed as the square root of -1.4709258215942544e-29. Fixing this requires implementing a more robust SVD algorithm (there are several candidates for that).

[AtharvaKhare]

@khinsen is Golub-Reinsch SVD algorithm good candidate? The GNU Scientific Library uses that.

[khinsen]

@AtharvaKhare Yes, that's one of the robust ones. And GSL is more generally a good source of inspiration, it's a very robust library.

One ingredient to Golub-Reinsch is LU decomposition, which would be nice to have in PolyMath for other reasons as well (there are other higher-level algorithms that build on it).

Question to @SergeStinckwich: what's the policy on API stability in PolyMath? I find it somewhat regrettable that SVD uses decompose, because there are many commonly used matrix decompositions. I'd prefer methods named decomposeLU, decomposeSVD etc., which tell the reader more precisely what is happening.

[SergeStinckwich]

@khinsen at the moment, we have no rules for API and I guess that not that much other projects depends on PolyMath. We can change the API before releasing version 1.0.

I agree with decomposeLU, decomposeSVD would be better name. Can you open a new issue about that ?

[khinsen]

Done: #116.