-
Notifications
You must be signed in to change notification settings - Fork 7
TridiagonalConjugation #194
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
Conversation
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #194 +/- ##
==========================================
+ Coverage 85.26% 87.05% +1.78%
==========================================
Files 13 11 -2
Lines 1595 1367 -228
==========================================
- Hits 1360 1190 -170
+ Misses 235 177 -58 ☔ View full report in Codecov by Sentry. |
The If there's always this better way to handle weighted derivatives (meaning differentiating |
OK that makes sense. It is possible to recover the other version as follows: julia> BidiagonalConjugation(R0', D0', R1', :L)'
ℵ₀×ℵ₀ Bidiagonal{Float64, InfiniteLinearAlgebra.BidiagonalConjugationBand{Float64}} with indices OneToInf()×OneToInf():
0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ …
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0 2.0 ⋅ …
⋮ ⋮ ⋮ ⋱ Though I wonder if the default should be swapped. |
The code is very slow though... will have to debug this: julia> n = 100_000; @time layout_getindex(B,1:n,1:n);
0.961080 seconds (16.18 M allocations: 521.940 MiB, 6.31% gc time)
julia> n = 100_000; @time layout_getindex(B,1:n,1:n);
0.001608 seconds (6 allocations: 1.526 MiB) |
Oh wait, that is |
Actually it was just slow because my EDIT: oh whoops it is still slow with a faster R0 |
If
Where would we avoid |
The problem with eg. |
Here's an example of the speed difference:
|
Ah I didn't think about the branching - I thought maybe branch prediction would help. My benchmarking shows around a 4x speedup with InfiniteLinearAlgebra.jl/src/banded/bidiagonalconjugation.jl Lines 67 to 75 in 0155e48
it could be better to not loop over each entry individually, using InfiniteLinearAlgebra.jl/src/banded/bidiagonalconjugation.jl Lines 39 to 61 in 0155e48
U and C and then loop over the allocated vectors?
|
Let me implement my own tridiagonal variant in a way that is fast and we can see how to revise the bidiagonal one accordingly. |
@DanielVandH I started looking at this but
BidiagonalConjugation
seems wrong to me. It is computinginv(U)XV
but I think we actually needU*X*inv(V)
.For example, consider differentiating Chebyshev U. We have banded conversions
T == U*R0
andU == C*R1
, and known derivativediff(T) == U*D0
. Thus we haveso we want to represent
R1*D0*inv(R0)
whereR1
andR0
are banded.Here is a code example:
Do you agree? Or there is a use case for
inv(U)XV
?