2
votes

I wish to expand matrix in two ways, but I can't figure out the code for it; I don't seek exact code but simply a pointer where to look in the literature. Problem below concerns linearization of commutator equation in quantum mechanics, I know theory from physics point of view, but I don't know the name of my problem in the programming world.

I have matrix A which is NxN matrix, for this example let N=3; So A is:

A=[a11 a12 a13
   a21 a22 a23
   a31 a32 a33] 

I wish to make a matrix B which is N^2xN^2 size. B has block form:

B=[A11 A12 A13
   A21 A22 A23
   A31 A32 A33]   

A11=[a11 0 0
     0 a11 0    ...
     0 0   a11]

each block element in B is NxN matrix and every block is linked by property A11=a11*eye(3,3), and so on (so (B(I,J)=A(i,j)*eye(3,3))).

My problem is that I have no idea how to construct such a loop, making temp 3x3 matrix from each element from A is easy, but I don't know how to stack these blocks in B (how to index positions in B).

Another matrix I'd like to form is also N^2xN^2 matrix which is much simpler:

C=[A 0 0
   0 A 0
   0 0 A]

it's block diagonal matrix made out of A. I think I can manage to find code for this (there is blkdiag function in matlab for it).

I'd like idea for solution in C++ and matlab. Problem is how to stack blocks into N^2xN^2 matrix and how many loops I need to use. I'd like the best computing approach for this.

My idea was to do something like this:

for i=1:N
    for j=1:N
        for k=1:N
          B(i*N+j-N,j+k*N-N)=A(i,k);
          C(i*N+j-N,k+i*N-N)=A(k,j);
        end
    end
end

This works perfectly (in C++ it would be without -N in B(?,?) and C(?,?) indexing)), but is it possible to do it with 2 for loops?

1

1 Answers

2
votes

I would suggest to use the "Kronecker" product, also known as "tensor product". A matlab function is already available "kron". You simply have to built the the matrix :

A=[a11 a12 a13;
   a21 a22 a23;
   a31 a32 a33] ;

Then do the kron product

B = kron(A,eye(3));

Therefore, "B" is the result you seek, and there is no need for any loop.

For the value of "C", you simply need to switch the 2 inputs of the "kron" function

C = kron(eye(3), A);