Двойное суммирование с векторизованными циклами в Matlab

Question

Я хочу векторизовать этот двойной цикл for, потому что это узкое место в моем коде. Поскольку Matlab является индексированным языком на основе одного, я должен создать дополнительный термин для M = 0.

R, r, лямбда - постоянные

Slm (L, M), Clm (L, M) - матрицы 70x70

Plm (L, M) - матрица 70x71

Cl (L), Pl (L) - векторы 70x1

% function dU_r
  s1 = 0;
  for L = 2:70
     s1 = s1 + ((R/r)^L)*(L+1)*Pl(L)*Cl(L);
     for m = 1:L
        s1 = s1 + ((R/r)^L)*(L+1)*Plm(L,M)*(Clm(L,M)*...
cos(M*lambda) + Slm(L,M)*sin(M*lambda));
     end
  end

  dU_r = -(mu_p/(r^2))*s1;


  % function dU_phi

  s2=0;
  for L = 2:70
      s2 = s2 + ((R/r)^L))*Plm(L,1)*Cl(L);
          for m = 1:l
          s2 = s2 + ((R/r)^L)*(Plm(L,M+1)-M*tan(phi)*Plm(L,M))*...
  (Clm(L,M)*cos(M*lambda) + Slm(L,M)*sin(M*lambda));
          end;
  end;
  dU_phi = (mu_p/r)*s2;

  % function dU_lambda

  s3=0;
      for L=2:70
          for m=1:L
          s3 = s3 + ((R/r)^L)*M*Plm(L,M)*(Slm(L,M)*cos(M*lambda)...
              - Clm(L,M)*sin(M*lambda));
          end;
  dU_lambda = (mu_p/r)*s3;

Answer 1

В этой конкретной ситуации для вас, вероятно, будет более эффективно , а не полностью векторизовать ваше решение. Как показывает Egon , можно заменить внутренний цикл путем векторизации, но замена внешнего цикла for также может потенциально замедлить вашего решения.

Причина?Если вы обратите внимание на то, как вы индексируете матрицы Plm, Clm и Slm в цикле, вы всегда будете использовать значения только из нижних треугольных частей из них.Все значения выше главной диагонали игнорируются.Путем полной векторизации вашего решения, так что в нем используются матричные операции вместо цикла, вы в итоге выполняете ненужные умножения с обнуленными верхними треугольными частями матриц.Это один из тех случаев, когда цикл for может быть лучшим вариантом.

Таким образом, вот уточненная и исправленная версия ответа Эгона , которая должна быть близка к оптимальной в отношенииколичество выполненных математических операций:

index = 2:70;
coeff = ((R/r).^index).*(index+1);
lambda = lambda.*(1:70).';  %'
cosVector = cos(lambda);
sinVector = sin(lambda);
s2 = coeff*(Pl(index).*Cl(index));
for L = index
  M = 1:L;
  s2 = s2 + coeff(L-1)*((Plm(L,M).*Clm(L,M))*cosVector(M) + ...
                        (Plm(L,M).*Slm(L,M))*sinVector(M));
end

Answer 2

Я думаю, что следующий код может быть частью решения (только частично векторизованного), для полной векторизации вы можете посмотреть на meshgrid, ndgrid, bsxfun и / или repmat, I предположим.

s2     = 0;
L      = 2:70;
PlCl   = Pl.*Cl;
PlmClm = Plm.*Clm;
Rrl    = (R/r).^(L).*(L+1);
COS    = cos((1:70)*lambda);
SIN    = sin((1:70)*lambda);

for l=L
    M  = 1:l;
    s1 = sum(Rrl(l-1)*PlmClm(l,M).*(COS(M) + Slm(l,M).*SIN(M)));
    s2 = s2 + s1;
end;
s2 = s2 + sum(Rrl(L).*PlCl(L));

Я не пытался запустить этот код, так что, если он может жаловаться на некоторые измерения. Вы должны быть в состоянии понять это (только некоторые переносы здесь или там могут сделать).

Просто примечание на стороне: старайтесь держаться подальше от коротких имен переменных (и, конечно, неоднозначных, таких как l (это может быть неправильно истолковано для 1, I или L). Кроме того, это может быть проще векторизовать ваш код, если у вас есть фактические (то есть не кодированные) выражения для начала.

edit : примененные исправления, предложенные gnovice

Answer 3

Взяв в качестве примера решение, данное Яном Симоном 1 , я написал следующий код для моей проблемы

Rr = R/r;
RrL = Rr;  
cosLambda = cos((1:70)* lambda);
sinLambda = sin((1:70)* lambda);
u1 = uint8(1);
s1 = 0;
s2 = 0;
s3 = 0;
for j = uint8(2):uint8(70)
   RrL = RrL * Rr;
   q1 = RrL * (double(j) + 1);
   t1 = Pl(j) * datos.Cl(j);
   q2 = RrL;
   t2 = Plm(j,1) * Cl(j);
   t3 = 0;
   for m = u1:j
      t1 = t1 + Plm(j,m) * ...
         (Clm(j, m) * cosLambda(m) + ...
          Slm(j, m) * sinLambda(m));
        t2 = t2 + (Plm(j,m+1)-double(m)*tan_phi*Plm(j,m))*...
    (Clm(j,m)*cosLambda(m) + Slm(j,m)*sinLambda(m));
        t3 = t3 + double(m)*Plm(j,m)*(Slm(j,m)*cosLambda(m)...
                - Clm(j,m)*sinLambda(m));
     end
     s1 = s1 + q1 * t1;
     s2 = s2 + q2 * t2;
     s3 = s3 + q2 * t3;
  end
dU_r = -(mu_p/(r^2))*s1;
dU_phi = (mu_p/r)*s2;
dU_lambda = (mu_p/r)*s3