Apply 3D Operator to Matrix and get new Matrix

Question

I know I can formulate an operator for a vector as a matrix, then apply that matrix to my vector to get a new vector. For example, if I define a left shift operator which shifts all elements left I can write such a matrix as follows (I will denote as $L$). I've assumed periodic boundary conditions.

 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0

Then I could apply this to some $8\times 1$ vector, say $[1, 2, 3, 4, 5, 6, 7, 8]^T$ and get $[2, 3, 4, 5, 6, 7, 8, 1]^T$.

I'd like to understand how I can extend this up one dimension. In other words, I want to apply a 3D operator to a matrix and get a new matrix out. I think I can form the 3D left shift operator by using the Kronecker product.

So then I have my operator, which is of size $64 \times 64$ and equal to $H = L \otimes L$.

Let's say I have a matrix $A$ of size $8 \times 8$ with $1, 2, 3, \ldots, 8$ on the diagonal.

 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  3.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  4.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  5.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  6.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  7.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  8.0

Applying a left shift operator that only impacts the column-wise direction should give (again, I think?)

 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 2.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  3.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  4.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  5.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  6.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  7.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  8.0  0.0

I am confused though how to perform this operation (assuming my operator was constructed correctly). In the matrix-vector case, it is clear to me that you contract over the common index and you get a new vector back.

In this case, though, what are the common indices that should be contracted so the result is a new matrix? I suppose I could form $H$ into a $512 \times 8$ matrix and do matrix-matrix multiplication with $A$ but the result would be the wrong size (my expectation is result should be a new $8 \times 8$ matrix).

I may be overthinking this - I'd appreciate any intuition or tips on how to approach this.

Answer 1

I am not sure why you call $H$ a "3D operator", so I am not completely sure that I understood the question right, but here is my attempt at explanation.

Consider the space $\mathbb{R}^{n \times n}$ of $n \times n$ matrices. It is a vector space, and if we denote by $E_{ij}$ the matrix which has $1$ in the $i$-th row and the $j$-th column, then all $E_{ij}$ form a basis of $\mathbb{R}^{n \times n}$. The coordinates of a matrix $A = (a_{ij})$ in this basis are exactly the matrix elements $a_{ij}$: $$A = \sum_{1 \leq i,j \leq n} a_{ij} E_{ij}$$ If we want, we can forget that $A$ is a matrix, and instead see it as a vector of $n^2$ coordinates. If we order $ij$ lexicographically, this vector will contain all rows of our matrix one after the other. Numerically, this is a noop - we can view the same array of length $n^2$ as a $n \times n$ matrix, or as a $n^2$ vector, and a function like numpy.reshape will translate between these views.

Now consider a linear map $H \colon \mathbb{R}^{n \times n} \to \mathbb{R}^{n \times n}$. In the basis $E_{ij}$ it is given by a $n^2 \times n^2$ matrix $(h_{ij,kl})$ and the application of $H$ can be translated into matrix multiplication: if $B = H(A)$ with $A = (a_{ij})$ and $B = (b_{ij})$, then $$b_{ij} = \sum_{1 \leq k,l\leq n} h_{ij,kl} \, a_{kl}$$

Instead of looking at $H$ as a $n^2 \times n^2$ matrix and $A$ and $B$ as $n^2$ vector, we can remember that $A$ and $B$ are matrices, and in this case we interpret $H$ as an order $4$ tensor $H = (h_{ijkl})$ of format $n \times n \times n \times n$. The preceding equation is then tensor contraction along two indices $$b_{ij} = \sum_{k = 1}^n \sum_{l = 1}^n h_{ijkl} \, a_{kl}$$ Again, this is only a question of interpretation - the equations in coordinates are essentially the same. If we use numpy, then the first interpretation (matrix-vector multiplication) can be done with

B = numpy.reshape(H @ numpy.reshape(A, n*n), (n, n)),

and the second (tensor contraction) as

B = numpy.tensordot(H, A, [[2,3],[0,1]])

in the first case H must be of format $n^2 \times n^2$, and in the second case — $n \times n \times n \times n$.

