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

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 .

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

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

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