SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QDiag.cpp
Go to the documentation of this file.
1 #ifdef HAVE_EIGEN3
2 
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 N = C.dims[0];
17  int T = C.dims[2];
18 
20  if (V0.num_rows == N && V0.num_cols == N)
21  {
22  V = V0.clone();
23  }
24  else
25  {
26  V = SGMatrix<float64_t>(N,N);
27 
28  for (int i = 0; i < N; i++)
29  {
30  for (int j = 0; j < N; j++)
31  V(i,j) = CMath::randn_double();
32  }
33  }
34 
35  std::vector<float64_t> p(T,1.0/T);
36 
37  MatrixXd C0 = MatrixXd::Identity(N,N);
38 
39  Map<MatrixXd> EV(V.matrix,N,N);
40  EV = EV * VectorXd::Ones(EV.rows()).cwiseQuotient((EV.transpose() * C0 * EV).cwiseSqrt().diagonal()).asDiagonal();
41 
42  MatrixXd P = MatrixXd::Zero(N,N);
43  for (int i = 0; i < N; i++)
44  P(i,N-1-i) = 1;
45 
46  std::vector<bool> issymmetric;
47  issymmetric.reserve(T);
48  for (int l = 0; l < T; l++)
49  {
50  Map<MatrixXd> Ci(C.get_matrix(l),N,N);
51 
52  Ci = P * Ci * P.transpose();
53 
54  if ( (Ci - Ci.transpose()).sum() > 1e-6 )
55  issymmetric[l] = false;
56  else
57  issymmetric[l] = true;
58  }
59 
60  C0 = P * C0 * P.transpose();
61  EV = P * EV;
62 
63  // initialisations for OKN^3
64  MatrixXd D = MatrixXd::Zero(N,N);
65  for (int t = 0; t < T; t++)
66  {
67  Map<MatrixXd> Ci(C.get_matrix(t),N,N);
68  MatrixXd M1 = Ci * EV;
69 
70  if (issymmetric[t])
71  {
72  D = D + 2*p[t] * M1 * M1.transpose();
73  }
74  else
75  {
76  MatrixXd M2 = Ci.transpose() * EV;
77  D = D + p[t] * (M1*M1.transpose() + M2*M2.transpose());
78  }
79  }
80 
81  int iter = 0;
82  float64_t deltacrit = 1.0;
83  std::vector<float64_t> crit;
84  while ( iter < itermax && deltacrit > eps )
85  {
86  float64_t delta_w = 0.0;
87 
88  for (int i = 0; i < N; i++)
89  {
90  VectorXd w = EV.col(i);
91 
92  for (int t = 0; t < T; t++)
93  {
94  Map<MatrixXd> Ci(C.get_matrix(t),N,N);
95  VectorXd m1 = Ci * w;
96 
97  if (issymmetric[t])
98  {
99  D = D - 2*p[t] * m1 * m1.transpose();
100  }
101  else
102  {
103  VectorXd m2 = Ci.transpose() * w;
104  D = D - p[t] * (m1*m1.transpose() + m2*m2.transpose());
105  }
106  }
107 
108  EigenSolver<MatrixXd> eig;
109  eig.compute(D);
110 
111  // sort eigenvectors
112  MatrixXd eigenvectors = eig.pseudoEigenvectors();
113  VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
114 
115  bool swap = false;
116  do
117  {
118  swap = false;
119  for (int j = 1; j < D.rows(); j++)
120  {
121  if( eigenvalues[j] > eigenvalues[j-1] )
122  {
123  std::swap(eigenvalues[j],eigenvalues[j-1]);
124  eigenvectors.col(j).swap(eigenvectors.col(j-1));
125  swap = true;
126  }
127  }
128 
129  } while(swap);
130 
131  VectorXd w_new = eigenvectors.col(N-1);
132  delta_w = std::max(delta_w, std::min(sqrt((w-w_new).cwiseAbs2().sum()), sqrt((w+w_new).cwiseAbs2().sum())));
133 
134  for (int t = 0; t < T; t++)
135  {
136  Map<MatrixXd> Ci(C.get_matrix(t),N,N);
137 
138  VectorXd m1 = Ci * w_new;
139  if (issymmetric[t])
140  {
141  D = D + 2*p[t] * m1 * m1.transpose();
142  }
143  else
144  {
145  VectorXd m2 = Ci.transpose() * w_new;
146  D = D + p[t] * (m1*m1.transpose() + m2*m2.transpose());
147  }
148  }
149  EV.col(i) = w_new;
150  }
151 
152  // err
153  crit.push_back(0.0);
154  EV = EV * (EV.transpose() * C0 * EV).diagonal().cwiseSqrt().asDiagonal().inverse();
155  for (int t = 0; t < T; t++)
156  {
157  Map<MatrixXd> Ci(C.get_matrix(t),N,N);
158  MatrixXd eD = EV.transpose() * Ci * EV;
159  eD.diagonal() = VectorXd::Zero(eD.rows());
160  crit.back() = crit.back() + p[t]*eD.cwiseAbs2().sum();
161  }
162  crit.back() = crit.back() / (N*N - N);
163 
164  if (iter > 1)
165  deltacrit = CMath::abs( crit[iter] - crit[iter-1] );
166 
167  iter++;
168  }
169 
170  EV = (P.transpose() * EV).transpose();
171 
172  if (iter == itermax)
173  SG_SERROR("Convergence not reached\n")
174 
175  return V;
176 }
177 #endif //HAVE_EIGEN3

SHOGUN Machine Learning Toolbox - Documentation