FM在特征組合中的應用

 特征組合

	x1年齡	x2北京	x3上海	x4深圳	x5男	x6女
用戶1	23	1	0	0	1	0
用戶2	31	0	0	1	0	1

 如上例特征X有6個維度，年齡是連續值，城市和性別用one-hot表示，假設我們用最簡單的線性擬合來預測y值。

$$\hat{y}=w_0+\sum_{i=1}^n{w_ix_i}$$

實際中“北京的男性用戶”、“上海的女性用戶”這種組合特征可能是有用的，即$x_i,x_j$（$x_i,x_j$都是one-hot特征）同時為1時可能是一個很有用的特征，這種組合特征是$x_i$和$x_j$的線性組合所無法表示的。這樣一來乘積$x_ix_j$ 就成一個新的特征。為了不錯過任何一個這種可能有用的組合特征，我們窮舉所有的i,j組合，把$x_ix_j, 1\le{i}\le{n}, i<j\le{n}$都加到特征里面去，即使其中某些$x_i$不是one-hot特征或者某些 $x_ix_j$ 不是有用的特征，都沒關系，經過大量樣本的訓練，模型會把那些無用的特征的系數訓練為0。

$$\hat{y}=w_0+\sum_{i=1}^n{w_ix_i}+\sum_i^n{\sum_{j=i+1}^n{w_{ij}x_ix_j}}$$

這只是組合了2個特征，同樣道理我們組合任意三個特征、四個特征，隨着階數的提高，樣本會顯得非常稀疏，而且額外引入的參數呈指數增長。

Factorization Machines 

由於二次項系數$w_{ij}$，我們額外引入$\frac{n^2}{2}$個參數需要訓練。有沒有什么辦法可以減少參數？再來觀察二次項系數矩陣$W_{n\times n}$，它是對稱的方陣$w_{ij}=w_{ji}$，同時它是稀疏的，因為絕大部分的組合特征都是無用的，所以其系數應該為0。可以對$W_{n\times n}$進行矩陣分解$W_{n\times n}=V_{n\times k}V_{n\times k}^T$，即$w_{i,j}=<v_i,v_j>$。其中$k\ll n$，本來需要訓練的$n\times n$個參數，現在只需要訓練$n\times k$個。

$$\hat{y}=w_0+\sum_{i=1}^n{w_ix_i}+\sum_i^n{\sum_{j=i+1}^n{<v_i,v_j>x_ix_j}}$$

$$<v_i,v_j>=\sum_{f=1}^k{v_{if}v_{jf}}$$

根據$x$計算$\hat{y}$的時間復雜度是$O(kn^2)$

 $\sum_{i=1}^n{\sum_{j=1}^n{<v_i,v_j>x_ix_j}}$構成一個完整的對稱矩陣，$\sum_{i=1}^n{\sum_{j=i+1}^n{<v_i,v_j>x_ix_j}}$是這個對稱矩陣的上三角部分（不包含對角線）,所以$\sum_{i=1}^n{\sum_{j=i+1}^n{<v_i,v_j>x_ix_j}}$等於 $\sum_{i=1}^n{\sum_{j=1}^n{<v_i,v_j>x_ix_j}}$減去對角線再除以2。

\begin{equation}
 \sum_{i=1}^n{\sum_{j=i+1}^n{<v_i,v_j>x_ix_j}} \\
= \frac{1}{2}\sum_{i=1}^n{\sum_{j=1}^n{<v_i,v_j>x_ix_j}}-\frac{1}{2}\sum_{i=1}^n{<v_i,v_i>x_ix_i}\\
 = \frac{1}{2}\left(\sum_{i=1}^n{\sum_{j=1}^n{\sum_{f=1}^k{v_{if}v_{jf}x_ix_j}}}-\sum_{i=1}^n{\sum_{f=1}^k{v_{if}v_{if}x_ix_i}}\right) \\
 = \frac{1}{2}\left(\sum_{f=1}^k{\sum_{i=1}^n{v_{if}x_i\sum_{j=1}^n{v_{jf}x_j}}}-\sum_{i=1}^n{\sum_{f=1}^k{v_{if}v_{if}x_ix_i}}\right)
