從零開始構建邏輯回歸模型

邏輯回歸模型是針對線性可分問題的一種易於實現而且性能優異的分類模型。我們將分別使用Numpy和TensorFlow實現邏輯回歸模型訓練和預測過程，並且探討在大規模分布式系統中的工程實現。

從零構建

首先，我們通過Numpy構建一個邏輯回歸模型。
我們定義shape如下：
\(X\):（n,m）
\(Y\)：(1,m)
\(w\):（n,1）
\(b\):（1）
其中\(n\)代表特征維數，\(m\)代表樣本個數。
對於邏輯回歸二分類模型，其損失函數如下：

\[J(\theta)=-\frac{1}{m}\sum_{i=1}^N{[y_i\log{h_\theta(x_i)}+(1-y_i)\log{(1-h_\theta(x_i))}]} \]

對\(\theta\)求導得\(\theta_j\)的更新方式是：

\[ \theta_j:= \theta_j-\alpha \frac{1}{m}\sum_{i=1}^m (h_\theta(x_i)-y_i)x_i^j \]

所以，在代碼中，\(\theta\)的更新方式為：

dw = np.dot(X,(A-Y).T)/m

各個函數作用如下：

• sigmoid(x):激活函數實現
• initialization(dim):零值初始化w以及b
• propagate(w,b,X,Y)：前向傳播得到梯度以及代價函數值
• optimize(w,b,X,Y,learning_rate,epochs,print_cost=False)：反向傳播更新參數
• predict(w,b,X)：得到預測值
• model(X_train,Y_train,X_test,Y_test,epochs=200,learning_rate=0.01,print_cost=False)：建模

import numpy as np
def sigmoid(x):
    return 1/(1+np.exp(-x))

def initialization(dim):
    w = np.zeros((dim,1))
    b = 0
    return w,b

def propagate(w,b,X,Y):
    m = X.shape[1]
    A = sigmoid(np.dot(w.T,X)+b)
    cost = -1*np.sum(Y*np.log(A)+(1-Y)*np.log(1-A))/m
    dw = np.dot(X,(A-Y).T)/m
    db = np.sum(A-Y)/m
    grads = {'dw':dw,'db':db}
    return grads,cost

def optimize(w,b,X,Y,learning_rate,epochs,print_cost=False):
    costs = []
    for epoch in range(epochs):
        grads,cost = propagate(w,b,X,Y)
        dw = grads['dw']
        db = grads['db']
        w -= learning_rate * dw
        b -= learning_rate * db
        if epochs % 100 == 0:
            costs.append(cost)
            if print_cost:
                print('epochs:%i;cost:%f'%(epoch,cost))
    params = {'w':w,'b':b}
    return params,costs

def predict(w,b,X):
    predictions = sigmoid(np.dot(w.T,X)+b)
    return (predictions>0.5).astype(int)

def model(X_train,Y_train,X_test,Y_test,epochs=200,learning_rate=0.01,print_cost=False):
    dim = X_train.shape[0]
    w,b = initialization(dim)
    params, costs = optimize(w,b,X_train,Y_train,learning_rate,epochs,print_cost)
    w,b = params['w'],params['b']
    Y_predictions = predict(w,b,X_test)
    print('Test Acc:{}%'.format(100-np.mean(abs(Y_predictions-Y_test))*100))

if __name__ == '__main__':
    n = 20 #特征維度
    m = 200 #樣本數目
    X_train = np.random.random((n,m))
    Y_train = np.random.randint(0,2,size=(1,m))
    X_test = np.random.random((n,10))
    Y_test = np.random.randint(0,2,size=(1,10))

    model(X_train,Y_train,X_test,Y_test,epochs=200,learning_rate=0.01,print_cost=False)

TensorFlow版本

下面，我們實現下 TensorFlow版本的邏輯回歸模型。
這里采用了mnist數據集，將每張圖片\(28*28\)像素作為特征，使用\(softmax\)作為激活函數，最終我們的損失函數定義只有一行代碼：

cost = tf.reduce_mean(-tf.reduce_sum(y*tf.log(pred)))

下面，是具體的構建過程：

import tensorflow as tf

from tensorflow.examples.tutorials.mnist import input_data
mnist = input_data.read_data_sets('/tmp/data',one_hot=True)

learning_rate = 0.1
training_steps = 100
display_steps = training_steps//10
batch_size =64

X = tf.placeholder(tf.float32,[None,28*28])
y = tf.placeholder(tf.float32,[None,10])

