Add BDC-SVD (CPU) and euclidean norm to linalg

src/shogun/mathematics/linalg/LinalgBackendBase.h

@@ -41,6 +41,7 @@
 #include <shogun/lib/config.h>
 #include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/linalg/GPUMemoryBase.h>
+#include <shogun/mathematics/linalg/LinalgEnums.h>
 #include <shogun/mathematics/linalg/internal/Block.h>
 
 namespace shogun
@@ -638,7 +639,7 @@ namespace shogun
 #define BACKEND_GENERIC_SVD(Type, Container)                                   \
 	virtual void svd(                                                          \
 	    const Container<Type>& A, SGVector<Type> s, SGMatrix<Type> U,          \
-	    bool thin_U) const                                                     \
+	    bool thin_U, linalg::SVDAlgorithm alg) const                           \
 	{                                                                          \
 		SG_SNOTIMPLEMENTED;                                                    \
 	}

src/shogun/mathematics/linalg/LinalgBackendEigen.h

@@ -37,6 +37,7 @@
 #include <shogun/mathematics/Math.h>
 #include <shogun/mathematics/eigen3.h>
 #include <shogun/mathematics/linalg/LinalgBackendBase.h>
+#include <shogun/mathematics/linalg/LinalgEnums.h>
 #include <shogun/mathematics/linalg/LinalgMacros.h>
 
 namespace shogun
@@ -316,7 +317,7 @@ namespace shogun
 #define BACKEND_GENERIC_SVD(Type, Container)                                   \
 	virtual void svd(                                                          \
 	    const Container<Type>& A, SGVector<Type> s, Container<Type> U,         \
-	    bool thin_U) const;
+	    bool thin_U, linalg::SVDAlgorithm alg) const;
 		DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_SVD, SGMatrix)
 #undef BACKEND_GENERIC_SVD
 
@@ -592,8 +593,8 @@ namespace shogun
 		/** Eigen3 compute svd method */
 		template <typename T>
 		void svd_impl(
-		    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U,
-		    bool thin_U) const;
+		    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U, bool thin_U,
+		    linalg::SVDAlgorithm alg) const;
 
 		/** Eigen3 compute trace method */
 		template <typename T>

src/shogun/mathematics/linalg/LinalgEnums.h

@@ -0,0 +1,59 @@
+/*
+ * Copyright (c) 2017, Shogun-Toolbox e.V. <shogun-team@shogun-toolbox.org>
+ * All rights reserved.
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions are met:
+ *
+ *  1. Redistributions of source code must retain the above copyright notice,
+ *     this list of conditions and the following disclaimer.
+ *
+ *  2. Redistributions in binary form must reproduce the above copyright
+ *     notice, this list of conditions and the following disclaimer in the
+ *     documentation and/or other materials provided with the distribution.
+ *
+ *  3. Neither the name of the copyright holder nor the names of its
+ *     contributors may be used to endorse or promote products derived from
+ *     this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+
+#ifndef LINALG_ENUMS_H__
+#define LINALG_ENUMS_H__
+
+namespace shogun
+{
+
+	namespace linalg
+	{
+
+		/**
+		 * Enum for choosing the algorithm used to calculate SVD.
+		 * The <em>bidiagonal divide and conquer</em> algorithm
+		 * is faster on large matrices but it's available with
+		 * Eigen >= 3.3, furthermore it may produce inaccurate
+		 * results when compiled with unsafe math optimization.
+		 * For more details see:
+		 * https://eigen.tuxfamily.org/dox/classEigen_1_1BDCSVD.html
+		 * https://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html
+		 */
+		enum class SVDAlgorithm
+		{
+			BidiagonalDivideConquer,
+			Jacobi
+		};
+	}
+}
+
+#endif // LINALG_ENUMS_H__

src/shogun/mathematics/linalg/LinalgNamespace.h

@@ -34,6 +34,7 @@
 #define LINALG_NAMESPACE_H_
 
 #include <shogun/mathematics/linalg/LinalgBackendBase.h>
+#include <shogun/mathematics/linalg/LinalgEnums.h>
 #include <shogun/mathematics/linalg/SGLinalg.h>
 
 namespace shogun
@@ -1145,6 +1146,19 @@ namespace shogun
 			return infer_backend(a)->mean(a);
 		}
 
