Solving $(H_1H_2+H_2H_1)x_2=H_2b_2-H_1b_1$, where $H$'s are Householder's matrixes. Solving $(H_1H_2+H_2H_1)x_2=H_2b_2-H_1b_1$. $H_1$ and $H_2$ are Householder's matrices, $x_2,b_1,b_2$ are vectors.
Is there an algorithm to calculate $x_2$ out from this equation with linear cost?
It seems that we have to calculate $(H_1H_2+H_2H_1)^{-1}$ But finding inverses of matrices usually has $n^2$ cost. Can we somehow use konwledge that we have Householder's matrices here?
 A: Using your efforts in Find algorithm solving system of equations involving Householder matrixes. you should be easily be able to find 
\begin{align}
H_1H_2+H_2H_1
&=2I-4(v_1v_1^T+v_2v_2^T)+4(v_1^Tv_2)(v_1v_2^T+v_2v_1^T)
\\
&=2(I-2VCV^T)
\end{align}
where $V=(v_1,v_2)$ and $C=\pmatrix{1&-v_1^Tv_2\\-v_1^Tv_2&1}$.
This makes this sum a rank-2 modifiation of a diagonal matrix that can be efficiently inverted, as I linked in the answer to that question, using the Sherman-Morrison-Woodbury matrix identity. The inverse has the form 
$$
\frac12(I+2VDV^T)\text{ with } D=(I-2CV^TV)^{-1}C
$$
Solving linear systems with this matrix only requires the inversion of a $2×2$ matrix, in this case it reduces to a scalar division.

Another approach is to rotate the basis so that $v_1=e_1$ and $v_2=c_{2}e_1+s_{2}e_2$. Then 
$$
H_1=\pmatrix{-1&0\\0&1},\;H_2=\pmatrix{-c&-s\\-s&c}\text{ and }H_1H_2+H_2H_1=2cI
$$
where $c=c_{2}^2-s_{2}^2=2c_{2}^2-1=2(v_1^Tv_2)^2-1$. This gives the solvabilitiy condition as $c\ne 0$ or $v_1^Tv_2\ne\pm\sqrt{\frac12}$, that is, the angle between $v_1$ and $v_2$ should not be $45°$ or $135°$.
A: This problem may be solved with elementary approach covered in linear algebra course. Assume (like in your first post) that $v_i$ are unit vectors and $|v_2^Tv_1| \neq \frac{\sqrt{2}}{2}$.
Observe that $v_1$ and $v_2$ are eigenvectors of $H_1H_2 + H_2H_1$ with non-zero eigenvalues. Indeed, since: 
$$
H_ix = x  - 2(v_i^Tx)v_i
$$
by direct computation (check it yourself) we obtain:
$$
(H_1H_2 + H_2H_1)v_1 = (2 - 4(v_2^Tv_1)^2)v_1
$$
$$
(H_1H_2 + H_2H_1)v_2 = (2 - 4(v_2^Tv_1)^2)v_2
$$
Since $|v_2^Tv_1| \neq \frac{\sqrt{2}}{2}$ we have $2 - 4(v_2^Tv_1)^2 \neq 0$. Let $\lambda = 2 - 4(v_2^Tv_1)^2$.
Consider plane $P = \text{lin}(v_1,v_2)$. Since $v_1$ and $v_2$ are eigenvectors with the same eigenvalue, all vectors $v \in P$ are eigenvectors with that eigenvalue. Using Gram-Schmidt algorithm we may find orthonormal basis for $P$ (cost for that is linear since you need to compute three scalar products and subtract some numbers): $u_1$, $u_2$.
Moreover, using Gram-Schmidt algorithm one may obtain orthonormal basis for $P^\bot$ (orthogonal complement of $P$). By definition of $H_i$, if $x$ is perpendicular to $v_1$ and $v_2$ one obtains $(H_1H_2 + H_2H_1)x = 2x$ so these vectors are eigenvectors with eigenvalue 2. It follows that any vector from $P^\bot$ is eigenvector with eigenvalue 2.
Therefore, any vector $x \in \mathbb{R}^n$ may be expressed as linear combination:
$$
x = c_1u_1 + c_2u_2 + u
$$
where $u_1$, $u_2$ and $u$ are orthogonal with respect to each other ($u \in P^\bot$). Applying $H_1H_2 + H_2H_1$ to this vector we obtain:
$$
H_1H_2 + H_2H_1(x) = c_1\lambda u_1 + c_2\lambda u_2 + 2u
$$
To solve the system $(H_1H_2 + H_2H_1)x = c$ we need to compute $c_1$, $c_2$ and $u$. By orthogonality:
$$
c_1 = <c,u_1>
$$
$$
c_2 = <c,u_2>
$$
with linear cost. Moreover, we can find $u$ by simple subtraction:
$$
u = \frac{1}{2}(c - c_1\lambda u_1 + c_2\lambda u_2)
$$
So to sum up, we need to find $u_1$ and $u_2$ (linear cost, Gram-Schmidt algorithm), $c_1$ and $c_2$ (two dot products, linear) and $u$ (simple operation, linear). 
