It must be observed that the system eq2 may be written in the equivalent form:

eq2 = {
u1 - z1 == 0,
u2 - z2 == 0,
u3 - z3 == -1,
a1 u1 + a2 u2 + a3 u3 == -a3,
b1 u1 + b2 u2 + b3 u3 == -b3
}

The matrix of this system of 5 linear equations and 6 variables (u1, u2, u3, z1, z2, z3) is:

{{1, 0, 0, -1, 0, 0}, {0, 1, 0, 0, -1, 0}, {0, 0, 1, 0, 0, -1}, {a1, 
   a2, a3, 0, 0, 0}, {b1, b2, b3, 0, 0, 0}} // MatrixForm

The rank of this matrix is equal to 3! Conclusion: there are no redundant equation in the system obtained after Eliminate!

Replacement in eq1 is not an equivalent transform with adding 3 independent equations (see 3 first equations in eq2).

So, all is correct!

Yes! a and b are regarded as the coefficients!

I started from initial interpretation that a and b are two given vectors, i.e. their coordinates are fixed. But, if we think about coordinates of a and b as variables, the equations are bi-linear that makes equations somewhat simpler for Eliminate[].

There are textbook references for Groebner bases, but even those are far from elementary. The one by Cox, Little, and O'Shea is perhaps most accessible, which still means it takes time and some amount of undergraduate math prerequisite. Possibly a decent place to look would be early articles by Bruno Buchberger (who named them for his thesis advisor).

In very brief terms, Groebner bases generalize both Gaussian elimination (the linear case) and the Euclidean algorithm (the univariate case). Given a set of polynomials they come up with another set generating the same "polynomial space" (ideal, in the language of algebra) but with a certain nice property: any element in the space can be "reduced" to zero, by a form of generalized division by the Groebner basis. This property will not in general hold for the original set. A set with this property gives a nice algorithmic way of determining membership (and that was the original purpose of Buchberger's work, 50+ years ago).

These bases are with respect to what are known as "term orders". It happens that certain term orders can be used for elimination of variables (meaning: the result is a Groebner basis for the subset of polynomials in the original set that do not contain the eliminated variables).

I realize this is probably too brief a description to be helpful. But if one wants to know what elimination of variables is doing, one has to deal with a bit of hard-core algebra (I just hope "deal with" is not "suffer with").

Hello Daniel, I am busily reading Cox, Little, and O'Shea. It is hard work, but well spent and not only "digestible" but even delicious. Here is what I found so far ( hopefully I understood everything o.k. and am not barking up the wrong tree):

Lets start with the equation from Eliminate (from above):

a3+a1 u1+a2 u2+a3 u3==0 &&
b3+b1 u1+b2 u2+b3 u3==0 &&
a1 (b3+b2 u2+b3 u3)==b1 (a3+a2 u2+a3 u3)

I claim that this set is linear dependent. As you pointed out, I wrongly claimed that the third equation is dependent on the first two. Because you said this, I stopped thinking (intimidation-by-expert-spell!) and it took me some time to get over it and realize that the set IS linear dependent. You can e.g. express the first one by the remaining ones:

eq= {a3+a1 u1+a2 u2+a3 u3,
b3+b1 u1+b2 u2+b3 u3,
a1 (b3+b2 u2+b3 u3)-(b1 (a3+a2 u2+a3 u3)) };

{a1/b1,-1/b1} . eq[[2;;3]] == eq[[1]] //Simplify
True

Well, why does "Eliminate " give linear dependent answers? Let's look what "Groebner Basis" has to tell. We start with eq2, rewritten from equation to polynomial form and calculate the Groebner basis (using lexicographical ordering. Remember z and u are variables, a and b parameters):

eq2a={
z1-u1,
z2-u2,
z3-1-u3,
a1 u1+a2 u2+a3 (1+u3),
b1 u1+b2 u2+b3 (1+u3)};

vars={z1,z2,z3,u1,u2,u3};
GroebnerBasis[eq2a,vars]

{-a3 b1+a1 b3-a2 b1 u2+a1 b2 u2-a3 b1 u3+a1 b3 u3,
b3+b1 u1+b2 u2+b3 u3,
a3+a1 u1+a2 u2+a3 u3,
-1-u3+z3,
-u2+z2,
-u1+z1}

I claim that this is not a reduced Groebner basis and therefor not minimal. To be a reduced basis, no monomial term of any basis polynomial can lie in the affine variety of the leading terms of the remaining polynomials. (Whow quite a mouth full!)

We need only look at the leading terms. The leading terms are:
{(-a2 b1+a1 b2) u2, b1 u1, a1 u1, z3, z2,z1}

it is quite obvious that the second and third term are dependent. I think this is the reason that "Eliminate" gives linear dependent solution sets.

Daniel, thanks a lot. The hint about "CoefficientDomain" really helps, I do not think I would have realized this myself.

Now the problem with "Eliminate". Is it too naive an approach to simply test the output of Eliminate for linear dependencies and eliminate if necessary. Or is the problem more subtle?

You are right, I must have made a mistake when checking for dependency. However, now it is even worse,the third equation is additional information. Where does it coming from?