SHOGUN  4.1.0
QDiag.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 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
Definition: SGMatrix.h:20
SGMatrix< T > clone()
Definition: SGMatrix.cpp:260
static float64_t randn_double()
Definition: Math.h:1132
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=CMath::MACHINE_EPSILON, int itermax=200)
Definition: QDiag.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
Matrix::Scalar max(Matrix m)
Definition: Redux.h:66
static T abs(T a)
Definition: Math.h:179

SHOGUN Machine Learning Toolbox - Documentation