more efficient way to write PPCA model?

[Lei Li]

I am experimenting with Tipping & Bishop's PPCA model. However, two versions I tried in Stan seem slow. Is there better way to write an efficient PPCA model?
Thanks!
Lei

PPCA version 1:

data {
  int<lower=1> N;              // num datapoints
  int<lower=1> D;              // num dimension
  int<lower=1> K;              // num basis
  real<lower=0> sigma1;
  //real<lower=0> sigma2;
  real<lower=0> sigma3;
  real<lower=0> sigma4;      
  real y[D,N];
}
transformed data {
  matrix[K,K] Sigma;
  matrix[D,D] Sigma2;
  vector<lower=0>[K] diag_elem;
  vector<lower=0>[K] zr_vec;
  vector<lower=0>[D] diag_elem2;
  vector<lower=0>[D] zr_vec2;
  for (k in 1:K) zr_vec[k] <- 0;
  for (k in 1:K) diag_elem[k] <- sigma3;
  Sigma <- diag_matrix(diag_elem);
  for (d in 1:D) zr_vec2[d] <- 0;
  for (d in 1:D) diag_elem2[d] <- sigma1;
  Sigma2 <- diag_matrix(diag_elem2);
}
parameters {
  matrix[K,D] A;            // basis
  vector[K] x[N];           // coefficients
}
model {
 
  //vector[K] x[N];
 
  for (i in 1:N)
      x[i] ~ multi_normal(zr_vec, Sigma);

  for (k in 1:K)
    A[k] ~ multi_normal(zr_vec2, Sigma2);
   
  for (i in 1:N)
    //y[i] ~ multi_normal(A * x[i], Sigma);
    for (d in 1:D)
      //y[d,i] ~ normal(dot_product(row(A, d), x[i]), sigma4);
      increment_log_prob(normal_log(y[d,i], dot_product(col(A, d), x[i]), sigma4));
}

PPCA version 2:

data {
  int<lower=1> N;              // num datapoints
  int<lower=1> D;              // num dimension
  int<lower=1> K;              // num basis
  real<lower=0> sigma1;
  //real<lower=0> sigma2;
  real<lower=0> sigma3;
  real<lower=0> sigma4;      
  real y[D,N];
}
transformed data {
  matrix[K,K] Sigma;
  vector<lower=0>[K] diag_elem;
  vector<lower=0>[K] zr_vec;
  for (k in 1:K) zr_vec[k] <- 0;
  for (k in 1:K) diag_elem[k] <- sigma3;
  Sigma <- diag_matrix(diag_elem);
}
parameters {
  matrix[D,K] A;            // basis
  vector[K] x[N];           // coefficients
}
model {  
  //vector[K] x[N];
 
  for (i in 1:N)
      x[i] ~ multi_normal(zr_vec, Sigma);
     
  for (d in 1:D)
    for (k in 1:K)
      A[d,k] ~ normal(0, sigma1);
 
  for (i in 1:N)
    //y[i] ~ multi_normal(A * x[i], Sigma);
    for (d in 1:D)
      //y[d,i] ~ normal(dot_product(row(A, d), x[i]), sigma4);
      increment_log_prob(normal_log(y[d,i], dot_product(row(A, d), x[i]), sigma4));
}

[Bob Carpenter]

Without worrying about what the model's doing, I can give
you some performance tips. The basic idea is that you
don't want to create a diagonal covariance matrix and then
sample from a multivariate normal.

Because Sigma = diag_matrix(diag_elem) and diag_elem is just copies of sigma3,
and zr_vec is all zeroes, this is much more efficiently done as:

to_vector(x) ~ normal(0,sigma3);

>
> for (k in 1:K)
> A[k] ~ multi_normal(zr_vec2, Sigma2);

Same again

>
> for (i in 1:N)
> //y[i] ~ multi_normal(A * x[i], Sigma);
> for (d in 1:D)
> //y[d,i] ~ normal(dot_product(row(A, d), x[i]), sigma4);
> increment_log_prob(normal_log(y[d,i], dot_product(col(A, d), x[i]), sigma4));

If you can reduce this to matrix algebra with A and x and then
vectorize, it'll be a lot faster.

The second version is similar in terms of Stan coding.

- Bob

[Herra Huu]

It's possible to analytically integrate out either:

a) x, as was done in the original article by Tipping & Bishop (see section 3.1 http://research.microsoft.com/pubs/67218/bishop-ppca-jrss.pdf)
b) A, see article: Probabilistic Non-linear Principal Component Analysis with
Gaussian Process Latent Variable Models by Neil Lawrence (link: http://jmlr.org/papers/volume6/lawrence05a/lawrence05a.pdf)

Both will give simple mv-Gaussian likelihood with certain mean&covariance matrix. By doing option b), we can easily see that PCA correspond to Gaussian process with linear covariance function (and by changing the covariance function we can get non-linear dimension reduction models etc)

[Bob Carpenter]

> Both will give simple mv-Gaussian likelihood with certain mean&covariance matrix. By doing option b), we can easily see that PCA correspond to Gaussian process with linear covariance function (and by changing the covariance functions we can get non-linear dimension reduction models etc, these are called Gaussian Process latent variable models, GP-LVM)

Cool --- thanks for sharing this. I'd love to hear how the direct
vs. integrated-out implementation works in Stan.

- Bob

[Lei Li]

Thanks for the pointer!
Yes. The integrated version is EM approach. Stan's algorithm is HMC.

Of course, one can even do PCA via SVD. That certainly has value. But I am particularly interested in how to make it efficient in current Stan's framework.

Lei

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/DYofzSAz1sk/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

You can integrate out parameters in Stan, too. In fact,
you have to for discrete parameters. You just need to
write down the appropriate density from which you want to sample.
You can then use the expectations of the parameters, which is
typically all you need for inference.

The bonus is that integrating out some of the parameters
almost always makes things faster (computing expectations
is almost always more efficient both computationally and
statistically than sampling --- it's the Rao-Blackwell theorem
in action).

- Bob

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.