Matrix factorization algorithms factorize a matrix D into two matrices P and Q, such that D ≈ PQ. By limiting the dimensionality of P and Q, PQ provides a lowrank approximation of D. While singular value decomposition (SVD) can also be used for this same task, the matrix factorization algorithms considered in this post accommodate missing data in matrix D, unlike SVD.
For an overview of matrix factorization, I recommend Albert Au Yeung’s tutorial. That post describes matrix factorization, motivates the problem with a ratings prediction task, derives the gradients used by stochastic gradient descent, and implements the algorithm in Python.
For exploratory work, it would be great if the algorithm could be implemented in such a way that the gradients could be automatically derived. With such an approach, gradients would not have to be rederived when e.g., a change is made to the loss function (either the error term and/or the regularization term). In general, automatically derived gradients for machine learning problems facilitate increased exploration of ideas by removing a timeconsuming step.
Theano is a Python library that allows users to specify their problem symbolically using NumPybased syntax. The expressions are compiled to run efficiently on actual data. Theano’s webpage provides documentation and a tutorial.
The following code includes a Theanobased implementation of matrix factorization using batch gradient descent. The parameters are similar to those in the quuxlabs matrix factorization implementation. D is a secondorder masked numpy.ndarray (e.g., a ratings matrix, where the mask indicates missing ratings), and P and Q are the initial matrix factors. The elements of P and Q are the parameters of the model, which are initialized by the function’s caller. The rank of the factorization is specified by the dimensions of P and Q. For a rankk factorization, P must be and Q must be (where D is an matrix). Additional parameters specify the number of iterations, the learning rate, and the regularization importance.
The code doesn’t contain any derived gradients. It specifies the loss function and the parameters that the loss function will be minimized with respect to. Theano figures out the rest!

import numpy as np 

import numpy.ma as ma 

import theano 

from theano import tensor as T 



floatX = theano.config.floatX 





def getmask(D): 

return ma.getmaskarray(D) if ma.isMA(D) else np.zeros(D.shape, dtype=bool) 





def matrix_factorization_bgd( 

D, P, Q, steps=5000, alpha=0.0002, beta=0.02): 

P = theano.shared(P.astype(floatX)) 

Q = theano.shared(Q.astype(floatX)) 

X = T.matrix() 

error = T.sum(T.sqr(~getmask(D) * (P.dot(Q) – X))) 

regularization = (beta/2.0) * (T.sum(T.sqr(P)) + T.sum(T.sqr(Q))) 

cost = error + regularization 

gp, gq = T.grad(cost=cost, wrt=[P, Q]) 

train = theano.function(inputs=[X], 

updates=[(P, P – gp * alpha), (Q, Q – gq * alpha)]) 

for _ in xrange(steps): 

train(D) 

return P.get_value(), Q.get_value() 