2
$\begingroup$

Update: you can now skip the "My attempted solution" part and jump to "My first attempt for improvement" in the end of the question.

My question is related to Best way to create symmetric matrices . The difference is that the (large) square matrix I want to populate has higher symmetry (and as in the mentioned question I want to use these symmetries to not evaluate every single element while populating the matrix, a Gram matrix in my case).

Besides being symmetric, my matrix has block structure: all the submatrices are also symmetric with zeros on their diagonals and bands along the diagonal with equal elements. Example of a 6x6=(2x3)x(2x3) matrix (consisting of submatrices of dimension 3):

enter image description here

Hence the set of all independent elements of any submatrix can be taken for example as the set of all elements on its first row (except for the first one, because it's a diagonal element and therefore zero). So, if $N_1$ is the number of submatrices in each row/column and $N_2$ the dimension of each (square) submatrix, instead of calculating $(N_1\times N_2)^2$ elements, I would like to:

A) evaluate $(N_2-1)N_1(N_1+1)/2$ elements (each of $N_1(N_1+1)$ independent submatrices has $N_2-1$ independent elements, see above)

B) efficiently populate the rest of the matrix using elements from A)


My attempted solution: (now that I added to the end of this question a much shorter, although definitely not yet optimal code, you may skip this lengthy, though working solution)

1) First, define matrix elements (which could be in principle used to determine elements of the whole matrix - later it will be a Gram matrix)

element[{i1_, i2_}, {j1_, j2_}] := Subscript[x, {i1, i2}, {j1, j2}]; 

where indexes $i_1,j_1$ run from 1 to $N_1$ (number of submatrices in each row/column of the large matrix) and $i_2,j_2$ run from 1 to $N_2$ (dimension of each submatrix).

2) Next, we populate the superdiagonal elements of the large matrix. We (a) populate each subdiagonal with zeros, (b) evaluate an independent set of $N_2-1$ elements (which determine the rest of the respective submatrix); if $i_1\neq j_1$, i.e. we are not on the main diagonal, we also symmetrize by occasion (c) populate the rest of the upper triangle of the respective submatrix (d) in the submatrices on the main diagonal only superdiagonal elements are calculated (for eventually we symmetrize the whole matrix)

Do[
  If[i1 != j1,
    Do[ g[{i1, i2}, {j1, i2}] = 0, {i2, 1, N2}
    ]; (*(a)*)
    Do[ g[{i1, 1}, {j1, j2}] = g[{i1, j2}, {j1, 1}] 
       = element[{i1, 1}, {j1, j2}], {j2, 2, N2}
    ]; (*(b)*)
    Do[ g[{i1, i2}, {j1, i2 + k}] = g[{i1, i2 + k}, {j1, i2}] 
       = g[{i1, 1}, {j1, k + 1}], {i2, 2, N2 - 1}, {k, 1, N2 - i2}
    ], (*(c)*)

    Do[ g[{i1, 1}, {j1, j2}] = element[{i1, 1}, {j1, j2}], {j2, 2, N2}
    ]; 
    Do[ g[{i1, i2}, {j1, i2 + k}] = g[{i1, 1}, {j1, k + 1}], {i2, 2, N2 - 1}, {k, 1, N2 - i2}
    ] (*(d)*)
 ],
{i1, 1, N1}, {j1, i1, N1}
];

3) That's all the superdiagonal terms, now we supplement the subdiagonal triangle (it again depends whether we are located in a submatrix on the main diagonal or not):

Do[
  If[i1 != j1,
   Do[
    g[{i1, i2}, {j1, j2}] = g[{j1, i2}, {i1, j2}], {i2, 1, N2}, {j2,1, N2}
   ],
   Do[ g[{i1, i2}, {j1, j2}] = g[{i1, j2}, {j1, i2}], {i2, 1, N2}, {j2,1,i2 - 1}
   ];
  ], {i1, 1, N1}, {j1, 1, i1}
];

