Vectorz's, Clatrix's svd return incompatible results

[mars0i]

Vectorz's svd returns a full singular-value decomposition, while the Clatrix core.matrix svd returns some kind of reduced singular-value decomposition. That means that the shapes of :U and :V* matrices can differ between the two versions of svd. The raw Clatrix interface allows you to get the full SVD with optional keyword arguments (:type :full). Maybe the optional arguments should be added to the core.matrix Clatrix wrapper? Or should Vectorz return a reduced representation? (I can't figure out where the core.matrix wrapper for Clatrix lives.)

user=> (run! pm (vals (svd (matrix :vectorz [[1 2 3][4 5 6]]))))
[[0.386  0.922]
 [0.922 -0.386]]
[9.508 0.773]
[[ 0.429  0.566 0.704]
 [-0.806 -0.112 0.581]
 [ 0.408 -0.816 0.408]]
nil
user=> (run! pm (vals (svd (matrix :clatrix [[1 2 3][4 5 6]]))))
[[-0.386 -0.922]
 [-0.922  0.386]]
[9.508 0.773]
[[-0.429 -0.566 -0.704]
 [ 0.806  0.112 -0.581]]
nil
;; Notice that the V* matrices differ above.
user=> (require '[clatrix.core :as cx])
nil
user=> (cx/svd (matrix :clatrix [[1 2 3][4 5 6]]))
{:left  A 2x2 matrix
 -------------
-3.86e-01 -9.22e-01
-9.22e-01  3.86e-01
, :right  A 3x2 matrix
 -------------
-4.29e-01  8.06e-01
-5.66e-01  1.12e-01
-7.04e-01 -5.81e-01
, :values (9.508032000695724 0.7728696356734847), :rank 2}
user=> (cx/svd (matrix :clatrix [[1 2 3][4 5 6]]) :type :full)
{:left  A 2x2 matrix
 -------------
-3.86e-01 -9.22e-01
-9.22e-01  3.86e-01
, :right  A 3x3 matrix
 -------------
-4.29e-01  8.06e-01  4.08e-01
-5.66e-01  1.12e-01 -8.16e-01
-7.04e-01 -5.81e-01  4.08e-01
, :values (9.508032000695724 0.7728696356734847), :rank 2}
user=> (run! pm (vals (cx/svd (matrix :clatrix [[1 2 3][4 5 6]]) :type :full)))
[[-0.386 -0.922]
 [-0.922  0.386]]
[[-0.429  0.806  0.408]
 [-0.566  0.112 -0.816]
 [-0.704 -0.581  0.408]]
[9.508 0.773]

I gather that it's OK for the U and V* matrices to be negated--i.e. the vectorz versions are -1 times the clatrix versions.

Given that there are various reduced representations that people use for SVD, perhaps the docstring for svd in core.matrix should explicitly specify the dimensions of the :U and :V* matrices, i.e. that the first is m x m and that the second is n x n, and not only that they are unitary. (This is helpful if you know you need to do something with svd but don't really know much about singular-value decompositions. Then you go to learn about them, and find out that there are different return values that you could get... and have to figure out that unitary implies square, and reduced representations are not necessarily square, so therefore svd must not be returning a reduced representation. Explicit dimensions short-cuts that reasoning a bit. :-) Related to #321.

[mikera]

Clatrix implementation is here: https://github.com/tel/clatrix/blob/master/src/clatrix/core.clj . This would be the best place to file issues with the Clatrix implementation.

I agree that results should be consistent, and the core.matrix docstring should describe the dimensions as you suggest. Probably Clatrix should return the full SVD by default.

[mars0i]

OK, good--will do.

I hadn't noticed mp/PSVDDecomposition (svd [m options] ... in that file. I'd only seen the (defn svd, which has a different interface.