# Accurately Computing a Positive Vector in the Nullspace of a Matrix

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.

• (1) If your matrices are singular by design, that "condition number" is meaningless; it will be just rounding errors. (2) Similarly, $\|Ax\|_\infty$ is useless because it is absolute (unscaled); plot $\|Ax\|_\infty / \|A\|_\infty\|x\|_\infty$ instead. (3) What do red and green points mean exactly? I can't understand the description. (4) What are you using as a threshold to distinguish zero and nonzero singular values? – Federico Poloni Jun 20 '20 at 8:43
• Anyway, that looks like a difficult problem; I think linear programming is going to be the best alternative in general, unless your problem has something special ($A$ isn't an M-matrix perchance?). But first you need to have a reliable way to determine the rank of $A$. – Federico Poloni Jun 20 '20 at 8:51
• Just to clarify, $max(|(Ax)_i|)$ is the maximal component of $Ax$. I am unsure if this is what you mean by $||Ax||_\inf$ but if so I will update the plots as you suggest. I have plots with the ratio of the maximum and minimal values of A instead of the singular value which I will include. Despite these being singular matrices, the determinant comes out as $\pm \inf$ using np.linalg.det, – cyfirx Jun 20 '20 at 8:56
• Yes, $\|y\|_{\infty}$ is a common notation for the maximum norm of a vector. Do not trust det; is useless for numerical computations; on moderately large matrices it will overflow. – Federico Poloni Jun 20 '20 at 9:00
• Yeah, that's what I'm thinking, it explains why there are so few red points at high rank. – cyfirx Jun 20 '20 at 12:03

• 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.
• 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).