Method fast_ica_parallel_compute

src/data_handler/fastica.cpp:53–98

}

Eigen::MatrixXd FastICA::fast_ica_parallel_compute (const Eigen::MatrixXd &X)
{
    int cols = (int)X.cols ();
    Eigen::MatrixXd W (num_components, num_components);
    random_normal (W);

    //  sW <- La.svd(W)
    Eigen::BDCSVD<Eigen::MatrixXd> sW (W, Eigen::ComputeThinU | Eigen::ComputeThinV);
    // W <- sW$u %*% Diag(1/sW$d) %*% t(sW$u) %*% W
    W = sW.matrixU () * (sW.singularValues ().array ().inverse ()).matrix ().asDiagonal () *
        sW.matrixU ().transpose () * W;
    Eigen::MatrixXd W1 = W;
    // lim <- rep(1000, maxit)
    std::vector<double> lim (max_it, 1000);
    // iteration counter
    int it = 0;

    while (lim[it] > tol && it < (max_it - 1))
    {
        //  wx <- W %*% X
        // gwx <- tanh(alpha * wx)
        // alpha = 1 , so ignore
        Eigen::MatrixXd gwx = (W * X).array ().tanh ().matrix ();
        // v1 <- gwx %*% t(X)/cols
        Eigen::MatrixXd v1 = gwx * (X.array () / cols).matrix ().transpose ();
        // g.wx <- alpha * (1 - (gwx)^2)
        // nb alpha == 1
        gwx = 1 - gwx.array ().square ();
        // v2 <- Diag(apply(g.wx, 1, FUN = mean)) %*% W
        Eigen::MatrixXd v2 = gwx.array ().rowwise ().mean ().matrix ().asDiagonal () * W;
        // W1 <- v1 - v2
        W1 = v1 - v2;
        // sW1 <- La.svd(W1)
        Eigen::BDCSVD<Eigen::MatrixXd> sW1 (W1, Eigen::ComputeThinU | Eigen::ComputeThinV);
        // W1 <- sW1$u %*% Diag(1/sW1$d) %*% t(sW1$u) %*% W1
        W1 = sW1.matrixU () * (sW1.singularValues ().array ().inverse ()).matrix ().asDiagonal () *
            sW1.matrixU ().transpose () * W1;
        // lim[it + 1] <- max( Mod(   Mod(  diag(W1 %*% t(W) )  )  - 1 ) )
        lim[it + 1] = ((W1 * W.transpose ()).diagonal ().array ().abs () - 1).abs ().maxCoeff ();
        // W <- W1
        W = W1;
        ++it;
    }

    return W;
}

void FastICA::scale (Eigen::Ref<Eigen::MatrixXd> M, bool center, bool normalize,
    bool ignore_invariants, std::vector<int> *zeros)