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));
52 EV = eigenvalues.cwiseAbs().cwiseSqrt().inverse() * eigenvectors.transpose();
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)
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]);
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=1e-12, int itermax=200)
all of classes and functions are contained in the shogun namespace
static SGMatrix< T > create_identity_matrix(index_t size, T scale)