SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
JADiagOrth.cpp
Go to the documentation of this file.
00001 #ifdef HAVE_EIGEN3
00002 
00003 #include <shogun/mathematics/ajd/JADiagOrth.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 float64_t givens_stack(float64_t *A, int M, int K, int p, int q);
00014 void left_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s);
00015 void right_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s);
00016 void left_rot_simple(float64_t *A, int m, int n, int p, int q, float64_t c, float64_t s);
00017 
00018 SGMatrix<float64_t> CJADiagOrth::diagonalize(SGNDArray<float64_t> C, SGMatrix<float64_t> V0,
00019                         double eps, int itermax)
00020 {
00021     int m = C.dims[0];
00022     int L = C.dims[2];
00023 
00024     SGMatrix<float64_t> V;
00025     if (V0.num_rows == m && V0.num_cols == m)
00026         V = V0.clone();
00027     else
00028         V = SGMatrix<float64_t>::create_identity_matrix(m,1);
00029 
00030     bool more = true;
00031     int rots = 0;
00032 
00033     while (more)
00034     {
00035         more = false;
00036 
00037         for (int p = 0; p < m; p++)
00038         {
00039             for (int q = p+1; q < m; q++)
00040             {
00041                 // computation of Givens angle
00042                 float64_t theta = givens_stack(C.array, m, L, p, q);
00043 
00044                 // Givens update
00045                 if (fabs(theta) > eps)
00046                 {
00047                     float64_t c = cos(theta);
00048                     float64_t s = sin(theta);
00049                     left_rot_stack (C.array, m, m, L, p, q, c, s);
00050                     right_rot_stack(C.array, m, m, L, p, q, c, s);
00051                     left_rot_simple(V.matrix, m, m, p, q, c, s);
00052                     rots++;
00053                     more = true;
00054                 }
00055             }
00056         }
00057     }
00058 
00059     return V;
00060 }
00061 
00062 /* Givens angle for the pair (p,q) of a stack of K M*M matrices */
00063 float64_t givens_stack(float64_t *A, int M, int K, int p, int q)
00064 {
00065     int k;
00066     float64_t diff_on, sum_off, ton, toff;
00067     float64_t *cm; // A cumulant matrix
00068     float64_t G11 = 0.0;
00069     float64_t G12 = 0.0;
00070     float64_t G22 = 0.0;
00071 
00072     int M2 = M*M;
00073     int pp = p+p*M;
00074     int pq = p+q*M;
00075     int qp = q+p*M;
00076     int qq = q+q*M;
00077 
00078     for (k=0, cm=A; k<K; k++, cm+=M2)
00079     {
00080         diff_on = cm[pp] - cm[qq];
00081         sum_off = cm[pq] + cm[qp];
00082 
00083         G11 += diff_on * diff_on;
00084         G22 += sum_off * sum_off;
00085         G12 += diff_on * sum_off;
00086     }
00087 
00088     ton  = G11 - G22;
00089     toff = 2.0 * G12;
00090 
00091     return -0.5 * atan2 ( toff , ton+sqrt(ton*ton+toff*toff) );
00092 }
00093 
00094 /*
00095    Ak(mxn) --> R * Ak(mxn) where R rotates the (p,q) rows R =[ c -s ; s c ]
00096    and Ak is the k-th matrix in the stack
00097 */
00098 void left_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s )
00099 {
00100     int k, ix, iy, cpt;
00101     int MN = M*N;
00102     int kMN;
00103     float64_t nx, ny;
00104 
00105     for (k=0, kMN=0; k<K; k++, kMN+=MN)
00106     {
00107         for (cpt=0, ix=p+kMN, iy=q+kMN; cpt<N; cpt++, ix+=M, iy+=M)
00108         {
00109             nx = A[ix];
00110             ny = A[iy];
00111             A[ix] = c*nx - s*ny;
00112             A[iy] = s*nx + c*ny;
00113         }
00114     }
00115 }
00116 
00117 /* Ak(mxn) --> Ak(mxn) x R where R rotates the (p,q) columns R =[ c s ; -s c ]
00118    and Ak is the k-th M*N matrix in the stack */
00119 void right_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s )
00120 {
00121     int k, ix, iy, cpt, kMN;
00122     int pM = p*M;
00123     int qM = q*M;
00124     float64_t nx, ny;
00125 
00126     for (k=0, kMN=0; k<K; k++, kMN+=M*N)
00127     {
00128         for (cpt=0, ix=pM+kMN, iy=qM+kMN; cpt<M; cpt++)
00129         {
00130             nx = A[ix];
00131             ny = A[iy];
00132             A[ix++] = c*nx - s*ny;
00133             A[iy++] = s*nx + c*ny;
00134         }
00135     }
00136 }
00137 
00138 /*
00139    A(mxn) --> R * A(mxn) where R=[ c -s ; s c ]   rotates the (p,q) rows of R
00140 */
00141 void left_rot_simple(float64_t *A, int m, int n, int p, int q, float64_t c, float64_t s)
00142 {
00143     int ix = p;
00144     int iy = q;
00145     float64_t nx, ny;
00146 
00147     for (int j = 0; j < n; j++, ix+=m, iy+=m)
00148     {
00149         nx = A[ix];
00150         ny = A[iy];
00151         A[ix] = c*nx - s*ny;
00152         A[iy] = s*nx + c*ny;
00153     }
00154 }
00155 #endif //HAVE_EIGEN3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation