Transformation invariant subspace of three-component array

Question

I have a large sparse array A with dimensions $m \times n \times q $. I also have two transformation matrices, B ($n \times n$) and C ($q \times q$) from which I act from the left on A. So, for instance, I might use

TensorTranspose[
  B . TensorTranspose[
        C . TensorTranspose[A, {2, 3, 1}],
        {2, 3, 1}
      ],
  {2, 3, 1}
]

where the final tensor transpose is there to put the indices back in the correct order at the end.

Now, A has the property such that $m$ 'rows' (i.e. if I flatten A into an m by n*q matrix) of the transformed matrix are always linear combinations of the original rows. I want to find the subspace of the space spanned by the m 'rows' of the original matrix which are actually invariant under the transformation which I have described above.

My procedure has been as follows:

1. Perform the transformation as described above

2. Flatten the original A and the transformed A into m by n*q matrices, and use LinearSolve to identify the image of each transformed row as a linear combination of the original rows.

3. Form an m by m matrix X using the results of 2, which acts on A from the left and captures the original combined action of B and C on A.

4. Take Z = NullSpace[X - IdentityMatrix[{m,m}]]; this acts as a projector so that Z.A gives the subspace I wanted.

The problem is that the use of LinearSolve in 2) gets extremely slow (these matrices are generally very large); I feel it should be possible to carry out this process directly using just the matrices B and C, since they themselves are sufficient to transform A, but I can't wrap my head around the required linear algebra.

Answer 1

I do not fully understand what you mean by "the m rows". Here is as far as I can follow you:

You want to calculate:

This can be done using Transpose and Do. As in MMA you should not use capitalized variable names, I am using ma, mb and mc. Note that the Dot product contracts the last index of the first tensor with the first index of the second tensor. Therefore, in a first step we exchange the first and second index of ma to multiply B with A. In a second step we then exchange index 1 and 3 of the result to get C times B times A. Here is an explicit example with m==n==q==2:

The matrices are:

(ma = Array[Subscript[a, #1, #2, #3] &, {2, 2, 2}]) // MatrixForm
(mb = Array[Subscript[b, #1, #2] &, {2, 2}]) // MatrixForm
(mc = Array[Subscript[c, #1, #2] &, {2, 2}]) // MatrixForm

Now we make the different multiplications:

(tmp = mb . Transpose[ma, {3, 2, 1}]) // MatrixForm
(tmp = mc . Transpose[tmp, {2, 1, 3}]) // MatrixForm

Finally we must restore the original order of the indices:

(tmp = Transpose[tmp, {3, 1, 2}]) // MatrixForm

If we now flatten the first level of ma and the result, we get:

Flatten[ma, 1]
Flatten[tmp, 1]

Or if we reshape them as:

ArrayReshape[ma, {2, 2}]
ArrayReshape[tmp, {2, 2}]

Unfortunately, I do not understand what should now be a linear combination of what. Please clarify this.

Answer 2

I found a solution to this in the end.

Transform A using B and C, and then flatten both A and the transformed A into m by (n*q) matrices. Applying NullSpace to this matrix gives a projector of A onto the invariant subspace.

To transform A, use

ATrans = TensorTranspose[C.TensorTranspose[B.TensorTranspose[A,{3,2,1}],{3,2,1}],{3,2,1}].

Then, flatten both A and ATrans in their second two arguments:

AFlat = ArrayReshape[A,{m,n*q}]
ATransFlat = ArrayReshape[ATrans,{m,n*q}]

Then take the null space of their difference to get the projector onto the invariant subspace.

projector = NullSpace[AFlat-ATransFlat]

Dot this with A to get the invariant subspace.

projector.A