I know that in general, multivariate polynomial satisfiability is equivalent to 3-SAT; however, I'm wondering if there are any good techniques in the quadratic case, specifically if there is a polynomial time solution.
I guess the more general question would be, are there any classes of multivariate polynomials for which the satisfiability problem is efficiently solvable?
You can decide if a quadratic polynomial $p: \mathbb{R}^n \rightarrow \mathbb{R}$ has real roots with some linear algebra. As you note, the general case should be hard.
Observe first that $p(x) \neq 0$ for all $x \in \mathbb{R}^n$ if either $p(x) > 0$ or $p(x) < 0$ for all $x$ (this follows by continuity). So it is enough to be able to decide if $p(x) > 0$ for all $x$. In general this is related to complexity-theoretic versions of Hilbert's 17th problem: a polynomial $p(x)$ is positive over the reals if and only if you can write $p$ as the sum of squares of rational functions and a positive constant $c$ (this is a theorem by Artin). Finding this decomposition or solving the decision problem is most likely hard in general, but the quadratic case is easy, because of the magic of the spectral theorem. For more information about the general case, look at Devanur,Lipton,Vishnoi, and Monique Laurent's survey.
Let us write $p(x) = p_2(x) + p_1(x) + c$ where $p_2$ is homogeneous of degree 2, $p_1$ is linear, and $c$ is a constant. Let us define $q(x_0,x) = p_2(x) + x_0p_1(x) + cx_0^2$ to be the homogenization of $p$, where $x_0$ is an additional variable.
Claim. $p(x) > 0 \Leftrightarrow \forall x_0 \neq 0: q(x_0,x) > 0$
The "if" direction is easy. In the non-trivial direction, assume $p(x) > 0$ for all $x$ and assume $x_0 \neq 0$: $$ q(x_0,x) = x_0^2q\left(1,\frac{x}{x_0}\right) = x_0^2 p\left(\frac{x}{x_0}\right) > 0. $$ QED
Notice also that, because $q(x_0, x)$ is continuous, if $p(x) > 0$ for all $x$ then $q(x_0,x) \geq 0$ for all $(x_0,x)$ (including $x_0 = 0$).
Since $q$ is homogeneous, we can write $q(x_0,x) = y^TQy$, where $Q$ is a symmetric matrix and $y = (x_0,x)$. By the above, if $p(x) > 0$ for all $x$, then $Q$ is positive semidefinite. Moreover, $q(x_0, x) > 0$ for all $x_0 \neq 0$ if and only if the kernel of $Q$ is a subset of the hyperplane $\{y = (0,x): x \in \mathbb{R}^n\}$. Both these conditions can be decided in polynomial time by computing the SVD of $Q$.
Here is an approach that will work for some polynomials, but is not guaranteed to work for all of them.
It is guaranteed to work for all multilinear quadratic polynomials. It is also guaranteed to work if there is any variable $x$ such that $p(\cdots,x,\cdots)$ does not contain a $x^2$ term. However, for polynomials of the form $p(x,\cdots,z) = c_x \cdot x^2 + \cdots + c_z \cdot z^2 + p'(x,\cdots,z)$ where $p'$ is a multilinear quadratic polynomial (i.e., polynomials where every variable appears squared), I am only able to present a heuristic. The heuristic might often work OK in practice, but I have no proof.
(I feel like there must be a cleaner solution lurking here somewhere....)
Let $\mathbb{F}$ be a field. Suppose $p(x,y,z,\cdots) \in \mathbb{F}[x,y,z,\cdots]$ is a multivariate quadratic polynomial with coefficients in the field $\mathbb{F}$. Factor out factors of $x$, to get
$$p(x,y,z,\cdots) = c \cdot x^2 + q(y,z,\cdots) \cdot x + r(y,z,\cdots).$$
Notice that $q$ must be linear/affine, $r$ must be quadratic, and $c$ must be a constant ($c \in \mathbb{F}$). We want to find an assignment of values to $x,y,z,\cdots$ that makes
$$c \cdot x^2 + q(y,z,\cdots) \cdot x + r(y,z,\cdots) = 0.$$
Now let $\Delta$ denote the discriminant of this quadratic equation (in $x$), i.e.,
$$\Delta(y,z,\cdots) = q(y,z,\cdots)^2 - 4c \cdot r(y,z,\cdots).$$
The original polynomial is satisfiable if and only if you can find an assignment of values to $y,z,\cdots$ that makes $\Delta(y,z,\cdots)$ be a square (i.e., a quadratic residue) in $\mathbb{F}$. Notice that $\Delta(y,z,\cdots)$ is itself a multivariate quadratic polynomial (since $q$ is affine and $r$ is quadratic).
At the point, our solution approach is going to branch, based upon whether $c=0$ (an easy case) or $c\ne 0$ (the harder case).
If $c=0$, solving this equation is easy: it is basically just a linear equation (it is linear in $x$, once you fix all the other variables).
As a result, there is a special class of polynomials where this problem is especially easy to solve: any polynomial $p$ where there exists a variable -- for simplicity, we will call it $x$ -- such that $x^2$ does not appear in $p$. For that class, you will have $c=0$ above, and then the following algorithm suffices to solve the problem:
Check whether $q$ has any assignment of values to $y,z,\cdots$ that makes $q(y,z,\cdots) \ne 0$. Since $q$ is linear, this is easy to check.
If yes, then choose any such assignment of values to $y,z,\cdots$. Next, set $x = - r(y,z,\cdots) q(y,z,\cdots)^{-1}$. By definition, $q(y,z,\cdots) \ne 0$, so $q(y,z,\cdots)$ has an inverse (since we're working in a field), so such a value of $x$ exists. This immediately yields an assignment of values to $x,y,z,\cdots$ that makes $p$ zero, and we're done.
If no, then we know $p(x,y,z,\cdots) = r(y,z,\cdots)$ for all possible assignments of values to the variables of $p$. Effectively, $p$ never depended on $x$ in the first place, so it was a mistake to think of it as a function of $x$. Or, put another way, we have eliminated one variable from $p$, and we get a new instance of the original problem. Recursively apply our methods to $r(y,z,\cdots)$. Exception: if $p$ was a function of a single variable (i.e., $p(x)$), then in this case you can conclude that $p$ is not satisfiable (it is identically zero for all possible values of its variables).
This algorithm handles the easy case: namely, polynomials where there exists at least one variable $x$ that appears in $p$ but where $x^2$ does not appear in $p$. For this case, the algorithm runs in polynomial time and determines whether there exists an assignment of values to the variables that makes the polynomial equal to zero.
Now back to the hard case. If it doesn't fall into the easy case, then we must have
$$p(x,\cdots,z) = c_x \cdot x^2 + \cdots + c_z \cdot z^2 + p'(x,\cdots,z).$$
where $p'(x,\cdots,z)$ has no squared-terms (it is multilinear). Let's look at what methods we can use for this case.
Well, here's one situation where we can solve this. Suppose there exist a pair of variables $x,y$ such that $-c_x/c_y$ is a square (a quadratic residue) in $\mathbb{F}$, say, $-c_x/c_y=\alpha^2$. Then we can apply the change of variables $y' = y+\alpha x$. Conveniently, we have
$(y')^2 = (y+\alpha x)^2 = y^2 + 2\alpha \cdot xy + \alpha^2 \cdot x^2 = y^2 + 2\alpha cdot xy - c_x/c_y \cdot x^2.$$
Plugging this change of variables into $p$, we get
$$p(x,y',\cdots) = c_x \cdot x^2 + c_y \cdot y^2 + 2\alpha c_y \cdot xy - c_y \cdot (c_x/c_y) \cdot x^2 + \dots,$$
i.e.,
$$p(x,y',\cdots) = c_y \cdot y^2 + 2\alpha c_y \cdot xy + \dots$$
where the omitted part does not contain a $x^2$ term (since the $c_x \cdot x^2$ and $\alpha^2 c_y \cdot x^2$ terms cancel). Effectively, we have eliminated the $x^2$ term, so now we can apply the method above to $p(x,y',\cdots)$. Notice that $p(x,y',\cdots)=0$ is satisfiable if and only if $p(x,y,\cdots)=0$ is. When we find a solution that makes this zero, we can back-solve for $y$, and we obtain an assignment to $x,y,z,\cdots$ that makes $p(x,y,z,\cdots)=0$.
Is it possible that no pair of variables allows us to eliminate a squared-term? If $-1$ is a quadratic residue (a square) in $\mathbb{F}$, and if we have at least 3 variables, then that is not possible: $-c_x/c_z = (-1) \times (-c_x/c_y) \times (-c_y/c_z)$, and since the product of a quadratic non-residue and a quadratic non-residue is a quadratic residue, we are guaranteed that at least one of $-c_x/c_z$, $-c_y/c_z$, or $-c_x/c_z$ is a quadratic residue so at least one of the squared terms can be cancelled. However, if $-1$ is a quadratic non-residue, or we have only 1 or 2 variables, then it might not be possible to eliminate a squared-term. That's the remaining difficult case. (The situation with only 1 or 2 variables is not difficult to handle, so really the difficult case is where $-1$ is a quadratic non-residue, and where each variable appears squared in $p$.) I have no general solution for this difficult case.
However, I can suggest a heuristic that may often work, even for this difficult case.
In any field $\mathbb{F}$, approximately half of all field elements are squares (quadratic residues); if you pick a random field element, it will be a square with probability about $1/2$. Therefore, if we pick values for $y,z,\ldots$ randomly, we can predict that (heuristically) $\Delta(y,z,\cdots)$ should be a square with probability $1/2$, if $\Delta(y,z,\cdots)$ acts like a random function. This suggests the following heuristic algorithm:
Randomly pick a variable to eliminate, say $x$.
Pick values for $y,z,\ldots$ randomly.
If $\Delta(y,z,\cdots)$ is a square in $\mathbb{F}$, then the equation $p(x,y,z,\cdots)=0$ has a solution for $x$, namely $x=(-q(y,z,\cdots) \pm \sqrt{\Delta(y,z,\cdots)})/(2c)$ (assuming $\mathbb{F}$ does not have characteristic 2). This gives us an assignment to the variables $x,y,z,\ldots$ that makes $p(x,y,z,\cdots)=0$, so we're done.
If $\Delta(y,z,\cdots)$ is not a square, go back to step 1.
If after many steps, you do not find any solution, then you might guess that the equation is unsatisfiable (but this is a heuristic, so your guess could be wrong).
This algorithm might work well for many of the remaining difficult polynomials, but I have no proof that it will always work. In particular, there is a risk that all the possible discriminant polynomials $\Delta(\cdots)$ might have the property that their values are always non-squares (quadratic non-residues), in which case the algorithm above will fail.
Your question does not specify the field you are working in. Sasho Nikolov has given an excellent answer if you are working over $\mathbb{R}$, but it does not extend to finite fields -- so I think it remains interesting to look at your problem, for the case of a finite field $\mathbb{F}$. I do not have a complete solution, but here is an algorithm will work for all multilinear quadratic polynomials except for one specific class of polynomials (which I do not know how to handle yet). It feels close: If we could find some extension to solve this problem for all bivariate polynomials, the techniques below could take care of the rest.
If $p$ is univariate, i.e., $p(x)=ax^2 + bx + c$, then the problem is easy: simply test whether $b^2-4ac$ is a square in $\mathbb{F}$.
If $p$ is bivariate, i.e., $p(x,y) = a x^2 + b y^2 + c xy + dx + ey + f$, then we can use the following procedure:
We will first make some transformations to "force" $a$ to zero, since this is particularly convenient for what follows:
If $a=0$ already, then there is nothing to do.
If $b=0$, we can swap the role of $x$ and $y$.
Alternatively, suppose $a\ne 0$ and $b \ne 0$. If $-b/a$ is a non-square in $\mathbb{F}$ I am screwed, and I don't have a solution (in this case I give up and go home). Otherwise, apply a change of variables $x' = \sqrt{-b/a} \cdot x$, $y' = x + y$ and define $p'(x,y)=p(x',y')$, noticing that $p'(x,y)$ has no $x^2$ term (it cancels after the change of variables). Notice that $p$ is equisatisfiable with $p'$ (any solution to $p'(x,y)$ can be immediately translated to a solution to $p(x,y)$ by applying the change of variables, and vice versa).
After this step, we have a polynomial with no $x^2$ term. So, without loss of generality we can assume $a=0$ ($p(x,y)$ has no $x^2$ term).
Since $p$ has no $x^2$ term, we can rewrite it as
$$p(x,y) = q(y) \cdot x + r(y).$$
Conveniently, this is a linear function of $x$. So, we can test whether there is a solution to $p(x,y)=0$, as follows:
Check whether there exists $y$ such $q(y) \ne 0$. Conveniently, degree considerations ensure that $q(y)$ is linear/affine, so this is easy to check.
If yes, then let $y$ be any value such that $q(y) \ne 0$. Now we can let $x = -r(y) q(y)^{-1}$. Since $q(y) \ne 0$, the inverse is sure to exist, so this gives an assignment $x,y$ that makes $p$ zero, and we're done.
If no, then we know that $p(x,y)=r(y)$. By degree considerations, $r$ must be a quadratic polynomial. Therefore, we can recursively apply our algorithm to $r$ (which has one fewer variable). There will be a solution to $p(x,y)=0$ if and only if there is a solution to $r(y)=0$.
This algorithm handles all bivariate polynomials except for the unfortunate case where $-b/a$ is a non-square.
Next, suppose we have a trivariate quadratic polynomial. Strangely, this case is actually easier.
$$p(x,y,z) = a x^2 + b y^2 + c z^2 + \cdots.$$
If all of $a,b,c$ are non-zero, we will eliminate the $x^2$ term using a cute trick. Find $\alpha,\beta \in \mathbb{F}$ such that $b \cdot \beta^2 + c \cdot \gamma^2 = -a$ (e.g., using the algorithm here; such $\beta,\gamma$ are always guaranteed to exist). Now apply the change-of-variable transformation $x'=x$, $y'=\beta x + y$, $z'=\gamma x + z$, yielding the polynomial $p'(x,y,z)=p(x',y',z')$. Notice that, thanks to the careful choice of the transformation, $p'(x,y,z)$ has no $x^2$ term.
It's smooth sailing from here. We can write
$$p'(x,y,z) = q(y,z) \cdot x + r(y,z),$$
where $q$ is linear and $r$ is quadratic. If $q$ is identically zero, we recurse on $r$. Otherwise, we can pick any values of $y,z$ that make $q(y,z) \ne 0$ and then immediately solve for $x$, as before.
(Isn't it interesting that the case of 2 variables seems to be harder than the case of 3 variables or 1 variable?)
I'm sure you can see how to extend these ideas to a multivariate quadratic polynomial with more than three variables.
This algorithm runs in polynomial time. Unfortunately, there is a class of unfortunate polynomials where I am screwed and I am forced to give up. Can anyone see how to handle the remaining class of bivariate quadratic polynomials?
This does not give a full answer for the question, but provides an interesting way to understand the space of potential solutions to an multivariate quadratic polynomial.
Let $ x^T A x + b^T x + c = 0$ be a multivariate quadratic polynomial equation. Note that $(Dx + h)^T(Dx + h) = x^T D^T D x + x^T D^T h + h^T Dx + h^T h$. If D is a diagonal matrix, this gives:
$$ (Dx + h)^T(Dx + h) = x^T D^2 x + 2 h^T Dx + h^T h $$
Returning to our original equation, $x^T A x + b^T x + c = 0$, note that without loss of generality we can assume that $A$ is symmetric, then it has a unique positive definite square root which may be computed in polynomial time. Call this root $D$, then we have:
$$ x^T A x + b^T x + c = 0 \implies x^T D^T D x + b^T x + \frac{b^T (D^{-1})^2 b }{4} = \frac{b^T (D^-1)^2 b }{4} - c $$
Therefore,
$$ (Dx + \frac{b^T D^{-1}}{2})^T (Dx + \frac{b^T D^{-1}}{2}) = \sqrt{\frac{b^T (D^-1)^2 b }{4} - c}$$
So the solutions $x$ may be obtained by applying an affine transformation to vectors of norm $\sqrt{\frac{b^T (D^-1)^2 b }{4} - c} = \sqrt{\frac{b^T A^{-1} b }{4} - c}$.
