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

4 views (last 30 days)
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.

Accepted Answer

Mohammad Abouali
Mohammad Abouali on 28 Sep 2014
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);
  7 Comments
SK
SK on 2 Oct 2014
Edited: SK on 2 Oct 2014
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
Jiasen on 24 Oct 2014
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.

Sign in to comment.

More Answers (2)

SK
SK on 2 Oct 2014
Edited: SK on 2 Oct 2014
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.
  2 Comments
Jiasen
Jiasen on 24 Oct 2014
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
SK on 24 Oct 2014
Edited: SK on 24 Oct 2014
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]

Sign in to comment.


Jan
Jan on 27 Sep 2014
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.
  1 Comment
Jiasen
Jiasen on 28 Sep 2014
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?

Sign in to comment.

Products

Community Treasure Hunt

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

Start Hunting!