JuliaSymbolics / Symbolics.jl

Symbolic programming for the next generation of numerical software
https://docs.sciml.ai/Symbolics/stable/
Other
1.35k stars 151 forks source link

Symbolic (symmetric-toepltiz)matrix-vector multiplication doesn't work for higher dimensions of matrix #939

Closed yewalenikhil65 closed 1 year ago

yewalenikhil65 commented 1 year ago

I am trying following , and this works okay!

using Symbolics, ToeplitzMatrices

N =250

a = Symbolics.variables(:a, 0:N);

b= [1; (1:N).*a[2:end]]

A  = SymmetricToeplitz(a)
F_expr = A*b;

However, the code fails if i increase the code to N=280


using Symbolics, ToeplitzMatrices

N =280

a = Symbolics.variables(:a, 0:N);

b= [1; (1:N).*a[2:end]]

A  = SymmetricToeplitz(a)
F_expr = A*b;

ERROR: MethodError: no method matching plan_fft!(::Vector{Complex{Num}}, ::UnitRange{Int64})

Closest candidates are:
  plan_fft!(::StridedArray{T, N}, ::Any; flags, timelimit, num_threads) where {T<:Union{ComplexF64, ComplexF32}, N}
   @ FFTW ~/.julia/packages/FFTW/HfEjB/src/fft.jl:786
  plan_fft!(::AbstractArray; kws...)
   @ AbstractFFTs ~/.julia/packages/AbstractFFTs/hJ0Fz/src/definitions.jl:65

Stacktrace:
 [1] plan_fft!(x::Vector{Complex{Num}}; kws::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ AbstractFFTs ~/.julia/packages/AbstractFFTs/hJ0Fz/src/definitions.jl:65
 [2] plan_fft!(x::Vector{Complex{Num}})
   @ AbstractFFTs ~/.julia/packages/AbstractFFTs/hJ0Fz/src/definitions.jl:65
 [3] factorize(A::SymmetricToeplitz{Num, Vector{Num}})
   @ ToeplitzMatrices ~/.julia/packages/ToeplitzMatrices/5gW1c/src/linearalgebra.jl:150
 [4] mul!(y::Vector{Num}, A::SymmetricToeplitz{Num, Vector{Num}}, x::Vector{Num}, α::Bool, β::Bool)
   @ ToeplitzMatrices ~/.julia/packages/ToeplitzMatrices/5gW1c/src/linearalgebra.jl:37
 [5] mul!
   @ ~/software/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/matmul.jl:276 [inlined]
 [6] *(A::SymmetricToeplitz{Num, Vector{Num}}, x::Vector{Num})
   @ LinearAlgebra ~/software/julia-1.9.2/share/julia/stdlib/v1.9/LinearAlgebra/src/matmul.jl:56
 [7] top-level scope
   @ REPL[36]:1
shashi commented 1 year ago

It seems the multiplication is switching to FFT algorithm for large convolutions. I'm not sure what we can do about that. What's your goal with this? We could potentially make fft a primitive in Symbolics. cc @ChrisRackauckas

xtalax commented 1 year ago

I fully support this idea, there's a lot of simplification you can do in Fourier space

ChrisRackauckas commented 1 year ago

A lot of pseudospectral PDE codes could be nicely represented with it. I think more generally the question is how to register array functions. If registration works well on that, we can start getting coverage.

xtalax commented 1 year ago

Definitely agree on this point, I've run up against this a few times. Could just get the registration methods to wrap on Union{Num, Arr}

yewalenikhil65 commented 1 year ago

It seems the multiplication is switching to FFT algorithm for large convolutions. I'm not sure what we can do about that. What's your goal with this? We could potentially make fft a primitive in Symbolics. cc @ChrisRackauckas

@shashi My aim was to create that F_expr vector in symbolic form as written in the code for arbitrary large N . It creates the set of quadratic equations I was interested in solving (they form system of equations related to deep gravity waves 🌊). Only thing is such matrix-vector multiplication wasn't possible due to error mentioned above for N >= 280 . I found a way around by helpful comment here by avoiding certain conditional statements

But it would be great if making fft a primitive in Symbolics could solve this issue. Certainly the said matrix-vector product is possible and user should be allowed to create them!

What looks overall from error trace is that plan_fft! method entertaining Vector{Num} should be included somewhere.

shashi commented 1 year ago

Array registration is a top priority for me right now.

@yewalenikhil65 even if we support fft node right now, you won't be able to access individual expressions after the convolution. It's best to request Toeplitz multiplication to add an option to multiplication where it does not take the path of calling fft.

ChrisRackauckas commented 1 year ago

Array registration is a top priority for me right now.

Yeah I think that's all you need. Once that exists, the community can register a whole standard library etc. and take it from there. But without that registration, these kinds of issues will just pop up. You don't need to worry about an fft node, just a registration process and I or anyone else will get to the fft node.

xtalax commented 1 year ago

This probably needs a broad request to package maintainers to register important functions in package