Accuracy Problem when solving linear equation system using: lu(S) vs. decomposition(S,'lu')

3 views (last 30 days)
Hello Matlab Community,
I want to solve an equation system B = S * X with a sparse coefficient matrix S by using LU factorisation.
There are two ways to do that:
  1. [l,u,p,q,d] = lu(S) (with the 2 triangular matrices l and u, the 2 permutation matrices p and q and a scaling matrix d)
  2. dS = decomposition(S,'lu') (with the decomposition object dS).
To solve the equation system I use:
  1. X = q * (u \ (l \ (p * (d \ B)))) (see Matlab documentation)
  2. X = dS \ B (where mldivide is a member function of dS).
I calculate the error of the two solutions by using norm(B - S * X). The following code evaluates the error for different sizes of S.
rng(42)
DeltaSize = 500;
M = DeltaSize:DeltaSize:10000;
ErrorLU = zeros(numel(M),1);
ErrorLUObject = zeros(numel(M),1);
for cnt = 1:numel(M)
disp(num2str(cnt))
MatrixDimension = M(cnt);
NumberOfNNZ = 0.1 * MatrixDimension^2;
% Build S and B
row = randi([1,MatrixDimension],NumberOfNNZ,1);
col = randi([1,MatrixDimension],NumberOfNNZ,1);
v = rand(NumberOfNNZ,1) + 1i * rand(NumberOfNNZ,1);
S = sparse(row,col,v,MatrixDimension,MatrixDimension) + speye(MatrixDimension);
B = rand(MatrixDimension,1) + 1i * rand(MatrixDimension,1);
% use LU
[l,u,p,q,d] = lu(S);
X = q * (u \ (l \ (p * (d \ B))));
ErrorLU(cnt) = norm(full(B - S * X));
% use decomposition object
dS = decomposition(S,'lu');
X = dS \ B;
ErrorLUObject(cnt) = norm(full(B - S * X));
end
LogErrorLU = log10(ErrorLU);
LogErrorLUObject = log10(ErrorLUObject);
figure(1)
plot(M,LogErrorLU,'.-','Displayname','lu(S)')
hold on
plot(M,LogErrorLUObject,'.-','Displayname','decomposition(S,lu)')
hold off
ylim([-17 -5])
legend show
grid on
When I run the code, I get a difference of more than three orders of magnitude...see the following figure.
Why is the error of lu(S) bigger than the error of dS? When I go into the decomposition.m file and from that into SparseLU.m I can see that the same matrices are computed in SparseLU like in lu(S) (Line 43 in SparseLU: only with: p -> Pleft, q -> Pright, d -> S).
The only difference I've spotted so far is, that when matlab computes X = dS \ B it uses the MuPAD (see: solve.m line 293).
Is this the only reason for the big difference between the two errors?
Thank you for your hlep.

Accepted Answer

Christine Tobler
Christine Tobler on 23 Dec 2021
The version in decomposition does some optional steps of iterative refinement: It uses the same solution you have above based on calling lu (up to minor round-off because things are stored in different format). However, it then checks the resulting residual and does another solve with just that residual as a right-hand side if necessary:
rng(42)
DeltaSize = 500;
M = DeltaSize:DeltaSize:2000; % Reduced length for run-time
ErrorLU = zeros(numel(M),1);
ErrorLUrefined = zeros(numel(M),1);
ErrorLUObject = zeros(numel(M),1);
for cnt = 1:numel(M)
disp(num2str(cnt))
MatrixDimension = M(cnt);
NumberOfNNZ = 0.1 * MatrixDimension^2;
% Build S and B
row = randi([1,MatrixDimension],NumberOfNNZ,1);
col = randi([1,MatrixDimension],NumberOfNNZ,1);
v = rand(NumberOfNNZ,1) + 1i * rand(NumberOfNNZ,1);
S = sparse(row,col,v,MatrixDimension,MatrixDimension) + speye(MatrixDimension);
B = rand(MatrixDimension,1) + 1i * rand(MatrixDimension,1);
% use LU
[l,u,p,q,d] = lu(S);
X = q * (u \ (l \ (p * (d \ B))));
ErrorLU(cnt) = norm(full(B - S * X));
% use LU with iterative refinement
R = B - S * X;
Xerr = q * (u \ (l \ (p * (d \ R))));
X = X + Xerr;
ErrorLUrefined(cnt) = norm(full(B - S * X));
% use decomposition object
dS = decomposition(S,'lu');
X = dS \ B;
ErrorLUObject(cnt) = norm(full(B - S * X));
end
1 2 3 4
LogErrorLU = log10(ErrorLU);
LogErrorLUObject = log10(ErrorLUObject);
figure(1)
semilogy(M,ErrorLU,'.-','Displayname','lu(S)')
hold on
semilogy(M,ErrorLUrefined,'x-','Displayname','lu(S) iterative refinement')
semilogy(M,ErrorLUObject,'.-','Displayname','decomposition(S,lu)')
hold off
ylim([1e-17 1e-5])
legend show
grid on;
ax = gca;
ax.YMinorGrid = 'off';
This is because decomposition does its best to match exactly what S\B does.

More Answers (0)

Products


Release

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by