4) We can add zeros to the main diagonal of the resulting matrix (but eventually I will add another large, i.e. $(N_1N_2)x(N_1N_2)$ diagonal matrix to it)

Do[g[{i, j}, {i, j}] = 0, {i, 1, N1}, {j, 1, N2}];

5) To operate with the result, I introduce a list with all possible combinations of indexes:

indexes = Flatten[Table[{i, j}, {i, 1, N1}, {j, 1, N2}], 1];
dimension = Length[indexes];

We notice that "dimension" is the dimension of the large matrix. For some reason I prefer to work with "indexes", but perhaps it is not completely necessary. Also, maybe I could somehow use this notation in the code.

6) The code is quite messy but gives the desired result (see the matrix above for $N_1=2,N_2=3$).

(result=Table[ g[indexes[[i2]], indexes[[j2]]], {i2, 1, dimension}, {j2, 1,dimension}]) // MatrixForm

So, as I said it works but it is very unclean and there is certainly a lot of space for improvement (how to make it more efficient, how to treat blocks on the main diagonal, maybe use Band[] ...) As a beginner (well, as you can see not only in Mathematica but programming in general), I am still not very familiar with all the @,@@,&,#,&@,#~/; :-) and other stuff which makes the code more compressed (I am sure somebody will write an answer with 3 lines of code) so it's not quite easy for me to follow and implement the advice from Best way to create symmetric matrices but shoot, I will try to digest it. Oh, and if you have suggestions how to improve my question so as it is more clear, which tags to add,... please let me know, I'll try to do that. Thank you all in advance, it's great to learn from such a great community!

To summarize - How to efficiently create a matrix with the above mentioned symmetries?

The next thing will be to calculate the eigenvalues and eigenvectors of this matrix. Is there an option of Eigensystem to take such symmetries into account?


Update: My first attempt for improvement based on the suggestion of @Mike Honeychurch to use SparseArray:

element[{i1_, i2_}, {j1_, j2_}] := Subscript[x, {i1, i2}, {j1, j2}];

Do[
  pattern = Join[Table[Band[{1, i}] -> element[{i1, 1}, {j1, i}], {i, 2, N2}], 
                 Table[Band[{i, 1}] -> element[{i1, 1}, {j1, i}], {i, 2, N2}]];
  sa[i1, j1] = SparseArray[pattern, {N2, N2}], 
{i1, 1, N1}, {j1, 1, N1}
];

Where "element..." is the definition of the matrix element and the Do loop creates the desired matrix (pattern tells SparseArray how to create each submatrix sa[i1,j1]). Finally, we put everything together:

ArrayFlatten[Array[sa, {N1, N1}], 2] // MatrixForm

Which gives for $N_1=2$ and $N_2=3$:

enter image description here

There are two issues (you can see that the matrix is different from the one above).

1) As you can see in the code, in the "pattern..." the definition of the matrix element is used two times (the submatrices are symmetric, so it is not necessary).

2) In the subdiagonal terms, the definition of the matrix element is used unnecessarily yet another two times (the whole matrix is symmetric).

Some solutions came to my mind, sure, but they seem to complicate the code too much - ideally, I would like to use the definition of the matrix element only once for each independent matrix element (note the mentioned symmetry) and simplify the code at the same time:-) Any suggestions how to make it elegantly, shrink the code, make it faster?

$\endgroup$
  • 1
    $\begingroup$ Have you had a look at the SparseArray docs? $\endgroup$ – Mike Honeychurch Dec 2 '14 at 1:50
  • $\begingroup$ @MikeHoneychurch Thanks for the suggestion, could you please check my try above? $\endgroup$ – wondering Dec 2 '14 at 23:46
  • $\begingroup$ I suggested SparseArray under the assumption that you had some sort of algorithm for determining the indices. $\endgroup$ – Mike Honeychurch Dec 3 '14 at 1:04
  • $\begingroup$ @MikeHoneychurch Well, this is a kind of algorithm how to go through the independent elements, isn't it? Or what is your point? Thanks! $\endgroup$ – wondering Dec 3 '14 at 2:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.