Подтвердить что ты не робот

MATLAB Оптимизация взвешенной грамотрицательной ортогонализации

У меня есть функция в MATLAB, которая выполняет Gram-Schmidt Orthogonalisation с очень важным весом, применяемым к внутренним продуктам (я не знаю Не думаю, что встроенная функция MATLAB поддерживает это). Эта функция работает хорошо, насколько я могу судить, однако она слишком медленная на больших матрицах. Какой был бы лучший способ улучшить это?

Я пробовал конвертировать в MEX файл, но я теряю параллелизацию с используемым компилятором, и поэтому он медленнее.

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

Может ли кто-нибудь прорисовать этот код или сделать его быстрее? Я не уверен, как это сделать элегантно...

Я знаю, что улов stackoverflow здесь потрясающий, рассмотрите эту задачу:)

Функция

function [Q, R] = Gram_Schmidt(A, w)
    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));
    v = zeros(n, 1);

    for j = 1:n
        v = A(:,j);
        for i = 1:j-1
            R(i,j) = sum(   v    .* conj( Q(:,i) ) .* w ) / ...
                     sum( Q(:,i) .* conj( Q(:,i) ) .* w );
            v = v - R(i,j) * Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
    end
end

где A - матрица m x n комплексных чисел, а w - вектор вещественных чисел m x 1.

Бутылка шея

Это выражение для R(i,j), которое является самой медленной частью функции (не на 100% уверенно, если нотация правильная): Weighted Inner-Product

где w - неотрицательная весовая функция. Весовой внутренний продукт упоминается на нескольких страницах Википедии, это один из весовых функций и это одно из ортогональных функций.

Воспроизведение

Вы можете создавать результаты, используя следующий script:

A = complex( rand(360000,100), rand(360000,100));
w = rand(360000, 1);
[Q, R] = Gram_Schmidt(A, w);

где A и w - входы.

Скорость и вычисления

Если вы используете вышеуказанный script, вы получите результаты профилировщика, синонимичные следующим образом: Profiler TimesProfiler Code Times

Результат тестирования

Вы можете проверить результаты, сравнив функцию с приведенной выше, используя следующий script:

A = complex( rand( 100, 10), rand( 100, 10));
w = rand( 100, 1);
[Q , R ] = Gram_Schmidt( A, w);
[Q2, R2] = Gram_Schmidt2( A, w);
zeros1 = norm( Q - Q2 );
zeros2 = norm( R - R2 );

где Gram_Schmidt - функция, описанная ранее, и Gram_Schmidt2 - альтернативная функция. Результаты zeros1 и zeros2 должны быть очень близки к нулю.

Примечание:

Я попытался ускорить вычисление R(i,j) следующим, но безрезультатно...

