Vectorization of non equal matrix elements

Hello

I got a bit of a headache with whis problem. I can't find a solution myself and I am not sure it is possible at all in Matlab.

I got my code vectorized down so instead of 1000 lines looped over millions of times, it is 100 lines that is looped over millions of times. But this still need a lot of power to do. (And repmat is not an option or any other memory heavy operations.) The last part is a problem with non-equal vectors that has to be added to repeated indices. So I if it is not possible to vectorize, there may be a faster way to set up the code without demanding tremendously amounts of memory and it is here I need help.

To give an idea of the problem, here is the last part of the code:

% Explenation of variables
N = elements which can be E.g. 10.000.000
omegaGrid = Grid of values which the calculations has to be projected to
(also maybe 1.000.000)
lowerLine = a vector of Nx1 with lower position in omegaGrid to start
upperLine = a vector of Nx1 with upper position in omegaGrid to end
%
%%%%%%%%%%%Code to vectorize start here %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This loop is still quite fast, but with 1.000.000 elements it still takes
% alot of time to do.
for i = 1:length(omegaGrid)
    boundIndexLower = plus(boundIndexLower, (lowerLine >= omegaGrid(i)) > 0);
    boundIndexUpper = plus(boundIndexUpper, (upperLine >= omegaGrid(i)) > 0);
end
%
% Main calculation loop, which is the big hurdle
for i = 1:N
    % Function that is only consistent of simple calculations that can be
    % vectorized easily, if the problem with omegaGrid changing boundaries
    % can be fixed.
    lineProfile = ProfileVoigt(Calc1, Calc2, Calc3, ...
    omegaGrid(boundIndexLower(i):boundIndexUpper(i)));
    %
    % Iterative addition to the boundary index points.
    absorptionCoefficient(boundIndexLower(i):boundIndexUpper(i)) = ...
    absorptionCoefficient(boundIndexLower(i):boundIndexUpper(i)) + ...
    factor / Calc4 * Calc5 .* Calc6 .* lineProfile(:);
end

