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?


  • 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.

  • $\begingroup$ Thank you for the answer. If the least squares solution is used, then the solution will be approximate or exact? $\endgroup$ – alamaranka Jul 8 '15 at 19:41
  • $\begingroup$ 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). $\endgroup$ – Budenn Jul 8 '15 at 20:00
  • $\begingroup$ 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? $\endgroup$ – alamaranka Jul 8 '15 at 22:47
  • $\begingroup$ 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. $\endgroup$ – Budenn Jul 9 '15 at 9:57
  • $\begingroup$ 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? $\endgroup$ – alamaranka Jul 10 '15 at 16:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.