SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
FFDiag.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation