SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
UWedge.cpp
Go to the documentation of this file.
00001 #ifdef HAVE_EIGEN3
00002 
00003 #include <shogun/mathematics/ajd/UWedge.h>
00004 
00005 #include <shogun/base/init.h>
00006 
00007 #include <shogun/mathematics/Math.h>
00008 #include <shogun/mathematics/eigen3.h>
00009 
00010 using namespace shogun;
00011 using namespace Eigen;
00012 
00013 SGMatrix<float64_t> CUWedge::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
00014                     double eps, int itermax)
00015 {
00016     int d = C.dims[0];
00017     int L = C.dims[2];
00018 
00019     SGMatrix<float64_t> V;
00020     if (V0.num_rows == d && V0.num_cols == d)
00021     {
00022         V = V0.clone();
00023     }
00024     else
00025     {
00026         Map<MatrixXd> C0(C.get_matrix(0),d,d);
00027         EigenSolver<MatrixXd> eig;
00028         eig.compute(C0);
00029 
00030         // sort eigenvectors
00031         MatrixXd eigenvectors = eig.pseudoEigenvectors();
00032         MatrixXd eigenvalues = eig.pseudoEigenvalueMatrix();
00033 
00034         bool swap = false;
00035         do
00036         {
00037             swap = false;
00038             for (int j = 1; j < d; j++)
00039             {
00040                 if ( eigenvalues(j,j) > eigenvalues(j-1,j-1) )
00041                 {
00042                     std::swap(eigenvalues(j,j),eigenvalues(j-1,j-1));
00043                     eigenvectors.col(j).swap(eigenvectors.col(j-1));
00044                     swap = true;
00045                 }
00046             }
00047 
00048         } while(swap);
00049 
00050         V = SGMatrix<float64_t>::create_identity_matrix(d,1);
00051         Map<MatrixXd> EV(V.matrix, d,d);
00052         EV = eigenvalues.cwiseAbs().cwiseSqrt().inverse() * eigenvectors.transpose();
00053     }
00054     Map<MatrixXd> EV(V.matrix, d,d);
00055 
00056     index_t * Cs_dims = SG_MALLOC(index_t, 3);
00057     Cs_dims[0] = d;
00058     Cs_dims[1] = d;
00059     Cs_dims[2] = L;
00060     SGNDArray<float64_t> Cs(Cs_dims,3);
00061     memcpy(Cs.array, C.array, Cs.dims[0]*Cs.dims[1]*Cs.dims[2]*sizeof(float64_t));
00062 
00063     MatrixXd Rs(d,L);
00064     std::vector<float64_t> crit;
00065     crit.push_back(0.0);
00066     for (int l = 0; l < L; l++)
00067     {
00068         Map<MatrixXd> Ci(C.get_matrix(l),d,d);
00069         Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
00070         Ci = 0.5 * (Ci + Ci.transpose());
00071         Csi = EV * Ci * EV.transpose();
00072         Rs.col(l) = Csi.diagonal();
00073         crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
00074     }
00075 
00076     float64_t iter = 0;
00077     float64_t improve = 10;
00078     while (improve > eps && iter < itermax)
00079     {
00080         MatrixXd B = Rs * Rs.transpose();
00081 
00082         MatrixXd C1 = MatrixXd::Zero(d,d);
00083         for (int id = 0; id < d; id++)
00084         {
00085             // rowSums
00086             for (int l = 0; l < L; l++)
00087             {
00088                 Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
00089                 C1.row(id) += Csi.row(id) * Rs(id,l);
00090             }
00091         }
00092 
00093         MatrixXd D0 = B.cwiseProduct(B.transpose()) - B.diagonal() * B.diagonal().transpose();
00094         MatrixXd A0 = MatrixXd::Identity(d,d) + (C1.cwiseProduct(B) - B.diagonal().asDiagonal() * C1.transpose()).cwiseQuotient(D0+MatrixXd::Identity(d,d));
00095         EV = A0.inverse() * EV;
00096 
00097         Map<MatrixXd> C0(C.get_matrix(0),d,d);
00098         MatrixXd Raux = EV * C0 * EV.transpose();
00099         MatrixXd aux = Raux.diagonal().cwiseAbs().cwiseSqrt().asDiagonal().inverse();
00100         EV = aux * EV;
00101 
00102         crit.push_back(0.0);
00103         for (int l = 0; l < L; l++)
00104         {
00105             Map<MatrixXd> Ci(C.get_matrix(l),d,d);
00106             Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
00107             Csi = EV * Ci * EV.transpose();
00108             Rs.col(l) = Csi.diagonal();
00109             crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
00110         }
00111 
00112         improve = CMath::abs(crit.back() - crit[iter]);
00113         iter++;
00114     }
00115 
00116     if (iter == itermax)
00117         SG_SERROR("Convergence not reached\n")
00118 
00119     return V;
00120 
00121 }
00122 
00123 #endif //HAVE_EIGEN3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation