Solve linear system with $A_{i,j} = \langle e_i, e_j\rangle^2$, edges of a triangle I have three vectors in $e_i\in\mathbb{R}^3$ that form a triangle. Let us consider now the linear equation system $Ax=b$ with
$$
A_{i,j} = \langle e_i, e_j\rangle^2,\\
b_i = \langle e_i, e_i\rangle.
$$
I can solve these problems numerically, but somehow I feel I'm missing out. Perhaps the solution can even be constructed explicitly from $e_i$?
Any ideas?
 A: Just so you know, it is possible to write $P^T A P = D$ diagonal with the entries of $P$ being rational expressions in the entries of $A.$ If desired, we can take a diagonal matrix $Q$ involving square roots of the absolute values of the entries of $D$ to get $Q P^T A P Q = W,$ where $W$ is diagonal and all entries are $0,1,-1.$ Note that there is no need for $A$ to be positive, just symmetric. 
Here is the idea in symbols; one may change the order of operations if convenient.
? m = [ a,f,e; f,b,d; e,d,c]
%1 = 
[a f e]

[f b d]

[e d c]

? m - mattranspose(m)
%2 = 
[0 0 0]

[0 0 0]

[0 0 0]


? p1 = [ 1, -f/a, -e/a; 0,1,0; 0,0,1] 
%6 = 
[1 -f/a -e/a]

[0 1 0]

[0 0 1]

? m1 = mattranspose(p1) * m * p1
%7 = 
[a 0 0]

[0 (b*a - f^2)/a (d*a - e*f)/a]

[0 (d*a - e*f)/a (c*a - e^2)/a]


? p2 = [ 1,0,0; 0,1, (a * d - e * f) / (f^2 - a * b); 0,0,1]
%9 = 
[1 0 0]

[0 1 (d*a - e*f)/(-b*a + f^2)]

[0 0 1]

? m2 = mattranspose(p2) * m1 * p2
%10 = 
[a 0 0]

[0 (b*a - f^2)/a 0]

[0 0 ((-c*b + d^2)*a + (c*f^2 - 2*d*e*f + b*e^2))/(-b*a + f^2)]

? m
%11 = 
[a f e]

[f b d]

[e d c]

? matdet(m)
%12 = (c*b - d^2)*a + (-c*f^2 + 2*d*e*f - b*e^2)
? 

A: Let's first assume that the triangle is in the $z=0$-plane. After some computation SymPy spits out
(b0*c0 + b1*c1)/(a0*a0*b1*c1 - a0*a1*b0*c1 - a0*a1*b1*c0 + a1**2*b0*c0)
(a0*c0 + a1*c1)/(a0*b1*b1*c0 - a0*b0*b1*c1 - a1*b0*b1*c0 + a1*b0*b0*c1)
(a0*b0 + a1*b1)/(a0*b0*c1*c1 - a0*b1*c0*c1 - a1*b0*c0*c1 + a1*b1*c0*c0)

as a solution (where (a0, a1), (b0, b1), (c0, c1) are the edge coordinates). Clearly the numerator is the dot-product of the two "other" edges; the denominator of the first component equals (a0*b1 - a1*b0) * (a0*c1 - a1*c0). Together, this can be written as
$$
x_1 = \frac{\langle e_2, e_3\rangle}{\langle e_1\times e_2, e_1\times e_3\rangle}
$$
for the solution in the first component. (The other components likewise.)
The result holds true if the triangle is rotated in space since
$$
\langle R a, R b\rangle = \langle a, b\rangle \text{ and}\\
(R a)\times (R b) = R(a\times b)
$$
for any two vectors $a, b$ and any rotation matrix $R$.
