Block Gauss -Seidel Iterative Method for Overdetermined Linear Systems

I am interested in solving a large linear system with block Gauss-Seidel Iterative Method. Suppose I have the following block matrix $A$ for $Ax=b$ linear system:

\begin{bmatrix} B & 0 & C & 0 \\ 0 & C & D & 0 \\ D & 0 & D & E \\ 0 & 0 & 0 & F \end{bmatrix}

where each block matrix is $n \times n$. Right hand side vector ,$b$ , is a zero vector. So far everything is fine because $A$ matrix is a square. However I have one more equation:

$\sum{x e} = 1$

where $e$ is a vector of ones. According to my best knowledge I can use Gauss-Seidel method if $A$ is a square matrix. In my case, the linear system is overdetermined, thus not square. My question is:

Can I use Gauss-Seidel method for this problem with some pre-work? If not, what do you recommend me to use?

Thanks.

• You can calculate the least squares solution, that is, find $\underset{x}{\arg \min} ||A x - b||$, where $A$ is the non-square matrix $$\begin{bmatrix} B & 0 & C & 0 \\ 0 & C & D & 0 \\ D & 0 & D & E \\ 0 & 0 & 0 & F\\ e^T & e^T & e^T & e^T \end{bmatrix} .$$ Its solution is (in overdetermined case) $x = (A^T A)^{-1} A^T b$ (where $b = (\underbrace{0, \ldots, 0}_{4 n}, 1)^T$), which can be calculated with Gauss-Seidel method. However, $A^T A$ is, in this case, less sparse than $A$.

• Another method (that would obtain the same result as above) is to calculate the singular value decomposition of $A$ and then calculate.

• If the original $A$ (without the bottom row) is not full-rank, there is a simpler way to find the solution.

1. Calculate the basis of null space of original $A$ (for example, by inverse iteration method).
2. For every vector $v_k$ of that basis, divide it by the sum of its elements. If there is only one basis vector (e.g. if the entire problem is known to have a unique solution), after the division it is the solution. Otherwise:
3. The convex hull (the set of $\sum_k \alpha_k w_k$ for non-negative $\alpha_k$ such that $\sum_k \alpha_k = 1$) of the divided vectors $w_i$ is the set of solutions.

Implementation note for the last method: In Matlab, the first step of the last method can be done by using null function. However, it involves calculating the SVD of original $A$ and doesn't work for sparse matrices, so using eigs to find eigenvectors corresponding to the eigenvalue $0$ may be faster, especially since $A$ is sparse. Instead of providing $A$ for eigs, one can use a handle to a function that calculates $A^{-1} y$ for given $y$, for example, a Gauss-Seidel solver.

• Thank you for the answer. If the least squares solution is used, then the solution will be approximate or exact? – alamaranka Jul 8 '15 at 19:41
• That depends on the method used in calculation of $(A^T A)^{-1} A^T b$. E.g. Gauss-Seidel method and other iterative methods inherently gives an approximate solution. If an exact method is used, then the solution is exact (up to round-off errors, of course). – Budenn Jul 8 '15 at 20:00
• OK, got it. One more thing: Right hand side of the equation you've provided is independent from $x$ s. So how to proceed the iteration? – alamaranka Jul 8 '15 at 22:47
• In the equation I wrote, $b = (0, \ldots, 0, 1)^T$ (with $4 n$ zeros), according to the original linear system and the additional equation. – Budenn Jul 9 '15 at 9:57
• Thanks again. I have a question about Matlab. I coded the algorithm as you describe above. In this problem vector $x$ represents probabilities. Thus all entities of this vector are supposed to be between 0 and 1. However, when I run the code, for large $A$ matrix (10202 $\times$ 10201), solution contains some negative values (all values are less than 1). How can I prevent this? Is there a way to force all elements to be non-negative? – alamaranka Jul 10 '15 at 16:07