I was looking to get some feedback on this ~20 line MLE loop. Since moving to the new Eigen 3 trunk, I had to switch from cwise to array and I didn't know if there was a better way to do this.

Code: Select all

bool solve(Eigen::VectorXf& x, Eigen::MatrixXf* A, Eigen::VectorXf* m)
{
  x = Eigen::VectorXf::Ones(num_cols);
  Eigen::VectorXf column_sums(num_cols);
  column_sums = A->colwise().sum();
 
  if ( (m->array() < 0).any() )
  {
     throw negative_measurement("negative measurements encountered.");
  }
  // estimated measurements ( = A*x ) during iterations
  Eigen::VectorXf curr_estimate(num_rows);
 
  // main iteration loop
  for (int iter = 0; iter < m_maxIter; ++iter)
  {
     if (!m_abort)
     {
        // calculate current reconstruction estimate
        curr_estimate = (*A) * x;
        // adjust measurements by it...
        for (int row = 0; row < num_rows; ++row)
        {
           if (curr_estimate(row) > 0 )
           {
              curr_estimate(row) = (*m)(row) / curr_estimate(row);
           }
        }
        // ...and reproject, normalize to update solution
        x = x.array() * ( ( A->transpose() * curr_estimate ).array() / column_sums.array() ); // used to be cwise
        // update iteration count
        if(m_ei)
        {
           m_ei->iterationsComplete(iter+1);
        }
     }
  }
}

I'm only running 20-50 iterations max, on 1-12GB "A" matricies

Thanks for your time.
JM

It is already pretty good. Here is a slightly optimized version:

Code: Select all

x = Eigen::VectorXf::Ones(num_cols);
Eigen::ArrayXf temp(num_cols);
Eigen::ArrayXf column_sums(num_cols);
column_sums = A->colwise().sum();
 
/* ... */

curr_estimate.noalias() = (*A) * x; // saves one useless temporary/copy per iteration

curr_estimate = (curr_estimate.array()>0).select(m->cwiseQuotient(curr_estimate),curr_estimate);
// ^^ will be vectorized for you in the future

temp.noalias() = A->transpose() * curr_estimate;
// ^^ the temporary object needed for the product will be allocated only once

x.array() *= temp / column_sums;

/* ... */


for my curiosity could you tell us what does this solver solve?

it solves for activity distributions


Thanks a ton. I need to keep improving my linear algebra skills. I feel like with a good mastery of linear algebra you can solve anything.

EDIT:
My compiler is telling me that noalias() is valid on Eigen::VectorXf's but not on Eigen::ArrayXf's

2> error C2039: 'noalias' : is not a member of 'Eigen::Array<_Scalar,_Rows,_Cols>'
2> with
2> [
2> _Scalar=float,
2> _Rows=-1,
2> _Cols=1
2> ]

Namely, this line
temp.noalias() = A->transpose() * curr_estimate;
// ^^ the temporary object needed for the product will be allocated only once

Should temp be a VectorXf and not an ArrayXf? I'm using Eigen trunk latest, I just pull'd 5 minutes ago.

I think temp should be a vector and then it should be
temp.noalias() = A->transpose() * curr_estimate;
x.array() *= temp.array() / column_sums;

Thoughts?

I think temp should be a vector and then it should be
temp.noalias() = A->transpose() * curr_estimate;
x.array() *= temp.array() / column_sums;

yes that's probably simpler.