Merge branch '8-bind-svd' into 'master'

Binded SVD in the python level and documented

Closes #8

See merge request !10

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;

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;
+
+}
+
+
+

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);

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 */
 };

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)
+  
+  
+

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