Matrix dimension limit?

Developers
1

I just got an error when compiling a model (with CmdStan) that makes me think there is limit on matrix dimension, namely that matrices must have 8 or fewer dimensions. Is that correct? The actual error was about the “to_vec” function and finding no matches that took 9 arguments, but matches for 8 and 7, etc. And the model compiles and runs once I reduce to an 8-dimensional matrix. I’m interpreting this as I said above but perhaps it originates in something I am trying to do with the matrix?

Is this difficult to change? I can work around it. The matrix in question is a post-stratification weight matrix for a model with several terms. But rather than use a multi-dimensional matrix, I could just use a long vector and indexes for each factor. That would make the json for the matrix simpler, at the cost of needing a lot of those new indexes. Is that the way people usually handle this?

Thanks!

2

Can you post an example of the code that’s failing?

3

Sure! It’s…not minimal. Or correctly formatted. I’m generating via a Haskell library so it’s harder to isolate just the bit that fails. I’m compiling using cmdStan 2.26.1, on OS X Big Sur with the built-in clang++.

Thanks!

data {
  int<lower=2> J_State;
  int<lower=1> N;
  int<lower=1> State[N];
  int<lower=2> J_Education;
  int<lower=1> Education[N];
  int<lower=2> J_Race;
  int<lower=1> Race[N];
  int<lower=2> J_Age;
  int<lower=1> Age[N];
  int<lower=2> J_WNH;
  int<lower=1> WNH[N];
  int<lower=2> J_Ethnicity;
  int<lower=1> Ethnicity[N];
  int<lower=2> J_CD;
  int<lower=1> CD[N];
  int<lower=2> J_WhiteNonGrad;
  int<lower=1> WhiteNonGrad[N];
  int<lower=2> J_Sex;
  int<lower=1> Sex[N];
  int<lower=0> T[N];
  int<lower=0> S[N];
  int K_CD;
  matrix[J_CD,K_CD] X_CD;
  int<lower=0> N_WI_ACS_WNH_State;
  real WI_ACS_WNH_State_wgts[N_WI_ACS_WNH_State,J_Education,J_Race,J_Age,J_WNH,J_Ethnicity,J_CD,J_WhiteNonGrad,J_Sex];
  int<lower=0> N_WI_ACS_NWNH_State;
  real WI_ACS_NWNH_State_wgts[N_WI_ACS_NWNH_State,J_Education,J_Race,J_Age,J_WNH,J_Ethnicity,J_CD,J_WhiteNonGrad,J_Sex];
  int<lower=0> N_NI_ACS_WNH_State;
  real NI_ACS_WNH_State_wgts[N_NI_ACS_WNH_State,J_Education,J_Race,J_Age,J_WNH,J_Ethnicity,J_CD,J_WhiteNonGrad,J_Sex];
  int<lower=0> N_NI_ACS_NWNH_State;
  real NI_ACS_NWNH_State_wgts[N_NI_ACS_NWNH_State,J_Education,J_Race,J_Age,J_WNH,J_Ethnicity,J_CD,J_WhiteNonGrad,J_Sex];
}
transformed data {
  vector[K_CD] mean_X_CD;
  matrix[J_CD,K_CD] centered_X_CD;
  for (k in 1:K_CD) {
    mean_X_CD[k] = mean(X_CD[,k]);
    centered_X_CD[,k] = X_CD[,k] - mean_X_CD[k];
  }
  matrix[J_CD,K_CD] Q_CD_ast = qr_thin_Q(centered_X_CD) * sqrt(J_CD - 1);
  matrix[K_CD,K_CD] R_CD_ast = qr_thin_R(centered_X_CD) / sqrt(J_CD - 1);
  matrix[K_CD,K_CD] R_CD_ast_inverse = inverse(R_CD_ast);
}
parameters {
  real alpha;
  vector[K_CD] thetaX_CD;
  real eps_Sex;
  real eps_WNH;
  real<lower=0> sigma_Race;
  vector[J_Race] beta_Race_raw;
  real eps_Ethnicity;
  real eps_Age;
  real eps_Education;
  real eps_WhiteNonGrad;
  real<lower=0> sigma_State;
  vector[J_State] beta_State_raw;
  real<lower=0> sigma_WNH_State;
  vector[J_State] eps_WNH_State_raw;
}
transformed parameters {
  vector[K_CD] betaX_CD;
  betaX_CD = R_CD_ast_inverse * thetaX_CD;
  vector[J_Race] beta_Race = sigma_Race * beta_Race_raw;
  vector[J_State] beta_State = sigma_State * beta_State_raw;
  vector[J_State] eps_WNH_State = sigma_WNH_State * eps_WNH_State_raw;
  vector[N] y_WNH_State;
  for (n in 1:N) {
    y_WNH_State[n] = {eps_WNH_State[State[n]],-eps_WNH_State[State[n]]}[WNH[n]];
  }
}
model {
  alpha ~ normal(0,2);
  thetaX_CD ~ normal(0,2);
  eps_Sex ~ normal(0,2);
  eps_WNH ~ normal(0,2);
  sigma_Race ~ normal(0,2);
  beta_Race_raw ~ normal(0,1);
  eps_Ethnicity ~ normal(0,2);
  eps_Age ~ normal(0,2);
  eps_Education ~ normal(0,2);
  eps_WhiteNonGrad ~ normal(0,2);
  sigma_State ~ normal(0,2);
  beta_State_raw ~ normal(0,1);
  sigma_WNH_State ~ normal(0,2);
  eps_WNH_State_raw ~ normal(0,1);
  S ~ binomial_logit(T,alpha + Q_CD_ast[CD] * thetaX_CD + to_vector({eps_Sex,-eps_Sex}[Sex]) + beta_Race[Race] + to_vector({eps_Ethnicity,-eps_Ethnicity}[Ethnicity]) + to_vector({eps_Age,-eps_Age}[Age]) + to_vector({eps_Education,-eps_Education}[Education]) + to_vector({eps_WhiteNonGrad,-eps_WhiteNonGrad}[WhiteNonGrad]) + beta_State[State] + y_WNH_State);
}
generated quantities {
  vector[N] log_lik;
  for (n in 1:N) {
    log_lik[n] = binomial_logit_lpmf(S[n] | T[n],alpha + Q_CD_ast[CD[n]] * thetaX_CD + {eps_Sex,-eps_Sex}[Sex[n]] + beta_Race[Race[n]] + {eps_Ethnicity,-eps_Ethnicity}[Ethnicity[n]] + {eps_Age,-eps_Age}[Age[n]] + {eps_Education,-eps_Education}[Education[n]] + {eps_WhiteNonGrad,-eps_WhiteNonGrad}[WhiteNonGrad[n]] + beta_State[State[n]] + {eps_WNH_State[State[n]],-eps_WNH_State[State[n]]}[WNH[n]]);
  }
  vector[N_WI_ACS_WNH_State] WI_ACS_WNH_State = rep_vector(0,N_WI_ACS_WNH_State);
  for (n in 1:N_WI_ACS_WNH_State) {
    real WI_ACS_WNH_State_WgtSum = 0;
    for (n_Education in 1:J_Education) {
      for (n_Race in 1:J_Race) {
        for (n_Age in 1:J_Age) {
          for (n_WNH in 1:J_WNH) {
            for (n_Ethnicity in 1:J_Ethnicity) {
              for (n_CD in 1:J_CD) {
                for (n_WhiteNonGrad in 1:J_WhiteNonGrad) {
                  for (n_Sex in 1:J_Sex) {
                    real pWI_ACS_WNH_State = inv_logit(alpha + Q_CD_ast[n_CD] * thetaX_CD + {eps_Sex,-eps_Sex}[n_Sex] + beta_Race[n_Race] + {eps_Ethnicity,-eps_Ethnicity}[n_Ethnicity] + {eps_Age,-eps_Age}[n_Age] + {eps_Education,-eps_Education}[n_Education] + {eps_WhiteNonGrad,-eps_WhiteNonGrad}[n_WhiteNonGrad] + beta_State[n] + {eps_WNH_State[n],-eps_WNH_State[n]}[n_WNH]);
                    WI_ACS_WNH_State_WgtSum += WI_ACS_WNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                    WI_ACS_WNH_State[n] += pWI_ACS_WNH_State * WI_ACS_WNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                  }
                }
              }
            }
          }
        }
      }
    }
    WI_ACS_WNH_State[n] /= WI_ACS_WNH_State_WgtSum;
  }
  vector[N_WI_ACS_NWNH_State] WI_ACS_NWNH_State = rep_vector(0,N_WI_ACS_NWNH_State);
  for (n in 1:N_WI_ACS_NWNH_State) {
    real WI_ACS_NWNH_State_WgtSum = 0;
    for (n_Education in 1:J_Education) {
      for (n_Race in 1:J_Race) {
        for (n_Age in 1:J_Age) {
          for (n_WNH in 1:J_WNH) {
            for (n_Ethnicity in 1:J_Ethnicity) {
              for (n_CD in 1:J_CD) {
                for (n_WhiteNonGrad in 1:J_WhiteNonGrad) {
                  for (n_Sex in 1:J_Sex) {
                    real pWI_ACS_NWNH_State = inv_logit(alpha + Q_CD_ast[n_CD] * thetaX_CD + {eps_Sex,-eps_Sex}[n_Sex] + beta_Race[n_Race] + {eps_Ethnicity,-eps_Ethnicity}[n_Ethnicity] + {eps_Age,-eps_Age}[n_Age] + {eps_Education,-eps_Education}[n_Education] + {eps_WhiteNonGrad,-eps_WhiteNonGrad}[n_WhiteNonGrad] + beta_State[n] + {eps_WNH_State[n],-eps_WNH_State[n]}[n_WNH]);
                    WI_ACS_NWNH_State_WgtSum += WI_ACS_NWNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                    WI_ACS_NWNH_State[n] += pWI_ACS_NWNH_State * WI_ACS_NWNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                  }
                }
              }
            }
          }
        }
      }
    }
    WI_ACS_NWNH_State[n] /= WI_ACS_NWNH_State_WgtSum;
  }
  vector[N_NI_ACS_WNH_State] NI_ACS_WNH_State = rep_vector(0,N_NI_ACS_WNH_State);
  for (n in 1:N_NI_ACS_WNH_State) {
    real NI_ACS_WNH_State_WgtSum = 0;
    for (n_Education in 1:J_Education) {
      for (n_Race in 1:J_Race) {
        for (n_Age in 1:J_Age) {
          for (n_WNH in 1:J_WNH) {
            for (n_Ethnicity in 1:J_Ethnicity) {
              for (n_CD in 1:J_CD) {
                for (n_WhiteNonGrad in 1:J_WhiteNonGrad) {
                  for (n_Sex in 1:J_Sex) {
                    real pNI_ACS_WNH_State = inv_logit(alpha + Q_CD_ast[n_CD] * thetaX_CD + {eps_Sex,-eps_Sex}[n_Sex] + beta_Race[n_Race] + {eps_Ethnicity,-eps_Ethnicity}[n_Ethnicity] + {eps_Age,-eps_Age}[n_Age] + {eps_Education,-eps_Education}[n_Education] + {eps_WhiteNonGrad,-eps_WhiteNonGrad}[n_WhiteNonGrad] + beta_State[n]);
                    NI_ACS_WNH_State_WgtSum += NI_ACS_WNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                    NI_ACS_WNH_State[n] += pNI_ACS_WNH_State * NI_ACS_WNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                  }
                }
              }
            }
          }
        }
      }
    }
    NI_ACS_WNH_State[n] /= NI_ACS_WNH_State_WgtSum;
  }
  vector[N_NI_ACS_NWNH_State] NI_ACS_NWNH_State = rep_vector(0,N_NI_ACS_NWNH_State);
  for (n in 1:N_NI_ACS_NWNH_State) {
    real NI_ACS_NWNH_State_WgtSum = 0;
    for (n_Education in 1:J_Education) {
      for (n_Race in 1:J_Race) {
        for (n_Age in 1:J_Age) {
          for (n_WNH in 1:J_WNH) {
            for (n_Ethnicity in 1:J_Ethnicity) {
              for (n_CD in 1:J_CD) {
                for (n_WhiteNonGrad in 1:J_WhiteNonGrad) {
                  for (n_Sex in 1:J_Sex) {
                    real pNI_ACS_NWNH_State = inv_logit(alpha + Q_CD_ast[n_CD] * thetaX_CD + {eps_Sex,-eps_Sex}[n_Sex] + beta_Race[n_Race] + {eps_Ethnicity,-eps_Ethnicity}[n_Ethnicity] + {eps_Age,-eps_Age}[n_Age] + {eps_Education,-eps_Education}[n_Education] + {eps_WhiteNonGrad,-eps_WhiteNonGrad}[n_WhiteNonGrad] + beta_State[n]);
                    NI_ACS_NWNH_State_WgtSum += NI_ACS_NWNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                    NI_ACS_NWNH_State[n] += pNI_ACS_NWNH_State * NI_ACS_NWNH_State_wgts[n][n_Education, n_Race, n_Age, n_WNH, n_Ethnicity, n_CD, n_WhiteNonGrad, n_Sex];
                  }
                }
              }
            }
          }
        }
      }
    }
    NI_ACS_NWNH_State[n] /= NI_ACS_NWNH_State_WgtSum;
  }
  vector[N_WI_ACS_WNH_State] rtDiffWI = WI_ACS_WNH_State - WI_ACS_NWNH_State;
  vector[N_WI_ACS_WNH_State] rtDiffNI = NI_ACS_WNH_State - NI_ACS_NWNH_State;
  vector[N_WI_ACS_WNH_State] rtDiffI = rtDiffWI - rtDiffNI;
}
4

