10 using namespace shogun;
11 using namespace Eigen;
14 double eps,
int itermax)
27 EigenSolver<MatrixXd> eig;
31 MatrixXd eigenvectors = eig.pseudoEigenvectors();
32 MatrixXd eigenvalues = eig.pseudoEigenvalueMatrix();
38 for (
int j = 1; j < d; j++)
40 if ( eigenvalues(j,j) > eigenvalues(j-1,j-1) )
42 std::swap(eigenvalues(j,j),eigenvalues(j-1,j-1));
43 eigenvectors.col(j).swap(eigenvectors.col(j-1));
51 Map<MatrixXd> EV(V.
matrix, d,d);
52 EV = eigenvalues.cwiseAbs().cwiseSqrt().inverse() * eigenvectors.transpose();
54 Map<MatrixXd> EV(V.
matrix, d,d);
64 std::vector<float64_t> crit;
66 for (
int l = 0; l < L; l++)
70 Ci = 0.5 * (Ci + Ci.transpose());
71 Csi = EV * Ci * EV.transpose();
72 Rs.col(l) = Csi.diagonal();
73 crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
78 while (improve > eps && iter < itermax)
80 MatrixXd B = Rs * Rs.transpose();
82 MatrixXd C1 = MatrixXd::Zero(d,d);
83 for (
int id = 0;
id < d;
id++)
86 for (
int l = 0; l < L; l++)
89 C1.row(
id) += Csi.row(
id) * Rs(
id,l);
93 MatrixXd D0 = B.cwiseProduct(B.transpose()) - B.diagonal() * B.diagonal().transpose();
94 MatrixXd A0 = MatrixXd::Identity(d,d) + (C1.cwiseProduct(B) - B.diagonal().asDiagonal() * C1.transpose()).cwiseQuotient(D0+MatrixXd::Identity(d,d));
95 EV = A0.inverse() * EV;
98 MatrixXd Raux = EV * C0 * EV.transpose();
99 MatrixXd aux = Raux.diagonal().cwiseAbs().cwiseSqrt().asDiagonal().inverse();
103 for (
int l = 0; l < L; l++)
107 Csi = EV * Ci * EV.transpose();
108 Rs.col(l) = Csi.diagonal();
109 crit.back() += Csi.cwiseAbs2().sum() - Rs.col(l).cwiseAbs2().sum();
112 improve =
CMath::abs(crit.back() - crit[iter]);