How to improve speed of interp1 in for loop?

Hi Matlab Community. For an extensive simulation I have to determine a probability density function based on a value and then use this to generate a pseudo-random value. That means I currently have a for-loop over each input value. Within the for-loop, the corresponding PDF is selected, the CDF is determined, and then a value is calculated using inverse transform sampling.

Since I usually have 1e6 input values per simulation run and my current solution takes about 40 seconds per simulation run on my powerful machine (MacBook with Apple M1 Pro), I wanted to know how I can solve my problem (much) faster because I need a lot of simulation runs. I want to achive one simulation run in less than a sec.

Here is the code example: In the matrix pdfMatrix a PDF is stored for a value range per value. The values for the rows are stored in the xFine vector. In this vector I first look for the indices of the values between which my input value falls. The PDF is then linearly interpolated using the input value and the two matrix rows.The CDF is then formed from the PDF and inverse transform sampling is carried out.

clear;
close all;
clc;

%% generate example inputValues
Num = 1e4; 
a = 0;
b = 2.999;
inputValues = a + (b-a).*rand(Num, 1);

%% corresponding values for pdfMatrix 
[XFine, YFine]     = meshgrid( ...
    linspace(0, 3, 161), ...
    linspace(0, 5, 353));

xFine = XFine(1,:);
yFine = YFine(:,1);

%% generate example pdfMatrix (here from a polynomial-fun + random)
% (my real pdfs vary much more depending on the value)
p = [-0.000164891028955014;
      0.00416401268241395; 
     -0.0428329447031803; 
      0.229754093919656;
     -0.682114138568206; 
      1.09686911465066; 
     -0.879541290571242;
      0.319778906564097;
      0];

pdfMatrix = polyval(p, YFine) + (0.03 * rand(1,161));

figure; 
mesh(xFine,yFine, pdfMatrix);
xlabel('Values');
ylabel('Corresponding pdf');

%% from here the calculation starts

tic;
% check if inputValue is lying in xFine range
err = inputValues < min(xFine) | inputValues > max(xFine);

% get indices from xFine (between which inputValue is lying)
[~, idx1] = max(~(xFine <= inputValues), [], 2);
idx1 = idx1-1;

[~, idx2] = max((xFine >= inputValues), [], 2);

% for each input value generate a pseudorandom value from the corresponding
% pdf from pdfMatrix
for i = 1:length(inputValues)
    if ~err(i)
 
        pdf = interp1(...
        xFine(idx1(i):idx2(i)), ...
        pdfMatrix(:,idx1(i):idx2(i))', ...
        inputValues(i));

        temp = cumsum(pdf); % numerical integration of PDF to get CDF

        cdf = temp/max(temp); % normalization of CDF

        % inverse transform sampling
        outputValue(i) = interp1(cdf, ...
            yFine, ...
            rand(1), ...
            'linear', ...
            'extrap'); 
    end
end
toc

% Num = 1e4: Elapsed time is ~ 0.5 seconds.
% Num = 1e5: Elapsed time is ~   5 seconds.
% Num = 1e6: Elapsed time is ~  50 seconds. <- must be ~100 times faster

One possibility would be to make my pdfMatrix extremely fine. This would not give me an exact result, but would eliminate an interpolation. Also it would be possible to pre-calculate the cdf (either with a extreme fine grid or interplote the cdf). But in the end the inverse transform sampling with interp1 has to be inside the for-loop, which is extremely slow. I already vectorized the index search, to reduce code inside the for-loop. par-for yielded no speed improve.

EDIT 1: I added example data to run the script. I hope now the problem is clearer.

EDIT 2: I modified my code with a interp1q and parfor loop. GriddedInterpolant for the first interp1q call was somehow slower. (maybe someone could explain this, as it should be faster)

clear;
close all;
clc;
parpool;

%% generate example inputValues
Num = 1e4; 
a = 0;
b = 3;
inputValues = a + (b-a).*rand(Num, 1);

%% corresponding values for pdfMatrix 
[XFine, YFine]     = meshgrid( ...
    linspace(a, b, 81), ...
    linspace(0, 5, 177));

xFine = XFine(1,:);
yFine = YFine(:,1);

%% generate example pdfMatrix (here from a polynomial-fun)
p = [-0.000164891028955014;
      0.00416401268241395; 
     -0.0428329447031803; 
      0.229754093919656;
     -0.682114138568206; 
      1.09686911465066; 
     -0.879541290571242;
      0.319778906564097;
      0];

pdfMatrix = polyval(p, YFine) + (0.03 * rand(1,81)); % cols = pdf

figure; 
mesh(xFine, yFine, pdfMatrix);
xlabel('Values');
ylabel('Corresponding pdf');

%% from here the calculation starts

tic;
% check if inputValue lies in xFine range
err = inputValues < min(xFine) | inputValues > max(xFine);

% get indices from xFine (between which inputValue is lying)
[~, idx1] = max(~(xFine <= inputValues), [], 2);
idx1 = idx1-1;

[~, idx2] = max((xFine >= inputValues), [], 2);

