JuliaApproximation / SingularIntegrals.jl

A Julia package for computing singular integrals
MIT License
6 stars 4 forks source link

Performance of RecurrenceArray #17

Open ioannisPApapadopoulos opened 1 year ago

ioannisPApapadopoulos commented 1 year ago

Sometimes RecurrenceArray is slower than getindex when actually assembling the values e.g.

using ClassicalOrthogonalPolynomials
using SingularIntegrals

import ClassicalOrthogonalPolynomials: recurrenecoefficients
import SingularIntegrals: RecurrenceArray

P = Jacobi(1/3,1/3)
xc = range(-1,1; length=10000)
n = 10000

# getindex baseline ~ 1.102867 seconds
@time pp = P[xc,1:n];
# Lazy is fast ~  0.000466 seconds
@time pr = RecurrenceArray(xc, recurrencecoefficients(P), P[xc, 1:2]');
# Actually getting the matrix is "slow" ~ 1.706700 seconds 
@time pr[1:n,1:length(xc)]';

pr ≈ pp
ioannisPApapadopoulos commented 1 year ago

@dlfivefifty just in case you missed this issue you asked me to raise :)

dlfivefifty commented 1 year ago

Your code doesn't run since is not defined

ioannisPApapadopoulos commented 1 year ago

ah sorry! fixed. It's particularly slow on the first run of @time pr[1:n,1:length(xc)]';

dlfivefifty commented 1 year ago

Note 0.6s is just making a second copy. That is, when you call pr[1:n,1:length(xc)] it first populates a matrix and then copies this data to a new matrix.

On the other hand, P[xc,1:n] populates the returned matrix directly.

So I don't think there's anything obvious that can be done about this. That is, we can't both cache and return a matrix without the overhead of a copy.

dlfivefifty commented 1 year ago

Unless we decide to get rid of caching and only do populating. That is, RecurrenceArray would become completely lazy and not cache any computations. But is it worth the time?

ioannisPApapadopoulos commented 1 year ago

For the paper, I've not used the RecurrenceArrays at all.. probably not worth it for now.