# $LU$ Factorization of a nonsingular matrix with a particular pattern

Consider $S\in\mathbb{R}^{n\times n}$ whose nonzero elements have the following pattern for $n = 8$: $$\begin{pmatrix} 1 & 0 & 0 & 0 & \mu_1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & \mu_2 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & \mu_3 & 0 & 0 & 0\\ 0 & 0 & 0 & \alpha & \beta & 0 & 0 & 0\\ 0 & 0 & 0 & \gamma & \delta & 0 & 0 & 0\\ 0 & 0 & 0 & \delta_1 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & \delta_2 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & \delta_3 & 0 & 0 & 0 & 1\\ \end{pmatrix}$$ The pattern generalizes to any $n$ easily. Assume that for any $n$, $S$ is a nonsingular matrix.

a.) We have considered several basic transformations (Gauss transforms, Gauss Jordan transforms, elementary permutations, Household reflectors) that can be used to compute the factorizations efficiently.

Using what ever combination of these transformations, describe an algorithm to compute stably a factorization of $S$ for any $n$ that can be used to solve $Sx = b$. Your algorithm should be designed to require as few computations as possible. Your solution must include a description of how you exploit the structure of the matrix and its factors.

Attempted solution - Right off the back to me this matrix looks ugly so I am thinking of putting it into a better form so essentially applying permutations to the matrix. In part of my professors solution he says: One way to approach this problem performs an priori symmetric permutation i.e., re ordering the rows and columns with the same permutation. Consider $P^T S P$ where $$P = \begin{pmatrix} e_1 & e_2 & e_3 & e_6 & e_7 & e_8 & e_4 & e_5\end{pmatrix}$$ I am wondering about the logic of how we gets $P$. Note that $e_i$ just denotes the position of the $1$. For example $e_1$ is just an $8$ by $1$ matrix with $1$ in the first entry followed by $0$'s in the remaining entries.

So that is question $1$ for me.

Now applying this we get $$\tilde{S} = P^T S P = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & \mu_1\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & \mu_2\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & \mu_3\\ 0 & 0 & 0 & 1 & 0 & 0 & \delta_1 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & \delta_2 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & \delta_3 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \alpha & \beta\\ 0 & 0 & 0 & 0 & 0 & 0 & \gamma & \delta\\ \end{pmatrix}$$ Alright, so now we have a better form to the matrix. Now I want to apply the Gauss transform to this matrix in order to get $U$ in the $LU$ factorization. Now since all the nonzero entries that we need to stress over is in the 2 last columns I believe we just need to find $M_6$ and $M_7$ to get my $U$ i.e., I am applying $$M_{n-1}\ldots M_1 A^{(n-1)} = U$$ I hope that is clear if not please let me know. Now once I do this correctly I can simply apply a forward and backward solve to compute the system $Sx = b$.

Any suggestions or comments in regards to this question are greatly appreciated.

Somebody has been studying the Spike algorithm for solving tridiagonal linear systems using two processors :)

Regardless, while permuting the matrix is a valuable theoretical tool, the operation is unnecessary in practice and requires additional time in the form of data movement.

Instead focus on the central 2 by 2 block which engages the 2 unknowns $x_4$ and $x_5$ and nothing else. You permutation reveals why this system is non singular, right? Once this system has been solved using a simple LU factorization, it is straight forward to recover $x_1, x_2$, and $x_3$ and $x_6, x_7$, and $x_8$. If $f$ denotes the right hand side of the (Spike) system $Sx = f$, then you have

$$x_j + \mu_j x_5 = f_j, \quad j=1,2,3,$$

and

$$\delta_j x_4 + x_{5 + j} = f_{5+j}, \quad j = 1,2,3.$$

• Which 2 by 2 block are you referring to? Yes my permutation does reveal that it's nonsingular. I am little confused by this your answer. It turns out that we actually haven't been studying the spike algorithm. What I was planning to do is use the permutation then that would cut down some computations since we can just hit this permitted matrix by the gauss transform matrix $M_6$ and $M_7$ then we will have the upper trapezoidal matrix $U$. Then again your the expert so if your method is easier I will try to keep re reading it – Wolfy Jun 2 '16 at 8:34
• I am refering to the submatrix which contains the numbers $\alpha, \beta, \gamma, \delta$. Ignore the reference to the Spike algorithm. Your professor is a friend of Ahmed Sameh who introduced the Spike algorithm. Your professor simply seized upon the opportunity to make a nice problem and he did not change the notation which made me smile. There is no harm in your approach and it is mathematically equivalent to mine. I was thinking of how it would be implemented and that is irrelevant at this point. Always concentrate on getting a reliable algorithm before you worry about speed. – Carl Christian Jun 2 '16 at 10:18
• I see, by the way one think that I was confused about in the question I posed was how he got the permutation matrix in the form that he defined? – Wolfy Jun 4 '16 at 10:25
• There is no trick. Placing, say, $e_6$ in the fourth column of $P$ simply has the effect that the fourth column of AP is equal to $a_6$, the sixth column of $A$. I recommend that you experiment a bit with different permutation matrices. – Carl Christian Jun 4 '16 at 11:12