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];
39 std::vector<float64_t> crit;
40 while (df > eps && inum < itermax)
49 int e =
CMath::ceil(log2(W.array().abs().rowwise().sum().maxCoeff()));
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]);
static float64_t ceil(float64_t d)
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)
T * get_matrix(index_t matIdx) const
void getW(float64_t *C, int *ptN, int *ptK, float64_t *W)
all of classes and functions are contained in the shogun namespace
Matrix::Scalar max(Matrix m)
static SGMatrix< T > create_identity_matrix(index_t size, T scale)