24 using namespace Eigen;
43 return m_cumulant_matrix;
60 VectorXd mean = (EX.rowwise().sum() / (
float64_t)T);
66 std::cout <<
"cov" << std::endl;
67 std::cout << cov << std::endl;
71 SelfAdjointEigenSolver<MatrixXd> eig;
75 std::cout <<
"eigenvectors" << std::endl;
76 std::cout << eig.eigenvectors() << std::endl;
78 std::cout <<
"eigenvalues" << std::endl;
79 std::cout << eig.eigenvalues().asDiagonal() << std::endl;
83 VectorXd scales = eig.eigenvalues().cwiseSqrt();
84 MatrixXd B = scales.cwiseInverse().asDiagonal() * eig.eigenvectors().transpose();
87 std::cout <<
"whitener" << std::endl;
88 std::cout << B << std::endl;
95 int dimsymm = (m * ( m + 1)) / 2;
101 VectorXd Xim = VectorXd::Zero(m);
102 VectorXd Xjm = VectorXd::Zero(m);
103 VectorXd Xijm = VectorXd::Zero(m);
106 for (
int im = 0; im < m; im++)
109 Xijm = Xim.cwiseProduct(Xim);
110 Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R - 2*R.col(im)*R.col(im).transpose();
111 CM.block(0,Range,m,m) = Qij;
113 for (
int jm = 0; jm < im; jm++)
116 Xijm = Xim.cwiseProduct(Xjm);
117 Qij = SPX.cwiseProduct(Xijm.replicate(1,m).transpose()) * SPX.transpose() / (float)T - R.col(im)*R.col(jm).transpose() - R.col(jm)*R.col(im).transpose();
118 CM.block(0,Range,m,m) = sqrt(2)*Qij;
124 std::cout <<
"cumulatant matrices" << std::endl;
125 std::cout << CM << std::endl;
135 for (
int i = 0; i < nbcm; i++)
138 EM = CM.block(0,i*m,m,m);
144 EQ = -1 * EQ.inverse();
147 std::cout <<
"diagonalizer" << std::endl;
148 std::cout << EQ << std::endl;
154 C = EQ.transpose() * B;
157 VectorXd A = (B.inverse()*EQ).cwiseAbs2().colwise().sum();
162 for (
int j = 1; j < n; j++)
166 std::swap(A(j),A(j-1));
167 C.col(j).swap(C.col(j-1));
174 for (
int j = 0; j < m/2; j++)
175 C.row(j).swap( C.row(m-1-j) );
178 VectorXd signs = VectorXd::Zero(m);
179 for (
int i = 0; i < m; i++)
186 C = signs.asDiagonal() * C;
189 std::cout <<
"un mixing matrix" << std::endl;
190 std::cout << C << std::endl;
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
virtual CFeatures * apply(CFeatures *features)
SGMatrix< float64_t > get_cumulant_matrix() const
all of classes and functions are contained in the shogun namespace
class ICAConverter Base class for ICA algorithms
The class Features is the base class of all feature objects.
SGMatrix< float64_t > m_mixing_matrix
static void inverse(SGMatrix< float64_t > matrix)
inverses square matrix in-place