[瞎玩儿系列] 使用SQL实现Logistic回归

本来想发在知乎专栏的，但是文章死活提交不了，我也是醉了，于是乎我就干脆提交到CNBLOGS了。

前言

前段时间我们介绍了Logistic的数学原理和C语言实现，而我呢？其实还是习惯使用Matlab进行计算的，而且是不带C的Matlab。（主要我们都用Windows）
那为什么要用SQL实现呢?（准确的说是PL/SQL）
因为我发现数据一次性加载进内存里面太大了，直接在SELECT的时候OutOfMemory了（其实数据是勉强能装进内存的，只是SELECT的时候产生的对象太多）
更主要的原因是因为我的电脑另有它用，留下的内存也不多了。
卧槽，为什么不用服务器算呢？
最近在重装系统，等我的小服务器安装好了，次回我会可能带来使用Hadoop/Spark的Logistic回归。
关于为什么我对Logistic回归这么着迷，并不是不会其它的模型，第一它简单，第二可解释性好，易于并行或者处理数据流。
为什么使用SQL呢？性能并不是其优势，反而是其软肋，但是可以把压力转嫁到服务器上，对于我残破不堪的工作电脑也是一个解脱，其次，对于特别大量的数据要做到随机梯度下降防止陷入局部极小，用SQL也算是一个解决方案。

诚然这个方案是不合适的，但是那又怎么样，我和SQL只是玩玩而已。

最近本人在找工作，希望找一个能让我做机器学习的岗位，我希望这家公司是一个脚踏实地的公司，有可持续的盈利模式，不会随便的喊出深度学习、人工智能和大数据之类的词汇，能从业务的角度来选择技术，那么我不会让你们失望。

准备数据

首先你有一张表，这张表列数不多，但是行数挺多的，其中一列是y，其余的是x，当然还可以有ID之类的一些其他信息。
我们这次的表结构是这样的：

CREATE TABLE public.jfeatures_cntf
(
    ---y
    cbuy integer,
    ---X
    cview integer,
    cadd integer,
    cdel integer,
    cstar integer,
    cclick integer
)

计算

首先你需要新建一个函数，该函数能做到从数据中随机取N行数据给你，因为数据量比较大，我们可能只有在Fine Tune的时候才会使用全部数据，平时的计算主要还是使用Radom Batch Gradient Descend/Ascend。
而Logistic的核心是：求偏导，我们也不需要什么都让SQL做，只要让SQL完成数据量最大的计算就行了。

随机取数据的函数

那我们首先构建一个PLSQL的函数：

CREATE OR REPLACE FUNCTION public.get_rand_x_record(x integer)
  RETURNS SETOF jfeatures_cntf AS
$BODY$
DECLARE
    N INTEGER;
BEGIN
    --- N<-Length-x
    SELECT count(*)-x INTO N FROM public.jfeatures_cntf;
    --- Random select
    RETURN QUERY SELECT * FROM public.jfeatures_cntf
        OFFSET floor(random()*N) LIMIT x;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE

这个函数可以每次从数据库中取出N条数据，比ORDER BY random()快不少。
由于在运算的时候我们也不会增删记录，所以可以预先获取数据大小N，随后使用这个函数：

CREATE OR REPLACE FUNCTION public.get_rand_x_record(
    x integer,
    n bigint)
  RETURNS SETOF jfeatures_cntf AS
$BODY$
DECLARE
BEGIN
    --- Random select
    RETURN QUERY SELECT * FROM public.jfeatures_cntf
        OFFSET floor(random()*N) LIMIT x;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE

Logistic回归求偏导的函数

得到数据以后，我们首先会求y，也就是1/(1+exp(1+b^Tx))，随后将(y-t)广播的乘到X上，最后求和就得到了结果。

CREATE OR REPLACE FUNCTION public.log_model_lr_random(
    batch_size integer,
    pview double precision,
    padd double precision,
    pdel double precision,
    pstar double precision,
    pclick double precision)
  RETURNS double precision[] AS
$BODY$
DECLARE
    OUT_VALUE RECORD;
BEGIN
    SELECT ---grad<-sum((y-t)*X)
        sum(log(cview+1)*D)/batch_size as pdview,
        sum(log(cadd+1)*D)/batch_size as pdadd,
        sum(log(cdel+1)*D)/batch_size as pddel,
        sum(log(cstar+1)*D)/batch_size as pdstar,
        sum(log(cclick+1)*D)/batch_size as pdclick
    INTO OUT_VALUE
    FROM
    (   ---get y-t from data
        SELECT *,
            (1.0/(exp(
            pview*log(cview+1) + 
            padd*log(cadd+1) +
            pdel*log(cdel+1) +
            pstar*log(cstar+1) +
            pclick*log(cclick+1) + 1.0
            )+1.0)
            ) - (case when cbuy>0 then 1.0 else 0.0 end)
            AS D
        FROM get_rand_x_record(batch_size)
    ) AS SUBS;
    return ARRAY[
        OUT_VALUE.pdview, 
        OUT_VALUE.pdadd,
        OUT_VALUE.pddel,
        OUT_VALUE.pdstar,
        OUT_VALUE.pdclick];
END;$BODY$
  LANGUAGE plpgsql VOLATILE

当然，这个函数也可以由N确定的版本（也就是如果你在计算过程中保证行数不变化的话可以使用的版本）,我最终使用的也就是这个版本。

这个就由大家自己写吧！

性能

关于性能方面，对3,000,000条数据求偏导需要1min！要知道，这在Matlab上（使用bsxfun做了并发）只需要0.5秒，这个性能差了100多倍（当然PostgreSQL在单次任务上不支持并行计算也是一个软肋），但是这个是有限定的，一个是内存计算，一个是外存计算，当数据量大到一定程度的时候，往往就需要外存算法。
Logistic是支持并行的，用SQL明显委屈他了，下次咱用Spark发挥出他最大的优势。