SHOGUN
v3.2.0
|
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