Creating the design matrices for a canonical correlation analysis
Here you want to move away from for loops for two reasons:


your number of for loops depends on your number of datasets, so using for loops prevents you from generalizing your code to any number of datasets.
many for loops will likely slow down your code execution when working with a larger number of datasets.


I think the key to vectorizing your for loops is to use the expand.grid function. If you have

n &lt;- 4


datasets, then you have

p &lt;- n * (n - 1) / 2    # 6


degrees of freedom (the number of for loops in your code, or the number of items on the lower triangle of each matrix). If for each of these you can pick among

w &lt;- seq(from = 0, to = 1, length.out = n)


then you can build the matrix of all possible combinations by doing:

W &lt;- as.matrix(expand.grid(rep(list(w), p)))


Here W is a big matrix with 4096 rows, each row representing a different combination of your (i1, i2, i3, i4, i5, i6) variables:

&gt; head(W)
          Var1      Var2 Var3 Var4 Var5 Var6
[1,] 0.0000000 0.0000000    0    0    0    0
[2,] 0.3333333 0.0000000    0    0    0    0
[3,] 0.6666667 0.0000000    0    0    0    0
[4,] 1.0000000 0.0000000    0    0    0    0
[5,] 0.0000000 0.3333333    0    0    0    0
[6,] 0.3333333 0.3333333    0    0    0    0


These 6 columns are only part of the n * n = 16 values needed in each matrix. We can expand using the following:

X &lt;- matrix(1:(n*n), n, n) # pattern matrix of indices
A &lt;- matrix(0, nrow(W), n * n)
A[,    X[lower.tri(X)]] &lt;- W
A[, t(X)[lower.tri(X)]] &lt;- W


A is similar to W in that it is a matrix with 4096 rows, but each row now has the n * n = 16 values of a symmetric matrix.

From there, you can reshape A into a 3D array:

dim(A) &lt;- c(nrow(W), n, n)


and your 4096 matrices can be accessed as follows:

A[1, , ]
#      [,1] [,2] [,3] [,4]
# [1,]    0    0    0    0
# [2,]    0    0    0    0
# [3,]    0    0    0    0
# [4,]    0    0    0    0

A[10, , ]
#           [,1]      [,2]      [,3] [,4]
# [1,] 0.0000000 0.3333333 0.6666667    0
# [2,] 0.3333333 0.0000000 0.0000000    0
# [3,] 0.6666667 0.0000000 0.0000000    0
# [4,] 0.0000000 0.0000000 0.0000000    0

A[4096, , ]
#      [,1] [,2] [,3] [,4]
# [1,]    0    1    1    1
# [2,]    1    0    1    1
# [3,]    1    1    0    1
# [4,]    1    1    1    0


If I were you I would probably stop here, i.e., keep the data in this form. A 3d array might allow you to continue writing vectorized code if the rest of your analysis allows for it. However if you absolutely want a list of matrices, you can do:

C_list2 &lt;- lapply(seq(nrow(A)), function(i) A[i, , ])


(note that the order of the matrices in my data and yours do not match. Let me know if this is a concern, it's probably a matter of reorganizing the rows and/or columns of the W matrix.)