There was a great question on the newsgroup this past week asking how to determine if a system of equations had a solution. The poster wasn't (at least yet) concerned with what the solution was.
The first recommendation was to use det, the determinant. While this is pedagogically correct for some cases, it is insufficient since it doesn't correctly account for non-singular systems which do have solutions. John D'Errico followed up with some examples and some smart math.
John points out that using det would not give the correct answer for this problem.
A = ones(2); b = [2;2];
Does a solution exist for A*x = b? It does despite the singularity of A. det(A) is zero of course.
ans = 0
Yet the answer is just x = [1;1]. Find it using pinv.
ans = 1 1
Using rank, check to see if the rank([A,b]) == rank(A)
rank([A,b]) == rank(A)
ans = 1
If the result is true, then a solution exists.
Let's try it for a problem that has no solution.
c = [1;2]; rank([A,c]) == rank(A)
ans = 0
I had a previous post about collinearity in which we tried a bunch of solutions, and the superior solution again was rank vs. det. Did you know that some singular systems had valid solutions? Have you had similar instances where what the textbooks sometimes recommend isn't ideal for the full spectrum of possible conditions. Post your thoughts here.
Get the MATLAB code
Published with MATLAB® 7.9
6 CommentsOldest to Newest
When the condition rank([A,c]) == rank(A) is satisfied, the linear system has a unique solution. When the condition is not satisfied, it is not that the linear system has no solution but has infinite solution (i.e. no unique solution), which is called ill posed.
If you want a solution of an ill posed linear system, you should add extra constraint. x=pinv(A)*b is always a solution of the linear system. Particularly, it is a solution that has minimal l2 norm.
correction to the previous post: solution of linear system Ax=b:
rank([A,b])~=rank(A): no solution.
rank([A,b])=rank(A) & size(A,1)size(A,2): infinite solutions.
when the linear system has solution, x=pinv(A)b gives a solution that has minimal l2 norm.
When the linear system has no solution. x=pinv(A)b gives a solution that has the least of square error.
In my opinion, the problem is not that textbook recommendations are not ideal for the full spectrum of possible conditions. Unless the textbook is written by someone not familiar with the field (in which case it shouldn’t be used in the first place), the problem is usually caused by not exercising enough care when applying its concepts. In this particular case the problem has been caused by reading the converse of a theorem as true when that is not the case. The logical statement relevant in this context is “if A is a non-singular matrix (i.e., the determinant of A is non-zero) then Ax = b is consistent and in fact the solution is unique.” Clearly, checking for consistency of a square system of linear equations by checking the determinant is using the converse of this theorem, which does not hold.
Another issue is the difference between a concept that is theoretically (“continuously”) true but might collapse when applied numerically (“discretely”). Ill conditioning for systems of linear equations or eigenvalue problems stand out here.
Another issue that comes to mind is computational efficiency (i.e., speed). An example would be solving a “large,” square, consistent, system of linear equations using Cramer’s rule (with a determinant to compute for the each component of the solution), which, in Gilbert Strang’s words, “is a disastrous way to go because to compute these determinants takes like, approximately, forever.”
I should just point out why rank works well here. The pair of calls to rank identify when the right hand side for the system lies in the column space of A. Thus, the code fragment
rank(A) == rank([A,b])
tests to see if b can be represented by some linear combination of the columns of A. This works regardless of the nature of A, singular or not, or even if the problem is under-determined.
For example, change the problem slightly so that
A = ones(2,3); b = ones(2,1);
Here, a solution exists. In fact, infinitely many solutions exist. We get one of them, the minimum norm solution, from pinv.
x = pinv(A)*b x = 0.33333 0.33333 0.33333
But or course, the true solution space is vastly larger than a single point for this problem. We can add any linear combination of the vectors spanning the null space of A to this solution to get an equally valid one.
V = null(A) V = -0.8165 -2.216e-16 0.40825 -0.70711 0.40825 0.70711 x = pinv(A)*b + V(:,1) x = -0.48316 0.74158 0.74158 A*x ans = 1 1
The problem is if the right hand side of the system fails to lie in the column space of A. The call to rank identifies this failure. And while pinv will return a minimum norm result, it will still not be a valid solution to the system.
c = [1;2]; rank([A,c]) == rank(A) ans = 0 x = pinv(A)*c x = 0.5 0.5 0.5 A*x ans = 1.5 1.5
An under-determined system does not assure infinitely many solutions. We may still have no solution at all.
As for the use of determinant to flag the singularity of a system, it is often SLOWER in MATLAB to compute the rank of a system, at least if the system is at all significant in size.
A = rand(100); tic,rank(A);toc Elapsed time is 0.016543 seconds. tic,det(A);toc Elapsed time is 0.003492 seconds.
The reason is that det in MATLAB uses an LU decomposition of A, whereas rank uses the SVD under the hood. The SVD takes considerably more work to do here. Of course, if one actually used the old formulas found in textbooks to compute det, that function would truly be incredibly slow. The problem in MATLAB is not really so much that det is slow as that it is terribly inaccurate as an indicator of singularity. Consider this example:
A = rand(100,99); A = [A,A*rand(99,1)]; rank(A) ans = 99 det(A) ans = 8.492e+10
Obviously A is singular, since I made the last column up as a linear combination of the first 99 columns. The call to rank shows this to be true. Even so, the determinant seems pretty large at first glance. One might almost look at that number and believe that A was not singular. On the other hand, if we wished to make A look singular, it takes only an inconsequential change to force det to yield a near zero result!
det(A/2) ans = 6.699e-20
For the case of multivariate Gaussian random variables, sometimes the solution (to, for example, an estimation problem via Kalman filtering) lies on a subset of the distribution (for the 2-dimensional case the covariance matrix would give the equation for a line segment rather than an ellipse). In the case of multivariate Gauss characteristic functions are much better to work with than distribution functions for the reason that they can solve the singular case.
If you’re working with medium-sized integer or rational arrays, check out the Fractions Toolbox http://www.mathworks.com/matlabcentral/fileexchange/24878-fractions-toolbox which does exact arithmetic on fractions. One useful feature is a 2-output version of mldivide (\) that can help diagnose indeterminate and inconsistent systems; also I’ve included the rank function in today’s update.