R(i,j) = ( w' * (   v    .* conj( Q(:,i) ) ) ) / ...
         ( w' * ( Q(:,i) .* conj( Q(:,i) ) ) );
4b9b3361

ОТВЕТЫ

Ответ 1

1)

Моя первая попытка векторизации:

function [Q, R] = Gram_Schmidt1(A, w)
    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));

    for j = 1:n
        v = A(:,j);
        QQ = Q(:,1:j-1);
        QQ = bsxfun(@rdivide, bsxfun(@times, w, conj(QQ)), w.' * abs(QQ).^2);
        for i = 1:j-1
            R(i,j) = (v.' * QQ(:,i));
            v = v - R(i,j) * Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
    end
end

К сожалению, он оказался медленнее исходной.

2)

Тогда я понял, что столбцы этой промежуточной матрицы QQ строятся постепенно и предыдущие не изменяются. Итак, вот моя вторая попытка:

function [Q, R] = Gram_Schmidt2(A, w)
    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));
    QQ = complex(zeros(m, n-1));

    for j = 1:n
        if j>1
            qj = Q(:,j-1);
            QQ(:,j-1) = (conj(qj) .* w) ./ (w.' * (qj.*conj(qj)));
        end
        v = A(:,j);
        for i = 1:j-1
            R(i,j) = (v.' * QQ(:,i));
            v = v - R(i,j) * Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
    end
end

Технически никакой основной векторизации не было сделано; Я только предварительно вычислил промежуточные результаты и переместил вычисление вне внутреннего цикла.

На основе быстрого теста эта новая версия, безусловно, быстрее:

% some random data
>> M = 10000; N = 100;
>> A = complex(rand(M,N), rand(M,N));
>> w = rand(M,1);

% time
>> timeit(@() Gram_Schmidt(A,w), 2)     % original version
ans =
    1.2444
>> timeit(@() Gram_Schmidt1(A,w), 2)    % first attempt (vectorized)
ans =
    2.0990
>> timeit(@() Gram_Schmidt2(A,w), 2)    % final version
ans =
    0.4698

% check results
>> [Q,R] = Gram_Schmidt(A,w);
>> [Q2,R2] = Gram_Schmidt2(A,w);
>> norm(Q-Q2)
ans =
   4.2796e-14
>> norm(R-R2)
ans =
   1.7782e-12

EDIT:

Следуя комментариям, мы можем переписать второе решение, чтобы избавиться от if-statmenet, переместив эту часть в конец внешнего цикла (т.е. сразу после вычисления нового столбца Q(:,j), мы вычислим и сохраним соответствует QQ(:,j)).

Функция идентична на выходе, и выбор времени тоже не отличается; код немного короче!

function [Q, R] = Gram_Schmidt3(A, w)
    [m, n] = size(A);
    Q = zeros(m, n, 'like',A);
    R = zeros(n, n, 'like',A);
    QQ = zeros(m, n, 'like',A);

    for j = 1:n
        v = A(:,j);
        for i = 1:j-1
            R(i,j) = (v.' * QQ(:,i));
            v = v - R(i,j) * Q(:,i);
        end
        R(j,j) = norm(v);
        Q(:,j) = v / R(j,j);
        QQ(:,j) = (conj(Q(:,j)) .* w) ./ (w.' * (Q(:,j).*conj(Q(:,j))));
    end
end

Обратите внимание, что я использовал синтаксис zeros(..., 'like',A) (новый в последних версиях MATLAB). Это позволяет нам запускать функцию без изменений на GPU (при условии, что у вас есть Parallel Computing Toolbox):

% CPU
[Q3,R3] = Gram_Schmidt3(A, w);

против.

% GPU
AA = gpuArray(A);
[Q3,R3] = Gram_Schmidt3(AA, w);

К сожалению, в моем случае это было не так быстро. На самом деле это было во много раз медленнее работать на GPU, чем на процессоре, но стоило сделать это:)

Ответ 2

Здесь долгая дискуссия, но, чтобы перейти к ответу. Вы взвешивали числитель и знаменатель вычисления R вектором w. Утяжеление происходит во внутренней петле и состоит из трехточечного продукта, точки Q dot w в числителе и Q точки Q точки w в знаменателе. Если вы сделаете одно изменение, я думаю, что код будет работать значительно быстрее. Напишите num = (dot sqrt (w)) dot (Q dot sqrt (w)) и напишите den = (Q dot sqrt (w)) dot (Q dot sqrt (w)). Это перемещает вычисления продукта (A dot sqrt (w)) и (Q dot sqrt (w)) из внутреннего цикла.

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

"Цели" ГСО в два раза. Во-первых, чтобы разрешить решение уравнения типа Ax = y, где A имеет гораздо больше строк, чем столбцы. Эта ситуация часто возникает при измерении данных, поскольку легко измерять больше данных, чем количество состояний. Подход к первой цели состоит в том, чтобы переписать A как QR так, что столбцы Q ортогональны и нормированы, а R - треугольная матрица. Я считаю, что алгоритм, который вы предоставили, достигает первой цели. Q представляет собой базисное пространство матрицы A, а R представляет собой амплитуду каждого базисного пространства, необходимого для генерации каждого столбца A.

Вторая цель GSO - ранжировать базисные векторы в порядке значимости. Это шаг, который вы еще не сделали. И, включив этот шаг, может увеличить время решения, результаты будут определять, какие элементы x важны, согласно данным, содержащимся в измерениях, представленных A.

Но, я думаю, с этой реализацией решение быстрее, чем предлагаемый вами подход.

Aij = Qij Rij, где Qj ортонормированы и Rij верхнетреугольный, Ri, j > я = 0. Qj - ортогональные базисные векторы для A, а Rij - участие каждого Qj для создания столбца в A. Итак,

A_j1 = Q_j1 * R_1,1
A_j2 = Q_j1 * R_1,2 + Q_j2 * R_2,2
A_j3 = Q_j1 * R_1,3 + Q_j2 * R_2,3 + Q_j3 * R_3,3

При проверке вы можете написать

A_j1 = ( A_j1 / | A_j1 | )  * | A_j1 | = Q_j1 * R_1,1

Затем вы проектируете Q_j1 из каждого другого столбца A, чтобы получить элементы R_1, j

R_1,2 = Q_j1 dot Aj2
R_1,3 = Q_j1 dot Aj3
...
R_1,j(j>1) = A_j dot Q_j1

Затем вы вычтите элементы проекта Q_j1 из столбцов A (это установило бы первый столбец равным нулю, поэтому вы можете игнорировать первый столбец

for j = 2,n
  A_j = A_j - R_1,j * Q_j1
end

Теперь был удален один столбец из A, первый ортонормированный базисный вектор, Q, j1, и был определен вклад первого базисного вектора в каждый столбец R_1, j, а вклад первого базисный вектор был вычтен из каждого столбца. Повторите этот процесс для остальных столбцов A, чтобы получить оставшиеся столбцы Q и строки R.

for i = 1,n
  R_ii = |A_i|                    A_i is the ith column of A, |A_i| is magnitude of A_i
  Q_i = A_i / R_ii                Q_i is the ith column of Q
  for j = i, n
    R_ij = | A_j dot Q_i | 
    A_j = A_j - R_ij * Q_i
  end
end

Вы пытаетесь весить строки A, с w. Вот один из подходов. Я бы нормализовал w и включил эффект в R. Вы "удалили" эффекты w умножением и делением на w. Альтернативой "удалению" эффекта является нормализация амплитуды w до единицы.

w = w / | w |
for i = 1,n
  R_ii = |A_i inner product w|        # A_i inner product w = A_i .* w                  
  Q_i = A_i / R_ii              
  for j = i, n
    R_ij = | (A_i inner product w) dot Q_i |      # A dot B = A' * B 
    A_j = A_j - R_ij * Q_i
  end
end

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

Использование следующего может помочь в ускорении вашего кода

A inner product B = A .* B
A dot w = A' w
(A B)' = B'A'
A' conj(A) = |A|^2

Вышеуказанное может легко быть векторизованным в Matlab, как написано.

Но вам не хватает второй части рейтинга A, которая сообщает вам, какие состояния (элементы x в x = y) значительно представлены в данных

Процедуру ранжирования легко описать, но я расскажу вам о деталях программирования. Вышеупомянутая процедура по существу предполагает, что столбцы A находятся в порядке значимости, а первый столбец вычитается из всех остальных столбцов, затем второй столбец вычитается из остальных столбцов и т.д. Первая строка R представляет собой вклад первый столбец Q для каждого столбца A. Если вы суммируете абсолютное значение первой строки вкладов R, вы получаете измерение вклада первого столбца Q в матрицу A. Итак, вы просто оцениваете каждый столбец A как первый (или следующий) столбец Q и определить оценку ранжирования вклада этого столбца Q в остальные столбцы A. Затем выберите столбец A, который имеет самый высокий ранг в качестве следующего столбца Q. Кодирование по существу сводится к предварительной оценке следующей строки R для каждого оставшегося столбца в A, чтобы определить, какая ранговая величина R имеет наибольшую амплитуду. Наличие вектора индекса, который представляет собой исходный порядок столбцов A, будет полезен. Посредством ранжирования базисных векторов вы получаете "основные" базовые векторы, которые представляют собой A, который обычно намного меньше по числу, чем число столбцов в A.

Кроме того, если вы оцениваете столбцы, нет необходимости вычислять каждый столбец R. Когда вы знаете, какие столбцы A не содержат никакой полезной информации, нет никакой реальной выгоды для сохранения этих столбцов.

В структурной динамике одним из подходов к уменьшению числа степеней свободы является вычисление собственных значений, если у вас есть репрезентативные значения для матрицы массы и жесткости. Если вы подумаете об этом, приведенный выше подход может быть использован для "вычисления" матриц M и K (и C) из измеренного отклика, а также определения "форм ответа на измерение", которые значительно представлены в данных. Они отличаются друг от друга и потенциально более важными, чем формы режима. Таким образом, вы можете решить очень сложные проблемы, т.е. Оценку матриц состояний и числа степеней свободы, представленных, от измеренного отклика, вышеупомянутым подходом. Если вы читаете N4SID, он сделал что-то подобное, за исключением того, что использовал SVD вместо GSO. Мне не нравится техническое описание N4SID, слишком большое внимание уделяется векторной проекции, которая является просто точечным продуктом.

В приведенной выше информации может быть одна или две ошибки, я пишу это с ног до головы, прежде чем спешить с работы. Итак, проверьте алгоритм/уравнения, как вы реализуете... Удача

Возвращаясь к вашему вопросу, как оптимизировать алгоритм при весе с w. Вот базовый алгоритм GSO без сортировки, написанный совместимый с вашей функцией.

Обратите внимание, что код ниже находится в октаве, а не в таблице. Есть некоторые незначительные отличия.

function [Q, R] = Gram_Schmidt_2(A, w)

    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));

    # Outer loop identifies the basis vectors    
    for j = 1:n
      aCol = A(:,j);
      # Subtract off the basis vector
      for i = 1:(j-1)
        R(i,j) = ctranspose(Q(:,j)) * aCol;
        aCol = aCol - R(i,j) * Q(:,j);
      end
      amp_A_col = norm(aCol);
      R(j,j) = amp_A_col;
      Q(:,j) = aCol / amp_A_col;
    end
end

Чтобы получить ваш алгоритм, измените только одну строку. Но вы теряете большую скорость, потому что "ctranspose (Q (:, j)) * aCol" является векторной операцией, но "sum (aCol. * Conj (Q (:, i)). * W)" - это строка операции.

function [Q, R] = Gram_Schmidt_2(A, w)

    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));

    # Outer loop identifies the basis vectors    
    for j = 1:n
      aCol = A(:,j);
      # Subtract off the basis vector
      for i = 1:(j-1)
        # R(i,j) = ctranspose(Q(:,j)) * aCol;
        R(i,j) = sum(   aCol    .* conj( Q(:,i) ) .* w ) / ...
                 sum( Q(:,i) .* conj( Q(:,i) ) .* w );
        aCol = aCol - R(i,j) * Q(:,j);
      end
      amp_A_col = norm(aCol);
      R(j,j) = amp_A_col;
      Q(:,j) = aCol / amp_A_col;
    end
