LR和FTRL-Proximal

看到别人使用过FM+FTRL的模型实现了一个CTR算法,印象很深。自己使用的是DNN+Embedding的方式做的一个算法模型,从理论上看embedding肯定比one-hot encoder的方式更加先进且能真实反馈特征数据的相关性,但是实际效果看对方的FM_FRTL得到的AUC比我高近一个百分点,而且可以在10G的数据上一个多小时跑完,而我的DNN+Embedding算法,因为没有GPU主机,跑一次需要十二个小时,严重影响了调参的积极性,这也是我非常想掌握FTRL的出发点。

更重要的是,考虑到刷竞赛与实际算法是否可工程化的的角度,FTRL结合LR或者FM是一个非常好的方向,比Top1开源的代码使用了PCA、NLP处理特征以及多个模型stacking的技巧,更具有学习或者借鉴的价值。

本文主要根据谷歌给出的FTRL理论论文,以及FTRL+LR的工程化实现论文,从理论到工程化实现LR+FTRL的开发,任一后端开发人员都能根据文末给出的python代码,简单的开发就能实现一个简单、高性能、高可靠的CTR预测模型。

LR

关键概念点:

1.logistic distribution
2.几率
3.对数几率
4.损失函数
5.参数估计
6.误差计算
7.随机梯度下降

可以看出LR也可以认为是用线性回归模型的预测结果来预测事件发生的对数几率。

LR有很多优点,比如:

  • 作为统计模型与机器学习的结合点,具有较好的预测结果以及可解释性。
  • 直接对分类的可能性建模,无需事先假设数据的分布,这就避免了假设分布不准确带来的影响。
  • 不仅预测得到分类,还有分类对应的概率,这对很多需要使用概率辅助决策的任务很有用。
  • sigmoid函数是高阶可导的凸函数,具有很好的数学性质,很多数值优化的算法都可以直接用于求解最优解。

参数估计

当样本数据里N很大的时候,通常采用的是随机梯度下降法,算法如下所示:

while {
    for i in range(0,m):
        w_j = w_j + a * g_j
}

随机梯度下降的好处是可以实现分布式并行化,具体计算流程是:

  1. 在每次迭代的时候,随机抽样一定比例的样本作为当前迭代的计算样本。
  2. 对计算样本中的每一个样本,分别计算不同特征的计算梯度。
  3. 通过聚合函数,对所有计算样本的特征的梯度进行累加,得到每一个特征的累积梯度以及损失。
  4. 最后根据最新的梯度以及之前的参数,对参数进行更新。
  5. 根据更新的参数计算损失函数误差值,如果损失函数误差值达到允许的范围,那么停止迭代,否则重复步骤1。

工程化实现思路(FTRL-Proximal)

主要是需要实现参数更新计算以及损失函数计算。
LR+FTRL工程化实现

FTRL-Proximal全称Followthe-Regularized-Leader Proximal,是谷歌公司提出的在线学习算法,在处理带非光滑正则项(例如l_1范数)的凸优化问题上表现出色。传统的基于batch的算法无法有效地处理大规模的数据和在线数据,而许多互联网应用,例如广告,数据是一条一条过来的,每过来一条样本数据,模型的参数需要根据这个样本进行迭代更新。面对这样的应用场景,谷歌提出了FTRL-Proximal算法,并给出了工程化实现。在线学习算法根据每一个样本更新模型参数时,由于梯度方向不是全局的,会存在误差,FTRL算法很好地解决了这个问题,在保证模型精度的同时还获得了更好的稀疏性(减轻了线上预测时的内存消耗和计算压力)。

原理

