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
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]