5

我正在使用本书中的 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);

有没有我不知道的更快的代码？我正在调用这个函数数百万次，这需要太多时间。

4

7 回答 7

9

MATLAB 有一整套内置的线性代数例程 - 类型help slash，help lu或者help chol开始使用一些常用方法在 MATLAB 中有效地求解线性方程组。

在幕后，这些函数通常调用优化LAPACK/BLAS库例程，这通常是在任何编程语言中进行线性代数的最快方法。与像 MATLAB 这样的“慢”语言相比，如果它们比 m 文件实现快几个数量级也就不足为奇了。

希望这可以帮助。

于 2012-01-27T04:02:55.530 回答
8

除非您特别想实现自己的，否则您应该使用 Matlab 的反斜杠运算符 ( mldivide)，或者，如果您想要这些因素，则应使用lu. 请注意，mldivide 可以做的不仅仅是高斯消除（例如，它会在适当的时候进行线性最小二乘法）。

使用的算法mldivide和lu来自 C 和 Fortran 库的算法，你自己在 Matlab 中的实现永远不会那么快。但是，如果您决定使用自己的实现并希望它更快，那么一种选择是寻找向量化您的实现的方法（也许从这里开始）。

另一件需要注意的事情：问题中的实现不做任何旋转，因此它的数值稳定性通常会比旋转的实现更差，甚至对于一些非奇异矩阵会失败。

存在不同的高斯消元变体，但它们都是 O( n 3 ) 算法。如果任何一种方法比另一种更好，取决于您的特定情况，并且您需要进行更多调查。

于 2012-01-27T04:04:58.147 回答
5 
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
于 2013-11-15T13:12:43.997 回答
1

假设 Ax=d 其中 A 和 d 是已知矩阵。我们希望使用嵌入在 matlab 中的“LU 分解”函数将“A”表示为“L U”，因此： LUx = d 这可以在 matlab 中完成，如下所示： [L,U] = lu(A) 术语返回一个上限U 中的三角矩阵和 L 中的置换下三角矩阵，使得 A = L U。返回值 L 是下三角矩阵和置换矩阵的乘积。（https://www.mathworks.com/help/matlab/ref/lu.html）

那么如果我们假设 Ly = d 其中 y=Ux。由于 x 是未知的，因此 y 也是未知的，通过知道 y 我们找到 x 如下： y=L\d; x=U\y

并且解决方案存储在 x 中。

这是求解线性方程组的最简单方法，前提是矩阵不是奇异的（即矩阵 A 和 d 的行列式不为零），否则，解的质量不会像预期的那样好，并且可能会产生错误结果。

如果矩阵是奇异的因此不能逆，则应使用另一种方法来求解线性方程组。

于 2017-08-01T19:19:28.347 回答
0 
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
于 2014-02-28T13:26:33.293 回答
0

对于 n x 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
于 2012-09-27T03:07:12.150 回答
0

我认为您可以使用 matlab 函数 rref：

[R,jb] = rref(A,tol)

它以简化的行梯形形式生成矩阵。就我而言，这不是最快的解决方案。在我的情况下，下面的解决方案快了大约 30%。

function C = gauss_elimination(A,B)
i = 1; % loop variable
X = [ A B ]; 
[ nX mX ] = size( X); % determining the size of matrix  

while i <= nX % start of loop 
    if X(i,i) == 0 % checking if the diagonal elements are zero or not
        disp('Diagonal element zero') % displaying the result if there exists zero 
        return
    end
    X = elimination(X,i,i); % proceeding forward if diagonal elements are non-zero
    i = i +1;
end

C = X(:,mX);

function X = elimination(X,i,j)
% Pivoting (i,j) element of matrix X and eliminating other column
% elements to zero
[ nX mX ] = size( X); 
a = X(i,j);
X(i,:) = X(i,:)/a;

for k =  1:nX % loop to find triangular form 
    if k == i
        continue
    end
    X(k,:) = X(k,:) - X(i,:)*X(k,j); 
end
于 2019-11-04T18:03:47.240 回答