How do you extract each column of a matrix as a column vector

General
1

Hi! Thank you for your help.Now I have some troubles in the matrix. For a matrix(p*n),I want to extract each column of the matrix as a column vector named as mi(i=1…n）, then calculate mi times transpose(mi), then we will get n matrices.After that ,I want to calculate the sum of n matrices.
Thank you for you help .

2

I am a novice with Stan so there may be far better ways to do this. If I understood you correctly, the R code below makes the calculation directly in R and then in Stan. I did the calculation in R so I could compare the results with what the Stan code produced. I included as comments in the Stan code the declaration and population of an array to hold each of the p x p matrices, if those are of interest.

library(rstan)
MAT <- matrix(1:20, nrow = 4)
print(MAT)
#        [,1] [,2] [,3] [,4] [,5]
# [1,]    1    5    9   13   17
# [2,]    2    6   10   14   18
# [3,]    3    7   11   15   19
# [4,]    4    8   12   16   20

#R code to do the same calcualtion
Tot <- matrix(rep(0, 16), nrow = 4)

for (i in 1:5){
  Tot <- Tot + (MAT[, i] %*% t(MAT[, i]))
}
print(Tot)
#       [,1] [,2] [,3] [,4]
# [1,]  565  610  655  700
# [2,]  610  660  710  760
# [3,]  655  710  765  820
# [4,]  700  760  820  880

MODEL <- "
data {
  int p; //number of rows
  int n; //number of columns
  matrix[p,n] M;
}
transformed data {
 // matrix[p,p] Mats[n];
  matrix[p,p] Total;
  vector[p] tmp;

  Total = rep_matrix(0, p, p);

  for ( i in 1:n) {
    for (j in 1:p) {
      tmp[j] = M[j,i];
    }
    //Mats[i] = tmp * tmp';
    Total += tmp * tmp';
  }

  print(Total);
}
model {

}
"
dataList <- list(p = nrow(MAT), n = ncol(MAT), M = MAT)
FIT <- stan(model_code = MODEL, data = dataList, iter = 1, chains = 1, 
            warmup = 0, algorithm = "Fixed_param")
# [[565,610,655,700],[610,660,710,760],[655,710,765,820],[700,760,820,880]]
3 Likes
3

Here is slightly improved Stan code using the col() function

MODEL <- "
data {
  int p; //number of rows
  int n; //number of columns
  matrix[p,n] M;
}
transformed data {
 // matrix[p,p] Mats[n];
  matrix[p,p] Total;
  vector[p] tmp;

  Total = rep_matrix(0, p, p);

  for ( i in 1:n) {
      tmp = col(M,i);
      //Mats[i] = tmp * tmp';
      Total += tmp * tmp';
  }

  print(Total);
}
"
2 Likes
4

Thank you very much ! It’s very useful for my working.