MATLAB: Solving a system of equation (Ax=b) with a block of zeros

dgesv

Dear All, I recently converted my FORTRAN code for contact-impact problems to MATLAB code. IN any contact problem, I have to solve a system of equation Ax=b for the contact node. I am using the Lagrange multiplier method (LMM) to enforce the contact constraints. Now, due to the way LMM method is constructed, a block of zeros is obtained in A which also leads to some of the diagonal entries of A becoming zero. I cannot use x=A\b since I have zeros on the diagonal. To overcome this, I use a perturbed Lagrangian approach where a small positive value (1e-11 or 1e-10) is added to the diagonal so that now the invention is possible. However, still, the condition number is very bad. When I was using FORTRAN version of the MATLAB code, I used DGESV without any problem. So, what should I do. I see there is https://in.mathworks.com/matlabcentral/fileexchange/8006-lapack-function-dgesv-in-matlab but I have no prior experience in mex file generation.

Any help in this regard would be highly appreciated.

PS: For matrix A the size is (ID+ROWS_Q1) x (ID+ROWS_Q1) where A(1:ID,1:ID) entires are all zeros.

Best Answer

This is simply not true that having zeros on the diagonal prevents you from using A\b.

A= rand(3);
A = A - diag(diag(A))
A =
            0      0.19566      0.10139
       0.1568            0      0.34752
      0.24236      0.54753            0

As you can see, the entire diagonal is zero, identically so. Yet A is full rank.

rank(A)
ans =
     3

The condition number is even quite low.

cond(A)
ans =
       5.9645
b = rand(3,1);
x = A\b
x =
       2.9739
    -0.010528
       1.1608

So A\b is quite well posed.

Likewise, having a block of zeros is not a problem, unless the block is seriously large.

A = rand(10,10);
A(1:3,1:3) = 0;
rank(A)
ans =
    10

So if we go overboard, then of course things will fail.

A(1:6,1:6) = 0;
rank(A)
ans =
     8

My guess is you need to use a tool like pinv. It is appropriate, since you are essentially trying to regularize the problam anyway. In that case, just use

x = pinv(A)*b;

There is no need to hack with A by adding tiny numbers to the diagonal.

Note that DGESV (see this link for doc) does nothing special. It just computes a partially pivoted LU decomposition to solve a linear system. So if you were getting results from a rank deficient matrix, they were probably garbage anyway.

Related Solutions

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

What do you mean by "tolerance"?

A QR solution has no need for a tolerance, nor does linsolve, as they are not iterative methods.

Anyway, you cannot specify the desired accuracy of a solution in a linear system like this, solved using QR. Just wanting a more accurate solution is not relevant. The solution exists, and is unique, or there are infinitely many solutions, all of which are equally good (or poor), or there is no exact solution, and that which qr/linsolve returns is entirely valid.

Perhaps you have a singular system, or one that is nearly singular. In that case, a tolerance is still not really relevant. Yes, you could decide which singular values to discard, in case of a pinv style solution. But that is not a "tolerance" in the way tolerances are traditionally used. And if your system is singular, then again, tolerances have no real meaning.

As far as your matrix being non-square, that just means (assuming you have an over-determined system) that you will never be able to achieve an exact solution. (Well, not entirely true. It IS possible that an over-determined system has an exact solution. But that is terribly rare.) So if no exact solution exists, then specifying a tolerance is again a meaningless goal.

If your problem is non-square in the sense that it is under-determined, then again, a tolerance is irrelevant. A solution may exist, or not. If one solution does exist, then there will be infinitely many solutions, all equally good. A tolerance is meaningless here.

I think that perhaps you are being confused, thinking that the solution to a linear system is like that of a nonlinear system, where you will employ convergence tolerances, etc. Or, perhaps you think that simply by cranking down on a tolerance, you can always gain an arbitrarily good solution. (Many people seem to think that. They are flat out wrong.) Or perhaps you are thinking that a solution like that provided by linsolve/qr is like an iterative solver, like that from lsqr.

If this has not cleared things up, then I will suggest you need to explain CLEARLY what you want to do, and why you think that a tolerance is a meaningful thing to ask for here. As well, you need to explain why it is that you think what I have said above does not apply.

MATLAB: Is det(x) better than rcond(x) in determining non-singularity here