W = tf.Variable(tf.zeros([784,10]))
b = tf.Variable(tf.zeros([10]))

pred = tf.nn.softmax(tf.add(tf.matmul(X,W),b))

cost = tf.reduce_mean(-tf.reduce_sum(y*tf.log(pred)))

optimizer = tf.train.GradientDescentOptimizer(learning_rate=learning_rate).minimize(cost)

init = tf.global_variables_initializer()

with tf.Session() as sess:
    sess.run(init)
    for epoch in range(training_steps):
        avg_cost = 0.
        total_batch = int(mnist.train.num_examples/batch_size)
        for i in range(total_batch):
            batch_x,batch_y = mnist.train.next_batch(batch_size)
            _,c = sess.run([optimizer,cost],feed_dict={X:batch_x,y:batch_y})
            avg_cost += c/total_batch
        if (epoch+1) % display_steps == 0:
            print('The epoch:%i,The cost:%9f'%(epoch+1,avg_cost))
    print('Finished')
    correct_prediction = tf.equal(tf.argmax(pred,1),tf.argmax(y,1))
    accuracy = tf.reduce_mean(tf.cast(correct_prediction,tf.float32))
    print('acc:',sess.run(accuracy,feed_dict={X: mnist.test.images[:3000], y: mnist.test.labels[:3000]}))

分布式實現

實際情況中，由於受到單機處理能力和效率的限制，在利用大規模樣本數據進行訓練的時候往往需要將求解LR問題的過程進行並行化，本文從並行化的角度討論LR的實現。
邏輯回歸問題的求解方法主要有三種：

• 梯度下降法
• 牛頓法(Newton Methods)：牛頓法是在當前參數\(\theta\)下，利用二次泰勒展開近似目標函數，然后利用該近似函數來求解目標函數的下降方向。該算法需要計算海森矩陣，因此算法需要花費大量的時間，迭代時間較長。牛頓法要求Hession矩陣是正定的，但在實際問題中，很難保證是正定的。
• 擬牛頓法(Quasi-Newton Methods)：使用近似算法，計算海森矩陣，從而降低算法每次迭代的時間，提高算法運行的效率。在擬牛頓算法中較為經典的算法有兩種：BFGS算法和L-BFGS算法。BFGS算法是利用原有的所有歷史計算結果，近似計算海森矩陣，雖然提高了整個算法的效率，但是由於需要保存大量歷史結果，因此該算法受到內存的大小的局限，限制了算法的應用范圍；而L-BFGS則是正是針對BFGS消耗內存較大的特點，為了解決空間復雜度的問題，只保存有限的計算結果，大大降低了算法對於內存的依賴。L-BFGS在特征量大時比BFGS實用，可以非常容易用map/reduce實現分布式求解，mapper求部分數據上的梯度，reducer求和並更新參數。它與梯度法實現復雜一點的地方在，它需要保存前幾次的模型，才能計算當前迭代的更新值。

由邏輯回歸問題的求解方法中可以看出，無論是梯度下降法、牛頓法、擬牛頓法，計算梯度都是其最基本的步驟，並且L-BFGS通過兩步循環計算牛頓方向的方法，避免了計算海森矩陣。因此邏輯回歸的並行化最主要的就是對目標函數梯度計算的並行化。從梯度更新公式中可以看出，目標函數的梯度向量計算中只需要進行向量間的點乘和相加，可以很容易將每個迭代過程拆分成相互獨立的計算步驟，由不同的節點進行獨立計算，然后歸並計算結果。

將M個樣本的標簽構成一個M維的標簽向量，M個N維特征向量構成一個M*N的樣本矩。其中特征矩陣每一行為一個特征向量（M行），列為特征維度（N列）。

如果將樣本矩陣按行划分，將樣本特征向量分布到不同的計算節點，由各計算節點完成自己所負責樣本的點乘與求和計算，然后將計算結果進行歸並，則實現了“按行並行的LR”。按行並行的LR解決了樣本數量的問題，但是實際情況中會存在針對高維特征向量進行邏輯回歸的場景（如廣告系統中的特征維度高達上億），僅僅按行進行並行處理，無法滿足這類場景的需求，因此還需要按列將高維的特征向量拆分成若干小的向量進行求解。

數據分割

假設所有計算節點排列成m行n列（m*n個計算節點），按行將樣本進行划分，每個計算節點分配M/m個樣本特征向量和分類標簽；按列對特征向量進行切分，每個節點上的特征向量分配N/n維特征。如圖所示，同一樣本的特征對應節點的行號相同，不同樣本相同維度的特征對應節點的列號相同。