+		/**
+		 * Method that computes the euclidean norm of a vector.
+		 *
+		 * @param a SGVector
+		 * @return The vector norm
+		 */
+		template <typename T>
+		T norm(const SGVector<T>& a)
+		{
+			REQUIRE(a.size() > 0, "Vector cannot be empty!\n");
+			return CMath::sqrt(dot(a, a));
+		}
+
 		/**
 		 * Solve the linear equations \f$Ax=b\f$ through the
 		 * QR decomposition of A.
@@ -1393,11 +1407,15 @@ namespace shogun
 		 * @param U The matrix that stores the resulting unitary matrix U
 		 * @param thin_U Whether to compute the full or thin matrix U
 		 * (default:thin)
+		 * @param alg Whether to compute the svd through bidiagonal divide
+		 * and conquer algorithm or Jacobi's algorithm (@see SVDAlgorithm)
+		 * (default: bidiagonal divide and conquer)
 		 */
 		template <typename T>
 		void
 		svd(const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U,
-		    bool thin_U = true)
+		    bool thin_U = true,
+		    SVDAlgorithm alg = SVDAlgorithm::BidiagonalDivideConquer)
 		{
 			auto r = CMath::min(A.num_cols, A.num_rows);
 			REQUIRE(
@@ -1425,7 +1443,7 @@ namespace shogun
 			                   "smaller dimension (%d).\n",
 			    s.vlen, r);
 
-			infer_backend(A)->svd(A, s, U, thin_U);
+			infer_backend(A)->svd(A, s, U, thin_U, alg);
 		}
 
 		/**

src/shogun/mathematics/linalg/backend/eigen/Decompositions.cpp

@@ -31,6 +31,7 @@
  */
 
 #include <shogun/mathematics/linalg/LinalgBackendEigen.h>
+#include <shogun/mathematics/linalg/LinalgEnums.h>
 #include <shogun/mathematics/linalg/LinalgMacros.h>
 
 using namespace shogun;
@@ -67,9 +68,9 @@ DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_EIGEN_SOLVER_SYMMETRIC, SGMatrix)
 #define BACKEND_GENERIC_SVD(Type, Container)                                   \
 	void LinalgBackendEigen::svd(                                              \
 	    const Container<Type>& A, SGVector<Type> s, Container<Type> U,         \
-	    bool thin_U) const                                                     \
+	    bool thin_U, linalg::SVDAlgorithm alg) const                           \
 	{                                                                          \
-		return svd_impl(A, s, U, thin_U);                                      \
+		return svd_impl(A, s, U, thin_U, alg);                                 \
 	}
 DEFINE_FOR_NON_INTEGER_PTYPE(BACKEND_GENERIC_SVD, SGMatrix)
 #undef BACKEND_GENERIC_SVD
@@ -185,14 +186,37 @@ void LinalgBackendEigen::eigen_solver_symmetric_impl(
 
 template <typename T>
 void LinalgBackendEigen::svd_impl(
-    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U, bool thin_U) const
+    const SGMatrix<T>& A, SGVector<T>& s, SGMatrix<T>& U, bool thin_U,
+    linalg::SVDAlgorithm alg) const
 {
 	typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
 	typename SGVector<T>::EigenVectorXtMap s_eig = s;
 	typename SGMatrix<T>::EigenMatrixXtMap U_eig = U;
 
-	auto svd_eig =
-	    A_eig.jacobiSvd(thin_U ? Eigen::ComputeThinU : Eigen::ComputeFullU);
-	s_eig = svd_eig.singularValues().template cast<T>();
-	U_eig = svd_eig.matrixU().template cast<T>();
+	switch (alg)
+	{
+	case linalg::SVDAlgorithm::BidiagonalDivideConquer:
+	{
+#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
+		auto svd_eig =
+		    A_eig.bdcSvd(thin_U ? Eigen::ComputeThinU : Eigen::ComputeFullU);
+		s_eig = svd_eig.singularValues().template cast<T>();
+		U_eig = svd_eig.matrixU().template cast<T>();
+		break;
+#else
+		SG_SWARNING(
+		    "At least Eigen 3.3 is required for BDC-SVD.\n"
+		    "Falling back on Jacobi-SVD.\n")
+#endif
+	}
+
+	case linalg::SVDAlgorithm::Jacobi:
+	{
+		auto svd_eig =
+		    A_eig.jacobiSvd(thin_U ? Eigen::ComputeThinU : Eigen::ComputeFullU);
+		s_eig = svd_eig.singularValues().template cast<T>();
+		U_eig = svd_eig.matrixU().template cast<T>();
+		break;
+	}
+	}
 }

tests/unit/mathematics/linalg/operations/Eigen3_operations_unittest.cc

@@ -881,6 +881,23 @@ TEST(LinalgBackendEigen, SGMatrix_multiply_by_rectified_linear_derivative)
 		EXPECT_NEAR(i * (A[i] != 0), B[i], 1e-15);
 }
 
+TEST(LinalgBackendEigen, SGVector_norm)
+{
+	const index_t n = 5;
+	SGVector<float64_t> v(n);
+	float64_t gt = 0;
+	for (index_t i = 0; i < n; ++i)
+	{
+		v[i] = i;
+		gt += i * i;
+	}
+	gt = CMath::sqrt(gt);
+
+	auto result = norm(v);
+
+	EXPECT_NEAR(result, gt, 1E-15);
+}
+
 TEST(LinalgBackendEigen, SGVector_qr_solver)
 {
 	const index_t n = 3;
@@ -1337,7 +1354,68 @@ TEST(LinalgBackendEigen, SGMatrix_block_rowwise_sum)
 	}
 }
 
