I have a simple (and small) linear homogeneous system $Ax=0$, where the entries of the $N\times M$ matrix $A$ are small integers. I do not need fancy methods which efficiently solve almost singular matrices and treat roundoff errors etc. But I need to know if a system is singular and also the general solution for underdetermined systems. It must be rock-solid. Any reference to such an algorithm is welcome. I am coding in C. Is SVD the way to go?
3 Answers
The Integer Matrix Library is a C library that claims to be able to compute the null space of an integer matrix. See also this answer on MathOverflow to the exact same question, which gives a list of libraries (including PARI, which also can be called from C and is still being updated). You could also take a look at LinBox, even though it's written in C++.
-
1$\begingroup$ The reference to IML was pretty useful to me. It does what I want. Also the question on MathOverflow is really the same. $\endgroup$ Commented Jan 6, 2013 at 15:06
The general solution of your underdetermined system is that $x$ is a member of the nullspace of $A$.
To find the nullspace of $A$, the most numerically stable method is to use an SVD. See Null-space of a rectangular dense matrix for further details.
The SVD will both (1) tell you what the null space of your matrix is (2) give you the pseudo-inverse of the matrix in the event that it is singular. The SVD is, generally speaking, a very handy matrix factorization and for your stated purposes should definitely get the job done.
If you're determined to write the code yourself, Trefethen and Bau's book has enough details to work from, either implementing the really bad way of finding the eigendecomposition of A*A or the proper way using bidiagonalization.
However, if you just happen to be working in C and don't feel the need to do it yourself, the Wikipedia page for the SVD has a long list of libraries that will implement it for you, chief among which is the GNU Scientific Library. Or, a google search of "SVD c source code" turns up a couple results as well.
N>M
andM>N
is both possible . Typical size could beN,M\sim 10
. I searched a bit further and found that 'reduced Row Echelon form' could be the solution. SinceA
has integer elements only rationals would appear and I could compute exact results. Anybody aware of an open implementation (any language). $\endgroup$NxM
matrices. I know that of Numerical recipies, but it is forNxN
. I understand now that SVD would work. But it seems a bit overkill to me with the additional advantage that it would use floats and probably be to long to adopt it to integers/rationals. $\endgroup$