end

Вы можете изменить его обратно на векторную операцию, взвешивая aCol и Q с помощью sqrt w.

function [Q, R] = Gram_Schmidt_3(A, w)

    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));
    Q_sw = complex(zeros(m, n));
    sw = w .^ 0.5;
    for j = 1:n
      aCol = A(:,j);
      aCol_sw = aCol .* sw;
      # Subtract off the basis vector
      for i = 1:(j-1)
        # R(i,j) = ctranspose(Q(:,i)) * aCol;
        numTerm = ctranspose( Q_sw(:,i) ) * aCol_sw;
        denTerm = ctranspose( Q_sw(:,i) ) * Q_sw(:,i);
        R(i,j) = numTerm / denTerm;
        aCol_sw = aCol_sw - R(i,j) * Q_sw(:,i);
      end
      aCol = aCol_sw ./ sw;
      amp_A_col = norm(aCol);
      R(j,j) = amp_A_col;
      Q(:,j) = aCol / amp_A_col;
      Q_sw(:,j) = Q(:,j) .* sw;
    end
end

Как указано JacobD, вышеуказанная функция работает не быстрее. Возможно, потребуется время для создания дополнительных массивов. Другой стратегией группировки тройного произведения является группа w с conj (Q). Надеюсь, что это быстрее...

