Scala-Spark digamma stackoverflow问题

这两天在用spark做点击率的贝叶斯平滑，参考雅虎的论文进行了一番尝试。

先上代码：

 1 # click_count, show_count # this method takes time
 2 def do_smooth(data_list):
 3     import scipy.special as sp
 4     a, b, i = 1.0, 1.0, 0
 5     da, db = a, b
 6     while i < 1000 and (da > 1.0E-10 or db > 1.0E-10):
 7         x1, y1, x2 = 0.0, 0.0, 0.0
 8         for lineList in data_list:
 9             x1 += sp.digamma((lineList[0]) + a) - sp.digamma(a)
10             y1 += sp.digamma((lineList[1]) + a + b) - sp.digamma(a + b)
11             x2 += sp.digamma((lineList[1]) - (lineList[0]) + b) - sp.digamma(b)
12         na, nb = a, b
13         a *= (x1 / y1)
14         b *= (x2 / y1)
15         da, db = abs(a - na), abs(b - nb)
16         i += 1
17     print i, a, b
18     return a, b

 

这是我之前用的python代码，改成scala也相当容易，digamma函数非常耗时，而且还要迭代1000次。最要命的是digamma在scala里面默认的实现会出现栈溢出！！！

var a, b, da, db: Double = 1.0
var index = 0
while (index < 1000 && (da > 1.0E-9 || db > 1.0E-9)) {
    var x1,x2,y1 = 0.0
    traindata.foreach(p => {
        x1 += MBlas.digamma(p(2) + a) - MBlas.digamma(a)
        y1 += MBlas.digamma(p(1) + a + b) - MBlas.digamma(a + b)
        x2 += MBlas.digamma(p(1) - p(2) + b) - MBlas.digamma(b)
        val na = a
        val nb = b
        a *= (x1 / y1)
        b *= (x2 / y1)
        da = Math.abs(a - na)
        db = Math.abs(b - nb)
    })
}

 

digamma 函数是个递归函数，问题就处在递归上了。

 1    public static double digamma(double x) {
 2         if (x > 0 && x <= S_LIMIT) {
 3             return -GAMMA - 1 / x;
 4         }
 5         if (x >= C_LIMIT) {
 6             double inv = 1 / (x * x);
 7             return FastMath.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252));
 8         }
 9         return digamma(x + 1) - 1 / x;
10     }

 

既然知道问题所在，是不是就可以重写递归为非递归呢？在Stack Overflow上找到了一个答案

 1 val GAMMA = 0.577215664901532860606512090082
 2 val GAMMA_MINX = 1.e-12
 3 val DIGAMMA_MINNEGX = -1250
 4 val C_LIMIT = 49
 5 val S_LIMIT = 1e-5
 6 var value = 0.0
 7 var x = input
 8 while (true) {
 9     if (x >= 0 && x < GAMMA_MINX) x = GAMMA_MINX
10     if (x < DIGAMMA_MINNEGX) {
11         x = DIGAMMA_MINNEGX + GAMMA_MINX
12     } else {
13         if (x > 0 && x <= S_LIMIT) return value + -GAMMA - 1 / x
14         if (x >= C_LIMIT) {
15             val inv = 1 / (x * x)
16             return value + Math.log(x) - 0.5 / x - inv * ((1.0 / 12) + inv * (1.0 / 120 - inv / 252))
17         }
18         value = value - 1.0 / x
19         x += 1
20     }
21 }

 

经测试，没看出什么问题，可以用了。
不过，上面的代码并没有解决慢的问题，当需要计算CTR的对象比较多时（几百万），仍然比较耗时。所以我决定用两个替代方法：

  1. 抽样，抽取能在可接受时间内出结果的样本数，得到α和β；
  2. 直接使用平均值作为α和β
  3. 使用平均值做迭代初值（推荐）

参考：
1. 雅虎专家的论文，如上
2. Stack Overflow 网友代码，如上

 

posted @ 2017-09-24 20:01  longwind09  阅读(624)  评论(4)  编辑  收藏  举报