Hi All,

I am kind of new to Eigen and i might have missed good number of opportunities where Eigen can optimize things (Using Expression templates, Lazy evaluation and vectorization etc) and reduce overall running time. I have pasted a piece of code with lots of computations (Using Eigen). I would be happy to see your comments/remarks regarding what i can remove/add in codes to reduce running time. Many thanks for your replies.

P.S: variables starting with m_ are matrix and variables starting with v_ are vector.

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;

bool ContinuousLBModel::EMRows()
{
  //Initialization
  m_Uil1_ = m_Dataij_*m_Rjl_;
  m_Uil2_ = m_Dataij2_*m_Rjl_;

  //Temporary variables
  MatrixReal m_sumik(nbSample_,nbRowCluster_),m_prodik(nbSample_,nbRowCluster_),m_trkl(nbRowCluster_,nbColCluster_);
  VectorReal v_sumikmax(nbSample_),v_sumi(nbSample_),Zsumk(nbRowCluster_);
  VectorReal::Index maxIndex;
  VectorReal Onesk = VectorReal::Ones(nbRowCluster_);
  VectorReal tmp1,tmp2;

  //EM algorithm
    for (int itr = 0; itr < IP::nbiterations_int_; ++itr) {
      //E-step
      tmp1 = 0.5*(((m_Sigma2kl_.array().log()+m_Mukl2_.array()/m_Sigma2kl_.array()).matrix())*v_Rl_);
      tmp2 = v_logPiek_;
      tmp2-=tmp1;
      m_sumik = VectorReal::Ones(nbSample_)*tmp2.transpose()
                -0.5*(m_Uil2_*((Array2DReal::Ones(nbRowCluster_,nbColCluster_)/m_Sigma2kl_.array()).matrix().transpose()))
                +(m_Uil1_*((m_Mukl_.array()/m_Sigma2kl_.array()).matrix().transpose()));

      Zsumk.setZero();
      for (int i = 0; i < nbSample_; ++i) {
        v_sumikmax(i) = m_sumik.row(i).maxCoeff(&maxIndex);
        Zsumk(maxIndex)+=1;
      }

      if((Zsumk.array()<0.00001).any())
      {
        empty_cluster = true;
        throw "Row clustering failed while running model.";
      }
      else
      {
        empty_cluster = false;
      }

      m_prodik=(m_sumik-v_sumikmax*Onesk.transpose()).array().exp();
      v_sumi = m_prodik.rowwise().sum();
      m_Tik_ = m_prodik.array()/(v_sumi*Onesk.transpose()).array();
      v_Tk_ = m_Tik_.colwise().sum();

      //M-step
      if(!IP::fixedproportions_) {
        v_logPiek_=(v_Tk_.array()/nbSample_).log();
      }

      m_Muklold2_ = m_Mukl_;
      m_trkl = (v_Tk_*v_Rl_.transpose());
      m_Mukl_ = (m_Tik_.transpose()*m_Uil1_).array()/m_trkl.array();
      m_Mukl2_ = m_Mukl_.array().square();
      m_Sigma2kl_ = (m_Tik_.transpose()*m_Uil2_).array()/m_trkl.array() - m_Mukl2_.array();

      //Termination check
      if((((m_Mukl_.array()-m_Muklold2_.array())/m_Mukl_.array()).abs().sum())<IP::epsilon_int_) {
        break;
      }

    }
}