Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inner product between bases on different grids, basis transforms #141

Open
jagot opened this issue Feb 4, 2023 · 1 comment
Open

Inner product between bases on different grids, basis transforms #141

jagot opened this issue Feb 4, 2023 · 1 comment

Comments

@jagot
Copy link
Member

jagot commented Feb 4, 2023

Say A isa Basis and B isa Basis, where axes(A,1) != axes(B,1), then we get the following error:

julia> A'B
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match

However, if I have a function $f$ that is expanded over the basis A, $|f\rangle = \sum_i c_i|a_i\rangle$, and I want to project it onto the basis B, $|f\rangle = \sum_j d_i|b_j\rangle$, then the above basis overlap is precisely the one I need, since $d_j = \sum_i c_i \langle b_j|a_i\rangle$, or equivalently $\mathbf{d}=A^HB\mathbf{c}$.

To me it is unimportant that the axes are different, i.e. I think it is up to the user to make sure that the end result makes sense (e.g. projecting onto a "larger" basis, or maybe I am interested in the subspace, etc).

I can of course override this on a case-by-case basis, but it would be convenient if this just worked, and the only thing you'd need to provide is the calculation of basis function overlaps $\langle b_j|a_i\rangle$. Then it would be very easy to transform between e.g. B-splines and finite-elements or -differences.

What do you think?

EDIT: It seems it will be difficult to work around this on a case-by-case basis, since I cannot circumvent the dimension check easily:

julia> ApplyQuasiArray(*, A', B)
ERROR: DimensionMismatch: Second axis of A, Inclusion(0.0..2.0), and first axis of B, Inclusion(0.0..2.5) must match
Stacktrace:
 [1] _check_mul_axes(A::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, B::QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false})
   @ ArrayLayouts ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:91
 [2] check_mul_axes
   @ ~/.julia/packages/ArrayLayouts/4IG3b/src/mul.jl:107 [inlined]
 [3] check_applied_axes
   @ ~/.julia/packages/LazyArrays/NYra8/src/linalg/mul.jl:22 [inlined]
 [4] instantiate
   @ ~/.julia/packages/LazyArrays/NYra8/src/lazyapplying.jl:63 [inlined]
 [5] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:53 [inlined]
 [6] ApplyQuasiArray
   @ ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:54 [inlined]
 [7] ApplyQuasiArray(M::LazyArrays.Applied{LazyArrays.MulStyle, typeof(*), Tuple{QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, QuasiArrays.SubQuasiArray{Float64, 2, FEDVR{Float64, Float64, FillArrays.Fill{Int64, 1, Tuple{Base.OneTo{Int64}}}}, Tuple{Inclusion{Float64, IntervalSets.ClosedInterval{Float64}}, UnitRange{Int64}}, false}}})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:55
 [8] ApplyQuasiArray(::Function, ::QuasiArrays.QuasiAdjoint{Float64, BSpline{Float64, Float64, LinearKnotSet{7, 7, 7, Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, Vector{Float64}, SparseArrays.SparseMatrixCSC{Float64, Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}}, ::Vararg{Any})
   @ QuasiArrays ~/.julia/packages/QuasiArrays/tF3vb/src/lazyquasiarrays.jl:59
@jagot
Copy link
Member Author

jagot commented Feb 6, 2023

For reference, A'B is implemented as basis_overlap(A,B) in JuliaApproximation/CompactBases.jl#63 and looks like this in action: https://juliaapproximation.github.io/CompactBases.jl/previews/PR63/basis_transforms/

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant