Loại bỏ Gauss trong Matlab

Question

Tôi đang sử dụng mã MATLAB từ cuốn sách này: http://books.google.com/books/about/Probability_Markov_chains_queues_and_sim.html?id=HdAQdzAjl60C Dưới đây là Code:

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); 

Có một mã nhanh hơn mà tôi không biết? Tôi gọi chức năng này hàng triệu lần và mất quá nhiều thời gian.

Answer 1

MATLAB có toàn bộ các phương thức đại số tuyến tính tích hợp - loại help slash, help lu hoặc help chol để bắt đầu với một số cách phổ biến để giải các phương trình tuyến tính hiệu quả trong MATLAB.

Theo dõi, các chức năng này thường được gọi là tối ưu hóa LAPACK/BLAS thói quen thư viện, thường là cách nhanh nhất để thực hiện đại số tuyến tính bằng bất kỳ ngôn ngữ lập trình nào. So với một ngôn ngữ "chậm" như MATLAB nó sẽ không được bất ngờ nếu họ được đơn đặt hàng của cường độ nhanh hơn so với một thực hiện m-file.

Hy vọng điều này sẽ hữu ích.

Answer 2

Trừ khi bạn đang muốn thực hiện riêng của mình, bạn nên sử dụng toán tử dấu gạch chéo ngược của Matlab (mldivide) hoặc, nếu bạn muốn các yếu tố, lu. Lưu ý rằng mldivide có thể làm được nhiều hơn loại bỏ Gaussian (ví dụ, nó có hình vuông nhỏ nhất tuyến tính, khi thích hợp).

Các thuật toán được sử dụng bởi mldivide và lu là từ thư viện C và Fortran và triển khai của riêng bạn trong Matlab sẽ không bao giờ nhanh như vậy. Tuy nhiên, nếu bạn quyết tâm sử dụng triển khai của riêng bạn và muốn nó được nhanh hơn, một lựa chọn là tìm cách để vector hóa việc triển khai của bạn (có thể bắt đầu here).

Một điều khác cần lưu ý: việc thực hiện từ câu hỏi không làm bất kỳ pivoting, vì vậy tính ổn định số nói chung sẽ tồi tệ hơn việc triển khai thực hiện xoay vòng và thậm chí sẽ thất bại đối với một số ma trận vô nghĩa.

Các biến thể khác nhau của loại bỏ Gaussian tồn tại, nhưng chúng là tất cả các thuật toán O (n). Nếu bất kỳ cách tiếp cận nào tốt hơn phương pháp khác phụ thuộc vào tình huống cụ thể của bạn và là điều bạn cần phải điều tra thêm.

Answer 3

Đối với cách tiếp cận ngây thơ (hay còn gọi là không trao đổi hàng) cho một n bởi ma trậ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 

Answer 4

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 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

Giả sử Ax = d đâu A và d được gọi ma trận. Chúng tôi muốn đại diện cho "A" là "L U" bằng cách sử dụng chức năng "LU phân hủy" được nhúng trong MATLAB: LUx = d Điều này có thể được thực hiện trong MATLAB sau: [L, U] = lu (A) trong đó về mặt trả về một ma trận tam giác trên trong U và ma trận hình tam giác thấp hơn được phép trong L sao cho A = L U. Giá trị trả về L là một sản phẩm của ma trận hình tam giác và hoán vị thấp hơn. (https://www.mathworks.com/help/matlab/ref/lu.html)

Sau đó, nếu chúng ta giả sử Ly = d trong đó y = Ux. Vì x là không xác định, do đó y cũng không biết, bằng cách biết y chúng ta tìm x như sau: y = L \ d; x = U \ y

và giải pháp được lưu trữ trong x. Đây là cách đơn giản nhất để giải quyết hệ phương trình tuyến tính cho rằng các ma trận không phải là số ít (tức là yếu tố quyết định của ma trận A và d không bằng 0), nếu không, chất lượng của dung dịch sẽ không tốt như mong đợi và có thể mang lại kết quả sai.

nếu ma trận là số ít, do đó không thể đảo ngược, một phương pháp khác nên được sử dụng để giải quyết hệ phương trình tuyến tính.