10 using namespace shogun;
11 using namespace Eigen;
16 double eps,
int itermax)
22 C_dims[0] = C0.
dims[0];
23 C_dims[1] = C0.
dims[1];
24 C_dims[2] = C0.
dims[2];
34 MatrixXd Id(n,n); Id.setIdentity();
35 Map<MatrixXd> EV(V.
matrix,n,n);
39 std::vector<float64_t> crit;
40 while (df > eps && inum < itermax)
42 MatrixXd W = MatrixXd::Zero(n,n);
49 int e =
CMath::ceil(log2(W.array().abs().rowwise().sum().maxCoeff()));
50 int s = std::max(0,e-1);
54 MatrixXd d = MatrixXd::Zero(EV.rows(),EV.cols());
55 d.diagonal() = VectorXd::Ones(EV.diagonalSize()).cwiseQuotient((EV * EV.transpose()).diagonal().cwiseSqrt());
58 for (
int i = 0; i < K; i++)
62 Ci = EV * C0i * EV.transpose();
66 for (
int i = 0; i < K; i++)
69 MatrixXd F = EV * C0i * EV.transpose();
70 f += (F.transpose() * F).diagonal().sum() - F.array().pow(2).matrix().diagonal().sum();
92 int auxij,auxji,auxii,auxjj;
96 for (
int i = 0; i < N; i++)
98 for (
int j = 0; j < N; j++)
105 for (
int i = 0; i < N; i++)
107 for (
int j = 0; j < N; j++)
109 for (
int k = 0; k < K; k++)
115 z[i][j] += C[auxii]*C[auxjj];
116 y[i][j] += 0.5*C[auxjj]*(C[auxij]+C[auxji]);
121 for (
int i = 0; i < N-1; i++)
123 for (
int j = i+1; j < N; j++)
127 W[auxij] = (z[j][i]*y[j][i] - z[i][i]*y[i][j])/(z[j][j]*z[i][i]-z[i][j]*z[i][j]);
128 W[auxji] = (z[i][j]*y[i][j] - z[j][j]*y[j][i])/(z[j][j]*z[i][i]-z[i][j]*z[i][j]);