一個樣本的特征向量被拆分到同一行不同列的節點中，即：

\[X_{r,k}=<X_{(r,1),k},X_{(r,2),k},...,X_{(r,n),k}> \]

其中\(X_{r,k}\)表示第\(r\)行的第\(k\) 個向量，\(X_{(r,c),k}\)表示\(X_{r,k}\)在第\(c\)列節點上的分量。同樣的，用\(W_c\)表示特征向量\(W\)在第\(c\)列節點上的分量，即：

\[W=<W_1,W_2,W_3...,W_n> \]

並行計算

觀察目標函數的梯度計算公式,其依賴於兩個計算結果：特征權重向量\(W_t\)和特征向量\(X_j\)的點乘，標量\((y_n-t_n)\)和特征向量\(X_j\)的相乘。其中\(y_n = \sigma(x)={1\over 1+e^{-XW}}\)可以將目標函數的梯度計算分成兩個並行化計算步驟和兩個結果歸並步驟：

• 各節點並行計算點乘，計算\(d_{(r,c),k,t}=W_{c,t}^T X_{(r,c),k}\)，其中\(k=1,2,\cdots,M/m\)，\(d_{(r,c),k,t}\)表示第\(t\)次迭代中節點\((r,c)\)上的第\(k\)個特征向量與特征權重分量的點乘，\(W_{c,t}\)為第\(t\)次迭代中特征權重向量在第\(c\)列節點上的分量。
• 對行號相同的節點歸並點乘結果：

\[d_{r,k,t}=W_t^TX_{r,k}=\sum_{c=1}^n d_{(r,c),k,t}=\sum_{c=1}^n W_{c,t}^T W_{(r,c),k} \]

計算得到的點乘結果需要返回到該行所有計算節點中，如圖所示。

• 各節點獨立算標量與特征向量相乘：

\[G_{(r,c),t}=\sum_{k=1}^{M/m}[\sigma(y_{r,k}d_{r,k,t})-1] y_{r,k} X_{(r,c),k} \]

\(G_{(r,c),t}\)可以理解為由第\(r\)行節點上部分樣本計算出的目標函數梯度向量在第\(c\)列節點上的分量。

• 對列號相同的節點進行歸並：

\[G_{c,t}=\sum_{r=1}^m G_{(r,c),t} \]

\(G_{c,t}\)就是目標函數的梯度向量\(G_t\)在第\(c\)列節點上的分量，對其進行歸並得到目標函數的梯度向量：

\[G_t=<G_{1,t},G_{2,t},...,G_{n,t}> \]

這個過程如圖所示。

綜合上述步驟，並行LR的計算流程如圖所示。並行LR實際上就是在求解損失函數最優解的過程中，針對尋找損失函數下降方向中的梯度方向計算作了並行化處理，而在利用梯度確定下降方向的過程中也可以采用並行化（如L-BFGS中的兩步循環法求牛頓方向）。

小結

1、按行並行。即將樣本拆分到不同的機器上去。其實很簡單，重新看梯度計算的公式：
\(\frac{1}{N}\sum_{n=1}^N(y_n-t_n)x_n\),其中\(y_n=\frac{1}{1+e^{-W^Tx_n}}\)。比如我們按行將其均分到兩台機器上去，則分布式的計算梯度，只不過是每台機器都計算出各自的梯度，然后 歸並求和再求其平均。為什么可以這么做呢？因為公式\(\frac{1}{N}\sum_{n=1}^N(y_n-t_n)x_n\)只與上一個時刻的以及當前樣本有關。所以就可以並行計算了。

2、按列並行。按列並行的意思就是將同一樣本的特征也分布到不同的機器中去。上面的公式為針對整個\(W\)，如果我們只是針對某個分量\(W_j\)，可得到對應的梯度計算公式\(\frac{1}{N}\sum_{n=1}^N (y_n-t_n)x_{n,j}\)即不再是乘以整個\(x_n\)，而是乘以\(x_n\)對應的分量\(x_{n,j}\)，此時可以發現，梯度計算公式僅與\(x_n\)中的特征有關系，我們就可以將特征分布到不同的計算上，分別計算\(W_j\)對應的梯度，最后歸並為整體的\(W\)，再按行歸並到整體的梯度更新。

參考：
【機器學習】LR的分布式（並行化）實現
LR如何並行化？