SHOGUN
v3.2.0
|
00001 #ifdef HAVE_EIGEN3 00002 00003 #include <shogun/mathematics/ajd/FFDiag.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 void getW(float64_t *C, int *ptN, int *ptK, float64_t *W); 00014 00015 SGMatrix<float64_t> CFFDiag::diagonalize(SGNDArray<float64_t> C0, SGMatrix<float64_t> V0, 00016 double eps, int itermax) 00017 { 00018 int n = C0.dims[0]; 00019 int K = C0.dims[2]; 00020 00021 index_t * C_dims = SG_MALLOC(index_t, 3); 00022 C_dims[0] = C0.dims[0]; 00023 C_dims[1] = C0.dims[1]; 00024 C_dims[2] = C0.dims[2]; 00025 SGNDArray<float64_t> C(C_dims,3); 00026 memcpy(C.array, C0.array, C0.dims[0]*C0.dims[1]*C0.dims[2]*sizeof(float64_t)); 00027 00028 SGMatrix<float64_t> V; 00029 if (V0.num_rows == n && V0.num_cols == n) 00030 V = V0.clone(); 00031 else 00032 V = SGMatrix<float64_t>::create_identity_matrix(n,1); 00033 00034 MatrixXd Id(n,n); Id.setIdentity(); 00035 Map<MatrixXd> EV(V.matrix,n,n); 00036 00037 float64_t inum = 0; 00038 float64_t df = 1; 00039 std::vector<float64_t> crit; 00040 while (df > eps && inum < itermax) 00041 { 00042 MatrixXd W = MatrixXd::Zero(n,n); 00043 00044 getW(C.get_matrix(0), 00045 &n, &K, 00046 W.data()); 00047 00048 W.transposeInPlace(); 00049 int e = CMath::ceil(log2(W.array().abs().rowwise().sum().maxCoeff())); 00050 int s = std::max(0,e-1); 00051 W /= pow(2,s); 00052 00053 EV = (Id+W) * EV; 00054 MatrixXd d = MatrixXd::Zero(EV.rows(),EV.cols()); 00055 d.diagonal() = VectorXd::Ones(EV.diagonalSize()).cwiseQuotient((EV * EV.transpose()).diagonal().cwiseSqrt()); 00056 EV = d * EV; 00057 00058 for (int i = 0; i < K; i++) 00059 { 00060 Map<MatrixXd> Ci(C.get_matrix(i), n, n); 00061 Map<MatrixXd> C0i(C0.get_matrix(i), n, n); 00062 Ci = EV * C0i * EV.transpose(); 00063 } 00064 00065 float64_t f = 0; 00066 for (int i = 0; i < K; i++) 00067 { 00068 Map<MatrixXd> C0i(C0.get_matrix(i), n, n); 00069 MatrixXd F = EV * C0i * EV.transpose(); 00070 f += (F.transpose() * F).diagonal().sum() - F.array().pow(2).matrix().diagonal().sum(); 00071 } 00072 00073 crit.push_back(f); 00074 00075 if (inum > 1) 00076 df = CMath::abs(crit[inum-1]-crit[inum]); 00077 00078 inum++; 00079 } 00080 00081 if (inum == itermax) 00082 SG_SERROR("Convergence not reached\n") 00083 00084 return V; 00085 00086 } 00087 00088 void getW(float64_t *C, int *ptN, int *ptK, float64_t *W) 00089 { 00090 int N=*ptN; 00091 int K=*ptK; 00092 int auxij,auxji,auxii,auxjj; 00093 float64_t z[N][N]; 00094 float64_t y[N][N]; 00095 00096 for (int i = 0; i < N; i++) 00097 { 00098 for (int j = 0; j < N; j++) 00099 { 00100 z[i][j] = 0; 00101 y[i][j] = 0; 00102 } 00103 } 00104 00105 for (int i = 0; i < N; i++) 00106 { 00107 for (int j = 0; j < N; j++) 00108 { 00109 for (int k = 0; k < K; k++) 00110 { 00111 auxij = N*N*k+N*i+j; 00112 auxji = N*N*k+N*j+i; 00113 auxii = N*N*k+N*i+i; 00114 auxjj = N*N*k+N*j+j; 00115 z[i][j] += C[auxii]*C[auxjj]; 00116 y[i][j] += 0.5*C[auxjj]*(C[auxij]+C[auxji]); 00117 } 00118 } 00119 } 00120 00121 for (int i = 0; i < N-1; i++) 00122 { 00123 for (int j = i+1; j < N; j++) 00124 { 00125 auxij = N*i+j; 00126 auxji = N*j+i; 00127 W[auxij] = (z[j][i]*y[j][i] - z[i][i]*y[i][j])/(z[j][j]*z[i][i]-z[i][j]*z[i][j]); 00128 W[auxji] = (z[i][j]*y[i][j] - z[j][j]*y[j][i])/(z[j][j]*z[i][i]-z[i][j]*z[i][j]); 00129 } 00130 } 00131 00132 return; 00133 } 00134 #endif //HAVE_EIGEN3