推荐系统（二）

 (原创文章，转载请注明出处！)

一、问题

实现对电影的推荐，数据集中有约1600部电影，有约900个用户对这些电影进行了评价。设每个电影有10个特征，根据推荐系统（一）描述的算法，每个用户也相应的有10的参数，那么总的参数个数 ≈ 1600 * 10 + 900 * 10 ≈ 16000 + 9000 ≈25000 个。电影评价的数据集可以访问链接：http://grouplens.org/datasets/movielens/

二、代码实现

用R编程实现推荐系统（一）描述的算法。在实现中使用R提供的一个求解最优化的函数，optim。推荐算法的优化目标是无约束最优化问题，使用optim时，选择使用共轭梯度算法求解无约束最优化问题。optim的参数需要两个函数：目标函数和目标函数的一阶导函数。

1. 目标函数的代码如下：

JThetaXLamda <- function(thetaX) {
    theta <- matrix(thetaX[1:(nUser*nMovieFeature)], nrow=nUser, ncol=nMovieFeature, byrow=TRUE)  
    X     <- matrix(thetaX[(nUser*nMovieFeature+1):length(thetaX)], nrow=nMovie, ncol=nMovieFeature, byrow=TRUE)
    
    return (  (1/2) * ( sum( (tcrossprod(X, theta) * R - Y)^2 )  
                        + lamda * ( sum(theta^2) + sum(X^2) )  ) )
}

2. 目标函数的一阶导函数代码如下：

JThetaXGrLamda <- function(thetaX) {
    theta <- matrix(thetaX[1:(nUser*nMovieFeature)], nrow=nUser, ncol=nMovieFeature, byrow=TRUE)  # 943 is the number of user,  each row is one user's parameter
    X     <- matrix(thetaX[(nUser*nMovieFeature+1):length(thetaX)], nrow=nMovie, ncol=nMovieFeature, byrow=TRUE)  # 1682 is the number of movie
     
    # gradient of theta
    grTheta <- matrix(numeric(),ncol=nMovieFeature)
    for (j in 1:nUser) { # jth user
        idx <- which(R[,j]==1)  # which movies are rated by jth user
        tmpX <- X[idx,] 
        tmpY <- matrix(Y[idx,j], nrow=1)
        grTheta <- rbind( grTheta, alpha * ( (tcrossprod(theta[j,], tmpX) - tmpY) %*% tmpX 
                                                       + lamda * theta[j,] )  )
    }

    # gradient of x
    # calculating gradient of each movie only need the users who rate it.
    grX <- matrix(numeric(),ncol=nMovieFeature)
    for (i in 1:nMovie) { # ith movie
        idx <- which(R[i,]==1)  # exclude the user who haven't rated movie i
        tmpTheta <- theta[idx,] 
        tmpY <- matrix(Y[i,idx], nrow=1) # tmpY is a vector including rating credit of movie i
        grX <- rbind( grX, alpha * ( (tcrossprod(X[i,], tmpTheta) - tmpY) %*% tmpTheta 
                                               + lamda * X[i,]) )
    }
    
    return (   c( t(grTheta), t(grX) )   )
}

3. 训练：

上面两个函数使用的变量：nUser、nMovie、nMovieFeature、Y、R是自由变量，在训练代码中定义和赋值。

推荐系统（一）中还提到对数据需要做进行均值正规化（Mean Normalization）, 代码如下：

## normalizae the training data
## Args :
##     Y  -  a matrix, rating of movies by users 
##     R  -  a matrix, flags, which movies are rated by user
## Returns :
##     a list contains two elements: Ynormalized - a matrix, normalized Y
##                                   Ymean - a vector contains mean of each movie
craEx8RS_meanNormalization <- function(Y, R)
{
    numOfMovies <- dim(Y)[1]
    l <- list( Ynormalized = Y, Ymean = numeric(length = numOfMovies) )
    for (i in 1:numOfMovies) {
        idx <- which(R[i,] == 1)
        meanOfOneMovie <- mean(Y[i,idx])
        l$Ynormalized[i,] <- Y[i,] - meanOfOneMovie
        l$Ymean[i] <- meanOfOneMovie
    }
    return( l )
}

训练代码如下：

## do the training
Y <- read.table(file="Y.txt")
R <- read.table(file="R.txt")
Y <- as.matrix(Y)
R <- as.matrix(R)
l <- meanNormalization(Y, R)  # normalize the data
Y <- l$Ynormalized
nMovie <- dim(Y)[1]
nUser  <- dim(Y)[2]
nMovieFeature <- 10 
lamda <- 0.01
alpha <- 0.01


## random initialize theta, X
theta <- matrix(runif(nUser*nMovieFeature,min=0,max=0.1), nrow=nUser, ncol=nMovieFeature)  # each row is parameter of one user
X     <- matrix(runif(nMovie*nMovieFeature,min=0,max=0.1), nrow=nMovie, ncol=nMovieFeature)  # each row is parameter of one movie
# convert the two initialize value matrix to vector
thetaX <- c( t(theta), t(X) ) 
# using conjugate gradient algorithm to train the parameters
optimCG <- optim( par=thetaX, fn=JThetaX, gr = JThetaXGr, method = "CG", control = list(maxit=500000, abstol=1e-6) )
optimCG$par # display the training reslut, i.e. best value has been found for parameters  

4. 预测：
通过训练得到参数值后，就可以通过推荐系统（一）中提到的计算公式:  θ^T*x + 平均值, 来计算每个用户对还没有评价电影的评分。通过比较计算得到评分，就能找到用户偏好的电影推荐给他。代码如下：

## Find the top-n un-rating movies of the jth user
## Args :
##     topn  -  a positive integer
##     theta_j  -  a vector, contains theta parameters of jth user
##     X  -  a matrix, contains X paramters
##     Y  -  a matrix, rating of movies by users 
##     R  -  a matrix, flags, which movies are rated by user
##     j  -  a interger, jth user
##     Ymean  -  a vector, contains mean of all movies
##     isAddMean  -  flag, is the mean added to rating credit, or not?
##                   TRUE  -  add the mean to rating credit
##                   FALSE  -  the mean will not be added to the rating credit
## Returns :
##     a list, contains the index of the top-n un-rating movies of jth user and the rating result
findTopN <- function(topn, theta_j, X, R, j, Ymean=NULL, isAddMean=FALSE)
{
    idx <- which(R[,j]==0)    
    Yvec <- rep(-1,length(R[,j]))
      
    if (isAddMean == TRUE) {
        if (is.null(Ymean)) {
            cat("Error: Ymean is NULL!!!
")
        }
        Yvec[idx] <- tcrossprod( X[idx,], matrix(theta_j, nrow=1) )   +   as.matrix( Ymean[idx], ncol=1)
    } else {
        Yvec[idx] <- tcrossprod( X[idx,], matrix(theta_j, nrow=1) )   
    }
    
    sortTating <- sort(Yvec, method = "sh", decreasing = TRUE, index.return = TRUE)
    objList <- list(ratingResult = Yvec, topnIndex = sortTating$ix[1:topn])
    return( objList )
}

（寻找N个数中最大的K个数的方法，可参考链接：http://zhangzhibiao02005.blog.163.com/blog/static/3736782020112123353695/）