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

SHOGUN Machine Learning Toolbox - Documentation