11 using namespace Eigen;
14 double eps,
int itermax)
28 for (
int i = 0; i < N; i++)
30 for (
int j = 0; j < N; j++)
35 std::vector<float64_t> p(T,1.0/T);
37 MatrixXd C0 = MatrixXd::Identity(N,N);
40 EV = EV * VectorXd::Ones(EV.rows()).cwiseQuotient((EV.transpose() * C0 * EV).cwiseSqrt().diagonal()).asDiagonal();
43 for (
int i = 0; i < N; i++)
46 std::vector<bool> issymmetric;
47 issymmetric.reserve(T);
48 for (
int l = 0; l < T; l++)
52 Ci = P * Ci * P.transpose();
54 if ( (Ci - Ci.transpose()).sum() > 1e-6 )
55 issymmetric[l] =
false;
57 issymmetric[l] =
true;
60 C0 = P * C0 * P.transpose();
65 for (
int t = 0; t < T; t++)
72 D = D + 2*p[t] * M1 * M1.transpose();
77 D = D + p[t] * (M1*M1.transpose() + M2*M2.transpose());
83 std::vector<float64_t> crit;
84 while ( iter < itermax && deltacrit > eps )
88 for (
int i = 0; i < N; i++)
90 VectorXd w = EV.col(i);
92 for (
int t = 0; t < T; t++)
99 D = D - 2*p[t] * m1 * m1.transpose();
103 VectorXd m2 = Ci.transpose() * w;
104 D = D - p[t] * (m1*m1.transpose() + m2*m2.transpose());
108 EigenSolver<MatrixXd> eig;
112 MatrixXd eigenvectors = eig.pseudoEigenvectors();
113 VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
119 for (
int j = 1; j < D.rows(); j++)
121 if( eigenvalues[j] > eigenvalues[j-1] )
123 std::swap(eigenvalues[j],eigenvalues[j-1]);
124 eigenvectors.col(j).swap(eigenvectors.col(j-1));
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())));
134 for (
int t = 0; t < T; t++)
138 VectorXd m1 = Ci * w_new;
141 D = D + 2*p[t] * m1 * m1.transpose();
145 VectorXd m2 = Ci.transpose() * w_new;
146 D = D + p[t] * (m1*m1.transpose() + m2*m2.transpose());
154 EV = EV * (EV.transpose() * C0 * EV).diagonal().cwiseSqrt().asDiagonal().inverse();
155 for (
int t = 0; t < T; t++)
158 MatrixXd eD = EV.transpose() * Ci * EV;
159 eD.diagonal() = VectorXd::Zero(eD.rows());
160 crit.back() = crit.back() + p[t]*eD.cwiseAbs2().sum();
162 crit.back() = crit.back() / (N*N - N);
165 deltacrit =
CMath::abs( crit[iter] - crit[iter-1] );
170 EV = (P.transpose() * EV).transpose();
static float64_t randn_double()
T * get_matrix(index_t matIdx) const
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)
all of classes and functions are contained in the shogun namespace
Matrix::Scalar max(Matrix m)