Можно ли векторизовать 3 вложенных для циклов (всего 3), когда внутренние параметры 2 для циклов меняются с каждой итерацией первого основного цикла for?

Question

Все примеры и вопросы, которые я нашел по векторизации, касались статических параметров.Я был успешным в векторизации таких циклов в прошлом, но я все еще новичок.В первую очередь получение опыта по программированию упражнений на машинном обучении Coursera в Стэнфордском университете.

Я мог бы использовать справочную векторизацию кода, где внутренние параметры цикла изменяются с каждой итерацией.Мой код зацикливается на двухмерном наборе данных.Область поиска определяется для каждой точки данных в наборе.Каждая точка данных в области поиска, которая находится в заданном радиусе, будет использоваться для совокупного расчета. Размер области поиска различен для каждой точки;следовательно, два внутренних цикла for будут иметь разные размеры для каждой итерации.

Каждая из голубых точек на этом изображении представляет точку данных, пересекаемую внешним циклом for.Пурпурный крест представляет текущую оцениваемую точку.Существует один массив для расположения точек и другой массив для данных точек.Черный прямоугольник - это область поиска для этой точки.Два внутренних цикла for проходят каждую точку сетки в области поиска.Точки сетки определяются галочками на графике, которые создают равномерную декартову сетку.Зеленый кружок выделяет текущую точку сетки, которая оценивается.Поскольку точка сетки находится внутри пурпурного радиуса, будет выполнена сверточная интерполяция с данными из синей точки (пурпурные перекрестия) и добавлена ​​к текущему значению в этой точке сетки (инициализируется в ноль).

В этом примере текущая точка сетки (красный круг) находится за пределами радиуса (пурпурный круг), поэтому ничего не делается.

Важно отметить, что в приведенных выше примерах каждая синяя точкабудет иметь только одну область поиска и использоваться только в одной итерации внешнего цикла for, но каждая точка сетки будет сравниваться со всеми окружающими синими точками, и интерполяция свертки будет добавляться к их значениям, когда они достаточно близки (внутри пурпурногокружок) к точке выборки (синяя точка).

Можно ли векторизовать этот код?Или это невозможно, а я просто теряю время?

Вот пример кода, в котором неактуальный код удален, а некоторые переменные изменены для простоты:

for loc_i = 1:length(dataLocations)

    xMin_i = floor(((xdata(loc_i)-searchRadius)));
    . . . condition check to set min here

    xMax_i = ceil(((xdata(loc_i)+searchRadius)));
    . . . condition check to set max here

    yMin_i = floor(((ydata(loc_i)-searchRadius)));
    . . . condition check to set min here

    yMax_i = ceil(((ydata(loc_i)+searchRadius)));
    . . . condition check to set max here
%
%--end search area boundary index determination-----------------------%

%--traverse the search area-------------------------------------------%
%
% nested for loops will traverse search area column by column starting
%   in the bottom left corner

% update output indexes to align with starting point of the current
%   search area
%   xMin_i * gridSize_n finds the current column
%   + yMin_i finds the current row
%

realOut_i = ((xMin_i) * gridSize_n + yMin_i); 
imagOut_i = ((xMin_i) * gridSize_n + yMin_i);

% update incrementer for output index to move to the starting point of
%   the next column within the search area
%
nxtCol_inc = gridSize_n - (yMax_i - yMin_i) - 1;

    % traverse columns of search area
    %
    for col_i = xMin_i:xMax_i

    % calculate x component of the distance vector between current
    %   trajectory point and current cartesian grid point
    %
    dx_n = col_i - (xTraj_d(loc_i ));


        % traverse rows of each search area
        %
        for row_i = yMin_i:yMax_i       

            % calculate y component of the distance vector between current
            %   trajectory point and current cartesian grid point
            %
            dy_n = row_i - (yTraj_d(loc_i));

            % calculate magnitude of the distance vector between the
            %   current trajectory point and the current cartesian grid 
            %   point
            %
            dk_n = sqrt(dx_n^2 + dy_n^2);

            if dk_n < searchRadius

                % determine where the grid point lies within the covolution
                %   radius, measured as the distance from the center
                %
                fractionKernel_i =...
                    (dk_n / (kernelRadius_n)*(kernTabPnts_n-1));

                % seperate the integer portion of the indexed distance to
                %   retrieve indexed values from the kernel table
                %
                integerKernel_i = floor(fractionKernel_i);

                % save the decimal portion of the indexed distance to
                %   weight the values retrieved from the kernel table
                % 
                %
                ratioKernel_n = fractionKernel_i - integerKernel_i;

                % determine kernel value for current grid point
                % value is linearly interpolated for precision
                % 
                %
                kernel = ...
                    kernelTable_d(integerKernel_i+1)*(1-ratioKernel_n)+...
                    kernelTable_d(integerKernel_i+2)*ratioKernel_n;

                % do convolution interpolation for real values
                %
                realOutput(realOut_i) = ...
                    realOutput(realOut_i) + ...
                    kernel*realSample_d(loc_i)*dcf(loc_i);

                % do convolution interpolation for imaginary values
                %
                imagOutput(imagOut_i) = ...
                    imagOutput(imagOut_i) +...
                    kernel*imagSample_d(loc_i)*dcf(loc_i);

            end % if dk_n < searchRadius

            % move to the next row in the search column
            %
            realOut_i = realOut_i + 1;
            imagOut_i = imagOut_i + 1;


        end % row_i = yMin_i:yMax_i

        % move to the starting point of the next column in the search area
        %
        realOut_i = realOut_i + nxtCol_inc;
        imagOut_i = imagOut_i + nxtCol_inc;

    end % col_i = xMin_i:xMax_i

    %-- finished traversing search area-----------------------------------%

end % for loc_i = 1:length(dataLocations)