-TEST(LinalgBackendEigen, SGMatrix_svd_thinU)
+TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_thinU)
+{
+	const index_t m = 5, n = 3;
+	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
+	                    0.30786772, 0.25503552, 0.34367041, 0.66491478,
+	                    0.20488809, 0.5734351,  0.87179189, 0.07139643,
+	                    0.28540373, 0.06264684, 0.56204061};
+	float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883};
+	float64_t result_U[] = {-0.60700926, -0.16647013, -0.56501385, -0.26696629,
+	                        -0.46186125, -0.69145782, 0.29548428,  0.5718984,
+	                        0.31771648,  -0.08101592, -0.27461424, 0.37170223,
+	                        -0.12681555, -0.53830325, 0.69323293};
+
+	SGMatrix<float64_t> A(data, m, n, false);
+	SGMatrix<float64_t> U(m, n);
+	SGVector<float64_t> s(n);
+
+	svd(A, s, U, true, SVDAlgorithm::Jacobi);
+
+	for (index_t i = 0; i < n; ++i)
+	{
+		auto c = CMath::sign(U[i * m] * result_U[i * m]);
+		for (index_t j = 0; j < m; ++j)
+			EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7);
+	}
+	for (index_t i = 0; i < (index_t)s.size(); ++i)
+		EXPECT_NEAR(s[i], result_s[i], 1e-7);
+}
+
+TEST(LinalgBackendEigen, SGMatrix_svd_jacobi_fullU)
+{
+	const index_t m = 5, n = 3;
+	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
+	                    0.30786772, 0.25503552, 0.34367041, 0.66491478,
+	                    0.20488809, 0.5734351,  0.87179189, 0.07139643,
+	                    0.28540373, 0.06264684, 0.56204061};
+	float64_t result_s[] = {1.75382524, 0.56351367, 0.41124883};
+	float64_t result_U[] = {
+	    -0.60700926, -0.16647013, -0.56501385, -0.26696629, -0.46186125,
+	    -0.69145782, 0.29548428,  0.5718984,   0.31771648,  -0.08101592,
+	    -0.27461424, 0.37170223,  -0.12681555, -0.53830325, 0.69323293,
+	    -0.27809756, -0.68975171, -0.11662812, 0.38274703,  0.53554354,
+	    0.025973184, 0.520631112, -0.56921636, 0.62571522,  0.11287970};
+
+	SGMatrix<float64_t> A(data, m, n, false);
+	SGMatrix<float64_t> U(m, m);
+	SGVector<float64_t> s(n);
+
+	svd(A, s, U, false, SVDAlgorithm::Jacobi);
+
+	for (index_t i = 0; i < n; ++i)
+	{
+		auto c = CMath::sign(U[i * m] * result_U[i * m]);
+		for (index_t j = 0; j < m; ++j)
+			EXPECT_NEAR(U[i * m + j], c * result_U[i * m + j], 1e-7);
+	}
+	for (index_t i = 0; i < (index_t)s.size(); ++i)
+		EXPECT_NEAR(s[i], result_s[i], 1e-7);
+}
+
+#if EIGEN_VERSION_AT_LEAST(3, 3, 0)
+TEST(LinalgBackendEigen, SGMatrix_svd_bdc_thinU)
 {
 	const index_t m = 5, n = 3;
 	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
@@ -1354,7 +1432,7 @@ TEST(LinalgBackendEigen, SGMatrix_svd_thinU)
 	SGMatrix<float64_t> U(m, n);
 	SGVector<float64_t> s(n);
 
-	svd(A, s, U, true);
+	svd(A, s, U, true, SVDAlgorithm::BidiagonalDivideConquer);
 
 	for (index_t i = 0; i < n; ++i)
 	{
@@ -1366,7 +1444,7 @@ TEST(LinalgBackendEigen, SGMatrix_svd_thinU)
 		EXPECT_NEAR(s[i], result_s[i], 1e-7);
 }
 
-TEST(LinalgBackendEigen, SGMatrix_svd_fullU)
+TEST(LinalgBackendEigen, SGMatrix_svd_bdc_fullU)
 {
 	const index_t m = 5, n = 3;
 	float64_t data[] = {0.68764958, 0.11456779, 0.75164207, 0.50436194,
@@ -1385,7 +1463,7 @@ TEST(LinalgBackendEigen, SGMatrix_svd_fullU)
 	SGMatrix<float64_t> U(m, m);
 	SGVector<float64_t> s(n);
 
-	svd(A, s, U, false);
+	svd(A, s, U, false, SVDAlgorithm::BidiagonalDivideConquer);
 
 	for (index_t i = 0; i < n; ++i)
 	{
@@ -1396,6 +1474,7 @@ TEST(LinalgBackendEigen, SGMatrix_svd_fullU)
 	for (index_t i = 0; i < (index_t)s.size(); ++i)
 		EXPECT_NEAR(s[i], result_s[i], 1e-7);
 }
+#endif
 
 TEST(LinalgBackendEigen, SGMatrix_trace)
 {