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

5 visualizzazioni (ultimi 30 giorni)
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.

Risposta accettata

Mohammad Abouali
Mohammad Abouali il 28 Set 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 Commenti
SK
SK il 2 Ott 2014
Modificato: SK il 2 Ott 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 il 24 Ott 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.

Accedi per commentare.

Più risposte (2)

SK
SK il 2 Ott 2014
Modificato: SK il 2 Ott 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 Commenti
Jiasen
Jiasen il 24 Ott 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 il 24 Ott 2014
Modificato: SK il 24 Ott 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]

Accedi per commentare.


Jan
Jan il 27 Set 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 Commento
Jiasen
Jiasen il 28 Set 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?

Accedi per commentare.

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by