77 sparseoptimization pattern discrepancy

[dimalvovs]

• add comments to chi-squared calc
• add 1 to denominator and numerator to avoid NaN in chi-squared calc

[drbergman]

I can't comment on how the 1+'s will affect calculations and their accuracy for a chisq calculation. I'll only ask: isn't the purpose of the SparseIterator that it will skip all 0 values? So dsq = get<1>(it) * get<1>(it) must be > 0?

Also, though it's clearly outside the scope of this PR, it feels like there's the possibility for computational savings by storing the results of the dot products inside the for loop to be re-used in the while loop. Similarly, I think the dsq values are constant throughout the CoGAPS run (they're just the element-wise square of the genes x sample matrix, right?). Could a static variable be used to store these?

[dimalvovs]

I can't comment on how the 1+'s will affect calculations and their accuracy for a chisq calculation. I'll only ask: isn't the purpose of the SparseIterator that it will skip all 0 values? So dsq = get<1>(it) * get<1>(it) must be > 0?

that is exactly right, I think the root cause is the implementation of the sparseIterator, so maybe this draft PR should be more of a proof that that the NaN comes from the denominator and demonstrate that it is curable by +1ing.

[dimalvovs]

additional to-dos from a call with @favorov:

• change comments for clarity: mDMatrix is data (D); mMatrix is A (left mult matrix) if D is samples x genes; and P is D is genes x samples; *mOtherMatrix contans const pointer to vis-a-vis of mMatrix;

CoGAPS/src/gibbs_sampler/DenseNormalModel.h

Line 56 in b06f01b

Matrix mDMatrix; // samples by genes for A, genes by samples for P

• same for DenseNormalModel;
• are mMatrix and mOtherMatrix stored in the same dimensions (so mOther is transposed - rows are selected for both below)?
  

  CoGAPS/src/gibbs_sampler/SparseNormalModel.cpp

  Line 47 in b06f01b

  float dot = gaps::dot(mMatrix.getRow(j), mOtherMatrix->getRow(i)); // dot = A * P

• where is the uncertainty matrix?
  

  CoGAPS/src/gibbs_sampler/SparseNormalModel.cpp

  Line 55 in 4334b81

  chisq += 1 + dot * (dot - 2 * get<1>(it) - dsq * dot) / dsq;

• why we multiply chisq() by mBeta?
  

  CoGAPS/src/gibbs_sampler/SparseNormalModel.cpp

  Line 59 in 4334b81

  return chisq * mBeta;