Non-standard quadratic matrix equation I have an equation that looks like the following:
$$
A\cdot\mathrm{diag}(x)\cdot x + B\cdot x + c = 0
$$
where $A, B, C \in \mathbb{R}^{n \times n}$ and $x, c \in \mathbb{R}^n$. $ x $ is unknown. The $\mathrm{diag}(x)\cdot x$ is a trick to get element-wise multiplication.
It is different from the Quadratic Matrix Equation (well-studied and solvable)
$$
Q(X) = AXX + BX + C = 0
$$
where $ A, B, C, X \in \mathbb{R}^{n\times n}$, $ X $ is unknown.
It's also different from the Quadratic Eigenvalue Problem (QEP) (also well-studied, solvable, and related to the Quadratic Matrix Equation).
$$
Q( \lambda ) x = (\lambda^2 A + \lambda B + C)\cdot x = 0
$$
where $ A, B, C \in \mathbb{R}^{n\times n}$, $x \in \mathbb{R}^n$ and $ \lambda \in \mathbb{R} $. $ \lambda$ and $ x $ are unknown.
I feel like there should be a reduction from my problem to either of the problems above, but I'm not well-versed enough in linear algebra to figure it out. Numerical approaches are fine too.
An alternative formulation is
$$
A\cdot\mathrm{diag}(x)\cdot \mathrm{diag}(x) + B\cdot \mathrm{diag}(x) + \mathrm{diag}(c) = 0
$$
This would be in the QME form, but I would have to choose a subset of solutions that fit $X = \mathrm{diag}(x)$.
From research, I think I could also use a Gröbner basis to solve this as a system of multivariate polynomial equations, but it seems non-trivial from first glance.
 A: Treat this as a nonlinear equation, i.e. $\,f(x)=0,\,$ and apply Newton's method.
$$\eqalign{
 f &= A(x\odot x) + Bx + c \\
df &= A(2x\odot dx) + B\;dx \\
   &= \Big(2A\cdot{\rm Diag}(x)+B\Big)\,dx \\
   &= \big(2AX+B\big)\,dx \\
\frac{\partial f}{\partial x} &= (2AX+B) \;\doteq\; J \qquad&\big({\rm Jacobian\,of\,}f\big) \\
x_+ &= x - J^{-1}f&\big({\rm Newton\;iteration}\big) \\
}$$
This iteration will converge to a solution near the initial guess $x=x_0$
Since this is a simple quadratic function, the Newton method should converge quickly.
If the rate of convergence isn't satisfactory, perform a line-search at each iteration.
$$\eqalign{
\lambda &= \arg\min_\lambda\;
  \Big\|\,f\Big(x-\lambda\,J^{-1}f(x)\Big)\,\Big\|_2^2 \\
x_+ &= x - \lambda\,J^{-1}f \\
}$$
A: If I understand correctly the question, then the equation is $A[x_1^2,\cdots,x_n^2]^T+B[x_1,\cdots,x_n]^T+c=0$. We assume that $A=[a_{i,j}],B=[b_{i,j}],c=[c_i]^T$ are generic matrices, that is, 
(*): $(a_{i,j}),(b_{i,j}),(c_{i})$ are assumed to be mutually transcendental over $\mathbb{Q}$, the rational numbers field. 
According to the Bézout's Theorem, this generic system ($n$ equations of degree $2$ in $n$ unknowns) has at most $2^n$ solutions. This equation seems to be simpler than the generic Riccati matrix equation $AX^2+BX+C=0$ where $X$ is a square matrix ; note that the generic Riccati equation is well studied but is not solvable by radicals except if $n=2$ and that it has ${2n}\choose{n}$ solutions that is much smaller than the bound $2^{n^2}$ given by Bezout's Theorem (but is greater than $2^n$). I think that, for our problem, the Bezout bound is reached.
Indeed, we can simulate the genericity by randomly chosing the matrices $A,B,c$. Then, using the Maple command Groebner, we obtain, when we solve our equation (with $n\leq 20$), exactly $2^n$ solutions.
