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