0

I have two large numpy matrices, X and r (appr. 5500x3500). I need to calculate digamma(X+r) and digamma(r), which I do with scipy.special.psi. The calculation happens in the gradient function that I am using for a fit by gradient descent. During optimization the gradient function is going to be called a lot of times, so it is imperative that it runs fast. Now, to the actual question:

  • digamma(r) takes ~11s to run
  • digamma(X+r) takes ~50s to run

out of curiosity I ran some other examples:

  • digamma(X.astype(np.float64)), takes ~1min to run
  • digamma(X), takes ~5-6min to run

why is that happening, and what can I do to improve runtimes?

  • X contains integers (numpy.int64) from 0 to 164.
  • r contains floats (numpy.float64) from 0. to ~4.6.
  • X is very sparse (about 75% of it is 0).
  • r is not sparse (about 10% of it is 0)

EDIT: corrected X.as_float

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
galicae
  • 101
  • 1
  • numpy arrays and matrices do not have an `as_float` attribute. Do you mean `X.astype(np.float64)`? – Warren Weckesser Jun 26 '17 at 14:40
  • yeah, I just didn't want to write it like that for the sake of brevity. Thanks for pointing it out, might as well edit. – galicae Jun 26 '17 at 16:33
  • So `digamma(X)` alone runs slower than `digamma(X+r)`? That's quite curious. – cadolphs Jun 26 '17 at 16:36
  • The docs for `digamma` hints that computation complexity depends on the variable value. In quick time tests `digamma(np.ones(1000))` is significantly slower than anything else. The wiki page for this function also shows considerable computational complexity. – hpaulj Jun 26 '17 at 17:38
  • it is quite puzzling... When I simulated matrices of similar sparsity to `X` and `r` with `scipy.sparse` and filled them with values from the corresponding ranges, they both needed close to 1 minute to run. vOv A glance at the source code of digamma confirms that it has different approaches for different value ranges. My working hypothesis is that the makeup of `X` and `r` are uniquely suited (or not) for digamma calculations. I just have to swallow the computational cost for now :( – galicae Jun 27 '17 at 08:03
  • I've run into the same problem and opened an [issue](https://github.com/scipy/scipy/issues/8127) about it. In the meantime, the scipy implementation of digamma is probably best avoided. For an alternative, [here](http://people.sc.fsu.edu/~jburkardt/py_src/asa103/asa103.html)'s a good place to start. – rmalouf Nov 06 '17 at 19:42
  • thanks a lot! Will definitely check it out! – galicae Nov 08 '17 at 10:31
  • Looks like it'll be improved in the next release: https://github.com/scipy/scipy/pull/8129 – rmalouf Nov 20 '17 at 22:13
  • yeah, I was following the thread on github. Glad to see that it is getting improved! If you or @hpaulj want to post this as an answer I can mark the question as closed :) – galicae Nov 22 '17 at 11:26

0 Answers0