[SciPy-Dev] Mass 1D interpolation, advice

[Jonathan Stickel]

On 4/15/15 09:24 , josef.pktd at gmail.com wrote:
> On Wed, Apr 15, 2015 at 11:04 AM,  <josef.pktd at gmail.com> wrote:
>> On Wed, Apr 15, 2015 at 10:33 AM, Jonathan Stickel <jjstickel at gmail.com> wrote:
>>> On 4/14/15 16:12 , Matthew Brett wrote:
>>>> Hi,
>>>>
>>>> In brain imaging (my world) - we need to do 1D interpolation of large
>>>> datasets, specifically, we have data in the order of (X, Y) = (200,
>>>> 4000), where we are interpolating at the same ~200 X locations for
>>>> each of the 4000 rows.
>>>>
>>>> We need to do this with extrapolation at the boundaries, so I think we
>>>> can't use scipy.interpolate.interp1d which is otherwise well suited
>>>> for the task.
>>>>
>>>> scipy.interpolate.InterpolatedUnivariateSpline does deal with
>>>> extrapolation, but only works with vectors.
>>>>
>>>> I think, in order to make speed reasonable, I will need to loop over
>>>> the rows in compiled code.
>>>>
>>>> My questions are:
>>>>
>>>> a) how easy would it be to try and link some external Cython code
>>>> against the fitpack code in scipy?
>>>> b) if that is not easy, what can I do that will also be of most use to
>>>> other scipy users?
>>>>
>>>> Cheers,
>>>>
>>>> Matthew
>>>
>>> I am not sure if I understand your problem statement. Most interpolation
>>> problems are posed as:  given data (x,y), determine yp at xp, where xp
>>> are different from x. This is fine for cases where you have a single y
>>> and possibly multiple sets xp. I had case awhile ago where I had a
>>> single xp and multiple y (at the same x). In this case, the problem
>>> statement is:  given known locations x and interpolant locations xp,
>>> what are the values for yp (corresponding to xp) given y (corresponding
>>> to xp). To solve this problem efficiently, I wrote the code I posted here:
>>>
>>> http://scipy-central.org/item/21/1/simple-piecewise-polynomial-interpolation
>>
>> nice,
>>
>> 3 questions
>>
>> Can you construct P = Xp*inv(X)*H  directly from the block diagonal
>> parts instead of going through sparse block-diagonal matrices?
>>

Perhaps. It has been 4 years since I worked on this, and so I don't 
remember all the details. What I did was use a loop (via list 
comprehension for speed) to calculate the inverse of all the small block 
matrices, then constructed the sparse Xinv. There may very well be a 
more elegant approach.

>> Why do you need monotonically increasing x? I would have thought it
>> works for any x
>>
>> (I'm not clear about this)
>>

I don't remember why x needs to be monotonic. I am usually careful to 
consider these things, so I am sure that there were at least nontrivial 
implementation problems if x is scattered (or repeated). Or maybe I just 
assumed that it would be easy for users (me!) to sort the data fist and 
didn't want to take the time to worry about it.

>> My guess is that replacing inv by pinv would be more stable. But,
>> would pinv also be able to handle the case when there are some
>> identical x (by implicitly taking an average)?
>

I am not so knowledgeable here. By taking the direct inverse of the 
individual small block matrices, I think the inverse method I used is 
quite stable.

> Something is wrong in the way I think about this.
> (I might be thinking about a different use case where I was recently
> using Kernel Ridge Regression.)
>
> Sorry for any noise

No worries. Nice to see some interest.

Jonathan

>
> Josef
>
>>
>> Josef
>>
>>>
>>> Although it is not true spline interpolation (rather a simpler local
>>> polynomial interpolation), it is reasonably fast and accurate.
>>> Extrapolation is allowed. However, I later found that
>>> scipy.interpolate.splrep/splev are so fast that they could be used
>>> instead with equivalent speed for my problem. By comparison,
>>> scipy.interpolate.interp1d is quite slow.
>>>
>>> HTH,
>>> Jonathan
>>> _______________________________________________
>>> SciPy-Dev mailing list
>>> SciPy-Dev at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-dev
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>