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.

[![!\[Maximum relative value of Ax vs rank of A][1]][1]
[![enter image description here][2]][2]

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]: https://i.sstatic.net/TTcc1.png
  [2]: https://i.sstatic.net/26lLv.png