Hi Toby, > So I spent some time today investigating a 'bug'. The prototype for nmf is: > > template <typename ScalarType> > void nmf(viennacl::matrix<ScalarType> const & V, > viennacl::matrix<ScalarType> & W, > viennacl::matrix<ScalarType> & H, > nmf_config const & conf) > > Now, if W or H is a zero matrix than prod(W, H) is also a zero > matrix. So prod(W, H) - V = -V. And ||V|| = ||-V||, so the difference > between the initial norm and the norm after the 1st iteration is 0, > which is always less than conf.tolerance(). So the routine quits. > > I think this is a bug in the algorithm: why should the behaviour depend > on the initial state of W and H, unless the initial W and H are supposed > explicitly to be taken as approximations? But if so, I don't think a > zero matrix approximation should cause the algorithm to quit!
hmm, on the one hand users want to have the option of providing a 'guess' for W and H, but on the other hand they should be warned if their guess is bad. > The following change is sufficient to fix this. I've committed it and a > change ensuring any ViennaCL objects created during the run are in the > right context. I don't know if it is ideal, since it means that in the > (rare) case where W and H in their initial state are adequate, two sets > of check_after_steps iterations are computed, where previously only one > set would have been. There are alternatives: of course, we could compute > norm_frobenius(prod(W, H)) and check that it is not 0 on the first > iteration, but that's much more expensive in every other case. A single matrix-matrix product shouldn't be too bad here, particularly as each iteration in nmf() takes several matrix-matrix products. Similar to the last fix you proposed in the follow-up email, what about checking norm_frobenius(W) and norm_frobenius(H) to be positive and to reinitialize them with all ones if not? This would fix the generic use case: matrix<T> W(N, K); matrix<T> H(K, M); nmf(A, W, H, conf); Best regards, Karli ------------------------------------------------------------------------------ Want fast and easy access to all the code in your enterprise? Index and search up to 200,000 lines of code with a free copy of Black Duck Code Sight - the same software that powers the world's largest code search on Ohloh, the Black Duck Open Hub! Try it now. http://p.sf.net/sfu/bds _______________________________________________ ViennaCL-devel mailing list ViennaCL-devel@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/viennacl-devel
