59、用户近邻模型

import random as rd
import math

def AB(A, B):
    return [[sum([A[i][k]*B[k][j] for k in range(len(A[0]))]) for j in range(len(B[0]))] for i in range(len(A))]

def Ax(A, x):
    return [sum([A[i][j]*x[j] for j in range(len(A[0]))]) for i in range(len(A))]

def xy(x, y):
    return sum([x[i]*y[i] for i in range(len(x))])

def mod(v):
    return math.sqrt(sum(v[i]*v[i] for i in range(len(v))))

def x_next(x, di, ai):
    return [x[i]+di[i]*ai for i in range(len(x))]

#评分矩阵生成，1-5分
def R_generate(m, n):

    R = []
    for i in range(m):
        Ei = rd.choice(range(1,6))
        if Ei == 1:
            R.append([Ei + rd.choice(range(0, 2)) for j in range(n)])
        elif Ei == 5:
            R.append([Ei + rd.choice(range(-1, 1)) for j in range(n)])
        else:
            R.append([Ei + rd.choice(range(-1, 2)) for j in range(n)])

    return R

def sim_cos(X1, X2):

    if mod(X1) == 0 or mod(X2) == 0:
        return 0

    return sum([X1[j]*X2[j] for j in range(len(X1))])/(mod(X1)*mod(X2))

def sim_pearson(X1, X2):

    n = len(X1)

    E1 = sum(X1)/n
    E2 = sum(X2)/n

    dX1 = [X1[i]-E1 for i in range(n)]
    dX2 = [X2[i]-E2 for i in range(n)]

    if mod(dX1) == 0 or mod(dX2) == 0:
        return 0

    return sum([dX1[i]*dX2[i] for i in range(n)])/(mod(dX1)*mod(dX2))

def Rk_generate(R, k):

    m = len(R)

    #给每个用户找前k个最相关的其他用户
    Rk = []
    for i in range(m):
        Ni = [[sim_cos(R[i], R[j]), j] for j in range(m) if j != i]
        Ni = sorted(Ni, key = lambda x:x[0], reverse = True)
        Rk.append([Ni[j][1] for j in range(k)])

    return Rk


def run_one(m, n, k):

    #评分矩阵生成
    R = R_generate(m, n)

    #这里仅以评分矩阵作为相似性的度量标准，如果有其他信息可加入该函数中
    Rk = Rk_generate(R, k)

    #整体均值E，用户均值Em与物品均值En
    Em = [sum(R[i])/n for i in range(m)]
    En = [sum([R[i][j] for i in range(m)])/m for j in range(n)]
    E = sum(Em)/m

    dEm = [Em[i]-E for i in range(m)]
    dEn = [En[j]-E for j in range(n)]

    #残差矩阵
    dR = [[R[i][j] - E - dEm[i] - dEn[j] for j in range(n)] for i in range(m)]
    H = [[0 for j in range(m*k)] for i in range(m*k)]
    for i in range(m):
        for t in range(k):
            for l in range(k):
                H[i*k+t][i*k+l] = sum([dR[Rk[i][t]][j]*dR[Rk[i][l]][j] for j in range(n)])


    #初始化权重
    x0 = [1 for i in range(m*k)]

    def f(x):
        return sum([math.pow(dR[i][j]-sum([dR[Rk[i][t]][j]*x[i*k+t] for t in range(k)]), 2) for i in range(m) for j in range(n)])

    # 确定下降方向
    # 负梯度方向
    def d_gradient(x, **kwargs):
        return [-sum([dR[Rk[i][t]][j]*(sum([dR[Rk[i][l]][j]*x[i*k+l] for l in range(k)]) - dR[i][j]) for j in range(n)]) for i in range(m) for t in range(k)]


    #共轭梯度法方向
    def d_conjugate(x, **kwargs):
        g = [sum([dR[Rk[i][t]][j]*(sum([dR[Rk[i][l]][j]*x[i*k+l] for l in range(k)]) - dR[i][j]) for j in range(n)]) for i in range(m) for t in range(k)]
        if kwargs["loop"] == 0:
            return [-g[i] for i in range(len(g))]
        di = kwargs["di"]
        b = xy(Ax(H, di), g) / xy(Ax(H, di), di)
        return [-g[i] + b * di[i] for i in range(len(g))]

    def a_conjugate(x, f, di, **kwargs):
        g = [sum([dR[Rk[i][t]][j]*(sum([dR[Rk[i][l]][j]*x[i*k+l] for l in range(k)]) - dR[i][j]) for j in range(n)]) for i in range(m) for t in range(k)]
        return -xy(g, di) / xy(Ax(H, di), di) if xy(Ax(H, di), di) > 0 else 0

    # 确定步长
    # 分割法步长
    def a_three(x, f, di, **kwargs):
        left = 0
        right = 1

        # 放缩法确定区间
        while f(x_next(x, di, right)) < f(x):
            right = right * (2 - kwargs['r'])

        # 三分法确定步长
        mid1 = left * 2 / 3 + right / 3
        mid2 = left / 3 + right * 2 / 3
        while abs(left - right) * mod(di) > kwargs['step_min']:
            if f(x_next(x, di, mid1)) < f(x_next(x, di, mid2)):
                right = mid2
            else:
                left = mid1
            mid1 = left * 2 / 3 + right / 3
            mid2 = left / 3 + right * 2 / 3

        return left

    # 牛顿法步长
    def a_newton(x, f, di, **kwargs):
        f1 = sum([2*(sum([dR[Rk[i][t]][j]*x[i*k+t] for t in range(k)])-dR[i][j])*sum([dR[Rk[i][t]][j]*di[i*k+t] for t in range(k)]) for i in range(m) for j in range(n)])
        f2 = sum([2*math.pow(sum([dR[Rk[i][t]][j]*di[i*k+t] for t in range(k)]),2) for i in range(m) for j in range(n)])
        return -f1/f2 if f2 > 0 else 0

    def opt(x0, f, d, a, **kwargs):

        x = x0

        print("\nf(x):")
        print("loop  value")
        loop = 0
        di = [0 for i in range(len(x))]

        while True:
            di = d(x, loop=loop, di=di, **kwargs)
            ai = a(x, f, di, **kwargs)
            if abs(ai) * mod(di) < kwargs['step_min']:
                break
            x = x_next(x, di, ai)
            loop += 1
            print(loop, ":  ", f(x))

        print("\nx:\n", x)
        print("\ndi:\n", di)
        print("\nai:\n", ai)
        print("\nmod(di):\n", mod(di))

        return x

    print("\nf(x0):\n", f(x0))
    x = opt(x0, f, d_conjugate, a_three, r = 0.8, step_min = 0.01)
    print("\nf(x):\n", f(x))


if __name__ == '__main__':

    #m个用户，n个电影
    m = 100
    n = 10
    k = 10
    run_one(m, n, k)