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

SHOGUN Machine Learning Toolbox - Documentation