% F = griddedInterpolant(XFine', YFine', pdfMatrix'); % is slower

% for each input value generate a pseudorandom value from the corresponding
% pdf from pdfMatrix
parfor i = 1:length(inputValues)
    if ~err(i)

        pdf = interp1q(...
        xFine(idx1(i):idx2(i)), ...
        pdfMatrix(:,idx1(i):idx2(i))', ...
        inputValues(i));
        
        % pdf = F(ones(177,1)*inputValues(i), yFine);  % is slower

        temp = cumsum(pdf); % numerical integration of PDF to get CDF

        cdf = temp/max(temp); % normalization of CDF

        % inverse transform sampling
        outputValue(i) = interp1q(cdf, ...
            yFine, ...
            rand(1)); 
    end
end
toc

% Num = 1e4: Elapsed time is ~ 0.2 seconds.
% Num = 1e5: Elapsed time is ~ 0.4 seconds.
% Num = 1e6: Elapsed time is ~  2 seconds. <- needs to be ~4 times faster

[Torsten]

For an extensive simulation I have to determine a probability density function based on a value and then use this to generate a pseudo-random value. That means I currently have a for-loop over each input value. Within the for-loop, the corresponding PDF is selected, the CDF is determined, and then a value is calculated using inverse transform sampling.

It must be wrong what you are doing here. You can't base a pdf and cdf on a single value.

[CL]

I have a matrix with predefined PDFs. Each row stands for a (evenly spaced) value. In this row is a pdf for this value stored.

[Stephen23]

https://www.mathworks.com/help/matlab/ref/interp1q.html

[CL]

@Stephen23 thanks, with interp1q its roughly 4 times faster. But for the number of samples my calculation still need over 10 sec

Answer 1

[Star Strider]

I’m not certain that I follow what you’re doing. However if you know the ‘inputValue’ and you want to find the closest indices to it, first do something like this:

L = numel(xFine); % Length Of Vector (Put Outside The Loop If This Does Not Change Within The Loop)

idx = find(diff(sign(xFine - inputValue))); % Closest Index Of ‘xFine’ To ‘inputValue’

idxrng = max(1,idx-1) : min(L,idx+1); % Index Range For Interpolation

pdfv = interp1(xFine(idxrng), pdfMatrix(:,idxrng), inputValue); % Interpolate

I don’t have a good feeling for your data (I don’t have any feeling at all for them, actually, since they’re not provided), however this is usually the easiest way to find the indices of the values in a curve (‘xFine’ here) that intersect a specific value (‘inputValue’ here). That might significantly shorten at least one of the loops, since it goes directly to the index-of-interest and then uses that for the interpolation. (The ‘idxrng’ vector contains safeguards that prevent it from addressing elements of ‘xFine’ that are beyond its index limits. If this will never be a problem, you can just use ‘idx-1 : idx+1’ for it.)

.

[CL]

Thanks so far. The search for the indices should increase with your approach. But there are still 2 interp1 calls in the for-loop. I will provide example data later today.

EDIT: i provided example data

[Star Strider]

I do not have any feel for what you are doing, so I have no idea if this does what you want.

This appears to be reasonably fast, however I have no idea if it gives the correct result —

% parpool;
%% generate example inputValues
Num = 1e6; 
a = 0;
b = 3;
inputValues = a + (b-a).*rand(Num, 1);
%% corresponding values for pdfMatrix 
[XFine, YFine]     = meshgrid( ...
    linspace(a, b, 81), ...
    linspace(0, 5, 177));
xFine = XFine(1,:);
yFine = YFine(:,1);
%% generate example pdfMatrix (here from a polynomial-fun)
p = [-0.000164891028955014;
      0.00416401268241395; 
     -0.0428329447031803; 
      0.229754093919656;
     -0.682114138568206; 
      1.09686911465066; 
     -0.879541290571242;
      0.319778906564097;
      0];
pdfMatrix = polyval(p, YFine) + (0.03 * rand(1,81)); % cols = pdf
figure; 
mesh(xFine, yFine, pdfMatrix);
xlabel('Values');
ylabel('Corresponding pdf');

%% from here the calculation starts
tic;
% check if inputValue lies in xFine range
err = inputValues < min(xFine) | inputValues > max(xFine);
% get indices from xFine (between which inputValue is lying)
[~, idx1] = max(~(xFine <= inputValues), [], 2);
idx1 = idx1-1;
[~, idx2] = max((xFine >= inputValues), [], 2);
% F = griddedInterpolant(XFine', YFine', pdfMatrix'); % is slower
% for each input value generate a pseudorandom value from the corresponding
% pdf from pdfMatrix
% parfor i = 1:2%length(inputValues)
for i = 1:length(inputValues)
    if ~err(i)
        % pdf = interp1q(...
        % xFine(idx1(i):idx2(i)), ...
        % pdfMatrix(:,idx1(i):idx2(i))', ...
        % inputValues(i));
        
        cidx = find(diff(sign(xFine-inputValues(i))));
        idxrng = max(1,cidx-1) : min(cidx+1,numel(xFine));
        % Q1 = xFine(idxrng)
        % Q2 = pdfMatrix(idxrng,:)
        pdf = interp1q(xFine(idxrng).', pdfMatrix(idxrng,:), inputValues(i));
        
        % pdf = F(ones(177,1)*inputValues(i), yFine);  % is slower
        temp = cumsum(pdf); % numerical integration of PDF to get CDF
        cdf = temp/max(temp); % normalization of CDF
        % inverse transform sampling
        outputValue(i) = interp1q(cdf, ...
            yFine, ...
            rand(1)); 
    end
end
toc

Elapsed time is 9.866645 seconds.

outputValue

outputValue = 1×1000000

0.4626 0.0596 1.2360 1.0639 2.1955 0.5270 1.1940 1.7477 1.0856 NaN 2.2201 0.5776 0.8550 NaN 1.8480 1.6365 1.1358 1.8638 0.9238 0.1815 0.0655 2.1606 1.1306 0.2366 0.8433 1.0835 0.2106 2.2149 1.6029 2.0146

Does this give the desired results?

.