Override constraint system solution

Question

I have real number n variables (I don’t know, I don’t really care), call them X[n] . I also have m >> n relationships between them, let's call them R[m] , of the form:

X[i] = alpha*X[j] , alpha is a nonzero positive real number, i and j different, but the pair (i, j) not necessarily unique (that is, between the same variables with different alpha -factor)

What I'm trying to do is find a set of alpha parameters that solve the redefined system in the least squares sense. The ideal solution would be to minimize the square of the sum of the differences between each parameter of the equation and the selected value, but I am satisfied with the following approximation:

If I turn m equations into an overridden system of n unknowns, any pseudo-reversible numerical solver will give me the obvious solution (all zeros). So, what I'm doing now is to add another equation to the mix, x[0] = 1 (virtually any constant will do) and solve the generated system in the sense of least squares using the Moore-Penrose pseudoinverse. Although this attempts to minimize the sum (x[0] - 1)^2 and the square sum x[i] - alpha*x[j] , I find this a good and numerically stable approximation to my problem. Here is an example:

 a = 1 a = 2*b b = 3*c a = 5*c 

in octave:

 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) 

What gives the values ​​for a , b , c : [0.99383; 0.51235; 0.19136] [0.99383; 0.51235; 0.19136] [0.99383; 0.51235; 0.19136] This gives me the following (reasonable) relationship:

 a = 1.9398*b b = 2.6774*c a = 5.1935*c 

So right now I need to implement this in C / C ++ / Java, and I have the following questions:

Is there a faster method to solve my problem, or am I on the right track with creating an overridden system and calculating a pseudo-inverse?

My current solution requires decomposing a singular value and three matrix multiplications, which depends a little on the fact that m can be 5000 or even 10000. Are there faster methods for calculating a pseudo-inverse (in fact, I only need the first one with its column, and not the whole matrix , given that B is zero, except for the first row), given the sparseness of the matrix (each row contains exactly two nonzero values, one of which is always one, and the other is always negative)

What math libraries do you propose to use for this? Is LAPACK OK?

I am also open to any other proposals, provided that they are numerically stable and asymptotically fast (say k*n^2 , where k can be large).

Answer 1

The SVD approach is numerically very stable, but not very fast. If you use SVD, LAPACK is a good library to use. If this is just a one-time calculation, then it is probably fast enough.

If you need a much faster algorithm, you may have to sacrifice stability. One possibility would be to use QR factorization. You will need to read about it to see the details, but some of the reasoning goes as follows. If AP = QR (where P is the permutation matrix, Q is the orthogonal matrix, and R is the triangular matrix) is the QR decomposition of economy A, then the equation AX = B becomes QRP ^ {- 1} X = B and the solution X = PR ^ {- 1} Q ^ T B. The following Octave code illustrates this using the same A and B as in your code.

 [Q,R,P] = qr(A,0) C(P) = R \ (Q' * B) 

The nice thing is that you can use sparse A by doing sparse QR decomposition. There is some explanation in the Octave help for the qr function, but it did not work for me right away.

Even faster (but still less stable) you should use the normal equations: if AX = B, then A ^ TAX = A ^ T B. The matrix A ^ TA is a square matrix (hopefully) of full rank, so you can use any solver for linear equations. Octave code:

 C = (A' * A) \ (A' * B) 

Again, sparseness can be used in this approach. There are many methods and libraries for solving rare linear systems; popular seems to be UMFPACK .

Added later: I do not know enough about this field for quantification. Entire books are written on this. Perhaps QR is about 3 or 5 faster than SVD and normal equations are twice as fast. The effect on numerical stability depends on your matrix A. Rare algorithms can be much faster (say, by a factor of m), but their computational cost and numerical stability are largely dependent on the problem, as it is sometimes not very clear.

In your use case, my recommendation would be to try to calculate the solution using SVD, see how long it takes, and if that is acceptable, just use this (I think it will be about a minute for n = 1000 and m = 10000) . If you want to study it further, try also QR and normal equations and see how fast they are and how accurate; if they give roughly the same solution as SVD, then you can be pretty sure that they are accurate enough for your purposes. Only if they are too slow, and you want to hurry for a while, look at sparse algorithms.

Answer 2

Your problem is incorrect. If you consider the problem as a function of n variables (the smallest squared difference), the function has exactly one global minimum hyperplane.

This global minimum will always contain zero unless you correct one of the variables other than zero, or reduce the scope of functions in any other way.

If you need a para-therapy of a hyperplane solution, you can get this from the Moore-Penrose pseudo-inversion http://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse and check the section on getting all the solutions.

(Note that I used the word “hyperplane” in a technically incorrect way. I mean a “flat” unlimited subset of your parameter space ... A row, a plane, that can be parameterized by one or more vectors. For some reason I can not find a common noun for such objects)