This tutorial is completely based on the explanation given by Simon Funk.
First the quick matrix version;
We have a matrix \(A\) which we are trying to approximate with the product of \(U\) and \(V\):
\(A \approx UV\)
We do this by minimising the sum squared error on all the elements in the matrix. If the row of \(U\) can be written \(u_i\) and the column of \(V\) can be written \(v_j\); the calculation for the \(i,j\) element of \(A\) will be the matrix product of these two vectors; \(u_i . v_j\). So we can write the error function as:
\(E = \sum_i \sum_j \left(A_{ij} - u_i v_j \right)^2\)
We want to minimise \(E\), with respect to these two vectors, \(U_i\) first:
\({{\delta E \over \delta u_i}} = 2 \sum_j \left(A_{ij} - u_i v_j \right) v_j\)
To detour slightly, and check this long-hand (i.e. avoid vector calculus):
The earlier expression for \(E\) can be written as:
\(E = \sum_i \sum_j \left(A_{ij} - \sum_f (U_{if} V_{fj}) \right)^2\)
Where \(f\) is the index of the 'feature'.
We can again differentiate, but this time with respect to \(U_{if}\):
\({{\delta E \over \delta U_{if}}} = 2 \sum_j \left(A_{ij} - \sum_f (U_{if} V_{fj}) \right) V_{fj}\)
Where we can simplify f (U{if} V_{fj}) as \(u_i v_j\) again, and note that if we apply the rule across all features \(f\) we will have the matrix-expression from earlier.
The same logic can be applied to find the derivative wrt \(v_j\), leaving us with:
\({{\delta E \over \delta v_j}} = 2 \sum_j \left(A_{ij} - u_i v_j \right) u_i\)
We can solve this numerically (which will be useful later) by simply following the gradient of the error function. Note also we can do this with a loop, over all j, (ignore the 2).
So one iteration looks like:
\(D = A - U.V\)
for every value of \(j\): \(\Delta u_i = r_{l} D_{ij} v_j\)
for every value of \(i\): \(\Delta v_j = r_{l} D_{ij} u_i\)
where \(r_l\) is the learning rate.
Properly, \(U\) and \(V\) should be updated once per iteration, but this will do for now (with a low learning rate).
The data from the last notebook;
import numpy as np
np.set_printoptions(precision=2)
np.set_printoptions(suppress=True)
A=np.array([[5,0,5,0],[4,1,3,0],[0,4,1,5],[5,1,3,1],[4,0,4,1],[1,3,0,4],[1,3,0,3],[3,2,4,1],[0,5,0,5],[0,4,1,4]])
A
W=np.array([[1,1,0,1],[0,1,1,0],[1,1,0,1],[1,0,1,1],[1,0,1,1],[1,0,0,1],[0,1,1,1],[0,1,1,0],[1,1,1,0],[1,1,1,1]])
W
The code to perform the algorithm described above,
#factorises matrix A, using features number of factors
def factoriseMatrix(A,features):
U = np.random.randn(A.shape[0],features) #initialising with random numbers
V = np.random.randn(features,A.shape[1]) #seems to be good, if non-deterministic
lr = 0.05 #the learning rate
for it in range(500): #could replace with a check on when err**2 has stopped changing much
lr = lr * 0.995 #simulated annealing sort of
err = A - np.dot(U,V) #how far we are from the actual matrix
deltaU = np.zeros(U.shape);
deltaV = np.zeros(V.shape);
for i in range(A.shape[0]): #for each user...
for j in range(A.shape[1]): #for each film...
deltaU[i,:] += err[i,j]*V[:,j]; #add the derivative to the vector
deltaV[:,j] += err[i,j]*U[i,:];
U += deltaU * lr;
V += deltaV * lr;
return U,V,np.sum(err**2)
#run the algorithm
U,V,e = factoriseMatrix(A,2)
print "%0.2f" % e
Compare the approximation with 2 factors with the actual data (with four columns):
print np.dot(U,V)
print A
We can see something of how the two feature vectors are organised by looking at \(U\) and \(V\). It is interesting that even these two rows aren't very independent, but if we reduce the number of features to one the estimate becomes very poor. This is because we've not subtracted the mean, explaining why two features are required.
V
If we subtract the mean from the two dimensions of A, and try with just one factor,
U,V,e = factoriseMatrix(A,1)
print "Error, with one factor, without demeaning: %0.2f" % e
meanA = A.mean(axis=1)
B = A - meanA[:, None]
meanB = B.mean(axis=0)
C = B - meanB[None, :]
U,V,e = factoriseMatrix(C,1)
print "Error, with one factor, *with* demeaning: %0.2f" % e
print np.dot(U,V)+meanA[:, None]+meanB[None, :]
print A
print "---"
print "The V matrix:"
print V
Clearly, to use one factor, we needed to subtract the mean. Hopefully it's fairly clear why this is the case (hint: If the mean isn't subtracted, one factor is effectively dedicated to adding the mean, while the other factor can deal with the differences between columns).
Finally, handling missing data. This is roughly the same function, but in this case only the non-zero elements of \(W\) are used. Note that \(A\) is really still a complete matrix, so we can compare the results from the function with the elements the function ignored, effectively using them as test data.
#W is a matrix of missing features (1=available, 0=missing)
#this function also subtracts the mean of A
def factoriseMatrixMissingData(A,W,features):
meanA = A.mean(axis=1)
A = A - meanA[:, None]
meanB = A.mean(axis=0)
A = A - meanB[None, :]
U = np.random.randn(A.shape[0],features)
V = np.random.randn(features,A.shape[1])
#U = np.ones([A.shape[0],2])*.0
#V = np.ones([2,A.shape[1]])*.0
lr = 0.05
for it in range(200):
lr = lr * 0.99
err = (A - np.dot(U,V))*W
deltaU = np.zeros(U.shape);
deltaV = np.zeros(V.shape);
for i in range(A.shape[0]): #for each user...
for j in range(A.shape[1]): #for each film...
deltaU[i,:] += err[i,j]*V[:,j];
deltaV[:,j] += err[i,j]*U[i,:];
U += deltaU * lr;
V += deltaV * lr;
guess = np.dot(U,V)+meanA[:, None]+meanB[None, :]
TrainErr =np.sum(((A - guess)*W)**2)
TestErr =np.sum(((A - guess)*(1-W))**2)
return U,V,TrainErr,TestErr, guess
Running the above function with different numbers of features (between 1 to 4), shows that the test data becomes less well fitted with more features (an example of over fitting). Not surprisingly, looking at the data, there's only one clear feature, which is why that does best.
errs = [];
for testit in range(10):
U,V,TrainErr,TestErr,guess = factoriseMatrixMissingData(A,W,1) #ONE FACTOR
errs.append(TestErr);
print errs
print guess
print A
print V
errs = [];
for testit in range(10):
U,V,TrainErr,TestErr,guess = factoriseMatrixMissingData(A,W,2) #TWO FACTORS
errs.append(TestErr);
print np.mean(errs)
errs = [];
for testit in range(10):
U,V,TrainErr,TestErr,guess = factoriseMatrixMissingData(A,W,3) #THREE FACTORS
errs.append(TestErr);
print np.mean(errs)
errs = [];
for testit in range(10):
U,V,TrainErr,TestErr,guess = factoriseMatrixMissingData(A,W,4) #FOUR FACTORS
errs.append(TestErr);
print np.mean(errs)
Next, applying this method to the movieLens data...