\end{equation}

因為$\sum_{i=1}^n{v_{if}x_i}$跟$j$沒有關系，$\sum_{j=1}^n{v_{jf}x_j}$跟$i$沒有關系，所以

$$\sum_{i=1}^n{v_{if}x_i\sum_{j=1}^n{v_{jf}x_j}}=\left(\sum_{i=1}^n{v_{if}x_i}\right)\left(\sum_{j=1}^n{v_{jf}x_j}\right)$$

\begin{equation}
\sum_{i=1}^n{\sum_{j=i+1}^n{<v_i,v_j>x_ix_j}}\\
= \frac{1}{2}\left(\sum_{f=1}^k{\left(\sum_{i=1}^n{v_{if}x_i}\right)\left(\sum_{j=1}^n{v_{jf}x_j}\right)}-\sum_{i=1}^n{\sum_{f=1}^k{v_{if}v_{if}x_ix_i}}\right)\\
= \frac{1}{2}\sum_{f=1}^k\left(\left(\sum_{i=1}^n{v_{if}x_i}\right)\left(\sum_{j=1}^n{v_{jf}x_j}\right)-\sum_{i=1}^n{v_{if}^2x_i^2}\right)\\
= \frac{1}{2}\sum_{f=1}^k\left(\left(\sum_{i=1}^n{v_{if}x_i}\right)^2-\sum_{i=1}^n{v_{if}^2x_i^2}\right)
\end{equation}

如此一來根據$x$求$\hat{y}$的時間復雜度就降為$O(kn)$。

用梯度下降法進行訓練時需要求$\hat{y}$對各個參數的偏導數：

$$\frac{\partial{\hat{y}}}{\partial{w_0}}=1$$

$$\frac{\partial{\hat{y}}}{\partial{w_i}}=x_i$$

$$\frac{\partial{\hat{y}}}{\partial{v_{if}}}=x_i\sum_{j=1}^n{v_{jf}x_j}-v_{if}x_i^2$$

在根據$x$計算$\hat{y}$的時候$\sum_{j=1}^n{v_{jf}x_j}$已經算好了，所以求$\frac{\partial{\hat{y}}}{\partial{v_{if}}}$的時間復雜度為O(1)，對所有參數求偏導的總的時間復雜度為$O(kn)$

代碼實現 

實際當中我們很少會直接用線性擬合來做預測，通常會再套一層sigmoid函數。

盡量避免使用for循環，盡量使用numpy的矩陣運算，因為numpy的矩陣運算做了並行處理。　　

# coding=utf-8
__author__ = 'orisun'

import numpy as np
np.random.seed(0)
import random


def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))


def sigmoid_prime(z):
    """
    sigmoid函數對z求一階偏導
    :param z:
    :return:
    """
    return sigmoid(z) * (1 - sigmoid(z))


class QuadraticCost(object):
    @staticmethod
    def fn(a, y):
        """
        平方誤差損失函數
        :param a: 預測值
        :param y: 真實值
        :return:
        """
        return 0.5 * np.linalg.norm(a - y) ** 2

    @staticmethod
    def delta(z, a, y):
        """
        損失函數對z求偏導
        :param z: x的線性函數
        :param a:
        :param y:
        :return:
        """
        return (a - y) * sigmoid_prime(z)


