5
$\begingroup$

I've started working with pymc3 over the past few days, and after getting a feel for the basics, I've tried implementing the Probabilistic Matrix Factorization model.

For validation, I use a subset of the Jester dataset. I take the first 100 users who rated all 100 jokes. I use the first 20 jokes and leave the ratings unchanged; they are in the range [-10, 10] for all jokes. For ease of reference, I've made this subset available here.

import pymc3 as pm
import numpy as np
import pandas as pd
import theano

data = pd.read_csv('jester-dense-subset-100x20.csv')    
n, m = data.shape
test_size = m / 10
train_size = m - test_size

train = data.copy()
train.ix[:,train_size:] = np.nan  # remove test set data
train[train.isnull()] = train.mean().mean()  # mean value imputation
train = train.values

test = data.copy()
test.ix[:,:train_size] = np.nan  # remove train set data
test = test.values    

# Low precision reflects uncertainty; prevents overfitting
alpha_u = alpha_v = 1/np.var(train)
alpha = np.ones((n,m)) * 2  # fixed precision for likelihood function
dim = 10  # dimensionality

# Specify the model.
with pm.Model() as pmf:
    pmf_U = pm.MvNormal('U', mu=0, tau=alpha_u * np.eye(dim),
                        shape=(n, dim))
    pmf_V = pm.MvNormal('V', mu=0, tau=alpha_v * np.eye(dim),
                        shape=(m, dim))
    pmf_R = pm.Normal('R', mu=theano.tensor.dot(pmf_U, pmf_V.T),
                      tau=alpha, observed=train)

    # Find mode of posterior using optimization
    start = pm.find_MAP()  # Find starting values by optimization

This all appears to work splendidly, but the values produced by find_MAP end up being all 0s for both U and V, as can be seen by running:

(start['U'] == 0).all()
(start['V'] == 0).all()

I am relatively new to both Bayesian modeling and pymc, so I could easily be missing something obvious here. Why am I getting all 0s?

$\endgroup$
4
$\begingroup$

I did two things to fix your code. One was to initialize the model away from zero, the other one was to use a non-gradient based optimizer:

import pymc3 as pm
import numpy as np
import pandas as pd
import theano
import scipy as sp

data = pd.read_csv('jester-dense-subset-100x20.csv')    
n, m = data.shape
test_size = m / 10
train_size = m - test_size

train = data.copy()
train.ix[:,train_size:] = np.nan  # remove test set data
train[train.isnull()] = train.mean().mean()  # mean value imputation
train = train.values

test = data.copy()
test.ix[:,:train_size] = np.nan  # remove train set data
test = test.values    

# Low precision reflects uncertainty; prevents overfitting
alpha_u = alpha_v = 1/np.var(train)
alpha = np.ones((n,m)) * 2  # fixed precision for likelihood function
dim = 10  # dimensionality

# Specify the model.
with pm.Model() as pmf:
    pmf_U = pm.MvNormal('U', mu=0, tau=alpha_u * np.eye(dim),
                        shape=(n, dim), testval=np.random.randn(n, dim)*.01)
    pmf_V = pm.MvNormal('V', mu=0, tau=alpha_v * np.eye(dim),
                        shape=(m, dim), testval=np.random.randn(m, dim)*.01)
    pmf_R = pm.Normal('R', mu=theano.tensor.dot(pmf_U, pmf_V.T),
                      tau=alpha, observed=train)

    # Find mode of posterior using optimization
    start = pm.find_MAP(fmin=sp.optimize.fmin_powell)  # Find starting values by optimization

    step = pm.NUTS(scaling=start)
    trace = pm.sample(500, step, start=start)

This is an interesting model that would make a great contribution. Please consider adding this, once you're certain it works as desired, to the examples folder and do a pull request.

$\endgroup$
  • $\begingroup$ Thanks for the fixes! Do you mind commenting on why these changes are necessary to avoid the 0-valued start matrices? $\endgroup$ – Mack Apr 16 '15 at 16:49
  • $\begingroup$ Does it produce sensible results that way? Unfortunately I don't have a good explanation. I kind of intuited that starting at 0 would be a very odd space in the likelihood with an ill-defined gradient (I'm not sure if that's actually the case). This still didn't work with the gradient based optimization so I tried Powell which usually is pretty robust. You could try still starting at 0 with Powell which might be good enough. $\endgroup$ – twiecki Apr 17 '15 at 13:21
  • $\begingroup$ Just going off of RMSE, it seems the MAP estimate is reasonable. The test RMSE for 5.0823, which is slightly better than what comes out from using a simple baseline. If I just replace all missing values by the global mean, the test RMSE is 5.0839. The best I can produce with libFM is 4.980. $\endgroup$ – Mack Apr 17 '15 at 17:44
  • $\begingroup$ Note that the 4.980 number mentioned above was incorporating bias terms, both global and per-item/per-joke. That model was also using 8 dimensions. If I leave out the bias terms and use 10 dimensions, as I do for the PMF model, then the test RMSE is around 5.2 using libFM. This model is similar to a regularized SVD, which PMF is supposed to improve upon. $\endgroup$ – Mack Apr 17 '15 at 18:27
  • $\begingroup$ Following the advice of Salakhutdinov and Mnih, p.886, I also monitored the Frobenius norms of the sampled U and V matrices. They converge after 7-10 samples, which I believe is quite rapid. This also leads me to believe the MAP estimate is reasonable. Interestingly, the MAP estimate actually produces a better RMSE than the sampler, which gives a test RMSE of 5.108 using 400 samples, after discarding 10 samples as burn-in. $\endgroup$ – Mack Apr 17 '15 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.