function [Q, R] = Gram_Schmidt_4(A, w)

    [m, n] = size(A);
    Q = complex(zeros(m, n));
    R = complex(zeros(n, n));

    for j = 1:n
      aCol = A(:,j);
      for i = 1:(j-1)
        cqw = conj(Q(:,i)) .* w;
        R(i,j) = ( transpose(  aCol )  * cqw) ...
               / (transpose( Q(:,i) ) * cqw);
        aCol = aCol - R(i,j) * Q(:,i);
      end
      amp_A_col = norm(aCol);
      R(j,j) = amp_A_col;
      Q(:,j) = aCol / amp_A_col;
    end
end

Здесь используется функция драйвера для разных версий.

function Gram_Schmidt_tester_2
      nSamples = 360000;
      nMeas = 100;
      nMeas = 15;

      A = complex( rand(nSamples,nMeas), rand(nSamples,nMeas));
      w = rand(nSamples, 1);

      profile on;
      [Q1, R1] = Gram_Schmidt_basic(A);
      profile off;
      data1 = profile ("info");
      tData1=data1.FunctionTable(1).TotalTime;
      approx_zero1 = A - Q1 * R1;
      max_value1 = max(max(abs(approx_zero1)));

      profile on;
      [Q2, R2] = Gram_Schmidt_w_Orig(A, w);
      profile off;
      data2 = profile ("info");
      tData2=data2.FunctionTable(1).TotalTime;
      approx_zero2 = A - Q2 * R2;
      max_value2 = max(max(abs(approx_zero2)));

      sw=w.^0.5;
      profile on;
      [Q3, R3] = Gram_Schmidt_sqrt_w(A, w);
      profile off;
      data3 = profile ("info");
      tData3=data3.FunctionTable(1).TotalTime;
      approx_zero3 = A - Q3 * R3;
      max_value3 = max(max(abs(approx_zero3)));

      profile on;
      [Q4, R4] = Gram_Schmidt_4(A, w);
      profile off;
      data4 = profile ("info");
      tData4=data4.FunctionTable(1).TotalTime;
      approx_zero4 = A - Q4 * R4;
      max_value4 = max(max(abs(approx_zero4)));

      profile on;
      [Q5, R5] = Gram_Schmidt_5(A, w);
      profile off;
      data5 = profile ("info");
      tData5=data5.FunctionTable(1).TotalTime;
      approx_zero5 = A - Q5 * R5;
      max_value5 = max(max(abs(approx_zero5)));


      profile on;
      [Q2a, R2a] = Gram_Schmidt2a(A, w);
      profile off;
      data2a = profile ("info");
      tData2a=data2a.FunctionTable(1).TotalTime;
      approx_zero2a = A - Q2a * R2a;
      max_value2a = max(max(abs(approx_zero2a)));


      profshow (data1, 6);
      profshow (data2, 6);
      profshow (data3, 6);
      profshow (data4, 6);
      profshow (data5, 6);
      profshow (data2a, 6);

      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data1.FunctionTable(1).FunctionName,
         data1.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value1)
      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data2.FunctionTable(1).FunctionName,
         data2.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value2)
      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data3.FunctionTable(1).FunctionName,
         data3.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value3)
      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data4.FunctionTable(1).FunctionName,
         data4.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value4)
      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data5.FunctionTable(1).FunctionName,
         data5.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value5)
      sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
         data2a.FunctionTable(1).FunctionName,
         data2a.FunctionTable(1).TotalTime,     
         nSamples, nMeas, max_value2a)

