我使用的是这本书中的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);
有没有我不知道的更快的代码?我已经调用这个函数数百万次了,而且花费了太多时间。
MATLAB有一整套内置的线性代数例程-类型为help slash
、help lu
或help chol
,以开始了解在MATLAB中有效求解线性方程的一些常见方法。
在后台,这些函数通常调用优化的LAPACK
/BLAS
库例程,这通常是任何编程语言中执行线性代数的最快方法。与MATLAB这样的"慢"语言相比,如果它们比m文件实现快几个数量级,也就不足为奇了。
希望这能有所帮助。
除非您特别希望实现自己的反斜杠运算符,否则您应该使用Matlab的反斜杠操作符(mldivide
),如果您需要因子,则应该使用lu
。请注意,mldive可以做的不仅仅是高斯消去(例如,在适当的情况下,它可以做线性最小二乘)。
mldivide
和lu
使用的算法来自C和Fortran库,您自己在Matlab中的实现速度永远不会那么快。然而,如果您决心使用自己的实现并希望它更快,那么一种选择是寻找向量化实现的方法(也许从这里开始)。
需要注意的另一件事是:该问题的实现不进行任何旋转,因此其数值稳定性通常会比进行旋转的实现差,甚至对于一些非奇异矩阵也会失败。
高斯消去存在不同的变体,但它们都是O(n3)算法。如果任何一种方法比另一种方法更好,这取决于你的特定情况,并且是你需要更多调查的事情。
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
让我们假设Ax=d其中A和d是已知矩阵。我们希望使用嵌入在matlab中的"LU分解"函数将"A"表示为"LU",因此:LUx=d这可以在matlab中完成如下操作:[L,U]=lu(A)其返回U中的上三角矩阵和L中的置换下三角矩阵,使得a=LU。返回值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的行列式不为零),否则,解的质量将不如预期,并可能产生错误的结果。
如果矩阵是奇异的,因此不能求逆,则应使用另一种方法来求解线性方程组。
对于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 / %dn", 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
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
我认为你可以使用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