# Solving large, sparse system of linear equations

I have a system of linear equations as follows:

$$(A+I)x=B$$

where $I$ is the $n\times n$ identity matrix,

$A$ is a $n\times n$ matrix such that the first and last rows are blank, and, for every other row $i$, with $2\le i\le(n-1)$, we choose one $k_i, 1\le k_i \le \min(i-1,n+1)$ such that, for a known value $q$,

$$a_{i,i-k_i}=-q, a_{i,i+k_i}=q-1, 0<q<1$$

and $B$ is an $n\times 1$ matrix:

$$B = \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}$$

What would be the most efficient way to calculate the solution $x$ if

1. I required the full solution vector $x$?
2. I only require the first few elements of the solution vector $x$ (in the case of large $n$)?

I realize that my question is a bit specific, so I would appreciate if the answers could point me towards computational techniques, rather than giving me an exact answer. What I have in mind would be some decomposition technique that makes the calculation less computationally intensive. I'd also appreciate if someone could provide a better representation of $A$ that more clearly captures the symmetry.

My current approach is to directly calculating the solution. While this is faster than calculating the full inverse of the matrix, it does not exploiting the inherent symmetry and sparseness of the matrix.

-
First I would suggest to use the sparse LU factorization (possibly with some reasonable ordering to preserve sparsity of the triangular factors. The matrix itself is an M-matrix (en.wikipedia.org/wiki/M-matrix), which are generally considered "nice" to be solved. By the way, it is always faster to solve such system by factorization methods than trying to form the inverse. In fact, the inverse of a sparse matrix is usually dense (e.g., the inverse of a generic tridiagonal matrix is completely full). Are you looking for some hints for Matlab or Octave or rather some software libraries? –  Algebraic Pavel Sep 9 '13 at 9:22
@PavelJiranek That'll be great. I'm currently using Matlab, and so far a naive approach of calculating the inverse has worked out fine since I'm only handling small matrices up to n=1000. Specifically, I now only require to calculate the average value of the solution vector $x$ - any idea what possible speed-ups I could get with Matlab? –  Vincent Tjeng Sep 11 '13 at 23:45
In MATLAB, you can store your matrix in the sparse storage format (help sparfun) which internally keeps only the nonzeros of the matrix column by column. Some linear algebra routines in MATLAB are overloaded for sparse matrices, so for instance x=A\b uses sparse LU factorization (MATLAB command lu) in the case if A is a sparse matrix. Regarding the problem with special (also sparse in fact) right-hand side, there could be some savings possible in using LU if you had some leading zeros. I'm not sure now how to exploit the trailing zeros. –  Algebraic Pavel Sep 12 '13 at 10:15
Thanks! I'm using x=A\b now and it's providing some speed-up ... guess I will have to investigate the different possibilities and time them properly. Will update you if I come up with anything concrete. –  Vincent Tjeng Sep 12 '13 at 13:16