Binded SVD in the python level and documented

Implemented the signal standarization. The signs of U and V is defined analysing the first element of U.
If U[0, 0] < 0 then U=-1*U and V=-1*V
With that we have a consistent SVD binds along different implementations (http://prod.sandia.gov/techlib/access-control.cgi/2007/076422.pdf)
parent c72f6c80
Pipeline #6794 passed with stages
in 5 minutes and 33 seconds
......@@ -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;
......
......@@ -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;
}
......@@ -8,3 +8,4 @@
#include <Python.h>
PyObject* py_gsvd(PyObject*, PyObject* args, PyObject* kwds);
PyObject* py_svd(PyObject*, PyObject* args, PyObject* kwds);
......@@ -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 */
};
......
......@@ -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)
......@@ -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
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment