added uwedge approximate joint diagonalizer and unit test

src/shogun/mathematics/ajd/UWedge.cpp

@@ -0,0 +1,109 @@
+#ifdef HAVE_EIGEN3
+
+#include <shogun/mathematics/ajd/UWedge.h>
+
+#include <shogun/base/init.h>
+#include <shogun/lib/common.h>
+
+#include <shogun/mathematics/Math.h>
+#include <shogun/mathematics/eigen3.h>
+
+using namespace Eigen;
+
+typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
+typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
+
+using namespace shogun;
+
+SGMatrix<float64_t> CUWedge::diagonalize(SGNDArray<float64_t> &C, SGMatrix<float64_t> V0,
+					double eps, int itermax)
+{
+	int d = C.dims[0];
+	int L = C.dims[2];
+
+	SGMatrix<float64_t> V;
+	if (V0.num_rows != 0)
+	{
+		V = V0.clone();
+	}
+	else
+	{	
+		Eigen::Map<EMatrix> C0(C.get_matrix(0),d,d);
+		EigenSolver<EMatrix> eig;
+		eig.compute(C0);					
+		V = SGMatrix<float64_t>::create_identity_matrix(d,1);
+		Eigen::Map<EMatrix> EV(V.matrix, d,d);
+		EV = eig.pseudoEigenvalueMatrix().cwiseAbs().cwiseSqrt().inverse() * eig.pseudoEigenvectors().transpose();
+	}
+	Eigen::Map<EMatrix> EV(V.matrix, d,d);	
+
+	index_t * Cs_dims = SG_MALLOC(index_t, 3);
+	Cs_dims[0] = d;
+	Cs_dims[1] = d;
+	Cs_dims[2] = L;
+	SGNDArray<float64_t> Cs(Cs_dims,3);
+	memcpy(Cs.array, C.array, Cs.dims[0]*Cs.dims[1]*Cs.dims[2]*sizeof(float64_t));
+
+	EMatrix Rs(d,L);
+	std::vector<double> crit;
+	crit.push_back(0.0);
+	for(int l = 0; l < L; l++)
+	{
+		Eigen::Map<EMatrix> Ci(C.get_matrix(l),d,d);
+		Eigen::Map<EMatrix> Csi(Cs.get_matrix(l),d,d);
+		Ci = 0.5 * (Ci + Ci.transpose());		
+		Csi = EV * Ci * EV.transpose();
+		Rs.col(l) = Csi.diagonal();
+		crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();	
+	}
+
+	double iter = 0;
+	double improve = 10;
+	while (improve > eps && iter < itermax)
+	{
+		EMatrix B = Rs * Rs.transpose();
+
+		EMatrix C1 = EMatrix::Zero(d,d);
+		for(int id = 0; id < d; id++)
+		{
+			// rowSums
+			for(int l = 0; l < L; l++)
+			{
+				Eigen::Map<EMatrix> Csi(Cs.get_matrix(l),d,d);
+				C1.row(id) += Csi.row(id) * Rs(id,l);
+			}		
+		}
+
+		EMatrix D0 = B.cwiseProduct(B.transpose()) - B.diagonal() * B.diagonal().transpose();	
+		EMatrix A0 = EMatrix::Identity(d,d) + (C1.cwiseProduct(B) - B.diagonal().asDiagonal() * C1.transpose()).cwiseQuotient(D0+EMatrix::Identity(d,d));
+		EV = A0.inverse() * EV;
+
+		Eigen::Map<EMatrix> C0(C.get_matrix(0),d,d);
+		EMatrix Raux = EV * C0 * EV.transpose();
+		EMatrix aux = Raux.diagonal().cwiseAbs().cwiseSqrt().asDiagonal().inverse();
+		EV = aux * EV;
+
+		crit.push_back(0.0);
+		for(int l = 0; l < L; l++)
+		{
+			Eigen::Map<EMatrix> Ci(C.get_matrix(l),d,d);
+			Eigen::Map<EMatrix> Csi(Cs.get_matrix(l),d,d);	
+			Csi = EV * Ci * EV.transpose();
+			Rs.col(l) = Csi.diagonal();
+			crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();	
+		}
+
+		improve = CMath::abs(crit.back() - crit[iter]);
+		iter++;
+	}
+
+	if(iter == itermax)
+	{
+		SG_SERROR("Convergence not reached\n")
+	}
+
+	return V;
+
+}
+
+#endif //HAVE_EIGEN3

src/shogun/mathematics/ajd/UWedge.h

