SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
JADiagOrth.cpp
Go to the documentation of this file.
2 
3 #ifdef HAVE_EIGEN3
4 
5 #include <shogun/base/init.h>
6 
9 
10 using namespace shogun;
11 using namespace Eigen;
12 
13 float64_t givens_stack(float64_t *A, int M, int K, int p, int q);
14 void left_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s);
15 void right_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s);
16 void left_rot_simple(float64_t *A, int m, int n, int p, int q, float64_t c, float64_t s);
17 
19  double eps, int itermax)
20 {
21  int m = C.dims[0];
22  int L = C.dims[2];
23 
25  if (V0.num_rows == m && V0.num_cols == m)
26  V = V0.clone();
27  else
29 
30  bool more = true;
31  int rots = 0;
32 
33  while (more)
34  {
35  more = false;
36 
37  for (int p = 0; p < m; p++)
38  {
39  for (int q = p+1; q < m; q++)
40  {
41  // computation of Givens angle
42  float64_t theta = givens_stack(C.array, m, L, p, q);
43 
44  // Givens update
45  if (fabs(theta) > eps)
46  {
47  float64_t c = cos(theta);
48  float64_t s = sin(theta);
49  left_rot_stack (C.array, m, m, L, p, q, c, s);
50  right_rot_stack(C.array, m, m, L, p, q, c, s);
51  left_rot_simple(V.matrix, m, m, p, q, c, s);
52  rots++;
53  more = true;
54  }
55  }
56  }
57  }
58 
59  return V;
60 }
61 
62 /* Givens angle for the pair (p,q) of a stack of K M*M matrices */
63 float64_t givens_stack(float64_t *A, int M, int K, int p, int q)
64 {
65  int k;
66  float64_t diff_on, sum_off, ton, toff;
67  float64_t *cm; // A cumulant matrix
68  float64_t G11 = 0.0;
69  float64_t G12 = 0.0;
70  float64_t G22 = 0.0;
71 
72  int M2 = M*M;
73  int pp = p+p*M;
74  int pq = p+q*M;
75  int qp = q+p*M;
76  int qq = q+q*M;
77 
78  for (k=0, cm=A; k<K; k++, cm+=M2)
79  {
80  diff_on = cm[pp] - cm[qq];
81  sum_off = cm[pq] + cm[qp];
82 
83  G11 += diff_on * diff_on;
84  G22 += sum_off * sum_off;
85  G12 += diff_on * sum_off;
86  }
87 
88  ton = G11 - G22;
89  toff = 2.0 * G12;
90 
91  return -0.5 * CMath::atan2 (toff, ton+sqrt(ton*ton+toff*toff));
92 }
93 
94 /*
95  Ak(mxn) --> R * Ak(mxn) where R rotates the (p,q) rows R =[ c -s ; s c ]
96  and Ak is the k-th matrix in the stack
97 */
98 void left_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s )
99 {
100  int k, ix, iy, cpt;
101  int MN = M*N;
102  int kMN;
103  float64_t nx, ny;
104 
105  for (k=0, kMN=0; k<K; k++, kMN+=MN)
106  {
107  for (cpt=0, ix=p+kMN, iy=q+kMN; cpt<N; cpt++, ix+=M, iy+=M)
108  {
109  nx = A[ix];
110  ny = A[iy];
111  A[ix] = c*nx - s*ny;
112  A[iy] = s*nx + c*ny;
113  }
114  }
115 }
116 
117 /* Ak(mxn) --> Ak(mxn) x R where R rotates the (p,q) columns R =[ c s ; -s c ]
118  and Ak is the k-th M*N matrix in the stack */
119 void right_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s )
120 {
121  int k, ix, iy, cpt, kMN;
122  int pM = p*M;
123  int qM = q*M;
124  float64_t nx, ny;
125 
126  for (k=0, kMN=0; k<K; k++, kMN+=M*N)
127  {
128  for (cpt=0, ix=pM+kMN, iy=qM+kMN; cpt<M; cpt++)
129  {
130  nx = A[ix];
131  ny = A[iy];
132  A[ix++] = c*nx - s*ny;
133  A[iy++] = s*nx + c*ny;
134  }
135  }
136 }
137 
138 /*
139  A(mxn) --> R * A(mxn) where R=[ c -s ; s c ] rotates the (p,q) rows of R
140 */
141 void left_rot_simple(float64_t *A, int m, int n, int p, int q, float64_t c, float64_t s)
142 {
143  int ix = p;
144  int iy = q;
145  float64_t nx, ny;
146 
147  for (int j = 0; j < n; j++, ix+=m, iy+=m)
148  {
149  nx = A[ix];
150  ny = A[iy];
151  A[ix] = c*nx - s*ny;
152  A[iy] = s*nx + c*ny;
153  }
154 }
155 #endif //HAVE_EIGEN3
static SGMatrix< float64_t > diagonalize(SGNDArray< float64_t > C, SGMatrix< float64_t > V0=SGMatrix< float64_t >(NULL, 0, 0, false), double eps=CMath::MACHINE_EPSILON, int itermax=200)
Definition: JADiagOrth.cpp:18
Definition: SGMatrix.h:20
SGMatrix< T > clone()
Definition: SGMatrix.cpp:260
static float64_t atan2(float64_t y, float64_t x)
atan(x), x being a complex128_t not implemented
Definition: Math.h:796
void left_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s)
Definition: JADiagOrth.cpp:98
index_t num_cols
Definition: SGMatrix.h:378
index_t num_rows
Definition: SGMatrix.h:376
float64_t givens_stack(float64_t *A, int M, int K, int p, int q)
Definition: JADiagOrth.cpp:63
double float64_t
Definition: common.h:50
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
index_t * dims
Definition: SGNDArray.h:177
void right_rot_stack(float64_t *A, int M, int N, int K, int p, int q, float64_t c, float64_t s)
Definition: JADiagOrth.cpp:119
void left_rot_simple(float64_t *A, int m, int n, int p, int q, float64_t c, float64_t s)
Definition: JADiagOrth.cpp:141
static SGMatrix< T > create_identity_matrix(index_t size, T scale)

SHOGUN Machine Learning Toolbox - Documentation