Why is solving linear equation more stable than directly computing matrix inverse? To compute $A^{-1}B$ using numerical packages like numpy, it is suggested to use solve(A,B) rather than inv(A).dot(B). I know that you can apply LU factorization to speed up in the former case. But why it is more stable than the second choice? How does they work internally differently?
Many thanks.
 A: If $x$ is computed by LU factorization, the residual can be bounded by
$$
\|b-Ax\|\leq cu\||L||U|\|\|x\|.
$$
Assuming for simplicity that $A^{-1}$ is computed exactly and that the only source of error is the matrix-vector multiplication, we have
$$
\|b-Ax\|\leq cu\||A||A^{-1}|\|\|b\|.
$$
Here, $c$ is a "moderate constant", $u$ is the unit round-off, and $\|\cdot\|$ is the $\infty$-norm. 
Often, $\||L||U|\|\approx\|A\|$ (e.g., with a suitable pivoting strategy in LU). The bounds indicate that you might get a much smaller residual with the LU factorization provided that 
$$\|x\|\approx\|A^{-1}b\|\ll\|A^{-1}\|\|b\|.$$
For illustration, consider this piece of Matlab code:
n = 32;
A = 2 * eye(n) - triu(ones(n));
[U, S, V] = svd(A);
b = U(:,1);

x_LU  = A \ b;
x_INV = inv(A) * b;

norm(b - A * x_LU,  inf)   % 5.551e-17
norm(b - A * x_INV, inf)   % 9.543e-10

The choice of $A$ and $b$ are quite extreme ($A$ highly ill-conditioned, $b$ in the direction of the first left singular vector) but they illustrate clearly the issue here (disregarding of course the computational efficiency).
