A question on broadcasting fft for vectors of matrices

Specific Domains Numerics
, ,
1

In the following I am creating a 3-vector of 2x2 matrices M:

M = [ rand(2,2) for i in 1:3]

Is it any easy compact way to extract/view (e.g., for the purpose of broadcasting), say, the vector of [1,1] elements of the three matrices ?

My solution was to build a 3d matrix and use view:

view(reshape(hcat(M...),2,2,3),1,1,:)

but I am afraid this is quite inefficient for large matrices.

This problem I encountered trying to apply the fft elementwise to each vector of corresponding entries, but


fft.(M)

produces the wrong result. I wonder actually what?

Many thanks in advance for any help.

2

Not sure whether maybe your MWE is too simple but:

julia> M = [ rand(2,2) for i in 1:3]
3-element Vector{Matrix{Float64}}:
 [0.798300477803938 0.7394167094533922; 0.020090132086653556 0.8809955596136937]
 [0.486959574050713 0.9602589802986382; 0.18396003854515464 0.8687833148070098]
 [0.5894287539978448 0.9398716351164952; 0.8376087279393974 0.23575124677607961]

julia> view(reshape(hcat(M...),2,2,3),1,1,:)
3-element view(::Array{Float64, 3}, 1, 1, :) with eltype Float64:
 0.798300477803938
 0.486959574050713
 0.5894287539978448

julia> first.(M)
3-element Vector{Float64}:
 0.798300477803938
 0.486959574050713
 0.5894287539978448
1 Like
3

Please check if this transformation suits you:

M = [rand(2,2) for i in 1:3]

using TensorCast
@cast T[i,j][k] := M[k][i,j]    # := returns a view

using FFTW
fft.(T)
Result 
julia> fft.(T)
2×2 Matrix{Vector{ComplexF64}}:
 [0.676581+0.0im, 0.266587+0.162225im, 0.266587-0.162225im]  …  [1.30509+0.0im, 0.728171-0.156387im, 0.728171+0.156387im]
 [1.45385+0.0im, -0.3766+0.38125im, -0.3766-0.38125im]          [2.19093+0.0im, 0.018722-0.186763im, 0.018722+0.186763im]
1 Like
5

Isn’t this just

getindex.(M, 1)

? (basically the same as first.(M))

For more general indices

getindex.(M, CartesianIndex(1,1))
1 Like
6

Actually I have to make similar selections for all entries. Thanks for the nice solution.

1 Like
7

Many thanks for this idea. And yes,

getindex.(M,i,j)

does exactly what I needed for any indices.

However, on Julia1 1.7.1, the following fails (not sure I can understand the reason)

 getindex.(M, CartesianIndex(1,1))
ERROR: iteration is deliberately unsupported for CartesianIndex. Use `I` rather than `I...`, or use `Tuple(I)...`
8

I will check if this is faster than cycling over all indices with
getindex.(M,i,j). Many thanks.

9

The original question was about how to broadcast fft() on a vector of matrices though, you may wanna change the post title accordingly.

1 Like
10

Try this:

getindex.(M, (CartesianIndex(1,1),))

NB:
And my comment on changing the title of the OP to match what was needed is that it seems to be about how to access a subset of a vector of matrices as specified by some indices, and not so much about broadcasting fft().

1 Like
11

That is surprising. Not that one cannot iterate, I knew that, but expected it to therefore to be treated as a scalar.

1 Like