Гауссово исключение в Matlab

Question

Я использую код Matlab из этой книги: http://books.google.com/books/about/Probability_Markov_chains_queues_and_sim.html?id=HdAQdzAjl60C Вот код:

    function [pi] = GE(Q)

    A = Q';
    n = size(A);
    for i=1:n-1
      for j=i+1:n
         A(j,i) = -A(j,i)/A(i,i);
      end
         for j =i+1:n
            for k=i+1:n
        A(j,k) = A(j,k)+ A(j,i) * A(i,k);
         end
      end
      end

     x(n) = 1;
      for i = n-1:-1:1
        for j= i+1:n
          x(i) = x(i) + A(i,j)*x(j);
        end
       x(i) = -x(i)/A(i,i);
      end

      pi = x/norm(x,1);

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

Answer 1

В MATLAB есть целый набор встроенных подпрограмм линейной алгебры - введите help slash, help lu или help chol, чтобы начать работу с несколькими распространенными способами эффективного решения линейных уравнений в MATLAB.

Под капотом эти функции обычно вызывают оптимизированные LAPACK / BLAS библиотечные процедуры, которые, как правило, являются самым быстрым способом выполнения линейной алгебры в любом языке программирования.По сравнению с «медленным» языком, таким как MATLAB, не было бы неожиданным, если бы они были на несколько порядков быстрее, чем реализация m-файла.

Надеюсь, это поможет.

Answer 2

Если вы специально не планируете реализовать свою собственную, вам следует использовать оператор обратной косой черты Matlab (mldivide) или, если вам нужны коэффициенты, lu.Обратите внимание, что mldivide может делать больше, чем исключение по Гауссу (например, он делает линейные наименьшие квадраты, когда это необходимо).

Алгоритмы, используемые mldivide и lu, взяты из библиотек C и Fortran и вашей собственной реализациив Matlab никогда не будет так быстро.Однако, если вы полны решимости использовать собственную реализацию и хотите, чтобы она была быстрее, одним из вариантов является поиск способов векторизации вашей реализации (возможно, начните здесь ).

Еще одинНа что следует обратить внимание: реализация из этого вопроса не выполняет никакого поворота , поэтому его числовая стабильность, как правило, будет хуже, чем у реализации, которая выполняет поворот, и даже для некоторых невырожденных матриц не получится.

Существуют разные варианты исключения Гаусса, но все они являются алгоритмами O ( n ³ ).Если какой-либо из подходов лучше, чем другой, зависит от вашей конкретной ситуации и вам нужно больше изучить.

Answer 3

function x = naiv_gauss(A,b);
n = length(b); x = zeros(n,1);
for k=1:n-1 % forward elimination
      for i=k+1:n
           xmult = A(i,k)/A(k,k);
           for j=k+1:n
             A(i,j) = A(i,j)-xmult*A(k,j);
           end
           b(i) = b(i)-xmult*b(k);
      end
end
 % back substitution
x(n) = b(n)/A(n,n);
for i=n-1:-1:1
   sum = b(i);
   for j=i+1:n
     sum = sum-A(i,j)*x(j);
   end
   x(i) = sum/A(i,i);
end
end

Answer 4

Предположим, Ax = d Где A и d - известные матрицы. Мы хотим представить «A» как «L U», используя функцию «LU-декомпозиция», встроенную в matlab, таким образом: LUx = d Это может быть сделано в Matlab следующим образом: [L, U] = Лу (А) который в терминах возвращает верхнюю треугольную матрицу в U и переставленную нижнюю треугольную матрицу в L так, что A = L U. Возвращаемое значение L является произведением нижних треугольных и перестановочных матриц. (https://www.mathworks.com/help/matlab/ref/lu.html)

Тогда, если мы примем Ly = d, где y = Ux. Поскольку x неизвестен, значит, y тоже неизвестен, зная y, мы находим x следующим образом: у = Ь \ д; х = U \ у

и решение сохраняется в x.

Это самый простой способ решения системы линейных уравнений при условии, что матрицы не являются единичными (т. Е. Определитель матрицы A и d не равен нулю), в противном случае качество решения будет не таким хорошим, как ожидалось, и может дать неправильные результаты.

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

Answer 5

function Sol = GaussianElimination(A,b)


[i,j] = size(A);

for j = 1:i-1


    for i = j+1:i


        Sol(i,j) = Sol(i,:) -( Sol(i,j)/(Sol(j,j)*Sol(j,:)));


    end


end


disp(Sol);


end

Answer 6

Для наивного подхода (он же без перестановки строк) для матрицы n на n:

function A = naiveGauss(A)

% find's the size

n = size(A);
n = n(1);

B = zeros(n,1);

% We have 3 steps for a 4x4 matrix so we have
% n-1 steps for an nxn matrix
for k = 1 : n-1     
    for i = k+1 : n
        % step 1: Create multiples that would make the top left 1
        % printf("multi = %d / %d\n", A(i,k), A(k,k), A(i,k)/A(k,k) )
        for j = k : n
            A(i,j) = A(i,j) - (A(i,k)/A(k,k)) *  A(k,j);
        end
        B(i) = B(i) - (A(i,k)/A(k,k))  * B(k);
    end
end