SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
JediDiag.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 void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber,
14  float64_t *s_max, float64_t *sh_max, float64_t *A);
15 
16 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm,
17  float64_t *s_max, float64_t *sh_max, float64_t *A);
18 
20  double eps, int itermax)
21 {
22  int d = C.dims[0];
23  int L = C.dims[2];
24 
26  if (V0.num_rows == d && V0.num_cols == d)
27  V = V0.clone();
28  else
30 
31  int iter = 0;
32  float64_t sh_max = 1;
33  float64_t s_max = 1;
34  while (iter < itermax && ( (sh_max>eps) || (s_max>eps) ))
35  {
36  sh_max = 0;
37  s_max = 0;
38  sweepjedi(C.array,
39  &d, &L,
40  &s_max, &sh_max,
41  V.matrix);
42  iter++;
43  }
44 
45  if (iter == itermax)
46  SG_SERROR("Convergence not reached\n")
47 
48  Map<MatrixXd> EV(V.matrix,d,d);
49  EV = EV.inverse();
50 
51  return V;
52 }
53 
54 void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber,
55  float64_t *s_max, float64_t *sh_max, float64_t *A)
56 {
57  for (int n=2;n<=*pMatSize;n++)
58  for (int m=1;m<n;m++)
59  iterJDI(C, pMatSize, pMatNumber, &n, &m, s_max, sh_max, A);
60 
61  int MS=*pMatSize;
62  int MN=*pMatNumber;
63  float64_t col_norm[MS];
64 
65  for (int i=0;i<MS;i++)
66  {
67  col_norm[i] = 0;
68 
69  for (int j=0;j<MS;j++)
70  col_norm[i] += A[i*MS+j]*A[i*MS+j];
71 
72  col_norm[i] = sqrt(col_norm[i]);
73  }
74 
75  float64_t daux=1;
76  float64_t d[MS];
77 
78  for (int i=0;i<MS;i++)
79  daux*=col_norm[i];
80 
81  daux=pow((float64_t)daux,(float64_t) 1/MS);
82 
83  for (int i=0;i<MS;i++)
84  d[i]=daux/col_norm[i];
85 
86  for (int i=0;i<MS;i++)
87  for (int j=0;j<MS;j++)
88  A[j*MS+i] *= d[j];
89 
90  for (int k=0;k<MN;k++)
91  {
92  for (int i=0;i<MS;i++)
93  {
94  for (int j=0;j<MS;j++)
95  C[k*MS*MS+i*MS+j] *= (1/d[i])*(1/d[j]);
96  }
97  }
98 }
99 
100 void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn,int *ptm,
101  float64_t *s_max, float64_t *sh_max, float64_t *A)
102 {
103  int n=*ptn-1;
104  int m=*ptm-1;
105  int MN=*pMatNumber;
106  int MS=*pMatSize;
107 
108  float64_t tmm[MN];
109  float64_t tnn[MN];
110  float64_t tmn[MN];
111  for (int i = 0, d3 = 0; i < MN; i++, d3+=MS*MS)
112  {
113  tmm[i] = C[d3+m*MS+m];
114  tnn[i] = C[d3+n*MS+n];
115  tmn[i] = C[d3+n*MS+m];
116  }
117 
118  // here we evd
119  float64_t G[MN][3];
120  float64_t evectors[9], evalues[3];
121  for (int i = 0; i < MN; i++)
122  {
123  G[i][0] = 0.5*(tmm[i]+tnn[i]);
124  G[i][1] = 0.5*(tmm[i]-tnn[i]);
125  G[i][2] = tmn[i];
126  }
127  float64_t GG[9];
128  for (int i = 0; i < 3; i++)
129  {
130  for (int j = 0; j <= i; j++)
131  {
132  GG[3*j+i] = 0;
133 
134  for (int k = 0; k < MN; k++)
135  GG[3*j+i] += G[k][i]*G[k][j];
136 
137  GG[3*i+j] = GG[3*j+i];
138  }
139  }
140 
141  for (int i = 0; i < 3; i++)
142  GG[3*i] *= -1;
143 
144  Map<MatrixXd> EGG(GG,3,3);
145 
146  EigenSolver<MatrixXd> eig;
147  eig.compute(EGG);
148 
149  MatrixXd eigenvectors = eig.pseudoEigenvectors();
150  VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
151 
152  eigenvectors.col(0) = eigenvectors.col(0).normalized();
153  eigenvectors.col(1) = eigenvectors.col(1).normalized();
154  eigenvectors.col(2) = eigenvectors.col(2).normalized();
155 
156  memcpy(evectors, eigenvectors.data(), 9*sizeof(float64_t));
157  memcpy(evalues, eigenvalues.data(), 3*sizeof(float64_t));
158 
159  float64_t tmp_evec[3],tmp_eval;
160  if (fabs(evalues[1])<fabs(evalues[2]))
161  {
162  tmp_eval = evalues[1];
163  evalues[1] = evalues[2];
164  evalues[2] = tmp_eval;
165  memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t));
166  memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t));
167  memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t));
168  }
169  if (fabs(evalues[0])<fabs(evalues[1]))
170  {
171  tmp_eval = evalues[0];
172  evalues[0] = evalues[1];
173  evalues[1] = tmp_eval;
174  memcpy(tmp_evec,evectors,3*sizeof(float64_t));
175  memcpy(evectors,&evectors[3],3*sizeof(float64_t));
176  memcpy(&evectors[3],tmp_evec,3*sizeof(float64_t));
177  }
178  if (fabs(evalues[1])<fabs(evalues[2]))
179  {
180  tmp_eval = evalues[1];
181  evalues[1] = evalues[2];
182  evalues[2] = tmp_eval;
183  memcpy(tmp_evec,&evectors[3],3*sizeof(float64_t));
184  memcpy(&evectors[3],&evectors[6],3*sizeof(float64_t));
185  memcpy(&evectors[6],tmp_evec,3*sizeof(float64_t));
186  }
187 
188  float64_t aux[9];
189  float64_t diagaux[3];
190  float64_t v[3];
191  for (int i = 0; i < 9; i++)
192  aux[i] = evectors[i];
193 
194  for (int i = 0; i < 3; i++)
195  aux[3*i] *= -1;
196 
197  for (int i = 0; i < 3; i++)
198  {
199  diagaux[i] = 0;
200 
201  for (int j = 0; j < 3; j++)
202  diagaux[i] += evectors[3*i+j] * aux[3*i+j];
203 
204  diagaux[i] = 1/sqrt(fabs(diagaux[i]));
205  }
206 
207  int ind_min=2;
208 
209  for (int i = 0; i < 3; i++)
210  v[i] = evectors[3*ind_min+i]*diagaux[ind_min];
211 
212  if (v[2]<0)
213  for (int i = 0; i < 3; i++)
214  v[i]*=-1;
215 
216  float64_t ch = sqrt( (1+sqrt(1+v[0]*v[0]))/2 );
217  float64_t sh = v[0]/(2*ch);
218  float64_t c = sqrt( ( 1 + v[2]/sqrt(1+v[0]*v[0]) )/2 );
219  float64_t s = -v[1]/( 2*c*sqrt( (1+v[0]*v[0]) ) );
220  *s_max = std::max(*s_max,fabs(s));
221  *sh_max = std::max(*sh_max,fabs(sh));
222 
223  float64_t rm11=c*ch - s*sh;
224  float64_t rm12=c*sh - s*ch;
225  float64_t rm21=c*sh + s*ch;
226  float64_t rm22=c*ch + s*sh;
227 
228  float64_t h_slice1,h_slice2;
229  float64_t buf[MS][MN];
230 
231  for (int i = 0; i < MS ;i++)
232  {
233  for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
234  {
235  h_slice1 = C[d3+i*MS+m];
236  h_slice2 = C[d3+i*MS+n];
237  buf[i][j] = rm11*h_slice1 + rm21*h_slice2;
238  C[d3+i*MS+n] = rm12*h_slice1 + rm22*h_slice2;
239  C[d3+i*MS+m] = buf[i][j];
240  }
241  }
242 
243  for (int i = 0; i < MS; i++)
244  {
245  for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
246  C[d3+m*MS+i] = buf[i][j];
247  }
248  for (int i = 0; i < MS; i++)
249  {
250  for (int j = 0, d3 = 0; j < MN; j++, d3+=MS*MS)
251  C[d3+n*MS+i] = C[d3+i*MS+n];
252  }
253  for (int i = 0; i < MN; i++)
254  {
255  C[i*MS*MS+m*MS+m] = (rm11*rm11)*tmm[i]+(2*rm11*rm21)*tmn[i]+(rm21*rm21)*tnn[i];
256  C[i*MS*MS+n*MS+m] = (rm11*rm12)*tmm[i]+(rm11*rm22+rm12*rm21)*tmn[i]+(rm21*rm22)*tnn[i];
257  C[i*MS*MS+m*MS+n] = C[i*MS*MS+n*MS+m];
258  C[i*MS*MS+n*MS+n] = (rm12*rm12)*tmm[i]+(2*rm12*rm22)*tmn[i]+(rm22*rm22)*tnn[i];
259  }
260 
261  float64_t col1;
262  float64_t col2;
263 
264  for (int i = 0; i < MS; i++)
265  {
266  col1 = A[m*MS+i];
267  col2 = A[n*MS+i];
268  A[m*MS+i] = rm22*col1 - rm12*col2;
269  A[n*MS+i] = -rm21*col1 + rm11*col2;
270  }
271 }
272 #endif //HAVE_EIGEN3
void sweepjedi(float64_t *C, int *pMatSize, int *pMatNumber, float64_t *s_max, float64_t *sh_max, float64_t *A)
Definition: JediDiag.cpp:54
Definition: SGMatrix.h:20
SGMatrix< T > clone()
Definition: SGMatrix.cpp:260
index_t num_cols
Definition: SGMatrix.h:378
index_t num_rows
Definition: SGMatrix.h:376
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
#define SG_SERROR(...)
Definition: SGIO.h:179
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: JediDiag.cpp:19
Matrix::Scalar max(Matrix m)
Definition: Redux.h:66
void iterJDI(float64_t *C, int *pMatSize, int *pMatNumber, int *ptn, int *ptm, float64_t *s_max, float64_t *sh_max, float64_t *A)
Definition: JediDiag.cpp:100
static SGMatrix< T > create_identity_matrix(index_t size, T scale)

SHOGUN Machine Learning Toolbox - Documentation