BPMF modificato in PyMC3 utilizzando priori `LKJCorr`: PositiveDefiniteError utilizzando` NUTS`

Question

In precedenza ho implementato il modello originale Bayesian Probabilistic Matrix Factorization (BPMF) in pymc3. See my previous question per riferimento, origine dati e impostazione problemi. Per la risposta a questa domanda di @twiecki, ho implementato una variante del modello usando i priori LKJCorr per le matrici di correlazione e i priori uniformi per le deviazioni standard. Nel modello originale, le matrici di covarianza sono tratte dalle distribuzioni di Wishart, ma a causa delle attuali limitazioni di pymc3, la distribuzione di Wishart non può essere campionata correttamente. This answer a una domanda liberamente correlata fornisce una spiegazione succinta per la scelta dei priori LKJCorr. Il nuovo modello è sotto.

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) 

L'obiettivo con questo reimplementazione era quello di produrre un modello che potrebbe essere stimato utilizzando il NUTS campionatore. Purtroppo, sto ancora ricevendo lo stesso errore all'ultima riga:

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

Ho fatto tutto il codice per PMF, BPMF, e questo BPMF modificato disponibile in this gist di rendere semplice per replicare l'errore. Tutto quello che devi fare è scaricare i dati (anch'essi referenziati nel significato).

Answer 1

Sembra che si sta passando la matrice completa precisione nella distribuzione normale:

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

Presumo che si desidera solo passare i valori diagonali:

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

Questo cambiamento mu_u e mu_v risolverlo per te?