Geändert BPMF in PyMC3 mit `LKJCorr` Prioren: PositiveDefiniteError mit` NUTS`

Question

Ich zuvor implementiert das Original Bayesian Probabilistic Matrix Factorization (BPMF) Modell in pymc3. See my previous question für Referenz, Datenquelle und Problemeinrichtung. Nach der Antwort auf diese Frage von @twiecki habe ich eine Variation des Modells implementiert, wobei LKJCorr Priors für die Korrelationsmatrizen und einheitliche Priors für die Standardabweichungen verwendet wurden. Im ursprünglichen Modell wurden die Kovarianzmatrizen aus Wishart-Verteilungen erstellt, aber aufgrund der aktuellen Einschränkungen von pymc3 kann die Wishart-Verteilung nicht ordnungsgemäß abgetastet werden. This answer zu einer lose verwandten Frage bietet eine knappe Erklärung für die Wahl von LKJCorr priors. Das neue Modell ist unten.

import pymc3 as pm 
import numpy as np 
import theano.tensor as t 


n, m = train.shape 
dim = 10 # dimensionality 
beta_0 = 1 # scaling factor for lambdas; unclear on its use 
alpha = 2 # fixed precision for likelihood function 
std = .05 # how much noise to use for model initialization 

# We will use separate priors for sigma and correlation matrix. 
# In order to convert the upper triangular correlation values to a 
# complete correlation matrix, we need to construct an index matrix: 
n_elem = dim * (dim - 1)/2 
tri_index = np.zeros([dim, dim], dtype=int) 
tri_index[np.triu_indices(dim, k=1)] = np.arange(n_elem) 
tri_index[np.triu_indices(dim, k=1)[::-1]] = np.arange(n_elem) 

logging.info('building the BPMF model') 
with pm.Model() as bpmf: 
    # Specify user feature matrix 
    sigma_u = pm.Uniform('sigma_u', shape=dim) 
    corr_triangle_u = pm.LKJCorr(
     'corr_u', n=1, p=dim, 
     testval=np.random.randn(n_elem) * std) 

    corr_matrix_u = corr_triangle_u[tri_index] 
    corr_matrix_u = t.fill_diagonal(corr_matrix_u, 1) 
    cov_matrix_u = t.diag(sigma_u).dot(corr_matrix_u.dot(t.diag(sigma_u))) 
    lambda_u = t.nlinalg.matrix_inverse(cov_matrix_u) 

    mu_u = pm.Normal(
     'mu_u', mu=0, tau=beta_0 * lambda_u, shape=dim, 
     testval=np.random.randn(dim) * std) 
    U = pm.MvNormal(
     'U', mu=mu_u, tau=lambda_u, 
     shape=(n, dim), testval=np.random.randn(n, dim) * std) 

    # Specify item feature matrix 
    sigma_v = pm.Uniform('sigma_v', shape=dim) 
    corr_triangle_v = pm.LKJCorr(
     'corr_v', n=1, p=dim, 
     testval=np.random.randn(n_elem) * std) 

    corr_matrix_v = corr_triangle_v[tri_index] 
    corr_matrix_v = t.fill_diagonal(corr_matrix_v, 1) 
    cov_matrix_v = t.diag(sigma_v).dot(corr_matrix_v.dot(t.diag(sigma_v))) 
    lambda_v = t.nlinalg.matrix_inverse(cov_matrix_v) 

    mu_v = pm.Normal(
     'mu_v', mu=0, tau=beta_0 * lambda_v, shape=dim, 
     testval=np.random.randn(dim) * std) 
    V = pm.MvNormal(
     'V', mu=mu_v, tau=lambda_v, 
     testval=np.random.randn(m, dim) * std) 

    # Specify rating likelihood function 
    R = pm.Normal(
     'R', mu=t.dot(U, V.T), tau=alpha * np.ones((n, m)), 
     observed=train) 

# `start` is the start dictionary obtained from running find_MAP for PMF. 
# See the previous post for PMF code. 
for key in bpmf.test_point: 
    if key not in start: 
     start[key] = bpmf.test_point[key] 

with bpmf: 
    step = pm.NUTS(scaling=start) 

Das Ziel dieser Neuimplementierung war es, ein Modell zu erzeugen, die mit dem NUTS Sampler geschätzt werden können. Leider bin immer, ich immer noch die gleichen Fehler in der letzten Zeile:

PositiveDefiniteError: Scaling is not positive definite. Simple check failed. Diagonal contains negatives. Check indexes [ 0 1 2 3 ... 1030 1031 1032 1033 1034 ] 

ich den gesamten Code für PMF gemacht haben, BPMF, und das modifizierte BPMF in this gist verfügbar zu machen es einfach, den Fehler zu reproduzieren. Alles, was Sie tun müssen, ist das Herunterladen der Daten (auch im Geiste referenziert).

Answer 1

Es sieht aus wie Sie die komplette Präzision Matrix in die Normalverteilung sind vorbei:

mu_u = pm.Normal(
    'mu_u', mu=0, tau=beta_0 * lambda_u, shape=dim, 
    testval=np.random.randn(dim) * std) 

Ich nehme an, Sie nur die diagonalen Werte eingeben möchten:

mu_u = pm.Normal(
    'mu_u', mu=0, tau=beta_0 * t.diag(lambda_u), shape=dim, 
    testval=np.random.randn(dim) * std) 

Hat diese Änderung mu_u und mu_v beheben Sie es für Sie?