Accurately Computing a Positive Vector in the Nullspace of a Matrix

Question

I'm sure this question has been asked before yet after many hours of searching I am unable to find a definitive answer.

The problem at hand is solving the linear system: $$A \mathbf{x} = \mathbf{0}$$ with the constraint: $$ x_i \geq 0 $$ At present I apply SVD to obtain a set of vectors which span the nullspace of the matrix. The threshold used for zero singular values is max(1e-12, 1e-15 x $\sigma_0$), where $\sigma_0$ is the greatest singular value. From these I search for a basis vector with the components all of the same sign. This is clearly suboptimal as any solutions which are linear combinations of the nullspace vectors are will be missed but I am unsure how to solve the resulting system of inequalities (possibly linear programming).

However a much more significant is issue dealing with rounding errors. The SVD method worked well for small matrices but as they increase in size (currently 64x64) it seems to result in vectors which are not in the nullspace, or with huge rounding errors when I come to verify they are solutions. Below I have plotted the maximum absolute value of the result of plugging a candidate solution back into the equation against the condition number and rank of the matrix. This issue becomes worse for lower ranks of A.

The value of $\mathbf{x}$ is used to compute a function of the matrix and solution which can also be obtained via a prohibitively expensive alternative to finding positive nullspace vector. This is possible because x is normalised to a probability vector which can be obtained via Monte Carlo simulation. These probabilities are then used to calculate a single value which is a function of $\mathbf{x}$ and $A$. Points in green are in agreement with the Monte Carlo simulation and red are not.

I am sure there must be some way of solving this problem for matrices of this size. In general my matrix is ~90% zeros so sparse matrix methods may be suitable, particularly as the sparsity will increase as I increase the size of the matrices.

I would like some advice on how to proceed with the problems presented here, primarily the issue of accurately computing nullspace vectors and secondarily, solving a linear system of inequalities and if sparse methods may be helpful.

Answer 1

Quick answer to summarize my comments.

• Keep in mind that a delicate point is the choice of the truncation threshold in the SVD (what is "numerically zero" and what is not). If you do not see a clear drop in the singular values, then it means your precision is insufficient to identify zeros.
• Since $\|Ax\|_\infty / \|A\|_\infty \|x\|_\infty$ is of the order of machine precision, your plots show that the vectors you computed are (numerically) in the kernel of $A$, so the numerical method seems to be working correctly.
• A possible source for the discrepancy that you observe is that this problem has multiple solutions. Are you sure that you are computing the correct one?
• If $A$ (or $-A$) is an M-matrix, then there might be better solutions to your problem: irreducible M-matrices have a kernel of dimension $1$ spanned by a positive vector; so you just need to reduce your matrix to irreducible components (i.e., block-triangularize it) and compute the vector in the kernel of each of the singular diagonal blocks (which is unique up to scaling).