https://en.wikipedia.org/wiki/FastICA https://arnauddelorme.com/ica_for_dummies/
| 6 | // https://en.wikipedia.org/wiki/FastICA |
| 7 | // https://arnauddelorme.com/ica_for_dummies/ |
| 8 | int FastICA::compute (Eigen::MatrixXd &X) |
| 9 | { |
| 10 | int rows = (int)X.rows (); |
| 11 | int cols = (int)X.cols (); |
| 12 | int min_rows_cols = rows < cols ? rows : cols; |
| 13 | if ((num_components < 2) || (max_it < 1) || (rows < 2) || (cols < 2) || |
| 14 | (num_components > min_rows_cols)) |
| 15 | { |
| 16 | return (int)BrainFlowExitCodes::INVALID_ARGUMENTS_ERROR; |
| 17 | } |
| 18 | |
| 19 | W.resize (num_components, num_components); |
| 20 | A.resize (rows, num_components); |
| 21 | K.resize (num_components, rows); |
| 22 | S.resize (num_components, cols); |
| 23 | |
| 24 | scale (X, true, row_norm); |
| 25 | |
| 26 | // Whitening |
| 27 | // X %*% t(X)/rows |
| 28 | Eigen::MatrixXd V = X * (X.array () / cols).matrix ().transpose (); |
| 29 | // s <- La.svd(V) |
| 30 | Eigen::BDCSVD<Eigen::MatrixXd> s (V, Eigen::ComputeThinU | Eigen::ComputeThinV); |
| 31 | // D <- diag(c(1/sqrt(s$d))) |
| 32 | Eigen::MatrixXd D = s.singularValues ().array ().sqrt ().inverse ().matrix ().asDiagonal (); |
| 33 | // K <- D %*% t(s$u) |
| 34 | Eigen::MatrixXd K_temp = D * s.matrixU ().transpose (); |
| 35 | // K <- matrix( K[1:rows.comp, ], rows.comp, cols) |
| 36 | Eigen::MatrixXd K_temp2 = K_temp.block (0, 0, num_components, rows); |
| 37 | // X1 <- K %*% X |
| 38 | Eigen::MatrixXd X1 = K_temp2 * X; |
| 39 | Eigen::MatrixXd a = fast_ica_parallel_compute (X1); |
| 40 | // w <- a %*% K |
| 41 | Eigen::MatrixXd w = a * K_temp2; |
| 42 | // S <- w %*% X |
| 43 | S = w * X; |
| 44 | // A <- t(w) %*% solve(w %*% t(w)) |
| 45 | A = w.transpose () * (w * w.transpose ()).inverse (); |
| 46 | A.transposeInPlace (); |
| 47 | K = K_temp2; |
| 48 | W = a; |
| 49 | |
| 50 | return (int)BrainFlowExitCodes::STATUS_OK; |
| 51 | } |
| 52 | |
| 53 | Eigen::MatrixXd FastICA::fast_ica_parallel_compute (const Eigen::MatrixXd &X) |
| 54 | { |