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

SHOGUN Machine Learning Toolbox - Documentation