SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
UWedge.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 
14  double eps, int itermax)
15 {
16  int d = C.dims[0];
17  int L = C.dims[2];
18 
20  if (V0.num_rows == d && V0.num_cols == d)
21  {
22  V = V0.clone();
23  }
24  else
25  {
26  Map<MatrixXd> C0(C.get_matrix(0),d,d);
27  EigenSolver<MatrixXd> eig;
28  eig.compute(C0);
29 
30  // sort eigenvectors
31  MatrixXd eigenvectors = eig.pseudoEigenvectors();
32  MatrixXd eigenvalues = eig.pseudoEigenvalueMatrix();
33 
34  bool swap = false;
35  do
36  {
37  swap = false;
38  for (int j = 1; j < d; j++)
39  {
40  if ( eigenvalues(j,j) > eigenvalues(j-1,j-1) )
41  {
42  std::swap(eigenvalues(j,j),eigenvalues(j-1,j-1));
43  eigenvectors.col(j).swap(eigenvectors.col(j-1));
44  swap = true;
45  }
46  }
47 
48  } while(swap);
49 
51  Map<MatrixXd> EV(V.matrix, d,d);
52  EV = eigenvalues.cwiseAbs().cwiseSqrt().inverse() * eigenvectors.transpose();
53  }
54  Map<MatrixXd> EV(V.matrix, d,d);
55 
56  index_t * Cs_dims = SG_MALLOC(index_t, 3);
57  Cs_dims[0] = d;
58  Cs_dims[1] = d;
59  Cs_dims[2] = L;
60  SGNDArray<float64_t> Cs(Cs_dims,3);
61  memcpy(Cs.array, C.array, Cs.dims[0]*Cs.dims[1]*Cs.dims[2]*sizeof(float64_t));
62 
63  MatrixXd Rs(d,L);
64  std::vector<float64_t> crit;
65  crit.push_back(0.0);
66  for (int l = 0; l < L; l++)
67  {
68  Map<MatrixXd> Ci(C.get_matrix(l),d,d);
69  Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
70  Ci = 0.5 * (Ci + Ci.transpose());
71  Csi = EV * Ci * EV.transpose();
72  Rs.col(l) = Csi.diagonal();
73  crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
74  }
75 
76  float64_t iter = 0;
77  float64_t improve = 10;
78  while (improve > eps && iter < itermax)
79  {
80  MatrixXd B = Rs * Rs.transpose();
81 
82  MatrixXd C1 = MatrixXd::Zero(d,d);
83  for (int id = 0; id < d; id++)
84  {
85  // rowSums
86  for (int l = 0; l < L; l++)
87  {
88  Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
89  C1.row(id) += Csi.row(id) * Rs(id,l);
90  }
91  }
92 
93  MatrixXd D0 = B.cwiseProduct(B.transpose()) - B.diagonal() * B.diagonal().transpose();
94  MatrixXd A0 = MatrixXd::Identity(d,d) + (C1.cwiseProduct(B) - B.diagonal().asDiagonal() * C1.transpose()).cwiseQuotient(D0+MatrixXd::Identity(d,d));
95  EV = A0.inverse() * EV;
96 
97  Map<MatrixXd> C0(C.get_matrix(0),d,d);
98  MatrixXd Raux = EV * C0 * EV.transpose();
99  MatrixXd aux = Raux.diagonal().cwiseAbs().cwiseSqrt().asDiagonal().inverse();
100  EV = aux * EV;
101 
102  crit.push_back(0.0);
103  for (int l = 0; l < L; l++)
104  {
105  Map<MatrixXd> Ci(C.get_matrix(l),d,d);
106  Map<MatrixXd> Csi(Cs.get_matrix(l),d,d);
107  Csi = EV * Ci * EV.transpose();
108  Rs.col(l) = Csi.diagonal();
109  crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
110  }
111 
112  improve = CMath::abs(crit.back() - crit[iter]);
113  iter++;
114  }
115 
116  if (iter == itermax)
117  SG_SERROR("Convergence not reached\n")
118 
119  return V;
120 
121 }
122 
123 #endif //HAVE_EIGEN3
int32_t index_t
Definition: common.h:62
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
T * get_matrix(index_t matIdx) const
Definition: SGNDArray.h:72
double float64_t
Definition: common.h:50
static SGMatrix< float64_t > diagonalize(SGNDArray< float64_t > C, SGMatrix< float64_t > V0=SGMatrix< float64_t >(NULL, 0, 0, false), double eps=1e-12, int itermax=200)
Definition: UWedge.cpp:13
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< T > create_identity_matrix(index_t size, T scale)
static T abs(T a)
Definition: Math.h:179

SHOGUN Machine Learning Toolbox - Documentation