[瞎玩兒系列] 使用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發揮出他最大的優勢。