numba tuned version of mf_numba.py, see: https://gist.github.com/jhemann/5584536 for previous versions and provenance. This has been updated to handle default arguments correctly.

mf_numba.py

# -*- coding: utf-8 -*-
#!/usr/bin/python
#
# Created by Albert Au Yeung (2010)
# http://www.quuxlabs.com/blog/2010/09/matrix-factorization-a-simple-tutorial-and-implementation-in-python/#source-code
#
# Modifications by Josh Hemann and Aron Ahmadia
#
# An implementation of matrix factorization
#

import numpy as np
from numba import typeof, double, int_
from numba.decorators import autojit, jit
"""
@INPUT:
    R     : a matrix to be factorized, dimension N x M
    P     : an initial matrix of dimension N x K
    Q     : an initial matrix of dimension M x K
    K     : the number of latent features
    steps : the maximum number of steps to perform the optimisation
    alpha : the learning rate
    beta  : the regularization parameter
@OUTPUT:
    the final matrices P and Q
"""

def matrix_factorization(R, P, Q, K, steps=5000, alpha=0.001, beta=0.015):
    return _matrix_factorization(R, P , Q, K, steps, alpha, beta)


@autojit
def _matrix_factorization(R, P, Q, K, steps, alpha, beta):
    half_beta = beta / 2.0
    N, M = R.shape

    for step in xrange(steps):
        for i in xrange(N):
            for j in xrange(M):
                if R[i,j] > 0:

                    eij = R[i,j]
                    for p in xrange(K):
                        eij -= P[i,p] * Q[j,p]

                    for k in xrange(K):
                        P[i,k] += alpha * (2 * eij * Q[j,k] - beta * P[i,k])
                        Q[j,k] += alpha * (2 * eij * P[i,k] - beta * Q[j,k])

        e = 0.0
        for i in xrange(N):
            for j in xrange(M):
                if R[i,j] > 0:

                    temp = R[i,j]
                    for p in xrange(K):
                        temp -= P[i,p] * Q[j,p]

                    e = e + temp*temp
                    for k in xrange(K):
                        e += half_beta * (P[i,k]*P[i,k] + Q[j,k]*Q[j,k])

        if e < 0.001:
            break

    return P, Q, step, e