# Equivalence of linear systems, solving one instead of the other

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?

• You need that $A^T$ acts bijectively on the range of $L$. – shuhalo Sep 21 '12 at 11:10
• @Martin Could you make a more elaborate answer accounting for the questions I primarily posted? Thanks – usero Sep 21 '12 at 11:13

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

$$Ly = b, \tag 0$$

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

$$LAx = b, \tag 1$$

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

$$A^{T}LAx = A^{T}b. \tag 2$$

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.

• Your elaboration addressed the important assumption on $A$ that make my proposal meaningful. Thanks. – usero Sep 22 '12 at 18:06
• Nice answer, but you fail to mention that in order to guarantee that $A^TLA$ is non-singular you have to assume (a) that $L$ is positive definite or (b) that the columns of $A$ are chosen in a sensible way. If $L$ is pos. def., then $\min_y \frac12 y^TLy -y^Tb$ nicely maps to $\min_x \frac12 x^TA^TLAx - x^TA^Tb$, justifying the request of the orthogonality of the residuum to $A^T$. – Stefano M Sep 22 '12 at 19:54
• You're right that if $L$ is positive definite, the residual assumption maps nicely to a least squares problem. However, in applications (for example, in combustion, among other others), $L$ may not be positive definite (or even symmetric), and these methods are still used anyway. "Sensible" choices for $A$ are still a matter of debate (and consternation) in parts of the model reduction community. Singularity of $A^{T}LA$ depends on the application; for ODEs, the matrix $I - A^{T}LA$ is more meaningful for implicit methods, which is one reason I sidestepped the issue. One can only say so much. – Geoff Oxberry Sep 22 '12 at 23:58

You'd look at simple cases of your statements before bothering others. Take $2\times 1$ matrices:

If $A=0$ and $b\ne 0$, (1) has no solution but (2) has infinitely many.

• Let me add another trivial example: $A=\begin{bmatrix} 1 \\ 0 \end{bmatrix}$ and $L = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$. – Stefano M Sep 21 '12 at 20:38

Suppose $\mathop{\rm rank}(A)=m$.

If $L$ is positive definite then the solution of $A^TLAx = A^Tb$ (2) is the least squares solution of $LAx = b$ (1), provided that the residual $r = LAx -b$ is minimized in the norm $\| r\|^2_* = r^T L^{-1}r$: $$\min_x \| LAx - b \|^2_* = \min_x (LAx-b)^TL^{-1}(LAx-b) \quad \Leftrightarrow \quad A^TLAx = A^Tb$$

If $L$ is not definite then no simple general relation holds between the solutions of (2) and (1).