# 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? 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$. 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, 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. 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. Jun 20 '20 at 12:03

• 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).
• Thanks for the summary. I have attempted a linear programming approach which improves the accuracy but fails for a significant number of cases and does not perform well. For this reason I have begun to investigate M-matrices as you have suggested and I think I may be able to prove my A is an M-matrix using the characterisation detailed in "R. Plemmons, Nonnegative Matrices in the Mathematical Science" p149, however I need time to check this and familiarise myself with some issues with the notation. Jun 22 '20 at 15:07
• Numerical investigations thus far have supported the assumption that they are indeed M-matrices but I am not sure how to implement the block-triangularisation (I assume this means to split the matrix into blocks and perform row operations on the entire matrix until each block is zero above or below it's diagonal) and then extract the nullspace vector. Will I find a library function to do this eg scipy or will I need to write my own? Jun 22 '20 at 15:18
• I am glad to hear that the M-matrix tip helped! No, that is a different kind of reduction: you reorder column and row indices symmetrically so that your matrix becomes block triangular with irreducible blocks. Essentially, you wish to split the directed graph associated with the matrix into strongly connected components. This is purely a symbolic operation, without actual floating-point computations in it. Jun 22 '20 at 15:29
• There is probably something inside the scipy graph libraries, but I am not sure exactly where to look. Also, it is possible that this comes from free from the structure of your problem. Jun 22 '20 at 15:31
• Following your advice I have implemented something similar to the GTH algorithm to solve my problem nicely. Thank you very much! Jun 27 '20 at 16:49