Now let's look closely at Kronecker products and how to apply things to columns and rows. If we want to apply the same transformation $L$ to all columns, then this is just matrix multiplication, and $B = LA$ can be written as $$b_{ij} = \sum_{k = 1}^n l_{ik} a_{kj}$$ In order to write this equation in the form we have seen before, we can define $h_{ijkl} = l_{ik} \delta_{jl}$. Then $$\sum_{k = 1}^n \sum_{l = 1}^n h_{ijkl} \, a_{kl} = \sum_{k = 1}^n \sum_{l = 1}^n l_{ik} \delta_{jl} a_{kl} = \sum_{k = 1}^n l_{ik} a_{kj}$$ If we see $H$ as a matrix, then $h_{ij,kl} = l_{ik} \delta_{jl}$ defines exactly the Kronecker product $L \otimes I$. Similarly, if we take $H = I \otimes L$, then this is the same as applying $L$ to every row (transposed, so as matrix multiplication it is $B = AL^\top$), and $H = L \otimes L$ corresponds to $B = LAL^\top$ — applying $L$ to each column and each row.

UPD: When using functions like reshape it is important to be mindful about the order of elements. Julia uses column-major order, which means that reshape will return vector with columns of the original matrix, not rows. So if you use the reshape approach, the rows and columns will swap. The solution is to either transpose before reshape, or use the opposite order in the Kronecker product.

Here is an example using Julia (with $4 \times 4$ matrices instead of $8 \times 8$. I also added a non-diagonal entry).

julia> n = 4
4

julia> A = [1. 0 0 7.; 0 2. 0 0; 0 0 3. 0; 0 0 0 4.]
4×4 Matrix{Float64}:
 1.0  0.0  0.0  7.0
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0

julia> L = [0 1. 0 0; 0 0 1. 0; 0 0 0 1.; 1. 0 0 0]
4×4 Matrix{Float64}:
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0
 1.0  0.0  0.0  0.0

julia> L * A
4×4 Matrix{Float64}:
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0
 1.0  0.0  0.0  7.0

julia> using LinearAlgebra

julia> H = kron(I(n), L)
16×16 Matrix{Float64}:
 0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0

julia> Avec = reshape(A, n*n)
16-element Vector{Float64}:
 1.0
 0.0
 0.0
 0.0
 0.0
 2.0
 0.0
 0.0
 0.0
 0.0
 3.0
 0.0
 7.0
 0.0
 0.0
 4.0

julia> Bvec = H * Avec
16-element Vector{Float64}:
 0.0
 0.0
 0.0
 1.0
 2.0
 0.0
 0.0
 0.0
 0.0
 3.0
 0.0
 0.0
 0.0
 0.0
 4.0
 7.0

julia> B = reshape(Bvec, n, n)
4×4 Matrix{Float64}:
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0
 1.0  0.0  0.0  7.0

julia> Htensor = reshape(H, n, n, n, n)
4×4×4×4 Array{Float64, 4}:
[:, :, 1, 1] =
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 1.0  0.0  0.0  0.0

(...)

julia> using TensorOperations

julia> tensorcontract(Htensor, [0,1,2,3], A, [2,3])
4×4 Matrix{Float64}:
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0
 1.0  0.0  0.0  7.0

julia> B = zeros(n, n)
4×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> @tensor B[i,j] = Htensor[i,j,k,l]*A[k,l]
4×4 Matrix{Float64}:
 0.0  2.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  0.0  0.0  4.0
 1.0  0.0  0.0  7.0

Here is an example using numpy. Note different order in kronecker product.

>>> import numpy as np
>>> n = 4
>>> A = np.diag([1.,2.,3.,4.])
>>> A
array([[1., 0., 0., 0.],
       [0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.]])
>>> L = np.array([[0,1.,0,0],[0,0,1.,0],[0,0,0,1.],[1.,0,0,0]])
>>> L
array([[0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.],
       [1., 0., 0., 0.]])
>>> L @ A
array([[0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.],
       [1., 0., 0., 0.]])
>>> H = np.kron(L, np.eye(n))
>>> H
array([[0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
       [1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
>>> Avec = np.reshape(A, n*n)
>>> Avec
array([1., 0., 0., 0., 0., 2., 0., 0., 0., 0., 3., 0., 0., 0., 0., 4.])
>>> Bvec = H @ Avec
>>> Bvec
array([0., 2., 0., 0., 0., 0., 3., 0., 0., 0., 0., 4., 1., 0., 0., 0.])
>>> B = np.reshape(Bvec, (n,n))
>>> B
array([[0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.],
       [1., 0., 0., 0.]])
>>> Htensor = np.reshape(H, (n,n,n,n))
>>> B = np.tensordot(Htensor, A, [[2,3],[0,1]])
>>> B
array([[0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.],
       [1., 0., 0., 0.]])
>>> B = np.einsum('ijkl,kl', Htensor, A)
>>> B
array([[0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.],
       [1., 0., 0., 0.]])