当内部2个for循环的参数随着第一个主for循环的每次迭代而改变时,是否可以对3个嵌套for循环(总共3个)进行矢量化?

我发现的所有矢量化的例子和问题都涉及静态参数。我已经成功地矢量化这样的for循环在过去,但仍然是一个初学者。主要获得斯坦福大学Coursera机器学习编程练习的经验。

Each of the blues dots in this image represent a data point being traversed by the outer for loop. The magenta cross represents the current point being evaluated. There is one array for the points location and another array for the points data. The black rectangle is the search area for this point. The two inner for loops traverse each grid point in the search area. The grid points are defined by the tick marks on the plot that create a uniform Cartesian grid. The green circle highlights the current grid point being evaluated. Because the grid point is inside the magenta radius a convolution interpolation will be performed with data from the blue point(magenta cross hairs) and added to the current value at that grid point(initialized to zero).

In this example the current grid point(red circle) is outside the radius(magenta circle) so nothing is done.

在上面的示例中需要注意的是,每个蓝点将只有一个搜索区域,并且仅用于外部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)