Let me suggest that you misunderstand det, and why it is a terribly poor tool to use to learn if a matrix is singular or not. Yes, we are taught to use det by teachers, in school. Hey it works there, so why not use it in practice? The problem is those same teachers stopped short of explaining the real issues. I've often said that the use of det is a bad meme, one that just keeps propagating forever. One person learns about it, so they teach others. That it is a bad idea? Nobody ever told them. So det lives on - IT'S ALIVE! (With my apologies to the movie where that line was used. Thanks Mel.)

a = -1.3589e-10;
b = -1.7108e-9;
c = 12.5893;
d = -1e11;
x = [a b; c d];

So, is the 2x2 matrix singular or not? It is NUMERICALLY singular, thus in floating point double precision arithmetic, it can be risky to use a numerically singular matrix to solve a system of linear equations reliably, even though an inverse exists in theory.

A big problem with det is that it is sensitive to scaling. What result would indicate to you that a matrix is singular based on the determinant? Remember that a numerically computed determinant will virtually NEVER be exactly zero. For example, consider this matrix, A:

A = rand(4,5);
A = [A;sum(A,1)];
A
A =
     0.031833      0.82346     0.034446       0.7952      0.64631
      0.27692      0.69483      0.43874      0.18687      0.70936
     0.046171       0.3171      0.38156      0.48976      0.75469
     0.097132      0.95022      0.76552      0.44559      0.27603
      0.45206       2.7856       1.6203       1.9174       2.3864
      
det(A)
ans =
   3.0549e-18

Is A singular? It darn tootin is! The last row of A is the sum of the other 4 rows. But wait a second, det(A) is not returned as zero. Wwwwhhhhaaaatttt happened?

And, if you claim that 3e-18 is close enough to zero, then just multiply A by 10000.

det(10000*A)
ans =
       5005.1

Gosh, 10000*A is surely not singular based on what det just told us. But I thought that we decided A was singular? Mathematically, A is indeed singular. We know that from theory. Theory would never lie, would it?

Or try this:

A = eye(1200);
det(A)
ans =
     1

Hey, great! Det knows that an identity matrix is not singular! But then, what about this?

det(2*A)
ans =
   Inf
det(A/2)
ans =
     0

Hmm. While I am sure you know that A is remarkably well-posed for solving ANY system of equations, or computing a matrix inverse, for that matter, the same applies to A/2 as well as 2*A. So why does det tell us such a confusing story on such simple matrices?

Or, how about this matrix? Lets see, I think I recall that magic squares of even order are all singular matrices. It is certainly true for magic(4), at least.

A = magic(4)
A =
    16     2     3    13
     5    11    10     8
     9     7     6    12
     4    14    15     1
     
rank(A)
ans =
     3
     
cond(A)
ans =
    8.148e+16

So rank and cond both think that A is singular. Composed of purely integers, the determinant must also be an exact integer. The symbolic toolbox can do it for us exactly.

det(sym(A))
ans =
0
det(A)
ans =
  -1.4495e-12

But, det fails here on A itself, since large scale determinants cannot be efficiently computed using any routine in a reasonable amount of time. So, what det actually does is to use linear algebra to compute the determinant. So, if we were to look at the internal code for det, it would look something like this:

[~,U] = lu(A);
prod(diag(U))
ans =
  -1.4495e-12

Here we can compute the determinant of a matrix using no more than a call to LU. And the LU factorization is an O(n^3) operation in terms of flops, whereas computing the determinant using traditional methods as you were taught in school is an O(factorial(n)) operation. So det probably uses the high school formula for a 2x2 or 3x3 determinant as special cases. But it should use LU for matrices of size 4x4 or larger. The end result is that det is not something you can even trust to compute the determinant you so very much desire - to see if your matrix is singular.

The problem is that det is NOT indeed a good tool to identify if a matrix is well-posed for such a problem. Yes, again, I know that your teacher said, that mathematically, if a matrix is singular, then the determinant is zero. But in fact, det is a bad tool to identify when a matrix is singular, because we cannot trust what it reports. Just because a tool works in theory, does not make it a good tool for practical use. And here, practical use involves floating point arithmetic. So just bcause a computed determinant is near zero or is not near zero, it still does not tell you if the matrix was actually singular or not.

Instead, we tell you rely on tools like cond, rcond, svd, and rank to identify singular matrices and avoid det like the plague on mathematics it rightfully is. (There are some circumstances where it is indeed arguably appropiate to use det. However testing for singularity is a problem where you should essentially NEVER use det.)

Best Answer

Related Solutions

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

MATLAB: Is det(x) better than rcond(x) in determining non-singularity here

Related Question