Efficient calculation of block matrices with only the Linear Map given

Specific Domains
1

Hi!

I have the situation that I construct block matrices where I want to avoid to explicitly have them.
For example, with Fourier transforms we can either use slow matrices or fast FFTs:

using FFTW
using LinearAlgebra

z = zeros((2,2))
# my true vectors are 1_000_000 elements long
x1 = randn((2,))
x2 = randn((2,))

# FFT plan
p = plan_fft(x1)
p_mat = fft(I(2), (1,))

# slow mat mul
M = [p_mat z; z p_mat]
res1 = M * [x1; x2] 

# fast FFT
res2 = [p * x1; p * x2]

#### 
julia> res1 == res2
true

The issue is that my matrices are not necessarily block diagonals, it can be any kind of diagonal matrix. Also my matrix is not quadratic.

And finally I would like to invert the matrix M which obviously does not work straightforwardly since I can’t access matrix elements but only the forward pass M * x.
But maybe we can calculate the pseudoinverse because we should be able to calculate M' *x?

Is there any kind of package allowing to generate the block matrices and efficient calculations?
For example, it probably is perfect if I could write: (@efficient_block_matrix [p 0; 0 p]) * [x1; x2]

Best,

Felix

Existing packages:

2

Would something like Home · LinearMaps.jl help? There are also some types which look suitable: Types and methods · LinearMaps.jl
(Maybe there is a more suitable package, but with LinearMaps you can define yourself how the linear map is defined and it plays along well with other packages.)

1 Like
3

Thanks!
Yeah, I forgot to add LinearMaps.jl here.

Yes, it kind of helps but I still would need to implement the matching of the blocks with the vectors, wouldn’t I? (not super hard).

1 Like
4

I’m not sure if I understand your intended application, but are you looking for something like this type:

5

No, you don’t, that’s all included.

1 Like