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

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.

Answer 1

[Christine Tobler]

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.

[Sebastian Raabe]

Thank you for your answer. :)