class FM(object):
    def __init__(self, train, valid, k, eta, maxecho, r2, cost=QuadraticCost):
        """
        構造函數
        :param train: 訓練數據
        :param valid: 驗證數據
        :param k: 矩陣V的第2維
        :param eta: 固定學習率
        :param maxecho: 最多迭代次數
        :param r2: R2小於該值后可停止迭代
        :param cost: 損失函數
        """
        self.train_x = train[:, :-1]
        self.train_y = train[:, -1:]
        self.valid_x = valid[:, :-1]
        self.valid_y = valid[:, -1:]
        self.var_y = np.var(self.valid_y)  # y的方差，在每輪迭代后計算R2時要用到
        self.k = k
        self.eta = float(eta)
        self.maxecho = maxecho
        self.r2 = r2
        self.cost = cost
        # 用正態分布隨機初始化參數W和V
        self.w0 = np.random.randn()
        self.w = np.random.randn(1, self.train_x.shape[1])
        self.v = np.random.randn(self.train_x.shape[1], self.k)

    def shuffle_data(self):
        """
        每輪訓練之前都隨機打亂樣本順序
        :return:
        """
        ids = range(len(self.train_x))
        random.shuffle(ids)
        self.train_x = self.train_x[ids]
        self.train_y = self.train_y[ids]

    def predict(self, x):
        """
        根據x求y
        :param x:
        :return:
        """
        z = self.w0 + np.dot(self.w, x.T).T + np.longlong(
            np.sum((np.dot(x, self.v) ** 2 - np.dot(x ** 2, self.v ** 2)),
                   axis=1).reshape(len(x), 1)) / 2.0

        return z, sigmoid(z)

    def evaluate(self):
        """
        在驗證集上計算R2
        :return:
        """
        _, y_hat = self.predict(self.valid_x)
        mse = np.sum((y_hat - self.valid_y) ** 2) / len(self.valid_y)
        r2 = 1.0 - mse / self.var_y
        print "r2={}".format(r2)
        return r2

    def update_mini_batch(self, x, y, eta):
        """
        平方誤差作為損失函數，梯度下降法更新參數
        :param x:
        :param y:
        :param eta: 學習率
        :return:
        """
        batch = len(x)
        step = eta / batch
        z, y_hat = self.predict(x)
        y_diff = self.cost.delta(z, y_hat, y)
        self.w0 -= step * np.sum(y_diff)
        self.w -= step * np.dot(y_diff.T, x)
        delta_v = np.zeros(self.v.shape)
        for i in xrange(batch):
            xi = x[i:i + 1, :]  # mini_batch中的第i個樣本。為保持shape不變，注意這里不能用x[i]
            delta_v += (np.outer(xi, np.dot(xi, self.v)) - xi.T ** 2 * self.v) * (y_diff[i])
        self.v -= step * delta_v

    def train(self, mini_batch=100):
        """
        采用批量梯度下降法訓練模型
        :param mini_batch:
        :return:
        """
        for itr in xrange(self.maxecho):
            print "iteration={}".format(itr)
            self.shuffle_data()
            n = len(self.train_x)
            for b in xrange(0, n, mini_batch):
                x = self.train_x[b:b + mini_batch]
                y = self.train_y[b:b + mini_batch]
                learn_rate = np.exp(-itr) * self.eta  # 學習率指數遞減
                self.update_mini_batch(x, y, learn_rate)

            if self.evaluate() > self.r2:
                break


def fake_data(sample, dim, k):
    """
    構造假數據
    :param sample:
    :param dim:
    :param k:
    :return:
    """
    w0 = np.random.randn()
    w = np.random.randn(1, dim)
    v = np.random.randn(dim, k)
    x = np.random.randn(sample, dim)
    z = w0 + np.dot(w, x.T).T + np.longlong(
        np.sum((np.dot(x, v) ** 2 - np.dot(x ** 2, v ** 2)),
               axis=1).reshape(len(x), 1)) / 2.0
    y = sigmoid(z)
    data = np.concatenate((x, y), axis=1)
    return z, data


if __name__ == "__main__":
    dim = 9  # 特征的維度
    k = dim / 3
    sample = 100
    z, data = fake_data(sample, dim, k)

    train_size = int(0.7 * sample)
    valid_size = int(0.2 * sample)
    train = data[:train_size]  # 訓練集
    valid = data[train_size:train_size + valid_size]  # 驗證集
    test = data[train_size + valid_size:]  # 測試集
    test_z = z[train_size + valid_size:]

    eta = 0.01  # 初始學習率
    maxecho = 200
    r2 = 0.9  # 擬合系數r2的最小值
    fm = FM(train, valid, k, eta, maxecho, r2)
    fm.train(mini_batch=50)

    test_x = test[:, :-1]
    test_y = test[:, -1:]
    print 'z=', test_z
    print "y=", test_y
    z_hat, y_hat = fm.predict(test_x)
    print "z_hat=", z_hat
    print "y_hat=", y_hat