System of linear equations, resulting from a weighted graph. How to solve this numerically?

We have a problem that leads to a system of linear equations which has to be solved numerically. There are thousands of algorithms to solve linear equations, but I haven't found any that fits our special requirements.

I've tried to describe these requirements in the pdf linked below: http://pdfcast.org/pdf/linalgproblem (Document has three pages)

Since our problem is probably a common one, we believe that there are algorithms specially designed for the kind of matrices we are facing.

Do you know someone who had the same problem? Do you know a fitting algorithm? Do you know an implementation in c++? That would help us a lot!

PS: The matrices to solve will have about 1000-10000 rows and 5000-20000 lines.

EDIT: I still hope to hear about a more specific recommendation than to use an algorithm for sparse matrices. I believe I'm not the first one having to solve this problem.

-
This is not MathOverflow, and it's not a bother. But solving systems of linear equations is the most basic problem in linear algebra, and the standard software packages (Maple, Mathematica, Matlab, ...) will generally do a good job of it. They will also do least-squares fits. Is there some reason you can't use those? – Robert Israel Aug 28 '11 at 23:57
Sorry, I had this question on math overflow before, until someone recommended me this site. We need to integrate the calculation into some software, that's why Maple is no good option. Still, I should propably check it out to get a feeling about what's possible. – Konstantin Aug 28 '11 at 23:59
Have you looked into OpenMaple? It is possible to call Maple code from a compiled C, Java or VisualBasic program. Of course, the Maple software must be available on the system running the program. – Robert Israel Aug 29 '11 at 0:11
That sounds like an option. Still I would prefer to know the algorithm that's designed for this problem (I'm pretty sure it exists), and I would like to implement it by my own, for I know what happens. – Konstantin Aug 29 '11 at 0:18
Looks like a large sparse problem. There's UMFPACK and a number of other libraries for sparse problems... – J. M. Aug 29 '11 at 1:04

1. As you note, your system is likely to be (at least partially) overdetermined, and the data is likely to have errors. This makes a big difference in the kind of algorithm you're going to need.
2. You mention that consecutive routes are likely to differ by one node, e.g. A, AB, ABC, ABCD, etc. Are the measurement errors for such routes likely to be correlated? If so, you'll also need to take this into account.

For the first point, you might try formulating your problem as something like

\begin{aligned} d(A,B) + d(B,E) &= 3.5 + e_1 \\ d(A,B) + d(B,C) + d(C,N) &= 8.6 + e_2 \end{aligned}

etc., and try to minimize the sum of the squares of the error terms $e_i$ (possibly with weighting factors to account for varying measurement accuracy) subject to the given constraints. I'm not familiar enough with these kinds of linear fitting problems to recommend a specific algorithm or implementation, but I believe such do exists.

As for the correlations, in the simple case where $E[e_{i+1}|e_1,\dotsc,e_i] = e_i$ within a given chain of routes, (i.e. where errors within a chain are simply cumulative), there is fortunately an easy way to decorrelate the data: just subtract the previous route from the next within each chain. This will also conveniently simplify your data, giving you lots of independent single-edge measurements instead of lots of overlapping multi-edge measurements.

-
Thanks for reading it :-) As you point out correctly, the errors for consecutive routes will be correlated. But also, the error depends a lot on the last node, thus the correlation is not perfect. But for a first approach, it may be helpful to pre-calculate the differneces for consecutive routes, yes.- – Konstantin Sep 1 '11 at 12:40
Your suggestion to handle the overdetermination is interesting.--But as you point out, there is probably somebody who had the problem already(who already dealed with the overdeterminition), and solved it. And I don't feel like implementing something that has already been done by sb with much more knowledge. :-) – Konstantin Sep 1 '11 at 12:48

I think I have solved the problem. (Though, we haven't tested it ye.)

The answer is Singular Value Decomposition. It gives us the possibility to discover column degenracies in our equation. We can remove the unsolvable variables from the equation and then calculate the least square solution for the rest of the unknowns.

When I look at the algorithm, I feel pretty sure it will also be possible optimize it for integer values.

I've found all this in the book "Numerical Recipies, Second Edition, Cambridge University Press."

-
Well, the problem with SVD is that if your underlying matrices are large and sparse, the orthogonal factors are large and dense... of course, there are sparse adaptations of SVD, but I haven't had experience with them. – J. M. Oct 25 '11 at 23:31