9 namespace ublas = boost::numeric::ublas;
16 void FastICA(rmat &input, rmat &source_estimate, rmat &mixing_matrix)
18 mixing_matrix = ublas::identity_matrix<double>(input.size1());
19 rmat old_mix(mixing_matrix);
20 rmat new_mix(mixing_matrix);
21 rmat white_mat(mixing_matrix);
22 cmat evec(input.size1(), input.size1());
23 cvec eval(input.size1());
24 PCA(input, evec, eval);
25 white_mat = real(
WhiteMat(evec, eval));
26 std::cout <<
"white_mat: " << white_mat << std::endl;
27 input = prod(white_mat, input);
28 const int maxiterations = 1000;
29 const double min_improv = 0.01;
30 double improvement = 10.0;
32 while (niter < maxiterations)
35 ublas::noalias(source_estimate) = ublas::prod(mixing_matrix, input);
36 rmat pow3mat = prod(trans(input), mixing_matrix);
38 for (
size_t i = 0; i < pow3mat.size1(); ++i)
40 for (
size_t j = 0; j < pow3mat.size2(); ++j)
41 pow3mat(i, j) = std::pow(pow3mat(i, j), 3);
43 new_mix = prod(input, pow3mat);
44 new_mix /= input.size2();
45 new_mix -= 3.0 * mixing_matrix;
46 mixing_matrix = new_mix;
48 rmat wwt = prod(mixing_matrix, trans(mixing_matrix));
49 mixing_matrix /= sqrt(norm_inf(wwt));
50 for (
int i = 0; i < maxiterations; ++i)
52 rmat neww(mixing_matrix);
53 wwt = prod(mixing_matrix, trans(mixing_matrix));
55 mixing_matrix = prod(wwt, mixing_matrix);
57 neww -= mixing_matrix;
65 noalias(source_estimate) = prod(trans(mixing_matrix), input);
66 mixing_matrix = prod(trans(mixing_matrix), white_mat);
gplib::cmat WhiteMat(gplib::cmat &evectors, gplib::cvec &evalues)
Calculate the Whitening Matrix.
void FastICA(rmat &input, rmat &source_estimate, rmat &mixing_matrix)
void PCA(const UblasMatrix &observations, gplib::cmat &evectors, gplib::cvec &evalues)
This template function calculates the principal component rotation matrix from a matrix of observatio...