3
$\begingroup$

Hopefully this question makes sense.

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 new $8 \times 8$ matrix).

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

$\endgroup$
1
  • $\begingroup$ I have not read the already posted answer, and it looks sound, however, in short: use Kronecker product of your 1D matrices to form a 3D left-shift operator. You can also have a look to the construction of the (two-point) Finite Difference operator in one and then higher dimensions, which is basically a combination of left- and right-shift. $\endgroup$
    – davidhigh
    Commented Dec 20, 2023 at 11:43

1 Answer 1

5
$\begingroup$

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_{kl}$$ 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.]])
$\endgroup$
3
  • 1
    $\begingroup$ Amazing, thank you for this thorough response. I will review this and let you know if I have any questions. For reference, I am using Julia. $\endgroup$
    – Nukesub
    Commented Dec 20, 2023 at 4:14
  • $\begingroup$ I added some Julia code and a warning about order of elements $\endgroup$ Commented Dec 20, 2023 at 8:01
  • $\begingroup$ Great answer! I suggest adding a remark to address a common source of subtle bugs: in the identity $\operatorname{vec}(AXB) = (B^\top \otimes A) \operatorname{vec}(X)$, which you reference in a special case, that transposition is truly a transposition even when working with complex matrices: it does not become a conjugate transpose, unlike almost every other transpose that appears in linear algebra. $\endgroup$ Commented Jan 6 at 10:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.