FTRL算法的梯度更新方式如式(1)所示,
\begin{equation}W^{t+1} = \mathop{\arg\min}_{W} \{G^{(1:t)}W+\lambda_1||W||_1+\frac{1}{2}\lambda_2||W||_2^2+\frac{1}{2}\sum_{s=1}^{t}\sigma^{(1:t)}||W-W^{(s)}||_2^2\}\tag{1}\end{equation}
其中,G^{(1:t)}=\sum_{s=1}^{t}G^{(s)}\sigma^{(s)}=\frac{1}{\eta^(s)}-\frac{1}{\eta^{(s-1)}}\sigma^{(1:t)}=\sum^{t}_{s=1}\sigma^{(s)}=\frac{1}{\eta^{(t)}},对式(1)展开,
\begin{equation}W^{t+1} = \mathop{\arg\min}_{W} \{(G^{(1:t)}-\sum_{t=1}^s\sigma^{(s)}W^{(s)})W+\lambda_1||W||_1+\frac{1}{2}(\lambda_2+\sum_{t=1}^s\sigma^{(s)})||W||_2^2+\frac{1}{2}\sum_{s=1}^{t}\sigma^{(s)}||W^{(s)}||_2^2\}\tag{2}\end{equation}
式(2)中的\frac{1}{2}\sum_{s=1}^{t}\sigma^{(s)}||W^{(s)}||_2^2可以看作是常数,令Z^{(t)}=G^{(1:t)}-\sum_{t=1}^s\sigma^{(s)}W^{(s)},式(2)写为
\begin{equation}W^{t+1} = \mathop{\arg\min}_{W} \{Z^{(t)}W+\lambda_1||W||_1+\frac{1}{2}(\lambda_2+\sum_{t=1}^s\sigma^{(s)})||W||_2^2\}\tag{3}\end{equation}
拆分到每一个维度,
\begin{equation}w^{t+1}_i = \mathop{\arg\min}_{w_i} \{z^{(t)}_i w_i+\lambda_1|w_i|_1+\frac{1}{2}(\lambda_2+\sum_{t=1}^s\sigma^{(s)})w_i^2\}\tag{4}\end{equation}
通过计算得到,
\begin{equation}w_i^{(t+1)}=\begin{cases}0, if \quad |z_i^{(t)}|<\lambda_1\\-(\lambda_2+\sum_{s=1}^{t}\sigma(s))(z_i^{(t)}-\lambda_1*sign(z_i^{(t)})),otherwise\end{cases}\tag{5}\end{equation}
每一个维度的学习率都是单独考虑的,
\begin{equation}\eta_{i}^{(t)}=\frac{\alpha}{\beta+\sqrt{\sum_{s=1}^t (g_i^{(s)})^2}}\tag{6}\end{equation}
其中g_i^{(s)}W的第i个维度w_i在时间步s的梯度。

FTRL的工业化实现伪代码如下,


image.png

根据前面的学习,可以使用LR作为基本学习器,使用FTRL作为在线最优化的方法来获取LR的权重系数,从而达到在不损失精度的前提下获得稀疏解的目标。

工程化实现的几个核心点是:

梯度计算代码
损失函数计算代码
权重更新计算代码

算法代码实现如下所示,这里只给出了核心部分代码,主要做了以下优化:

  • 修复原代码更新梯度时候的逻辑错误。
  • 支持pypy3加速,5.6G的训练数据,使用Mac单机可以9分钟跑完一个模型,一个epoch之后的logloss结果为0.3916;这对于刷竞赛或者实际工程中模型训练已经是比较理想的性能了。
  • 使用16位浮点数保存权重数据,降低模型文件大小。
  • 使用json保存模型非0参数w,z,n,进一步压缩线上使用的时候占据的内存空间,可以进一步考虑压缩模型文件大小,比如使用bson编码+gzip压缩后保存等。
from datetime import datetime
from csv import DictReader
from math import exp, log, sqrt
import gzip
import random
import json
import argparse


