0
$\begingroup$

I want to quickly solve linear least squares problem for $x \in \mathbb{R}^n$

$$min_x \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 michg 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.

$\endgroup$
0
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
0
$\begingroup$

For software which does the LU factorization, then solve, look into the banded solvers in LINPACK or the newer LAPACK. They are available as free software, usually installed on Linux workstations. The LAPACK solvers are often included in ioptimi9zed libraries, like Intel's MKL or Cray's SCILIB.

Since it is positive definite, look to see if there are banded Cholesky factorizations since this should cut the operation count by a factor of 2.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.