@@ -0,0 +1,80 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Written (W) 2013 Kevin Hughes
+ *
+ * Thanks to Andreas Ziehe and Cedric Gouy-Pailler
+ */
+
+#ifndef UWEDGE_H_
+#define UWEDGE_H_
+
+#ifdef HAVE_EIGEN3
+
+#include <shogun/mathematics/ajd/ApproxJointDiagonalizer.h>
+
+#include <shogun/lib/common.h>
+#include <shogun/base/SGObject.h>
+#include <shogun/lib/SGMatrix.h>
+#include <shogun/lib/SGNDArray.h>
+
+#include <shogun/mathematics/Math.h>
+
+namespace shogun
+{
+
+/** @brief Class FFDiag
+ *
+ * An Approximate Joint Diagonalizer (AJD) Implementation
+ */
+class CUWedge : public CApproxJointDiagonalizer
+{
+	public:
+
+		/** constructor */
+		CUWedge()
+		{
+		};
+
+		/** destructor */
+		virtual ~CUWedge()
+		{
+		}
+
+		/** Computes the matrix V that best diagonalizes C 
+		 * @param C the set of matrices to be diagonalized
+		 * @param V0 an estimate of the matrix V
+		 * @param eps machine epsilon or desired epsilon
+		 * @param itermax maximum number of iterations
+		 * @return V the matrix the best diagonalizes C 
+		 */
+		static SGMatrix<float64_t> diagonalize(SGNDArray<float64_t> &C,
+       							SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
+							double eps=CMath::MACHINE_EPSILON,
+							int itermax=200);
+
+		/** Computes the matrix V that best diagonalizes C 
+		 * @param C the set of matrices to be diagonalized
+		 * @param V0 an estimate of the matrix V
+		 * @param eps machine epsilon or desired epsilon
+		 * @param itermax maximum number of iterations
+		 * @return V the matrix the best diagonalizes C 
+		 */
+		virtual SGMatrix<float64_t> compute(SGNDArray<float64_t> &C,
+						   SGMatrix<float64_t> V0= SGMatrix<float64_t>(NULL,0,0),
+						   double eps=CMath::MACHINE_EPSILON,
+						   int itermax=200)
+		{
+			m_V = diagonalize(C,V0,eps,itermax);
+			return m_V;	
+		}
+
+		/** @return object name */
+		virtual const char* get_name() const { return "UWedge"; }
+};
+}
+#endif //HAVE_EIGEN3
+#endif //UWEDGE_H_ 

tests/unit/mathematics/ajd/UWedge_unittest.cc

@@ -0,0 +1,75 @@
+#include <shogun/base/init.h>
+#include <shogun/lib/common.h>
+#include <gtest/gtest.h>
+
+#ifdef HAVE_EIGEN3
+
+#include <shogun/lib/SGVector.h>
+#include <shogun/lib/SGMatrix.h>
+#include <shogun/lib/SGNDArray.h>
+
+#include <shogun/mathematics/Math.h>
+#include <shogun/mathematics/eigen3.h>
+#include <shogun/mathematics/ajd/UWedge.h>
+
+#include <shogun/evaluation/ica/PermutationMatrix.h>
+
+using namespace Eigen;
+
+typedef Matrix< float64_t, Dynamic, Dynamic, ColMajor > EMatrix;
+typedef Matrix< float64_t, Dynamic, 1, ColMajor > EVector;
+
+using namespace shogun;
+
+TEST(CUWedge, diagonalize)
+{
+	// Generating diagonal matrices
+	index_t * C_dims = SG_MALLOC(index_t, 3);
+	C_dims[0] = 10;
+	C_dims[1] = 10;
+	C_dims[2] = 30;
+	SGNDArray< float64_t > C(C_dims, 3);
+
+	CMath::init_random(17);
+
+	for (int i = 0; i < C_dims[2]; i++)
+	{
+		Eigen::Map<EMatrix> tmp(C.get_matrix(i),C_dims[0], C_dims[1]);
+		tmp.setIdentity();
+
+		for (int j = 0; j < C_dims[0]; j++)
+		{
+			tmp(j,j) *= CMath::abs(CMath::random(1,5)); 
+		}
+	}
+
+	// Mixing and demixing matrices
+	EMatrix B(C_dims[0], C_dims[1]); 
+	B.setRandom();
+	EMatrix A = B.inverse();
+
+	for (int i = 0; i < C_dims[2]; i++)
+	{
+		Eigen::Map<EMatrix> Ci(C.get_matrix(i),C_dims[0], C_dims[1]);
+		Ci = A * Ci * A.transpose();
+	}	
+
+	/** Diagonalize **/
+	SGMatrix<float64_t> V = CUWedge::diagonalize(C);
+
+	// Test output size
+	EXPECT_EQ(V.num_rows, C_dims[0]);
+	EXPECT_EQ(V.num_cols, C_dims[1]);
+
+	// Close to a permutation matrix (with random scales)
+	Eigen::Map<EMatrix> EV(V.matrix,C_dims[0], C_dims[1]);
+	SGMatrix<float64_t> P(C_dims[0],C_dims[1]);
+	Eigen::Map<EMatrix> EP(P.matrix,C_dims[0], C_dims[1]);
+	EP = EV * A;
+
+	// Test if output is correct
+	bool isperm = is_permutation_matrix(P);
+	EXPECT_EQ(isperm,true);
+}
+
+#endif //HAVE_EIGEN3