class FTRLProximal(object):
    """
    FTRL Proximal engineer project with logistic regression
    Reference:
    https://static.googleusercontent.com/media/research.google.com/zh-CN//pubs/archive/41159.pdf

    """

    def __init__(self, alpha, beta, L1, L2, D,
                 interaction=False, dropout=1.0,
                 dayfeature=True,
                 device_counters=False):

        # parameters
        self.alpha = alpha
        self.beta = beta
        self.L1 = L1
        self.L2 = L2
        self.dayfeature = dayfeature
        self.device_counters = device_counters

        # feature related parameters
        self.D = D
        self.interaction = interaction
        self.dropout = dropout

        # model
        self.n = [0.] * D
        self.z = [0.] * D
        self.w = [0.] * D

    def _indices(self, x):
        '''
        A helper generator that yields the indices in x
        The purpose of this generator is to make the following
        code a bit cleaner when doing feature interaction.
        '''

        for i in x:
            yield i

        if self.interaction:
            D = self.D
            L = len(x)
            for i in range(1, L):  # skip bias term, so we start at 1
                for j in range(i + 1, L):
                    # one-hot encode interactions with hash trick
                    yield abs(hash(str(x[i]) + '_' + str(x[j]))) % D

    def predict(self, x, dropped=None):
        """
        use x and computed weight to get predict
        :param x:
        :param dropped:
        :return:
        """
        # wTx is the inner product of w and x
        wTx = 0.
        for j, i in enumerate(self._indices(x)):

            if dropped is not None and dropped[j]:
                continue

            wTx += self.w[i]

        if dropped is not None:
            wTx /= self.dropout

        # bounded sigmoid function, this is the probability estimation
        return 1. / (1. + exp(-max(min(wTx, 35.), -35.)))

    def update(self, x, y):
        """
        update weight and coordinate learning rate based on x and y
        :param x:
        :param y:
        :return:
        """

        ind = [i for i in self._indices(x)]

        if self.dropout == 1:
            dropped = None
        else:
            dropped = [random.random() > self.dropout for i in range(0, len(ind))]

        p = self.predict(x, dropped)

        # gradient under logloss
        g = p - y

        # update z and n
        for j, i in enumerate(ind):

            # implement dropout as overfitting prevention
            if dropped is not None and dropped[j]:
                continue

            g_i = g * i
            sigma = (sqrt(self.n[i] + g_i * g_i) - sqrt(self.n[i])) / self.alpha
            self.z[i] += g_i - sigma * self.w[i]
            self.n[i] += g_i * g_i

            sign = -1. if self.z[i] < 0 else 1.  # get sign of z[i]

            # build w on the fly using z and n, hence the name - lazy weights -
            if sign * self.z[i] <= self.L1:
                # w[i] vanishes due to L1 regularization
                self.w[i] = 0.
            else:
                # apply prediction time L1, L2 regularization to z and get
                self.w[i] = (sign * self.L1 - self.z[i]) \
                            / ((self.beta + sqrt(self.n[i])) / self.alpha + self.L2)

    def save_model(self, save_file):
        """
        保存weight数据到本地
        :param save_file:
        :return:
        """
        with open(save_file, "w") as f:
            w = {k: v for k, v in enumerate(self.w) if v != 0}
            z = {k: v for k, v in enumerate(self.z) if v != 0}
            n = {k: v for k, v in enumerate(self.n) if v != 0}
            data = {
                'w': w,
                'z': z,
                'n': n
            }
            json.dump(data, f)

    def load_weight(self, model_file, D):
        """
        loada weight data
        :param model_file:
        :return:
        """
        with open(model_file, "r") as f:
            data = json.load(f)
            self.w = data.get('w', [0.] * D)
            self.z = data.get('z', [0.] * D)
            self.n = data.get('n', [0.] * D)

    @staticmethod
    def loss(y, y_pred):
        """
        log loss for LR model
        :param y:
        :param y_pred:
        :return:
        """
        p = max(min(y_pred, 1. - 10e-15), 10e-15)
        return -log(p) if y == 1. else -log(1. - p)


def data(f_train, D, dayfilter=None, dayfeature=True, counters=False):
    ''' GENERATOR: Apply hash-trick to the original csv row
                   and for simplicity, we one-hot-encode everything

        INPUT:
            path: path to training or testing file
            D: the max index that we can hash to

        YIELDS:
            ID: id of the instance, mainly useless
            x: a list of hashed and one-hot-encoded 'indices'
               we only need the index since all values are either 0 or 1
            y: y = 1 if we have a click, else we have y = 0
    '''

    device_ip_counter = {}
    device_id_counter = {}

    for t, row in enumerate(DictReader(f_train)):
        # process id
        ID = row['id']
        del row['id']

        # process clicks
        y = 0.
        if 'click' in row:
            if row['click'] == '1':
                y = 1.
            del row['click']

        # turn hour really into hour, it was originally YYMMDDHH

        date = row['hour'][0:6]
        row['hour'] = row['hour'][6:]

        if dayfilter != None and not date in dayfilter:
            continue

        if dayfeature:
            # extract date
            row['wd'] = str(int(date) % 7)
            row['wd_hour'] = "%s_%s" % (row['wd'], row['hour'])

        if counters:
            d_ip = row['device_ip']
            d_id = row["device_id"]
            try:
                device_ip_counter[d_ip] += 1
                device_id_counter[d_id] += 1
            except KeyError:
                device_ip_counter[d_ip] = 1
                device_id_counter[d_id] = 1
            row["ipc"] = str(min(device_ip_counter[d_ip], 8))
            row["idc"] = str(min(device_id_counter[d_id], 8))

        # build x
        x = [0]  # 0 is the index of the bias term
        for key in row:
            value = row[key]
            # one-hot encode everything with hash trick
            index = abs(hash(key + '_' + value)) % D
            x.append(index)
        yield t, ID, x, y
参考文献
  • Follow-the-Regularized-Leader and Mirror Descent:Equivalence Theorems and L1 Regularization(Google, FTRL原理论文)
  • Ad Click Prediction: a View from the Trenches(Google,FTRL工程化文档)
  • 在线最优化求解(冯杨,讲在线最优化算法非常好的一篇文档)
  • 统计学习方法(李航)
  • spark MLlib机器学习(黄美玲)
  • 机器学习(周航)
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,163评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,301评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,089评论 0 352
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,093评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,110评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,079评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,005评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,840评论 0 273
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,278评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,497评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,667评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,394评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,980评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,628评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,796评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,649评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,548评论 2 352

推荐阅读更多精彩内容