Hi All,

After profiling my codes i figured out that the function GaussianLogsumRows (shown below in code) is being called several times (inside EM algorithm) and is where most of the time is being consumed.

After reading Eigen docs, i find that writing this way could give Eigen more opportunities for optimization using Expression templates and Lazy evaluation. My question is if there is yet better way of writing such expressions which could perhaps give Eigen even more opportunities for optimization.

Any advice is much appreciated !

Code: Select all

//Matrix containers
typedef Eigen::MatrixXf MatrixReal;
typedef Eigen::MatrixXi MatrixInteger;
typedef Eigen::Matrix<bool,Eigen::Dynamic,Eigen::Dynamic> MatrixBinary;
//Vector Containers
typedef Eigen::VectorXf VectorReal;
typedef Eigen::VectorXi VectorInteger;
typedef Eigen::Matrix<bool,Eigen::Dynamic,1> VectorBinary;
//2D array containers
typedef Eigen::ArrayXXf Array2DReal;

void GaussianLogsumRows(MatrixReal & m_Sumik)
{
  m_Sumik = MatrixReal::Ones(nbSample_,1)*v_logPiek_.transpose()
            - 0.5*(MatrixReal::Ones(nbSample_,1)*v_Rl_.transpose()*((m_Sigma2kl_.array().log()+m_Mukl2_.array()/m_Sigma2kl_.array()).matrix()).transpose()
            + m_Uil2_*((Array2DReal::Ones(nbRowCluster_,nbColCluster_)/m_Sigma2kl_.array()).matrix().transpose())
            - 2*(m_Uil1_*((m_Mukl_.array()/m_Sigma2kl_.array()).matrix().transpose())));
}


Sincerely,
Parmeet

Something like this should be faster:

m_Sumik = (v_logPiek_-0.5*v_Rl_.transpose()*((m_Sigma2kl_.array().log() + m_Mukl2_.array()/m_Sigma2kl_.array()).matrix())).eval().transpose().replicate(nbSample_,1)
- (0.5*m_Uil2_)*(m_Sigma2kl_.array().inverse()).matrix().transpose()
+ m_Uil1_*((m_Mukl_.array()/m_Sigma2kl_.array()).matrix().transpose());


Also makes sure to use Vector* type whenever that's possible instead of a Matrix* type with 1 column or 1 row. Knowing the types of v_Rl_, m_Mukl2_, m_Uil1_, etc. would help to further refine the expr. You can also decompose this expression with a few temporaries that are created anyway:

tmp1 = m_Sigma2kl_.array().log() + m_Mukl2_.array()/m_Sigma2kl_.array();
tmp2 = v_logPiek_;
tmp2 -= 0.5*v_Rl_.transpose()*tmp1;
m_Sumik = tmp2.transpose().colwise().replicate(nbSample_)
- (0.5*m_Uil2_)*(m_Sigma2kl_.array().inverse()).matrix().transpose()
+ m_Uil1_*((m_Mukl_.array()/m_Sigma2kl_.array()).matrix().transpose());

Thanks for your answer. Yes , i am able to gain good speed by making my expression removing unnecessary/repeated operations.

To my surprise, "replicate" is performing worst than multiplying vector by appropriate matrix/vector that is to say v_sumikmax*VectorReal::Ones(nbRowCluster_).transpose() is much better than v_sumikmax.replicate(1,nbRowCluster_)

where v_sumikmax is a vector.

Any ideas?

very strange.... is it nested into a more complex expression??

The expressions are shown below in the code.The second expression is pretty costly in terms of computation cost.

m_prodik=(m_sumik-v_sumikmax*VectorReal::Ones(nbRowCluster_).transpose()).array().exp();

m_prodik=(m_sumik-v_sumikmax.replicate(1,nbRowCluster_)).array().exp();

indeed, the fix is to use:

v_sumikmax.rowwise().replicate(nbRowCluster_) then it becomes faster than the product. This version avoids a very expensive modulo when addressing the elements.

Ye, its better than previous version but still performing worst than product version..

Oh indeed, the replicate expression is not vectorized yet. This has to be fixed.