------Question update------

Thank you for your answers. However, suppose that we have to solve $n$ times a system $Ax = b$, where $A$ is always the same matrix. Consider that

-$A$ is full, and thus $A^{-1}$ requires the same memory storage than $A$.

-The condition number of $A$ is small, hence $A^{-1}$ can be computed with accuracy.

In that case, would not it be more efficient to compute $A^{-1}$ rather than to use a LU decomposition? For example, I have tried this Matlab code:

%Set A and b:
M           = 1000; 
A           = rand(M,M);
A           = real(expm(1i*(A+A.')));
b           = rand(M,1);

%Times we solve the system:
n           = 3000;

%Performing LU decomposition:
disp('Performing LU decomposition')
tic
[L,U,P]     = lu(A);
toc
fprintf('\n')

%Solving the system n times with LU decomposition:
optsL.LT    = true;   %Options for linsolve
optsU.UT    = true;
disp('Solving the system n times using LU decomposition')
tic
for ii=1:n
    x1      = linsolve(U, linsolve(L,P*b,optsL) , optsU);
end
toc
fprintf('\n')

%Computing inverse of A:
disp('Computing inverse of A')
tic
Ainv        = inv(A);
toc
fprintf('\n')

%Solving the system n times with Ainv:
disp('Solving the system n times with A inv')
tic
for ii=1:n
    x2  = Ainv*b;
end
toc
fprintf('\n')

disp('Residuals')
disp(norm(b-A*x1))
disp(norm(b-A*x2))

disp('Condition number of A')
disp(cond(A))

For a matrix with condition number about 450, the residuals are $O(10^{-11})$ in both cases, but it takes 19 seconds for solving the system n times using the LU decomposition, whereas using the inverse of A it only takes 9 seconds.

I am aware about that inverting a matrix to solve a linear system is not a good idea, since it is not as accurate and as efficient as directly solving the system or using LU, Cholesky or QR decomposition.

However, I have not been able to check this with a practical example. I have tried this code (in MATLAB)

M = 500;

A = rand(M,M);

A = real(expm(1i*(A+A.')));

b = rand(M,1);

x1 = A\b;

x2 = inv(A)*b;

disp(norm(b-A*x1))

disp(norm(b-A*x2))

and the residuals are always of the same order (10^-13).

Could someone provide a practical example in which inv(A)*b is much less inaccurate than A\b?

Thanks, Manu

I am aware about that inverting a matrix to solve a linear system is not a good idea, since it is not as accurate and as efficient as directly solving the system or using LU, Cholesky or QR decomposition.

However, I have not been able to check this with a practical example. I have tried this code (in MATLAB)

M   = 500;    
A   = rand(M,M);
A   = real(expm(1i*(A+A.')));
b   = rand(M,1);

x1  = A\b;
x2  = inv(A)*b;

disp(norm(b-A*x1))
disp(norm(b-A*x2))

and the residuals are always of the same order (10^-13).

Could someone provide a practical example in which inv(A)*b is much less inaccurate than A\b?