Why using mldivide to solve high dimension linear equations would cause computational error

I tried to use A\b to solve the linear equation Ax=b, where A is a n*n matrix shown as followed:

when n is relatively small (like 10 or 20), it provided a good solution. But when n is greater than approximately 60, the solution provided is drifted from the real answer. The error between the provided solution and the actual solution becomes greater and greater when n increase.

I was wondering why this happened. A hint is said that the gaussian elimination for the matrix A is unstable, and coincidentally, the mldivide operation \ is using gaussian elimination. However, I still cannot find a obvious relation between the hint and the phenomenon.

Answer 1

[Mohammad Abouali]

I am guessing the condition number of your matrix.

The error in matrix inversion get's worse as the condition number of a matrix increases. When your matrix is 10x10 the condition number of your specific matrix is around 4. while when your matrix is 200x200 it's condition number is around 90. So this means that small perturbation will be amplified until the error will dominate the solution.

Although I have to say that these condition number are not that high. Actually they are small enough. But anyway, one thing that you can use is using some sort of conditioner matrix. So, instead of solving Ax=b you would solve PAX=Pb, where P is your conditioner which hopefully (PA) has better condition number.

before going to use conditioner, just make sure that you are constructing your A matrix properly and there isn't code error.

something like this should do it:

A=-tril(ones(n,n),-1)+eye(n); A(:,end)=1;

By the way, in your case, it appears that you can analytically invert these matrices. The inverted matrix can be constructed as follow:

n=200;
Ainv=eye(n,n)/2;
for i=1:n-1
  Ainv(i, (i+1:n-1) )= -2.^-(2:(n-i));
end
Ainv(n,1:n-1)=2.^-(1:n-1);
Ainv(1:n-1,n)=Ainv(1:n-1,n-1);
Ainv(n-1,n)=-0.5;
Ainv(n,n)=Ainv(n,n-1);

[Jiasen]

Thanks for your comment. But I think it is not condition number that matters. Please look at my comment to the first answer. There is actually another part in the question. The perturbed matrix Aeps has the same condition number as the original A, but Aeps could provide a much more accurate solution that A do.

[Mohammad Abouali]

Well, yes, the small perturbation that you are adding is not going to affect the condition number. Because you are not changing the size of the matrix and the elements are not changed that much really. The condition number however changes with increasing n, i.e. the matrix size.

By the way, the condition number amplifies the perturbation in b vector in Ax=b and consequently affects the solution of x.

Anyway, your condition number is not that high in my opinion. when n=200 the condition number is about 90.

Did you try the analytic inverse matrix? (The code at the bottom of my previous post).

[Jiasen]

Wow, the analytic inverse matrix is amazing! Learned something new :D

But I was still wondering why perturbation can improve the accuracy of A\b.

[Mohammad Abouali]

I am not sure what you mean by "perturbation improves the accuracy of A\b". To my understanding it is going to the revers. So, if there is a small perturbation in b it is going to be amplified by A\b. The higher the condition number of A is the more amplified that small perturbation in b would be in A\b. It becomes a threshold that your solution is actually dominated by the amplified perturbation of b and your answer would be all wrong. So, actually it is not a good thing in this case.

[Jiasen]

Anyway, thanks for your patients.

[SK]

The perturbation is added to the matrix A, and not to b, so the poster is asking why the solution is accurate with the perturbation added. For example, let A be 100 x 100 and let A1 = A + 0.0001*rand(100). Then, you can check that

1. The "perturbed A" solution is quite close to the actual solution (with A as the lhs) when b = rand(n,1).

2. The condition number of A1 is around the same as the condition number of A (around 45).

[Jiasen]

In this case it's not the condition number that matters. The true reason for the inaccuracy is that Matlab uses Gaussian Elimination to solve this but encounter ties when pivoting, which result in extremely large elements. It's perturbation that breaks the tie and make the Gaussian Elimination stable.

Answer 2

[SK]

To the OP:

Try taking different values of the RHS b. I'm pretty sure that for some particular vectors b, even the solution with perturbed A will not be accurate. There is no other explanation, since the condition numbers of A and A+perturbation are equally large.

[Jiasen]

In this case it's not the condition number that matters. The true reason for the inaccuracy is that Matlab uses Gaussian Elimination to solve this but encounter ties when pivoting, which result in extremely large elements. It's perturbation that breaks the tie and make the Gaussian Elimination stable.

[SK]

Thanks for taking the trouble to reply.

Yes, I think you are indeed right. In fact a condition number of 40 or 50 is actually quite small and should not lead to such large errors particularly when we are doing things upto 15 digits of precision. As further evidence that it is indeed instability in Gaussian elimination, take the following matrix:

B = A(end:-1:1, end:-1:1)

From the point of view of system of equations, this is exactly equivalent to the system obtained from A. Moreover condition number of A is the same as condition number of B. However taking b = rand(n,1) and solving:

y = B\b

the answer for y is very accurate, which shows that it has nothing to do with condition number and everything to do with the algorithm used for solution.

It would be good if you could post an answer and mark it answered yourself so that others could see the real reason and not be misled.

Regards.

[added later: Oh, I see that there is already an accepted answer so perhaps you could just add another answer]

Answer 3

[Jan]

The hint say: The Gauss-elimination is instable. It is applied and you observe an instable behavior of the results. The relation is not only obvious but trivial. Therefore it's not clear to me, what you are asking for.

[Jiasen]

Thanks. I'm sorry I forgot to post it, but actually there is anther part in the question: if you add a tiny random perturbation in the matrix A (eps * randn(n) + A, where eps is extremely small), the gaussian elimination applied to the perturbed matrix Aeps would become stable again. How to explain this?