# Solving system of linear equations that consists of a few matrices

I have a task to write a Matlab program that get LU-decomposition of a given matrix
$$A$$ and then solves a system of linear equations using the obtained decomposition.

$$Mz = f$$, where $$M = \left(\begin{array}[c c] - I & A\\ A^T & 0 \end{array}\right)$$

Before I start with coding, I have a problem with the theoretical part of this task, mainly how to approach such an equation that is made of a few matrices.

• I had made a number typographical errors in my answer. They have now been fixed. Upon futher reflection, I added some additional material. – Carl Christian Apr 16 at 21:34

Let $$2m$$ denote the dimension of $$M$$. You are familiar with the case of $$m=1$$. Here $$Mz: = \begin{bmatrix} 1 & a \\ a & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} =:f$$ The symbols ":=" and "=:" indicate that we are defining objects, i.e., $$M$$ and $$z$$ as well as $$f$$. This linear system is equivalent to \begin{align} z_1 + az_2 &= f_1, \\ a z_1 + 0 \cdot z_2 &=f_2. \end{align} In the case of $$m>1$$, the vectors $$f$$ and $$z$$ are still split in two, but the "components" are now vectors in $$\mathbb{R}^m$$. The linear system reads $$Mz: = \begin{bmatrix} I & A \\ A^T & 0 \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} =:f$$ where $$z_i \in \mathbb{R}^{m}$$ and $$f_i \in \mathbb{R}^{m}$$. This linear system is equivalent to \begin{align} z_1 + Az_2 &= f_1, \\ A^T z_1 + 0 \cdot z_2 &=f_2. \end{align} Your matrix $$M$$ is special case of a block matrix/partitioned matrix. A more general discussion of how to multiply block/partitioned matrices can be found in this answer to this question.
Standard Gaussian elimination is an example of a scalar algorithm, because it operates on individual real numbers, i.e., scalars. There is a block variant of Gaussian elimination which is structurally identical to the scalar algorithm but operates on blocks. Let us consider your example. It is clear that \begin{align} z_1 + Az_2 &= f_1, \\ A^T z_1 + 0 \cdot z_2 &=f_2. \end{align} implies \begin{align} A^T z_1 + A^TA z_2 &= A^T f_1, \\ A^T z_1 + 0 \cdot z_2 &=f_2. \end{align} We can now eliminate $$z_1$$ and find that $$- A^T A z_2 = f_2 - A^T f_1.$$ Once this system has been solved with respect to $$z_2$$ we can compute $$z_1$$ using $$z_1 = f_1 - Az_2.$$ There are a couple of points to mention. Above the horizontal line I treated $$A$$ as a square matrix. I deliberately sacrificed generality for the sake of clarity. Below the horizontal line, I am allowing $$A$$ to be an $$m$$ by $$n$$ matrix. This implies that $$z_1 \in \mathbb{R}^m$$ and $$z_2 \in \mathbb{R}^n.$$ I need $$A^TA$$ to be nonsingular, so $$A$$ must have full rank, i.e., I must have $$m \ge n$$. In other words, $$A$$ must be a tall matrix. Moreover, we should not form the matrix $$A^TA$$ because this matrix can be very ill-conditioned. Instead we can compute a QR factorization of $$A$$, i.e., $$A = QR$$, where $$Q$$ is orthgonal and $$R$$ is upper triangular. In reality, this is just the Gram-Schmidt algorithm applied to $$A$$. Then $$A^T A = R^T Q^T Q R = R^T R$$. We see that the Cholesky factorization of $$A^TA = LL^T$$ is immediately available to us with $$L = R^T$$.