end

На моем старом домашнем ноутбуке в Octave результаты

ans = Time for Gram_Schmidt_basic is 0.889 sec with 360000 samples and 15 meas, max value is 1.57009e-16
ans = Time for Gram_Schmidt_w_Orig is 0.952 sec with 360000 samples and 15 meas, max value is 6.36717e-16
ans = Time for Gram_Schmidt_sqrt_w is 0.390 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_4 is 0.452 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_5 is 2.636 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt2a is 0.905 sec with 360000 samples and 15 meas, max value is 6.68443e-16

Эти результаты показывают, что самым быстрым алгоритмом является алгоритм sqrt_w выше на 0,39 с, а затем группировка conj (Q) с w (выше) на 0,452 с, затем версия 2 решения Amro в 0,905 секунды, затем исходный алгоритм в вопросе на 0,952, затем в версии 5, которая меняет местами строки/столбцы, чтобы увидеть, было ли хранилище строк (код не включен) в 2.636 сек. Эти результаты показывают, что разделение sqrt (w) между A и Q является самым быстрым решением. Но эти результаты не согласуются с комментарием JacobD о том, что sqrt (w) не выполняется быстрее.

Ответ 3

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

function [Q, R] = Gram_Schmidt5(A, w)
Q = A;
n_dimensions = size(A, 2);
R = zeros(n_dimensions);
R(1, 1) = norm(Q(:, 1));
Q(:, 1) = Q(:, 1) ./ R(1, 1);
for i = 2 : n_dimensions
    Qw = (Q(:, i - 1) .* w)' * Q(:, (i - 1) : end);
    R(i - 1, i : end) = Qw(2:end) / Qw(1);
    %% Surprisingly this loop beats the matrix multiply
    for j = i : n_dimensions
        Q(:, j) = Q(:, j) - Q(:, i - 1) * R(i - 1, j);
    end
    %% This multiply is slower than above
    %    Q(:, i : end) = ...
    %     Q(:, i : end) - ...
    %     Q(:, i - 1) * R(i - 1, i : end);
    R(i, i) = norm(Q(:,i));
    Q(:, i) = Q(:, i) ./ R(i, i);
end