10 using namespace Eigen;
13 double eps,
int itermax)
27 for (
int i = 0; i < N; i++)
29 for (
int j = 0; j < N; j++)
34 std::vector<float64_t> p(T,1.0/T);
36 MatrixXd C0 = MatrixXd::Identity(N,N);
39 EV = EV * VectorXd::Ones(EV.rows()).cwiseQuotient((EV.transpose() * C0 * EV).cwiseSqrt().diagonal()).asDiagonal();
42 for (
int i = 0; i < N; i++)
45 std::vector<bool> issymmetric(T);
46 for (
int l = 0; l < T; l++)
50 Ci = P * Ci * P.transpose();
52 if ( (Ci - Ci.transpose()).sum() > 1e-6 )
53 issymmetric[l] =
false;
55 issymmetric[l] =
true;
58 C0 = P * C0 * P.transpose();
63 for (
int t = 0; t < T; t++)
70 D = D + 2*p[t] * M1 * M1.transpose();
75 D = D + p[t] * (M1*M1.transpose() + M2*M2.transpose());
81 std::vector<float64_t> crit;
82 while ( iter < itermax && deltacrit > eps )
86 for (
int i = 0; i < N; i++)
88 VectorXd w = EV.col(i);
90 for (
int t = 0; t < T; t++)
97 D = D - 2*p[t] * m1 * m1.transpose();
101 VectorXd m2 = Ci.transpose() * w;
102 D = D - p[t] * (m1*m1.transpose() + m2*m2.transpose());
106 EigenSolver<MatrixXd> eig;
110 MatrixXd eigenvectors = eig.pseudoEigenvectors();
111 VectorXd eigenvalues = eig.pseudoEigenvalueMatrix().diagonal();
117 for (
int j = 1; j < D.rows(); j++)
119 if( eigenvalues[j] > eigenvalues[j-1] )
121 std::swap(eigenvalues[j],eigenvalues[j-1]);
122 eigenvectors.col(j).swap(eigenvectors.col(j-1));
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())));
132 for (
int t = 0; t < T; t++)
136 VectorXd m1 = Ci * w_new;
139 D = D + 2*p[t] * m1 * m1.transpose();
143 VectorXd m2 = Ci.transpose() * w_new;
144 D = D + p[t] * (m1*m1.transpose() + m2*m2.transpose());
152 EV = EV * (EV.transpose() * C0 * EV).diagonal().cwiseSqrt().asDiagonal().inverse();
153 for (
int t = 0; t < T; t++)
156 MatrixXd eD = EV.transpose() * Ci * EV;
157 eD.diagonal() = VectorXd::Zero(eD.rows());
158 crit.back() = crit.back() + p[t]*eD.cwiseAbs2().sum();
160 crit.back() = crit.back() / (N*N - N);
163 deltacrit =
CMath::abs( crit[iter] - crit[iter-1] );
168 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)