Merge pull request #558 from karlnapf/master

convienience method for eigenvalues and check for correct results in example
shogun-toolbox · May 29, 2012 · 4b6f6f9 · 4b6f6f9
2 parents fb4d081 + 4c20d8a
commit 4b6f6f9
Show file tree

Hide file tree

Showing 3 changed files with 91 additions and 7 deletions.
diff --git a/examples/undocumented/libshogun/mathematics_lapack.cpp b/examples/undocumented/libshogun/mathematics_lapack.cpp
@@ -5,19 +5,67 @@
  * (at your option) any later version.
  *
  * Written (W) 2011 Sergey Lisitsyn
+ * Written (W) 2012 Heiko Strathmann
  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
  */
 
 #include <shogun/base/init.h>
 #include <shogun/lib/config.h>
+#include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/lapack.h>
 
 using namespace shogun;
 
-int main(int argc, char** argv)
+bool is_equal(float64_t a, float64_t b, float64_t eps)
 {
-	init_shogun();
+	return CMath::abs(a-b)<=eps;
+}
+
+void test_ev()
+{
+#ifdef HAVE_LAPACK
+
+	SGMatrix<float64_t> A(3,3);
+	A(0,0)=0;
+	A(0,1)=1;
+	A(0,2)=0;
+	A(1,0)=1;
+	A(1,1)=0;
+	A(1,2)=1;
+	A(1,0)=0;
+	A(2,1)=1;
+	A(2,2)=0;
+
+	SGVector<float64_t> ev=CMath::compute_eigenvectors(A);
+	CMath::display_matrix(A.matrix, A.num_rows, A.num_cols, "A");
+	CMath::display_vector(ev.vector, ev.vlen, "eigenvalues");
+
+	float64_t sqrt22=CMath::sqrt(2.0)/2.0;
+	float64_t eps=10E-16;
+
+	/* check for correct eigenvectors */
+	ASSERT(is_equal(A(0,0), 0.5, eps));
+	ASSERT(is_equal(A(0,1), -sqrt22, eps));
+	ASSERT(is_equal(A(0,2), 0.5, eps));
+
+	ASSERT(is_equal(A(1,0), -sqrt22, eps));
+	ASSERT(is_equal(A(1,1), 0, eps));
+	ASSERT(is_equal(A(1,2), sqrt22, eps));
+
+	ASSERT(is_equal(A(2,0), 0.5, eps));
+	ASSERT(is_equal(A(2,1), sqrt22, eps));
+	ASSERT(is_equal(A(2,2), 0.5, eps));
 
+	/* check for correct eigenvalues */
+	ASSERT(is_equal(ev[0], -sqrt22*2, eps));
+	ASSERT(is_equal(ev[1], 0, eps));
+	ASSERT(is_equal(ev[2], sqrt22*2, eps));
+
+#endif /* HAVE_LAPACK */
+}
+
+void test_lapack()
+{
 #ifdef HAVE_LAPACK
 	// size of square matrix
 	int N = 100;
@@ -41,13 +89,13 @@ int main(int argc, char** argv)
 			double_matrix[i*N+j] = ((double)(i-j))/(i+j+1);
 
 		double_matrix[i*N+i] += 100;
-	}	
+	}
 	status = 0;
 	wrap_dsygvx(1,'V','U',N,double_matrix,N,double_matrix,N,1,3,double_eigenvalues,double_eigenvectors,&status);
 	if (status!=0)
 	{
 		printf("DSYGVX/SSYGVX failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_eigenvectors;
 
@@ -59,7 +107,7 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DGEQRF/DORGQR failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_tau;
 
@@ -74,7 +122,7 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DGESVD failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_s;
 	delete[] double_U;
@@ -86,13 +134,22 @@ int main(int argc, char** argv)
 	if (status!=0)
 	{
 		printf("DSYEV failed with code %d\n",status);
-		return -1;
+		ASSERT(false);
 	}
 	delete[] double_eigenvalues;
 	delete[] double_matrix;
 
 	#endif // HAVE_LAPACK
+}
+
+int main(int argc, char** argv)
+{
+	init_shogun();
+
+	test_lapack();
+	test_ev();
 
 	exit_shogun();
 	return 0;
 }
+
diff --git a/src/shogun/mathematics/Math.cpp b/src/shogun/mathematics/Math.cpp
@@ -546,6 +546,22 @@ void CMath::inverse(SGMatrix<float64_t> matrix)
 	SG_FREE(ipiv);
 }
 
+SGVector<float64_t> CMath::compute_eigenvectors(SGMatrix<float64_t> matrix)
+{
+	if (matrix.num_rows!=matrix.num_rows)
+	{
+		SG_SERROR("CMath::compute_eigenvectors(SGMatrix<float64_t>): matrix"
+				" rows and columns are not equal!\n");
+	}
+
+	/* use reference counting for SGVector */
+	SGVector<float64_t> result(NULL, 0, true);
+	result.vlen=matrix.num_rows;
+	result.vector=compute_eigenvectors(matrix.matrix, matrix.num_rows,
+			matrix.num_rows);
+	return result;
+}
+
 double* CMath::compute_eigenvectors(double* matrix, int n, int m)
 {
 	ASSERT(n == m);

diff --git a/src/shogun/mathematics/Math.h b/src/shogun/mathematics/Math.h
@@ -494,6 +494,17 @@ class CMath : public CSGObject
 		}
 
 #ifdef HAVE_LAPACK
+		/** compute eigenvalues and eigenvectors of symmetric matrix using
+		 * LAPACK
+		 *
+		 * @param matrix symmetric matrix to compute eigenproblem. Is
+		 * overwritten and contains orthonormal eigenvectors afterwards
+		 * @return eigenvalues vector with eigenvalues equal to number of rows
+		 * in matrix
+		 * */
+		static SGVector<float64_t> compute_eigenvectors(
+				SGMatrix<float64_t> matrix);
+
 		/** compute eigenvalues and eigenvectors of symmetric matrix
 		 *
 		 * @param matrix  overwritten and contains n orthonormal eigenvectors