icic so yeah in the Stan compiler right now we limit it to 8 I’m pretty sure (fyi what you have above is an array of 8 dimensions). Are you getting an error when the stanc compiler compiles this or when the C++ compiles?

I’m pretty sure 8 was just an arbitrary number in the compiler. To lookup whether a function the user calls is accepted by the C++ we generate a big table of all the supported function signatures and I’m pretty sure the limit there is 8 just because the table gets to be pretty big after that. But if there’s a need for it I think we could increase that to 9 or 10.

I understand another compiler is generting this Stan code, though it may be good to talk to them as I think a more efficient form of this would just be a matrix with N_* rows and J_* columns

5

Thanks! I wrote the Haskell → Stan compiler–not nearly a complete one at this point, but helpful for managing the data and model bits all together. Anyway, I can re-write it to do as you suggest. Do you mean an index matrix, one with a row for every post-stratification cell and columns for the value of each group for that cell? So, in my case N x 8? Then a model term, e.g., for the State intercept, would look something like beta_State[PS_Indexes[n, Col_State]] ?

Is there some reason that would preferable to just generating a 1-D index array for each group the way I do for the modeled data? Thus leading to something more like beta_State[PS_State[n]] ?
I’d have to add each index to the data block but I wouldn’t need variables like Col_State to specify which column of the index matrix corresponds to each group.