Vectorization interpolation matlab

Question

I have a set of N measurements taken at S time points (time points are different for different measurements). I have two matrices:

V is the NxS matrix representing the measured values

T is the NxS matrix that represents the measurement time

I would like to generate a matrix VI representing linearly interpolated measurements at TI times. The following is a non-vector version of the code:

tic; VI = zeros([size(V,1), size(TI,2)]); for j = 1:size(V,1) VI(j,:) = interp1(T(j,:),V(j,:),TI); end toc; 

I would like to rewrite this code to exclude for-loop so that it is implemented using matrix operations and functions. Is it possible to vectorize it?

Answer 1

It's hard to say anything without data and start the profiler, but if your data is sorted, you can use interp1q instead of interp , which does not perform any data checks.

Taken from Matlab help:

For interp1q to work correctly, x must be a monotonically increasing column vector. Y must be a column vector or matrix with a length (x) of rows. xi must be a column vector

Answer 2

Matlab packs data into columns rather than rows, so I suspect you will see a speed improvement simply because the loop changes from iterating over rows into columns:

 [N, S] = size(V); VT = V'; % Value series in columns TT = T'; % Time series in columns VIT = zeros(length(TI), N); % Interpolated value series in columns for j = 1:N VIT(:, j) = interp1(VT(:, j), TT(:, j), TI); end VI = VIT'; % Interpolated value series in rows 

Unfortunately, I don’t think that too much can be done to further improve the performance of this code, since different series of values ​​do not have related time series. If they had the same times so that we could fold T into a vector with length(T) = S , then it would be easy to do this:

 VIT = interp1(VT, T, TI); % (see VIT and VT from above) 

Answer 3

If you have, as you say, a different measurement time for all measurements (T is a matrix, not a vector), you can do what you want with one call to arrayfun as follows:

 VI = arrayfun(@(x)(interp1(T(x,:),V(x,:),TI)), 1:size(V, 1), 'UniformOutput', false); VI = cell2mat(VI'); 

arrayfun is like a loop, but since it is a matlab internal function, it can be faster. It returns a cell of vectors, so the second line ensures that you have a matrix as output. You may not need this - it depends on what you later do with the data.

If, on the other hand, the measurements were performed at the same time for different values ​​of N (T is a vector of size S, not a matrix, or, in other words, all rows from T are equal), you can interpolate interp1 in a single call.

 VI = interp1(T(1,:), V', TI) 

Here you need to transpose V, as interp1 interpolates inside the columns. This is because MATLAB stores matrices in columns (columns are contiguous in memory). If you pass V as an SxN matrix, this potentially allows you to parallelize interp1 more efficiently, since all processors can access memory more efficiently. Therefore, I would suggest that you transfer your matrices to all of your code, unless, of course, you rely on this exact data format somewhere else for performance reasons.

Edit Due to the layout of the matrix columns, your source code can be improved by moving the matrices and the working column. The next version on my computer is about 20% higher for N = 1000, S = 10000 and TI 10000 elements. This is likely to grow with system size due to more efficient use of memory bandwidth.

 tic; VI = zeros(size(TI,2), size(V,2)); for j = 1:size(V,2) VI(:,j) = interp1(T(:,j),V(:,j),TI); end toc; 

Answer 4

I work and therefore do not have time to get to interp1 (I have never used it before), therefore, if the following answer is inappropriate, I apologize in advance.

Having said that, since the iterations of your loop are independent of each other, it is possible to vectorize. It seems to me that you can get rid of the explicit loop using mat2cell and cellfun .

A simple example of what I mean is the following:

 NumRow = 4; NumCol = 3; V = randn(NumRow, NumCol); VCell = mat2cell(V, ones(NumRow, 1), NumCol); A = cellfun(@sum, VCell); 

Of course, what I did is equivalent to sum(V, 2) . But I think that the method I used applies to your situation. The mat2cell function converts V to the column vector of the cell, where each cell contains the row V. Then the call to cellfun applies the sum function to each cell in VCell and returns the result in A You can simply replace @sum with @interp1 and, of course, adjust the inputs accordingly with cellfun .

Let me know if you cannot get it to work, and I will try to collect something more explicit as soon as I get home. In addition, if you earn it, but it does not speed up much, it would be very interesting for me to find out, so please post the result here.