SHOGUN
v3.2.0
|
00001 #ifdef HAVE_EIGEN3 00002 00003 #include <shogun/mathematics/ajd/QDiag.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> CQDiag::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0, 00014 double eps, int itermax) 00015 { 00016 int N = C.dims[0]; 00017 int T = C.dims[2]; 00018 00019 SGMatrix<float64_t> V; 00020 if (V0.num_rows == N && V0.num_cols == N) 00021 { 00022 V = V0.clone(); 00023 } 00024 else 00025 { 00026 V = SGMatrix<float64_t>(N,N); 00027 00028 for (int i = 0; i < N; i++) 00029 { 00030 for (int j = 0; j < N; j++) 00031 V(i,j) = CMath::randn_double(); 00032 } 00033 } 00034 00035 std::vector<float64_t> p(T,1.0/T); 00036 00037 MatrixXd C0 = MatrixXd::Identity(N,N); 00038 00039 Map<MatrixXd> EV(V.matrix,N,N); 00040 EV = EV * VectorXd::Ones(EV.rows()).cwiseQuotient((EV.transpose() * C0 * EV).cwiseSqrt().diagonal()).asDiagonal(); 00041 00042 MatrixXd P = MatrixXd::Zero(N,N); 00043 for (int i = 0; i < N; i++) 00044 P(i,N-1-i) = 1; 00045 00046 std::vector<bool> issymmetric; 00047 issymmetric.reserve(T); 00048 for (int l = 0; l < T; l++) 00049 { 00050 Map<MatrixXd> Ci(C.get_matrix(l),N,N); 00051 00052 Ci = P * Ci * P.transpose(); 00053 00054 if ( (Ci - Ci.transpose()).sum() > 1e-6 ) 00055 issymmetric[l] = false; 00056 else 00057 issymmetric[l] = true; 00058 } 00059 00060 C0 = P * C0 * P.transpose(); 00061 EV = P * EV; 00062 00063 // initialisations for OKN^3 00064 MatrixXd D = MatrixXd::Zero(N,N); 00065 for (int t = 0; t < T; t++) 00066 { 00067 Map<MatrixXd> Ci(C.get_matrix(t),N,N); 00068 MatrixXd M1 = Ci * EV; 00069 00070 if (issymmetric[t]) 00071 { 00072 D = D + 2*p[t] * M1 * M1.transpose(); 00073 } 00074 else 00075 { 00076 MatrixXd M2 = Ci.transpose() * EV; 00077 D = D + p[t] * (M1*M1.transpose() + M2*M2.transpose()); 00078 } 00079 } 00080 00081 int iter = 0; 00082 float64_t deltacrit = 1.0; 00083 std::vector<float64_t> crit; 00084 while ( iter < itermax && deltacrit > eps ) 00085 { 00086 float64_t delta_w = 0.0; 00087 00088 for (int i = 0; i < N; i++) 00089 { 00090 VectorXd w = EV.col(i); 00091 00092 for (int t = 0; t < T; t++) 00093 { 00094 Map<MatrixXd> Ci(C.get_matrix(t),N,N); 00095 VectorXd m1 = Ci * w; 00096 00097 if (issymmetric[t]) 00098 { 00099 D = D - 2*p[t] * m1 * m1.transpose(); 00100 } 00101 else 00102 { 00103 VectorXd m2 = Ci.transpose() * w; 00104 D = D - p[t] * (m1*m1.transpose() + m2*m2.transpose()); 00105 } 00106 } 00107 00108 EigenSolver<MatrixXd> eig; 00109 eig.compute(D); 00110 00111 // sort eigenvectors 00112 MatrixXd eigenvectors = eig.pseudoEigenvectors(); 00113 VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal(); 00114 00115 bool swap = false; 00116 do 00117 { 00118 swap = false; 00119 for (int j = 1; j < D.rows(); j++) 00120 { 00121 if( eigenvalues[j] > eigenvalues[j-1] ) 00122 { 00123 std::swap(eigenvalues[j],eigenvalues[j-1]); 00124 eigenvectors.col(j).swap(eigenvectors.col(j-1)); 00125 swap = true; 00126 } 00127 } 00128 00129 } while(swap); 00130 00131 VectorXd w_new = eigenvectors.col(N-1); 00132 delta_w = std::max(delta_w, std::min(sqrt((w-w_new).cwiseAbs2().sum()), sqrt((w+w_new).cwiseAbs2().sum()))); 00133 00134 for (int t = 0; t < T; t++) 00135 { 00136 Map<MatrixXd> Ci(C.get_matrix(t),N,N); 00137 00138 VectorXd m1 = Ci * w_new; 00139 if (issymmetric[t]) 00140 { 00141 D = D + 2*p[t] * m1 * m1.transpose(); 00142 } 00143 else 00144 { 00145 VectorXd m2 = Ci.transpose() * w_new; 00146 D = D + p[t] * (m1*m1.transpose() + m2*m2.transpose()); 00147 } 00148 } 00149 EV.col(i) = w_new; 00150 } 00151 00152 // err 00153 crit.push_back(0.0); 00154 EV = EV * (EV.transpose() * C0 * EV).diagonal().cwiseSqrt().asDiagonal().inverse(); 00155 for (int t = 0; t < T; t++) 00156 { 00157 Map<MatrixXd> Ci(C.get_matrix(t),N,N); 00158 MatrixXd eD = EV.transpose() * Ci * EV; 00159 eD.diagonal() = VectorXd::Zero(eD.rows()); 00160 crit.back() = crit.back() + p[t]*eD.cwiseAbs2().sum(); 00161 } 00162 crit.back() = crit.back() / (N*N - N); 00163 00164 if (iter > 1) 00165 deltacrit = CMath::abs( crit[iter] - crit[iter-1] ); 00166 00167 iter++; 00168 } 00169 00170 EV = (P.transpose() * EV).transpose(); 00171 00172 if (iter == itermax) 00173 SG_SERROR("Convergence not reached\n") 00174 00175 return V; 00176 } 00177 #endif //HAVE_EIGEN3