SHOGUN
v3.2.0
|
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