LU decomposition for cyclic tridiagonal matrices 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.
 A: 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}
$$
