# Solving sparse linear least squares or a positive definite 5-band matrix system fast

I want to quickly solve the following linear least-squares problem

$$\min_{x \in \mathbb{R}^n} \left\| A x - b \right\|_2^2$$

with a special sparse structure where each row in $$A$$ has only up to $$4$$ consecutive non-zero entries. This makes its normal matrix

$$C = A^T A$$

a positive definite $$7$$-band matrix with a condition number between $$8^2$$ and $$400^2$$. So, these condition numbers aren't too bad that solving the Gauß normal equation system

$$C x = A^T b$$

instead wouldn't get me into trouble numerically, I think. What are my options? I could try conjugate gradient methods but I would prefer direct solvers that can deal with these kinds of special cases in $$O(n)$$ time independent of the condition. I'm aware of algorithms for the tridiagonal case and I guess I could try to adapt them for $$5$$ bands (?). But before reinventing the wheel and/or testing many different algorithms I wanted to ask you about what approach might be the most efficient in terms of time because I have lots of these problems (millions) with values for $$n$$ of around $$5000$$ where $$A$$ has about $$4n$$ rows.

The LU factorization of $C$ along with forward and backward substitution works well in this case. The factorization can still be done completely in-place. So, there is no need to touch or create other off-band elements.