diff --git a/bob/math/cpp/svd.cpp b/bob/math/cpp/svd.cpp
index 1fbf99cacf78f210b2ab5bdee97476d2aa695a45..854cd47328f134cabab3c4f28040c879243ab409 100644
--- a/bob/math/cpp/svd.cpp
+++ b/bob/math/cpp/svd.cpp
@@ -79,6 +79,19 @@ static void svd_lapack( const char jobz, const int M, const int N,
if (info != 0)
throw std::runtime_error("The LAPACK dgesdd function returned a non-zero value. You may consider using LAPACK dgsevd instead (see #171) by enabling the 'safe' option.");
}
+
+ // Defining the sign of the eigenvectors
+ // Approch extracted from page 8 - http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf
+ if(U[0] < 0){
+ int ucol=0; ucol= (jobz=='A')? M : std::min(M,N);
+ for (int i=0; i<ldu*ucol; i++){
+ U[i] = -1*U[i];
+ }
+
+ for (int i=0; i<ldvt*N; i++){
+ VT[i] = -1*VT[i];
+ }
+ }
}
void bob::math::svd(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
@@ -148,7 +161,8 @@ void bob::math::svd_(const blitz::Array<double,2>& A, blitz::Array<double,2>& U,
// Call the LAPACK function
svd_lapack(jobz, N, M, A_lapack, lda, S_lapack, U_lapack, ldu,
- VT_lapack, ldvt, safe);
+ VT_lapack, ldvt, safe);
+
// Copy singular vectors back to U, V and sigma if required
if (!U_direct_use) Vt = U_blitz_lapack;
diff --git a/bob/math/gsvd.cpp b/bob/math/gsvd.cpp
index aa20a8156b06b67c1e39b8fa278d5230640ab81d..8d80eafce2bb89ab4ebec2f45aef45312011a990 100644
--- a/bob/math/gsvd.cpp
+++ b/bob/math/gsvd.cpp
@@ -9,6 +9,8 @@
#include <bob.blitz/cppapi.h>
#include <bob.blitz/cleanup.h>
#include <bob.math/gsvd.h>
+#include <bob.math/svd.h>
+#include <bob.math/linear.h>
PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
@@ -41,9 +43,9 @@ PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
auto B_bz = PyBlitzArrayCxx_AsBlitz<double,2>(B);
- const int M = A_bz->extent(0);
- const int N = A_bz->extent(1);
- const int P = B_bz->extent(0);
+ const int M = A_bz->extent(0);
+ const int N = A_bz->extent(1);
+ const int P = B_bz->extent(0);
// Creating the output matrices
blitz::Array<double,2> U(M, M); U=0;
@@ -75,3 +77,73 @@ PyObject* py_gsvd (PyObject*, PyObject* args, PyObject* kwds) {
return 0;
}
+
+
+PyObject* py_svd (PyObject*, PyObject* args, PyObject* kwds) {
+
+ /* Parses input arguments in a single shot */
+ static const char* const_kwlist[] = { "A", 0 /* Sentinel */ };
+ static char** kwlist = const_cast<char**>(const_kwlist);
+
+ PyBlitzArrayObject* A = 0;
+
+ if (!PyArg_ParseTupleAndKeywords(args, kwds, "O&", kwlist,
+ &PyBlitzArray_Converter, &A
+ ))
+ return 0;
+
+ auto A_ = make_safe(A);
+
+ if (A->ndim != 2 || A->type_num != NPY_FLOAT64) {
+ PyErr_Format(PyExc_TypeError, "`A` matrix only supports 2D 64-bit float array");
+ return 0;
+ }
+
+
+ auto A_bz = PyBlitzArrayCxx_AsBlitz<double,2>(A);
+
+ int M = A_bz->extent(0);
+ int N = A_bz->extent(1);
+
+ // Creating the output matrices
+ // Prepares to call LAPACK function:
+ // We recall that FORTRAN/LAPACK is column-major order whereas blitz arrays
+ // are row-major order by default.
+ // If A = U.S.V^T, then A^T = V.S.U^T
+
+ blitz::Array<double,2> V(M, M); V=0;
+ blitz::Array<double,1> S(std::min(M,N)); S=0;
+ blitz::Array<double,2> U(N, N); U=0;
+
+
+ try {
+ bob::math::svd(*A_bz,V,S,U);
+
+ // S for the python output.
+ // LAPACK returns an 1d matrix of size n.
+ // In order to nicelly fit A = U S V,
+ // I need to create a matrix S_output of size MxN and fits the diagonal of S (NxN) on it.
+ blitz::Array<double,2> S_output(M, N); S_output=0;
+ blitz::Array<double,2> S_diag(N, N); S_diag=0;
+ bob::math::diag(S, S_diag);
+ S_output(blitz::Range(0,N-1), blitz::Range(0,N-1)) = S_diag;
+ //bob::math::diag(S, (reinterpret_cast<blitz::Array<double,2>>)S_output(blitz::Range(0,N-1), blitz::Range(0,N-1)));
+ return Py_BuildValue("NNN",
+ PyBlitzArrayCxx_AsConstNumpy(U),
+ PyBlitzArrayCxx_AsConstNumpy(S_output),
+ PyBlitzArrayCxx_AsConstNumpy(V));
+ }
+ catch (std::exception& e) {
+ PyErr_SetString(PyExc_RuntimeError, e.what());
+ }
+ catch (...) {
+ PyErr_SetString(PyExc_RuntimeError, "norminv failed: unknown exception caught");
+ }
+
+
+ return 0;
+
+}
+
+
+
diff --git a/bob/math/gsvd.h b/bob/math/gsvd.h
index cc38f9fc38360e041536a83cca82c2e6ad3f33e8..1465400ec974636041915128499d0eb9c22de182 100644
--- a/bob/math/gsvd.h
+++ b/bob/math/gsvd.h
@@ -8,3 +8,4 @@
#include <Python.h>
PyObject* py_gsvd(PyObject*, PyObject* args, PyObject* kwds);
+PyObject* py_svd(PyObject*, PyObject* args, PyObject* kwds);
diff --git a/bob/math/main.cpp b/bob/math/main.cpp
index e78fe28de0135d7fdf1e9a78f94f803175404202..97886907b22fade0c9afd98dd54cd842071db4a5 100644
--- a/bob/math/main.cpp
+++ b/bob/math/main.cpp
@@ -376,6 +376,22 @@ static bob::extension::FunctionDoc s_gsvd = bob::extension::FunctionDoc(
;
+static bob::extension::FunctionDoc s_svd = bob::extension::FunctionDoc(
+ "svd",
+ "Computes the SVD",
+ "Computes the SVD (Singular Value Decomposition).\n"
+ "[U,S,V] = svd(A) returns :math:`U`, :math:`S` and :math:`V` such that ` \n\n"
+ ".. math:: A=U S V\n"
+ )
+ .add_prototype("A", "")
+ .add_parameter("A", "[array_like (float, 2D)]", "Must be :math:`m \\times n`")
+ .add_return("U", "[array_like (float, 2D)]", "The :math:`U` matrix of left singular vectors (size :math:`m \\times m`)")
+ .add_return("S", "[array_like (float, 2D)]", "The matrix of singular values :math:`S` of size :math:`m \\times n`")
+ .add_return("V", "[array_like (float, 2D)]", "The :math:`V^{T}` matrix of right singular vectors (size :math:`n \\times n`)")
+;
+
+
+
static PyMethodDef module_methods[] = {
{
s_histogram_intersection.name(),
@@ -497,6 +513,12 @@ static PyMethodDef module_methods[] = {
METH_VARARGS|METH_KEYWORDS,
s_gsvd.doc()
},
+ {
+ s_svd.name(),
+ (PyCFunction)py_svd,
+ METH_VARARGS|METH_KEYWORDS,
+ s_svd.doc()
+ },
{0} /* Sentinel */
};
diff --git a/bob/math/test_gsvd.py b/bob/math/test_gsvd.py
index 804203e0740ea8c58ed0a37057ac17d54bebff15..eac76a011776c66bbc7c3f743bb73550aca251eb 100644
--- a/bob/math/test_gsvd.py
+++ b/bob/math/test_gsvd.py
@@ -20,6 +20,7 @@ import bob.math
import numpy
import nose.tools
numpy.random.seed(10)
+import pkg_resources
def gsvd_relations(A,B):
@@ -39,12 +40,17 @@ def gsvd_relations(A,B):
nose.tools.eq_( (abs(B-B_check) < 1e-10).all(), True )
+def svd_relations(A):
+
+ [U, S, V] = bob.math.svd(A)
+ A_check = numpy.dot(numpy.dot(V,S), U)
+ nose.tools.eq_( (abs(A-A_check) < 1e-10).all(), True )
+
def test_first_case():
- """
- Testing the first scenario of gsvd:
- M-K-L >= 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
- """
+
+ #Testing the first scenario of gsvd:
+ #M-K-L >= 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
A = numpy.random.rand(10,10)
B = numpy.random.rand(790,10)
@@ -53,10 +59,10 @@ def test_first_case():
def test_second_case():
- """
- Testing the second scenario of gsvd:
- M-K-L < 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
- """
+
+ #Testing the second scenario of gsvd:
+ #M-K-L < 0 (check http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga4a187519e5c71da3b3f67c85e9baf0f2.html#ga4a187519e5c71da3b3f67c85e9baf0f2)
+
A = numpy.random.rand(4,5)
B = numpy.random.rand(11,5)
@@ -65,10 +71,9 @@ def test_second_case():
def test_corner_case():
- """
- Testing when P <= N.
-
- """
+
+ #Testing when P <= N.
+
A = numpy.random.rand(25, 25)
B = numpy.random.rand(25, 25)
@@ -76,3 +81,56 @@ def test_corner_case():
+def test_svd_relation():
+
+ ##Testing SVD
+
+ A = numpy.random.rand(25, 25)
+ svd_relations(A)
+
+ #A = numpy.random.rand(20, 25)
+ #svd_relations(A)
+
+ A = numpy.random.rand(30, 25)
+ svd_relations(A)
+
+
+
+def test_svd_signal():
+
+ ##Testing SVD signal
+ ##This test was imported from bob.learn.linear
+
+ A = numpy.array([[3,-3,100], [4,-4,50], [3.5,-3.5,-50], [3.8,-3.7,-100]], dtype='float64')
+
+ U_ref = numpy.array([[ 2.20825004e-03, -1.80819459e-03, -9.99995927e-01],
+ [ 7.09549949e-01, -7.04649416e-01, 2.84101853e-03],
+ [ -7.04651683e-01, -7.09553332e-01, -2.73037723e-04]])
+
+ [U,S,V] = bob.math.svd(A)
+ nose.tools.eq_((abs(U-U_ref) < 1e-8).all(), True)
+ svd_relations(A)
+
+
+
+def test_svd_signal_book_example():
+
+ ## Reference copied from here http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf
+
+ A = numpy.array([[2.5, 63.5, 40.1, 78, 61.1],
+ [0.9, 58.0, 25.1, 78, 94.1],
+ [1.7, 46.0, 65.0, 78, 106.4],
+ [1.2, 15.7, 102.1, 78, 173.0],
+ [1.5, 12.2, 100.0, 77, 199.7],
+ [2.0, 8.9, 87.8, 76, 176.0],
+ [3.8, 2.7, 17.1, 69, 373.6],
+ [1.0, 1.7, 140.0, 73, 283.7],
+ [2.1, 1.0, 55.0, 79, 34.7],
+ [0.8, 0.2, 50.4, 73, 36.4]])
+
+ [U,S,V] = bob.math.svd(A)
+ assert U[0,0] > 0
+ svd_relations(A)
+
+
+
diff --git a/doc/py_api.rst b/doc/py_api.rst
index 252233d3e0a009b637dd39cf2de4cc17448f1c99..7e63d97495da91b843e56fa886879b8909c0c054 100644
--- a/doc/py_api.rst
+++ b/doc/py_api.rst
@@ -16,6 +16,7 @@ Summary
bob.math.LPInteriorPointShortstep
bob.math.LPInteriorPointLongstep
bob.math.chi_square
+ bob.math.svd
bob.math.gsvd
bob.math.histogram_intersection
bob.math.kullback_leibler