Analog of Vandermonde determinant for fitting a quadratic form?

1D interpolation: finding a polynomial satisfying $\forall_i\ p(x_i)=y_i$ can be written as a system of linear equations, having well known Vandermonde determinant: $\det=\prod_{i<j} (x_i-x_j)$. Hence, the interpolation problem is well defined as long as the system of equations is determined ($\det\neq 0$), that is equivalent with condition of having no two repeating $x$-s.

I need something analogous for quadratic form in $n$ dimensions: we would like to find symmetric matrix $A$ satisfying $\forall_k\ f(x^k)=y_k$, where $f(x)=x^T A x$, this time $x^k$ are vectors.

We get a system of linear equations: $$ \forall_k\ \sum_i A_{ii} (x^k_i)^2 + 2\sum_{i<j} A_{ij} x_i^k x_j^k =y_k$$ for $D=n(n+1)/2$ coefficients of symmetric $A$.

In analogy to interpolation problem, having values in $D$ points, we would like to find $A$. However, it requires that $\det\neq 0$ for the above set of linear equations.

Is there known a compact form for this determinant? (in analogy to Vandermonde)

If not, are there some known conditions ensuring it is nonzero - making fitting quadratic form well defined? These conditions need to contain e.g. that no two points are in one line $(x^k=a\cdot x^l)$. In what I need we can assume that all points lie on a sphere.

Specifically, my motivation is that looking at eigenspaces of adjacency matrix, we can convert the graph isomorphism problem into a question if two sets of points differ only by rotation (page 9-11 here. For strongly regular graphs these points are on a sphere and form a very regular polyhedron. Hence, I wanted to use an affine space of quadratic forms defining "wobbling" ellipsoids, such that they intersect only in our set - then we could use characteristic polynomial to test if they differ only by rotation. The crucial question is if e.g. $\{x: x^T A x=1 \textrm{ for all }A=A_0 + a\cdot A_1\}$ doesn't add too many extra points to the description. Geometrically: if "wobbling" ellipsoids with fixed some points, what extra fixed points would their intersection have?

Here is example of 2D situation: describing two points as intersection of ellipses/hyperbolas. Intersection only adds symmetric ($-x$) points, the question is when it is true, also in higher dimensions:

I am not certain that I fully answer your question. Tell me if such is the case.

Let me introduce the following notation, different from yours : let us set

$$Y_{ij}:=X_i^TAX_j \in \mathbb{R}.$$

With this notation, do we agree that the system of constraints can be written :

$$\tag{1}X^TAX=Y$$

where $X=[X_1|X_2|\cdots|X_n]$ ? (the columns of $X$ are the $X_i$s).

But (1) can be written, using Kronecker (=tensor) product $\otimes$ and operator "vec" (which means a conversion of a $n \times n$ matrix into a $n^2 \times 1$ column vector). See for that "Matrix equations" in (https://en.wikipedia.org/wiki/Kronecker_product):

$$\tag{2}\underbrace{(X^T \otimes X^T)}_{\text{known}}\underbrace{vec(A)}_{\text{unknown}}=\underbrace{vec(Y)}_{\text{known}}$$

where "the compact expression" $X^T \otimes X^T$ is a $n^2 \times n^2$ matrix .

Thus we have transformed the issue into the resolution of a linear system.

As

$$\det(X^T \otimes X^T)= \det(X^T)^{2n}, $$ (Determinant of the Kronecker Product of Two Matrices)

the condition for this system to have a unique solution is still $\det(X^T) \ne 0 \iff \det(X) \ne 0$.

A remark: one can wonder if the constraint of symmetry is taken into account: this is the case because in (2), $vec(A)$ and $vec(Y)$ correspond to symmetrical unknowns (for $A$) and symmetrical data (for $Y$).