Crthuky

I want to solve the system $Ax = b$ where $A in mathbb R^n times n$ and $b in mathbb R^n$ with $n approx 10^6$. If $A$ would be a fully dense matrix this would be hopeless of course but luckily the matrix is sparse and highly structured.

The matrix $A$ can be written as

beginalign
A = beginpmatrix
B_0 & B_1 & & & \
& C_1 & B_2 & & \
& & C_2 & B_3 & \
& & & C_3 & B_4 \
D & & & & C_4 \
endpmatrix
endalign

where $B_i$ are upper triangular matrices, $C_i$ are diagonal with each entry being equal to $-1$, and $D$ is a square matrix. The figure below provides a schematic overview of $A$ where every possibly nonzero entry is green and all zeroes are gray.

The upper triangular matrices $B_i$ have approximately the same structure but do not necessarily have the same entry values. The submatrices $B_i$, $C_i$, and $D$ typically have dimension $mtimes m$ with $m approx 10^5$.

I don't know whether it is useful, but we also know for all non-zeroes with $i neq j$ that $a_ij in (0, 1]$ and the entries on the diagonal equal $-1$ except for the part that overlaps with $B_0$, i.e., we have $C_i = -I$.

Note that the corresponding system $(A + I)y = b$, where $I$ is an identity matrix, is relatively easy to solve as the submatrices $C_i$ become null matrices. Then we can first solve $D hat x = hat b$ where $hat x$ and $hat b$ are the last entries of $x$ and $b$. Next, we can iteratively solve the triangular matrices one by one starting with the one closest to the bottom. Maybe we can use the solution $y$ (or a transformation of $y$) as an initial solution for $x$ for an iterative approach.

Is there a better way to solve this system compared to naively feeding the system to a solver such as Matlab or LAPACK?

EDIT 1: Green entries indicate a possibly positive value. However, for the square at the bottom, around half of the entries are expected to be zero.

EDIT 2: Given the enormous size of this problem, approximate methods are also more than welcome.

The key is the fact that $C_i = - I$, i.e., the appropriately size identity matrix is trivially nonsingular and can be used to nullify the blocks above each copy. We have the following row equivalence
$$
A =
beginpmatrix
B_0 & B_1 & & & \
& -I & B_2 & & \
& & -I & B_3 & \
& & & -I & B_4 \
D & & & & -I \
endpmatrix
sim
beginpmatrix
B_0 & B_1 & & & \
& -I & B_2 & & \
& & -I & B_3 & \
B_4D & & & -I & \
D & & & & -I \
endpmatrix sim
beginpmatrix
B_0 & B_1 & & & \
& -I & B_2 & & \
B_3B_4D & & -I & & \
B_4D & & & -I & \
D & & & & -I \
endpmatrixsim
beginpmatrix
B_0 & B_1 & & & \
B_2B_3B_4D & -I & & & \
B_3B_4D & & -I & & \
B_4D & & & -I & \
D & & & & -I \
endpmatrixsim
beginpmatrix
B_0+B_1B_2B_3B_4D & & & & \
B_2B_3B_4D & -I & & & \
B_3B_4D & & -I & & \
B_4D & & & -I & \
D & & & & -I \
endpmatrix
$$
You can now factor the block $B_0 + B_1B_2B_3B_4D$ and compute $x_0$. The rest will be just be forward substitution, but compute them in descending order $x_4$, $x_3$, $x_2$, $x_1$ so that you do not have to regenerate the spike on the left.

You could try formulating it as a feasibility optimization problem such as
$$beginequation
beginarrayrrclcl
displaystyle min_x &0 \
textrms.t. & A x & = & b \
endarray
endequation$$

and feed it into a linear programming solver based on a branch-and-price algorithm such as GCG. Such algorithms are usually designed to exploit structured variables via a Dantzig-Wolfe decomposition. In fact, GCG is a branch-price-and-cut solver, meaning it also adds cutting planes in addition to the column generation approach via DW-decomposition. Very often such structured problems like yours can be solved fast via such a decomposition algorithm.

In the above problem formulation you can use any objective function as you are only looking for a feasible solution of your system $Ax=b$. Technically, you could replace the 0-objective function with any other objective function but why not just keep it simple?

The easiest quite fast way is probably to use Carl-Christians algebraic solution.

If you are interested in low-level calculational tricks for these matrices you can keep reading below line.

Given the matrix structure there are some key matrix elements at the tips of the triangles. If you let the solution diffuse from there as in forward-substitution-like scheme or implicitly with CG you will gain a lot.

Maybe you can get some inspiration from another question of mine, if I can find it. Any triangular matrix you can factor into sparse di-diagonal matrix "in-place" multiplication. This saves calculations. But I don't remember where I wrote this question right now so I will try explain.

You can factor:
$$B_i = F_1F_2cdots F_m$$

Where each $F_i$ is an identity matrix except for the entries (i,i) and (i,i+1) which need to be determined.

Whooops $m$ matrices you think, but $m$ is huge, but since we only need store 2 values for each, it will only need $2m$ in total for each $B_i$, compared to $m^2/2$ which the tridiagonal $B_i$ ones need.

In practice for $m=10^5$ you could reduce memory fingerprint from tens gigabytes to singles of megabytes.