SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
JediDiag.cpp
Go to the documentation of this file.
00001 #ifdef HAVE_EIGEN3
00002 
00003 #include <shogun/mathematics/ajd/JediDiag.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 sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber,
00014                 float64_t *s_max, float64_t *sh_max, float64_t *A);
00015 
00016 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm,
00017             float64_t *s_max, float64_t *sh_max, float64_t *A);
00018 
00019 SGMatrix<float64_t> CJediDiag::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
00020                         double eps, int itermax)
00021 {
00022     int d = C.dims[0];
00023     int L = C.dims[2];
00024 
00025     SGMatrix<float64_t> V;
00026     if (V0.num_rows == d && V0.num_cols == d)
00027         V = V0.clone();
00028     else
00029         V = SGMatrix<float64_t>::create_identity_matrix(d,1);
00030 
00031     int iter = 0;
00032     float64_t sh_max = 1;
00033     float64_t s_max = 1;
00034     while (iter < itermax && ( (sh_max>eps) || (s_max>eps) ))
00035     {
00036         sh_max = 0;
00037         s_max = 0;
00038         sweepjedi(C.array,
00039                     &d, &L,
00040                     &s_max, &sh_max,
00041                     V.matrix);
00042         iter++;
00043     }
00044 
00045     if (iter == itermax)
00046         SG_SERROR("Convergence not reached\n")
00047 
00048     Map<MatrixXd> EV(V.matrix,d,d);
00049     EV = EV.inverse();
00050 
00051     return V;
00052 }
00053 
00054 void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber,
00055                 float64_t *s_max, float64_t *sh_max, float64_t *A)
00056 {
00057     for (int n=2;n<=*pMatSize;n++)
00058         for (int m=1;m<n;m++)
00059             iterJDI(C, pMatSize, pMatNumber, &n, &m, s_max, sh_max, A);
00060 
00061     int MS=*pMatSize;
00062     int MN=*pMatNumber;
00063     float64_t col_norm[MS];
00064 
00065     for (int i=0;i<MS;i++)
00066     {
00067         col_norm[i] = 0;
00068 
00069         for (int j=0;j<MS;j++)
00070             col_norm[i] += A[i*MS+j]*A[i*MS+j];
00071 
00072         col_norm[i] = sqrt(col_norm[i]);
00073     }
00074 
00075     float64_t daux=1;
00076     float64_t d[MS];
00077 
00078     for (int i=0;i<MS;i++)
00079         daux*=col_norm[i];
00080 
00081     daux=pow((float64_t)daux,(float64_t) 1/MS);
00082 
00083     for (int i=0;i<MS;i++)
00084         d[i]=daux/col_norm[i];
00085 
00086     for (int i=0;i<MS;i++)
00087         for (int j=0;j<MS;j++)
00088             A[j*MS+i] *= d[j];
00089 
00090     for (int k=0;k<MN;k++)
00091     {
00092         for (int i=0;i<MS;i++)
00093         {
00094             for (int j=0;j<MS;j++)
00095                 C[k*MS*MS+i*MS+j] *= (1/d[i])*(1/d[j]);
00096         }
00097     }
00098 }
00099 
00100 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm,
00101             float64_t *s_max, float64_t *sh_max, float64_t *A)
00102 {
00103     int n=*ptn-1;
00104     int m=*ptm-1;
00105     int MN=*pMatNumber;
00106     int MS=*pMatSize;
00107 
00108     float64_t tmm[MN];
00109     float64_t tnn[MN];
00110     float64_t tmn[MN];
00111     for (int i = 0, d3 = 0; i < MN; i++, d3+=MS*MS)
00112     {
00113         tmm[i] = C[d3+m*MS+m];
00114         tnn[i] = C[d3+n*MS+n];
00115         tmn[i] = C[d3+n*MS+m];
00116     }
00117 
00118     // here we evd
00119     float64_t G[MN][3];
00120     float64_t evectors[9], evalues[3];
00121     for (int i = 0; i < MN; i++)
00122     {
00123         G[i][0] = 0.5*(tmm[i]+tnn[i]);
00124         G[i][1] = 0.5*(tmm[i]-tnn[i]);
00125         G[i][2] = tmn[i];
00126     }
00127     float64_t GG[9];
00128     for (int i = 0; i < 3; i++)
00129     {
00130         for (int j = 0; j <= i; j++)
00131         {
00132             GG[3*j+i] = 0;
00133 
00134             for (int k = 0; k < MN; k++)
00135                 GG[3*j+i] += G[k][i]*G[k][j];
00136 
00137             GG[3*i+j] = GG[3*j+i];
00138         }
00139     }
00140 
00141     for (int i = 0; i < 3; i++)
00142         GG[3*i] *= -1;
00143 
00144     Map<MatrixXd> EGG(GG,3,3);
00145 
00146     EigenSolver<MatrixXd> eig;
00147     eig.compute(EGG);
00148 
00149     MatrixXd eigenvectors = eig.pseudoEigenvectors();
00150     VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
00151 
00152     eigenvectors.col(0) = eigenvectors.col(0).normalized();
00153     eigenvectors.col(1) = eigenvectors.col(1).normalized();
00154     eigenvectors.col(2) = eigenvectors.col(2).normalized();
00155 
00156     memcpy(evectors, eigenvectors.data(), 9*sizeof(float64_t));
00157     memcpy(evalues, eigenvalues.data(), 3*sizeof(float64_t));
00158 
00159     float64_t tmp_evec[3],tmp_eval;
00160     if (fabs(evalues[1])<fabs(evalues[2]))
00161     {
00162         tmp_eval = evalues[1];
00163         evalues[1] = evalues[2];
00164         evalues[2] = tmp_eval;
00165         memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t));
00166         memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t));
00167         memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t));
00168     }
00169     if (fabs(evalues[0])<fabs(evalues[1]))
00170     {
00171         tmp_eval = evalues[0];
00172         evalues[0] = evalues[1];
00173         evalues[1] = tmp_eval;
00174         memcpy(tmp_evec,evectors,3*sizeof(float64_t));
00175         memcpy(evectors,&evectors[3],3*sizeof(float64_t));
00176         memcpy(&evectors[3],tmp_evec,3*sizeof(float64_t));
00177     }
00178     if (fabs(evalues[1])<fabs(evalues[2]))
00179     {
00180         tmp_eval = evalues[1];
00181         evalues[1] = evalues[2];
00182         evalues[2] = tmp_eval;
00183         memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t));
00184         memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t));
00185         memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t));
00186     }
00187 
00188     float64_t aux[9];
00189     float64_t diagaux[3];
00190     float64_t v[3];
00191     for (int i = 0; i < 9; i++)
00192         aux[i] = evectors[i];
00193 
00194     for (int i = 0; i < 3; i++)
00195         aux[3*i] *= -1;
00196 
00197     for (int i = 0; i < 3; i++)
00198     {
00199         diagaux[i]  = 0;
00200 
00201         for (int j = 0; j < 3; j++)
00202             diagaux[i] += evectors[3*i+j] * aux[3*i+j];
00203 
00204         diagaux[i] = 1/sqrt(fabs(diagaux[i]));
00205     }
00206 
00207     int ind_min=2;
00208 
00209     for (int i = 0; i < 3; i++)
00210         v[i] = evectors[3*ind_min+i]*diagaux[ind_min];
00211 
00212     if (v[2]<0)
00213         for (int i = 0; i < 3; i++)
00214             v[i]*=-1;
00215 
00216     float64_t ch = sqrt( (1+sqrt(1+v[0]*v[0]))/2 );
00217     float64_t sh = v[0]/(2*ch);
00218     float64_t c = sqrt( ( 1 + v[2]/sqrt(1+v[0]*v[0]) )/2 );
00219     float64_t s = -v[1]/( 2*c*sqrt( (1+v[0]*v[0]) ) );
00220     *s_max = std::max(*s_max,fabs(s));
00221     *sh_max = std::max(*sh_max,fabs(sh));
00222 
00223     float64_t rm11=c*ch - s*sh;
00224     float64_t rm12=c*sh - s*ch;
00225     float64_t rm21=c*sh + s*ch;
00226     float64_t rm22=c*ch + s*sh;
00227 
00228     float64_t h_slice1,h_slice2;
00229     float64_t buf[MS][MN];
00230 
00231     for (int i = 0; i < MS ;i++)
00232     {
00233         for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
00234         {
00235             h_slice1 = C[d3+i*MS+m];
00236             h_slice2 = C[d3+i*MS+n];
00237             buf[i][j] = rm11*h_slice1 + rm21*h_slice2;
00238             C[d3+i*MS+n] = rm12*h_slice1 + rm22*h_slice2;
00239             C[d3+i*MS+m] = buf[i][j];
00240         }
00241     }
00242 
00243     for (int i = 0; i < MS; i++)
00244     {
00245         for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
00246             C[d3+m*MS+i] = buf[i][j];
00247     }
00248     for (int i = 0; i < MS; i++)
00249     {
00250         for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
00251             C[d3+n*MS+i] = C[d3+i*MS+n];
00252     }
00253     for (int i = 0; i < MN; i++)
00254     {
00255         C[i*MS*MS+m*MS+m] = (rm11*rm11)*tmm[i]+(2*rm11*rm21)*tmn[i]+(rm21*rm21)*tnn[i];
00256         C[i*MS*MS+n*MS+m] = (rm11*rm12)*tmm[i]+(rm11*rm22+rm12*rm21)*tmn[i]+(rm21*rm22)*tnn[i];
00257         C[i*MS*MS+m*MS+n] = C[i*MS*MS+n*MS+m];
00258         C[i*MS*MS+n*MS+n] = (rm12*rm12)*tmm[i]+(2*rm12*rm22)*tmn[i]+(rm22*rm22)*tnn[i];
00259     }
00260 
00261     float64_t col1;
00262     float64_t col2;
00263 
00264     for (int i = 0; i < MS; i++)
00265     {
00266         col1 = A[m*MS+i];
00267         col2 = A[n*MS+i];
00268         A[m*MS+i] = rm22*col1 - rm12*col2;
00269         A[n*MS+i] = -rm21*col1 + rm11*col2;
00270     }
00271 }
00272 #endif //HAVE_EIGEN3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation