This question is related to recently posted one, but I guess it deserves a separate attention.Suppose a symmetric matrix $L\in\mathbb{R}^{n\times n}$ is given, and a rectangular matrix $A\in\mathbb{R}^{n\times m}$, $m<n$. A solution to the system $$LAx=b, \tag 1$$ is sought, for known $b\in\mathbb{R}^n$ and unknown $x\in\mathbb{R}^m$. Since $LA\in\mathbb{R}^{n\times m}$, a solution $x$ should be calculated in a least squares sense. Now, suppose that one premultiplies the above system by $A^T$, ie, $$A^TLAx=A^Tb \tag 2$$Could one replace solving (1) by solving (2), or vice versa? Why couldn't one replace solving (2) by solving (1)? If a solution to (1) exists, then it is also a solution to (2). Could it be that the solution to (1) does not exits (ie, exists only in LS sense), but the exact solution to (2) exists?
Equivalence of linear systems, solving one instead of the other
linear algebra;matrices;linear solver
Weirdly enough, this sort of scenario can arise in orthogonal model reduction methods. In that situation, $A$ is unitary, and therefore, $A$ can take the interpretation of a basis set used for model reduction, and $A^{T}$ is a pseudoinverse (to be precise, a $\{1,2\}$-inverse, in the language of generalized inverses) of $A$.As Arnold Neumaier points out, the two sets of solutions are, in general, not the same. If $x \in \mathbb{R}^{m}$ is a solution to (1), then it will also be a solution to (2). However, the  converse of that statement is not true in the general case.Nevertheless, for the case where $A$ is unitary, such a transformation can be quite useful. Suppose you want to solve\begin{equation}Ly = b, \tag 0\end{equation}but for whatever reason (for instance, memory, CPU time), it is impractical to solve this equation using standard dense or sparse (direct or iterative) linear algebra methods. One approach is to assume that $y \in \mathbb{R}^{n}$ is a linear combination of some set of orthonormal basis vectors $A \in \mathbb{R}^{n}$. In other words, assuming that $y = Ax$, replace $y$ with $Ax$ so that you can solve for $x \in \mathbb{R}^{m}$, which has the advantage of solving for fewer unknowns. So now you have\begin{equation}LAx = b, \tag 1\end{equation}as in equation (1), but this system is nonsquare, and underdetermined. It may not have a solution, because $b$ may not be in the range of $A$. A classic assumption in the model reduction literature made to obtain a square system is to assume that the quantity $LAx - b$ (you can think of it as the residual of a value of $x$ substituted into equation (1)) is orthogonal to each column of $A$, in which case, $A^{T}(LAx - b) = 0$ or\begin{equation}A^{T}LAx = A^{T}b. \tag 2\end{equation}This assumption may or may not be a good one. If the solution $y$ to (0) is in the range of $A$, then solving (2) does not induce any approximation error in the computed solution to (0). Solving (0) can yield substantial computational savings compared to solving (2) directly, provided $m$ is small enough (there is an overhead of two matrix multiplies to reap savings in the linear solve).However, if the solution to (0) is not in the range of $A$, then solving (2) instead of (0) will only yield an approximation of the solution to (0), and the error associated with such an approximation may be substantial. Approximate solutions can be acceptable for many real-world applications (for example, problem parameters may only be known to limited precision).If $A$ isn't unitary, all bets are off. The model reduction interpretation no longer holds, although it is still true that if $x$ is a solution to (1), it is also a solution to (2). It's possible to generalize the model reduction argument for general $A$, but then $A^{T}$ is typically replaced in that argument by another matrix, so the analogy to the question above no longer holds.