[Question] in the implementation of function hpbcpp

[ping-at-git]

In STMCfuncs.cpp, function hpbcpp, line 128, the following:

EB.each_row() %= arma::trans(sqrt(doc_cts))/sum(EB,0);

If i'm not mistaken, the above implementation is modifying the EB matrix in place, which also means sum(EB,0) should have a different value for each iteration. So my question is: is this the intended behavior? or was the intention to have the denominator as the constant version of the column sum of EB without any modification?

[bstewart]

I've checked this implementation against several R implementations (as I know R much better than C++) so I'd be really surprised if there is an error here. These other implementations were slower but much clearer because I didn't have to worry as much about speed.

When you say sum(EB,0) should have a different value 'for each iteration' you mean each row right? I think what is happening is it is precomputing the entire right hand side (for both speed and accuracy reasons) and then doing the element-wise multiplication. That said, I see what you mean and the documentation in armadillo isn't super clear on this point. Have you run tests suggesting this is wrong or just having a look?

[ping-at-git]

I was indeed first just having a look, and I did not know exactly what is the expected behaviour thus left the comment. It is always tricky when the iteratively modifying the object, as suggested, I just run a quick test:

the following function is implemented:
arma::mat testEB(arma::rowvec input_vec){ arma::mat mat(input_vec.n_elem, input_vec.n_elem, arma::fill::ones); mat.each_row() %= input_vec/sum(mat,0); return mat; }
and the following test is done:
`

testEB(c(1,1,1,1,1))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.2 0.2 0.2 0.2 0.2
[2,] 0.2 0.2 0.2 0.2 0.2
[3,] 0.2 0.2 0.2 0.2 0.2
[4,] 0.2 0.2 0.2 0.2 0.2
[5,] 0.2 0.2 0.2 0.2 0.2
`
I hope this is the expected behaviour, in which case there is no problem at all.

[bstewart]

Yeah just to clarify- expected behavior here is that the vector is precomputed and the multiplied elementwise by row. In other words, my intention was that each row is multiplied by the same vector (not that it would be updated line by line as one would get if one looped over the rows recomputing the vector each time). It looks like it is doing the right thing though- which is great.

If you see an equivalently fast and memory efficient way to do this but which is clearer, please let me know.