As it is seen there is a hurdle in the non consistensy of "(boundIndexLower(i):boundIndexUpper(i))" as they may not be the same length at each iteration. And at the same it is not possible to create a large matrix. (maybe sparse, but I can't see how to do it properly)

This may be impossible or a really complex system, but if there is anybody who can help me I'll grateful.

¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤

¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤

¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤

¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤ EDIT ¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤¤

To show the problem by code I have inculded a section of code that demostrate the problem. It can run on Matlab 9.1.0.441655 (R2016b)

% Initializing the script
clearvars
clc

% Setting up constants for simulation of the problem
numberOfRows   = 1e6; % Data resolution used from a big database
resolution     = 1;
lineShift      = 50;
factor         = 1; % This is a factor used for unit conversion
wavelengthMax  = 5e3;

% This is the number of points that the data has to be projected on
numberOfPoints = wavelengthMax/resolution; 

% Initializing realistic vectors as they would look in the real code
lineVector     = sort(rand(numberOfRows,1).*wavelengthMax);
omegaGrid      = sort(rand(numberOfPoints,1).*wavelengthMax);
abdundances    = rand(numberOfRows,1);
lineIntensity  = rand(numberOfRows,1).*1e-8;
shift0         = rand(numberOfRows,1).*1e1;
gammaD         = rand(numberOfRows,1).*1e1;
gamma0         = rand(numberOfRows,1).*1e1;

% Initializing the boundary index
boundIndexLower       = zeros(numberOfRows,1);
boundIndexUpper       = zeros(numberOfRows,1);
absorptionCoefficient = zeros(numberOfPoints,1);

% Finding the boundaries for the line vector
lowerLine = lineVector - lineShift;
upperLine = lineVector + lineShift;

% Determinding the points where the lower and upper index are located for
% all row elements. This was the fastest method I could think of/find, as 
% the grid always have the lowest number of points
for i = 1:numberOfPoints
    boundIndexLower = boundIndexLower + double((lowerLine>=omegaGrid(i))>0);
    boundIndexUpper = boundIndexUpper + double((upperLine>=omegaGrid(i))>0);
end

% making sure that the 0th value is the first index
boundIndexLower = boundIndexLower + 1;

% Looping thorugh the line centers
for i = 1:numberOfRows

      % Lineshape calculations
      lineProfile = ProfileVoigt(lineVector(i)+shift0(i),gammaD(i),gamma0(i),...
                         omegaGrid(boundIndexLower(i):boundIndexUpper(i)));

      % Saving the absorption coefficient
      absorptionCoefficient(boundIndexLower(i):boundIndexUpper(i)) = ...
          absorptionCoefficient(boundIndexLower(i):boundIndexUpper(i)) + ...
          factor/abdundances(i)*abdundances(i)*lineIntensity(i).*lineProfile;

end

function lineProfile = ProfileVoigt(lineIn, gammaD, gamma0, omegaGrid)
% This is to give an idea of what type of calculation are in this function

% Variable assigning
sg    = omegaGrid; % Vector that can be different sizes each loop
gamD  = gammaD;    % Single value taken out of a column vector
gam0  = gamma0;    % Single value taken out of a column vector
sg0   = lineIn;    % Single value taken out of a column vector

% These are not used in this function the way it is called
eta    = 0;
anuVC  = 0; 
shift0 = 0;
shift2 = 0;
gam2   = 0;

% Calculations
numberOfPoints = length(sg);

% Initialization
ATermGlobal = zeros(1,numberOfPoints);
BTermGlobal = zeros(1,numberOfPoints);

% Setting up complex numbers
cte = sqrt(log(2))/gamD;
rpi = sqrt(pi);
c0  = gam0 + 1i*shift0;
c2  = gam2 + 1i*shift2;
c0t = (1 - eta)*(c0 - 1.5*c2) + anuVC;
c2t = (1 - eta)*c2;

%_________________________________________________________________________%
% ------------------------------------------------------------------------%
% first part of the calculations

% The if statement can change, but this is just an example on the function
if abs(c2t) == 0 
    Z1  = (1i*(sg0 - sg) + c0t) * cte;
    xZ1 = -imag(Z1);
    yZ1 = real(Z1);

      t    = yZ1-1i*xZ1;
      cerf = 1/sqrt(pi)*t./(0.5+t.^2);

      % Finding the real and imaginary parts
      WR1 = real(cerf);
      WI1 = imag(cerf);

      ATermGlobal = rpi * cte * WR1 + 1i*WI1;

      indexZ1 = length(Z1(abs(Z1) <= 4e3));

      if indexZ1 > 0
          BTermGlobal = rpi*cte*((1 - Z1.^2).*(WR1 + 1i*WI1) + Z1./rpi);
      elseif indexZ1 == 0
          BTermGlobal = cte*(rpi*(WR1 + 1i*WI1) + 0.5./Z1 - 0.75./Z1.^3);
      end
  else
      X = (1i * (sg0 - sg) + c0t)./c2t;
      Y = 1/(2*cte*c2t)^2;

      csqrtY = (gam2 - 1i*shift2)/(2*cte*(1-eta)*(gam2^2+shift2^2));

      % Finding index parameters to determined which parts that has to be
      % calculated
      indexPart2 = logical(X.*(abs(X) <= 3e-8*abs(Y)));
      indexPart3 = logical(Y.*(abs(Y) <= 1e-15*abs(X) & ~indexPart2));
      indexPart4 = ~(indexPart2 | indexPart3);

      % This section continue with some other parts that are not relevant for
      % the question and therefore I cut the code here.
  end

LSpCqSDHC = (1/pi) * (ATermGlobal./(1 -(anuVC-eta*(c0-1.5*c2))...
             .*ATermGlobal + eta*c2*BTermGlobal));

lineProfile = real(LSpCqSDHC);

end

[Eng. Fredius Magige]

Hi Not yet clear at all; You might mix the two/three options, while vs for or vice versa. Thanks

[Benjamin]

I did not have time yesterday to make a simulation code of the problem, but there is included working code now. This could clear some misunderstanding

Answer 1

[Benjamin]

After a couple of weeks i found a solution to increase the speed of the boundary index loop:

for i = 1:numberOfPoints
    boundIndexLower = boundIndexLower + double((lowerLine>=omegaGrid(i))>0);
    boundIndexUpper = boundIndexUpper + double((upperLine>=omegaGrid(i))>0);
end

I did this by applying arithmetics on the "lowerLine" and "upperLine" vectors.

 gradientOmegaGrid = omegaGrid(2)-omegaGrid(1); 

 boundIndexLower = ceil(lowerLine/(gradientOmegaGrid(1)) - ...
                                      min(omegaGrid)/gradientOmegaGrid(1));                           
 boundIndexUpper = ceil(upperLine/(gradientOmegaGrid(1)) - ...
                                      min(omegaGrid)/gradientOmegaGrid(1));

This code generally calculate the index value based on the omegaGrid step size and the start value. By doing so it calculate the index values that the vectors refere to in the loop. I have not used tic-toc on it, but it do not show on the profiler compared to the loop - huge speed optimization!

For the last loop that is looped millions of times I have chosen to compile it. Dependt on the size of omegaGrid it optimize the speed between 2.5-10x. (in absolute values the 2.5x optimization gave 100 seconds of the calculation time for normally ~270 seconds)

Answer 2

[Jan]

I cannot imagine the purpose of:

for i = 1:length(omegaGrid)
  boundIndexLower = plus(boundIndexLower, (lowerLine >= omegaGrid(i)) > 0);
end

What is the use of the "> 0"? (lowerLine >= omegaGrid(i)) is a logical vector, therefore I assume the "> 0" wastes time only. What does this code do at all? For each element of "lowerLine" it is counted how many elements of "omegaGrid" are larger - correct? This should be much faster with an interpolation method, which replies the next larger element (not a linear interpolation): Collect for each value of "omegaGrid" how many elements are larger (this should be a simple sorting and unsing the index vector). Then find the next larger element for each value of lowerLine. This should avoid this loop.

This might contain several vector operations, when the scalar operations are processed from left to right:

X = X + factor / Calc4 * Calc5 .* Calc6 .* lineProfile(:)

Use parenthesis:

X = X + ((factor / Calc4 * Calc5 .* Calc6) .* lineProfile(:));

If Calc5 and Calc6 are vectors (please do not let us guess but post equivalent code, which we can run by copy and paste), than process this constants before the loop.

[Benjamin]

Hello and thanks for taking your time for answering

I am sorry that I did not have time to make a simulation code of the problem, but I made an edit with a working example of the problem that demonstrate all parts of the optimization problem. This should help on some of the misunderstanding.