Solving special (simple?) system of polynomial equations (only up to second degree)

I have the following system of equations, \begin{align} \mathrm{(I)}&&\quad \mathbf{A} f_1(\mathbf{x})&=\mathbf{b}_1,\\ \mathrm{(II)}&&\quad\mathbf{A} f_2(\mathbf{x})&=\mathbf{b}_2,\\ \mathrm{(III)}&&\quad\mathbf{A} f_3(\mathbf{x})&=\mathbf{b}_3, \end{align} where I only need to solve any two equations (for the same $\mathbf{x}$) and $\mathbf{A} \in \mathbb{R}^{9\times N}$, $\mathbf{x} \in \mathbb{R}^{N}$ , $\mathbf{b}_i \in \mathbb{R}^{9}$ and $N \in \mathbb{N}$ (probably necessary to be at least bigger or equal to $18$). $\mathbf{A} f_i(\mathbf{x})$ is a matrix-vector product in the following sense: The functions $f_i(x)$ are given by \begin{align} f_1(x) &= \frac{1}{2}\left(x^2-x+\frac{1}{4}\right),\\ f_2(x) &= -x^2+\frac{3}{4},\\ f_3(x) &= \frac{1}{2}\left(x^2+x+\frac{1}{4}\right) \end{align} and act component-wise on vectors, e.g.: $$f_i(\mathbf{x})=\begin{pmatrix} f_i(x) \\ f_i(y) \\ f_i(z) \end{pmatrix}.$$ The system can therefore also be written as (just to provide some options) \begin{align} \mathrm{(I')}&&\quad \mathbf{A} (\mathbf{x}^2-\mathbf{x}) = \mathbf{A} \mathbf{x}^2-\mathbf{A}\mathbf{x}&=\mathbf{b}_1',\\ \mathrm{(II')}&&\mathbf{A} \mathbf{x}^2&=\mathbf{b}_2',\\ \mathrm{(III')}&&\quad\mathbf{A} (\mathbf{x}^2+\mathbf{x}) = \mathbf{A} \mathbf{x}^2+\mathbf{A}\mathbf{x}&=\mathbf{b}_3', \end{align} where then, e.g., $$\mathbf{x}^2 = \begin{pmatrix} x^2 \\ y^2 \\ z^2 \end{pmatrix}.$$ I need to solve any two of the equations (I) - (III). I do not care which two exactly. My problem is that the functions $f_i$ are linearly independent, which makes any system of two equations non-linear. I tried the whole day and lots of things. For example \begin{align} f_1(x) = f_3(x) -x \Rightarrow \end{align} \begin{align} \mathrm{(I)} &&\quad\mathbf{A} f_3(\mathbf{x}) - \mathbf{A}\mathbf{x} &= \mathbf{b}_1\\ \mathrm{(III)} &&\quad\mathbf{A} f_3(\mathbf{x}) &= \mathbf{b}_3 \end{align} which nearly gives a system of linear equations (in $f_3(\mathbf{x})$) if it wouldn't be for that pesky $-\mathbf{A}\mathbf{x}$.

I need to find an efficient solution which I can do numerically in a fast way (that's why I tried to get a system of linear equations). Is there any nice mathematical relation I am missing? I looked into regular chains but got scared. My system is not that complicated, or is it?

I'd also be grateful if anybody knows a way to determine if there is a solution at all. This is not a priority, as, if there is an efficient algorithm, I'd just let the computer try and possibly fail while moving on to a different set of matrices.

EDIT: Thanks to the comment by @AlexRavsky I looked into the Moore-Penrose inverse $\mathbf{A}^+ \in \mathbb{R}^{N\times 9}$ and can now rewrite the problem. The solutions to \begin{align} \mathrm{(I')}&&\mathbf{A} \mathbf{x}^2-\mathbf{A}\mathbf{x}&=\mathbf{b}_1'\quad\mathrm{and}\\ \mathrm{(II')}&&\mathbf{A} \mathbf{x}^2&=\mathbf{b}_2' \end{align} are therefore given by \begin{align} \mathrm{(I')}&&\mathbf{x}&=\mathbf{A}^+(\mathbf{b}_2'-\mathbf{b}_1')+\left[\mathbf{I}-\mathbf{A}^+\mathbf{A}\right]\mathbf{k}_1 \quad \mathrm{and}\\ \mathrm{(II')}&&\mathbf{x}^2&=\mathbf{A}^+\mathbf{b}_2'+\left[\mathbf{I}-\mathbf{A}^+\mathbf{A}\right]\mathbf{k}_2, \end{align} where $\mathbf{I}\in \mathbb{R}^{N\times N}$ is the identity matrix and $\mathbf{k}_1, \mathbf{k}_2 \in \mathbb{R}^N$ can be chosen arbitrarily. Keep in mind that in $\mathbf{x}^2$ the square function acts element-wise. These are $2 N$ equations for my $N$ unknowns. Now I need help solving these equations for $\mathbf{x}$.

• This description of all solutions of a linear system in terms of Moore-Penrose inverse $A^+$ of the matrix $A$ may be helpful. – Alex Ravsky Jan 22 '18 at 14:41
• @AlexRavsky Thank you very much for your input. I looked into the details and applied them to my problem. It changes things, but I don't really know if for the better ;) – Nils_M Jan 22 '18 at 18:56
• Wouldn't it be cleaner to consider solving $\mathbf{A}\mathbf{x} = \mathbf{a}$ and $\mathbf{A}\mathbf{x^2} = \mathbf{b}$? You could consider the more general setting of having $\mathbf{A}\mathbf{y} = \mathbf{b}$ and solving a larger system of $\mathbf{B}\mathbf{z} = \mathbf{c}$ where $\mathbf{B}$ is block-diagonal with two blocks of $\mathbf{A}$, $\mathbf{z}$ consists of $\mathbf{x}$ and $\mathbf{y}$ concatenated and so on. Then the condition $\mathbf{y} = \mathbf{x^2}$ can be imposed later to restrict the solution set further. – Alexander Vlasev Jan 28 '18 at 9:44
• @AlexanderVlasev Thank you for your answer! What is the difference between your approach and using Eqs. (I') and (II')? Aren't these exactly the solutions you're refering to? The larger system consisting of the block diagonal matrix and the concatenated $\mathbf{x}$ and $\mathbf{b}'_i$ has the same solutions as just solving the two equations on their own, or am I mistaken? – Nils_M Feb 2 '18 at 18:29