It is known that a tridiagonal matrix $$ A = \begin{pmatrix} b_1 & c_1 & 0 & 0 & \dots & 0\\ a_2 & b_2 & c_2 & 0 & \dots & 0\\ 0 & a_3 & b_3 & c_3 & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & a_{n-1} & b_{n-1} & c_{n-1}\\ 0 & \dots & 0 & 0 & a_{n} & b_{n}\\ \end{pmatrix} $$ has a very simple LU decomposition, that is $$ A = \begin{pmatrix} 1 & 0 & 0 & 0 & \dots & 0\\ l_2 & 1 & 0 & 0 & \dots & 0\\ 0 & l_3 & 1 & 0 & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & l_{n-1} & 1 & 0\\ 0 & \dots & 0 & 0 & l_{n} & 1\\ \end{pmatrix} \begin{pmatrix} u_{1} & v_{1} & 0 & 0 & \dots & 0\\ 0 & u_{2} & v_{2} & 0 & \dots & 0\\ 0 & 0 & u_{3} & v_{3} & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & 0 & u_{n-1} & v_{n-1}\\ 0 & \dots & 0 & 0 & 0 & u_{n}\\ \end{pmatrix} $$ which can be computed easily in $O(n)$ operations via regular LU decomposition algorithm.
I'm interested in the LU decomposition of a cyclic tridiagonal matrix, that is $$ A' = \begin{pmatrix} b_1 & c_1 & 0 & 0 & \dots & a_1\\ a_2 & b_2 & c_2 & 0 & \dots & 0\\ 0 & a_3 & b_3 & c_3 & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & a_{n-1} & b_{n-1} & c_{n-1}\\ c_n & \dots & 0 & 0 & a_{n} & b_{n}\\ \end{pmatrix} $$ with $a_1$ and $c_n$ wrapped. I know that there are efficient $O(n)$ algorithms to solve a system with that kind of matrix. But I'm interested in a LU-like decomposition. I found that this matrix can be represented as $$ A' = \begin{pmatrix} 1 & 0 & 0 & 0 & \dots & l_1\\ l_2 & 1 & 0 & 0 & \dots & 0\\ 0 & l_3 & 1 & 0 & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & l_{n-1} & 1 & 0\\ 0 & \dots & 0 & 0 & l_{n} & 1\\ \end{pmatrix} \begin{pmatrix} u_{1} & v_{1} & 0 & 0 & \dots & 0\\ 0 & u_{2} & v_{2} & 0 & \dots & 0\\ 0 & 0 & u_{3} & v_{3} & \dots & 0\\ 0 & 0 & \ddots &\ddots &\ddots & 0\\ 0 & \dots & 0 & 0 & u_{n-1} & v_{n-1}\\ v_{n} & \dots & 0 & 0 & 0 & u_{n}\\ \end{pmatrix}. $$ If I denote $A' = A + a_1 Z + c_n Z^\top$ then $$ (L + l_1 Z)(U + v_{n} Z^\top) = \underbrace{LU + l_1 v_n ZZ^\top}_{\text{tridiagonal}} + l_1 u_n Z + v_n Z^\top $$ Thus $$ v_n = c_n\\ l_1 u_n = a_1\\ LU = A - l_1 v_n ZZ^\top = A - \frac{a_1 c_n}{u_n} ZZ^\top. $$ But I would like to find an effective algorithm to perform such factorization. I tried performing symbolical decomposition of $A - \alpha ZZ^\top$ with $\alpha$ being unknown, but that uses $\Omega(n^2)$ memory, thus is not an $O(n)$ algorithm.
Actually, I was mistaken. There is a way to perform a symbolical decomposition of $B = A - \frac{a_1 c_n}{u_n} ZZ^\top$ with $u_n$ being a variable effectively in $O(n)$ operations.
One can note, that the first diagonal element of the matrix $B$ that is $b_1 - \frac{a_1 c_n}{u_n}$ has the following form $$\frac{b_1 u_n - a_1 c_n}{u_n} = p_1\frac{u_n - q_2}{u_n - q_1}$$ with $p_1 = b_1, q_1 = 0, q_2 = \frac{a_1 c_n}{b_1}$. The elements of $L,U$ are rational functions of $u_n$.
Considering rank-1 update for the LU decomposition procedure $$ \begin{pmatrix} p_i\frac{u_n - q_{i+1}}{u_n - q_i} & c_i & \\ a_{i+1} & b_{i+1} & \ddots \\ & \ddots & \ddots \end{pmatrix} = \begin{pmatrix} 1 & 0 & \\ l_{i+1} & L_{i+1} & \ddots \\ & \ddots & \ddots \end{pmatrix} \begin{pmatrix} p_i \frac{u_n - q_{i+1}}{u_n - q_i} & v_i & \\ 0 & U_{i+1} & \ddots \\ & \ddots & \ddots \end{pmatrix} $$ one can obtain a recurrent formula to compute $p_i, q_i$: $$ v_i = c_i\\ l_{i+1} = \frac{a_{i+1}}{p_i}\frac{u_n - q_i}{u_n - q_{i+1}}\\ u_{i+1} = p_{i+1} \frac{u_n - q_{i+2}}{u_n - q_{i+1}} = b_{i+1} - l_{i+1} v_i = b_{i+1} - \frac{u_n - q_i}{u_n - q_{i+1}} \frac{a_{i+1} c_i}{p_i} =\\=\frac{u_n - q_{i+1}}{u_n - q_{i+1}}b_{i+1} - \frac{u_n - q_i}{u_n -q_{i+1}} \frac{a_{i+1} c_i}{p_i} = \frac{\left(b_{i+1} - \frac{a_{i+1} c_i}{p_i}\right)u_n - \left( q_{i+1}b_{i+1} - q_i\frac{a_{i+1} c_i}{p_i} \right)}{u_n - q_{i+1}}\\ p_{i+1} = b_{i+1} - \frac{a_{i+1} c_i}{p_i}\\ q_{i+2} = b_{i+1}\frac{q_{i+1}}{p_{i+1}} - \frac{a_{i+1} c_i}{p_i} \frac{q_i}{p_{i+1}} $$ Finally $u_n$ can be determined from $$ u_n = p_n \frac{u_n - q_{n+1}}{u_n - q_{n}}\\ u_n^2 - \left( q_n + p_n \right) u_n + p_n q_{n+1} = 0\\ u_n = \frac{(q_n + p_n) \pm \sqrt{(q_n + p_n)^2 - 4 p_n q_{n+1}}}{2}. $$ It seems that this type of decomposition is not unique, since the last equation has two roots. I'm curious which one of them gives a numerically stable decomposition and whether this algorithm is stable at all.
Applied to sample matrix $$ A = \begin{pmatrix} 3 & 1 & & & & 1 \\ 1 & 1 & 1 & & & \\ & 1 & 0 & 1 & & \\ & & 1 & 1 & 1 & \\ & & & 1 & 3 & 1 \\ 1 & & & & 1 & 4 \end{pmatrix} $$ this method yields two decompositions $$ A = \begin{pmatrix} 1 & & & & & \frac{1}{3} \left(4+\sqrt{10}\right) \\ 1+\sqrt{\frac{2}{5}} & 1 &&&&\\ & -\sqrt{\frac{5}{2}} & 1 &&&\\ && \sqrt{\frac{2}{5}} & 1 &&\\ &&& \frac{1}{3}\left(5+\sqrt{10}\right) & 1 &\\ &&&&2+\sqrt{\frac{5}{2}} & 1 \end{pmatrix} \cdot \qquad \\ \qquad \cdot \begin{pmatrix} \frac{1}{3} \left(5-\sqrt{10}\right) & 1 & &&& \\ & -\sqrt{\frac{2}{5}} & 1 && &\\ & & \sqrt{\frac{5}{2}} & 1 & & \\ & & & 1-\sqrt{\frac{2}{5}} & 1 & \\ & & & & \frac{1}{3} \left(4-\sqrt{10}\right) & 1 \\ 1 & & & & & 2-\sqrt{\frac{5}{2}} \end{pmatrix} = \\ = \begin{pmatrix} 1 & & & & & \frac{1}{2+\sqrt{\frac{5}{2}}} \\ 1-\sqrt{\frac{2}{5}} & 1 & & & & \\ & \sqrt{\frac{5}{2}} & 1 & & & \\ & & -\sqrt{\frac{2}{5}} & 1 & & \\ & & & \frac{1}{3} \left(5-\sqrt{10}\right) & 1 & \\ & & & & 2-\sqrt{\frac{5}{2}} & 1 \end{pmatrix} \cdot \qquad \\ \qquad \cdot \begin{pmatrix} \frac{1}{3} \left(5+\sqrt{10}\right) & 1 & & & & \\ & \sqrt{\frac{2}{5}} & 1 & & & \\ & & -\sqrt{\frac{5}{2}} & 1 & & \\ & & & 1+\sqrt{\frac{2}{5}} & 1 & \\ & & & & \frac{1}{3} \left(4+\sqrt{10}\right) & 1 \\ 1 & & & & & 2+\sqrt{\frac{5}{2}} \end{pmatrix} $$