title: Solving an overdetermined constraint system
body: <p>I have <code>n</code> real number variables (don't know, don't really care), let's call them <code>X[n]</code>.
  I also have <code>m &gt;&gt; n</code> relationships between them let's call them <code>R[m]</code>, of the form:</p>
  
  <p><code>X[i] = alpha*X[j]</code>, <code>alpha</code> is a nonzero positive real number, <code>i</code> and <code>j</code> are distinct but the <code>(i, j)</code> pair is not necessarily unique (i.e. there can be two relationships between the same variables with a different alpha factor)</p>
  
  <p>What I'm trying to do is find a set of <code>alpha</code> parameters that solve the overdetermined system in some least squares sense. The ideal solution would be to minimize the squared sum of differences between each equation parameter and it's chosen value, but I'm satisfied with the following approximation:</p>
  
  <p>If I turn the m equations into an overdetermined system of n unknowns, any pseudo-inverse based numeric solver will give me the obvious solution (all zeroes). So what I currently do is add another equation into the mix, <code>x[0] = 1</code> (actually any constant will do) and solve the generated system in the least squares sense using the Moore-Penrose pseudo-inverse. While this tries to minimize the sum of <code>(x[0] - 1)^2</code> and the square sum of <code>x[i] - alpha*x[j]</code>, I find it a good and numerically stable approximation to my problem. Here is an example:</p>
  
  <pre><code>a = 1
  a = 2*b
  b = 3*c
  a = 5*c
  </code></pre>
  
  <p>in Octave:</p>
  
  <pre><code>A = [
    1  0  0;
    1 -2  0;
    0  1 -3;
    1  0 -5;
  ]
  
  B = [1; 0; 0; 0]
  
  C = pinv(A) * B or better yet:
  C = pinv(A)(:,1)
  </code></pre>
  
  <p>Which yields the values for <code>a</code>, <code>b</code>, <code>c</code>: <code>[0.99383; 0.51235; 0.19136]</code>
  Which gives me the following (reasonable) relationships:</p>
  
  <pre><code>a = 1.9398*b
  b = 2.6774*c
  a = 5.1935*c
  </code></pre>
  
  <p>So right now I need to implement this in C / C++ / Java, and I have the following questions:</p>
  
  <p>Is there a faster method to solve my problem, or am I on the right track with generating the overdetermined system and computing the pseudo-inverse?</p>
  
  <p>My current solution requires a singular value decomposition and three matrix multiplications, which is a little much considering <code>m</code> can be 5000 or even 10000. Are there faster ways to compute the pseudo-inverse (actually, I only need the first column of it, not the entire matrix given that B is zero except for the first row) given the sparsity of the matrix (each row contains exactly two non-zero values, one of which is always one and the other is always negative)</p>
  
  <p>What math libraries would you suggest to use for this? Is LAPACK ok?</p>
  
  <p>I'm also open to any other suggestions, provided that they are numerically stable and asymptotically fast (let's say <code>k*n^2</code>, where <code>k</code> can be large).</p>