What to do when solving a system with a sparse matrix and a dense matrix results in a singular matrix, even after normalization?

The system I have consists of a matrix A of size 401*401, sparse doubles (tridiagonal), along with a dense matrix B of size 401*1. Initially, without normalising parameters of the system (which were in ranges of 10e-7 to 10e4), when I got singluar matrix as a result of x = A\B (mldivide), I thought normalising would make it function correctly. But it doesn't. Also, the B matrix is actually derived from multiplying another sparse matrix (tridiagonal) of size 401*401 with a dense matrix. Can anyone suggest what options do I have?

If it helps, I am trying to make a 1D Finite Element (FEM) program for drift-diffusion model of semiconductor. The exact error message is: "Matrix is singular to working precision."

[Torsten]

Check your FEM discretization for errors.

[Chaitanya Singh]

@Torsten Yes, I am doing that right now.

Answer 1

[John D'Errico]

Normalizing the matrix in some way will not help. In fact, it might invalidate the FEM solution, depending on what you did. Since you were to vague for us to know what you tried to do, it is irelevant, since normalization is not going to help anyway.

You want to solve the problem A*x=B. That A is sparse and B is dense is also irrelevant. It is litlerally meaningless to describe a vector as sparse or dense anyway. And how you derived B is again, COMPLETELY irrelevant, since the singularity is entirely in A. B IS IRRELEVANT.

Your issue is that A is singular. Nothing more, nothing less. Is A only numerically singular, and if you could use infinite precision, it would not be? That we cannot know yet. However, far more likely is that A is singular by design. In fact, you probably made it that way, without realizing you built that singularity into A.

And that A is sparse and only of size 401x401 is silly, since you will not gain anything there. A 401x401 solve is going to be blazingly fast anyway, unless you are doing it many thousands of times.

What should you do to diagnose why A is singular? I would look very carefully at the singular values of A. Are all but one of them significantly non-zero, and then there is one singular value of A that is zero to within floating point trash? (Divide it by the maximum singular value. If the ratio is less than roughly 1e-14, then the matrix will be deemed to be numerically singular.) If that is the case, with just one essentially zero singular value, then it suggests that you made a mistake in formulating the problem.

So what might you have done wrong? A common use of finite element modeling is to build physical models of things like a bridge truss. It will allow you to compute the stresses and strains of the truss, but if you fail to pin down one end of the truss, then the entire truss can translate in any direction without any energy penalty. The result will be a singular matrix, with one or perhaps several zero singular values. This is a common error made by those who are using FEM, and it is an easy mistake to make.

So first, make A a full matrix. Then compute the singular values. (SVD requires your matrix is full, so to compute the singular values, you need it that way.) Look carefully at what comes out. Then go back and loook at how you formulated the FEM problem. Have you sufficiently constrained the problem so that there are no effectively zero cost translations or rotations available to the system?

[Chaitanya Singh]

Wow. Thanks for your reply. I can't believe how naive I was. I will certainly look into the issues you mentioned.