The scipy.linalg.solve_banded() function in Python is used to solve banded linear systems. Banded linear systems are characterized by a matrix in which the non-zero elements are confined to a diagonal band. This can greatly reduce the computational cost and memory requirements compared to standard linear solvers, which do not exploit this structure.

The algorithm that is typically used for solving banded linear systems is a variation of the LU decomposition, which is adapted for banded matrices. Here is a general outline of the algorithm:

1. Factorization: The banded matrix `A` is factored into the product of a lower banded matrix `L` and an upper banded matrix `U`. In contrast to the standard LU decomposition, in the banded LU decomposition, the factors L and U also have a banded structure, similar to the original matrix.

2. Forward substitution: The lower banded matrix `L` is used to solve for an intermediate vector `y` in the equation `L * y = b`. This is done through forward substitution.

3. Backward substitution: The upper banded matrix `U` is used to solve for the desired solution vector `x` in the equation `U * x = y`. This is done through backward substitution.

This algorithm is particularly efficient for banded matrices since each step takes advantage of the limited number of non-zero elements. Specifically, the factorization process only computes the non-zero bands in `L` and `U`, and the forward/backward substitutions only involve the non-zero elements of `L` and `U`, respectively.

The scipy.linalg.solve_banded() function takes as an input the banded matrix in a condensed form, similar to what you have in `Ac`, and the right-hand-side vector `b`, and performs the steps outlined above.

Behind the scenes, `scipy.linalg.solve_banded()` is likely making use of LAPACK routines (Linear Algebra PACKage), which are highly optimized for performance. The specific LAPACK routine for solving banded linear systems is called `dgbsv` for double-precision real matrices.

Note that `scipy` is an open-source library, so you can also look directly at its source code on GitHub if you're interested in the specific implementation details.

Dear Mirtunjay Kumar,

Thank you for detailed answer with steps. I was able to factorize a square spare matrix like A by LU decomposition and then solve the Ax= b using the factors of A.

However, when A is stored into a condensed form like Ac, I could not find an algorithm to decompose Ac directly into L and U. I know one thing that L and U of Ac will also be condensed and banded like Ac .

Assuming Ac is banded matrix structure, a possible algorithm to perform LU Decomposition on Ac could look like this:

*** Assuming that no pivoting is required. If pivoting is required, the situation gets significantly more complicated for a banded matrix.

Dear Rana Hamza Shakil

My question is not about solving a system of equations without inverting a matrix in general. It is more specific towards solving a system of linear equations when the coefficient matrix is stored into a particular way i.e. into a banded_storage form.

Please go through the description of my question. Thank you

Regards,

Rana Hamza Shakil

Thomas algorithm is a simpler version of Gauss elimination, which is used to solve a system of equations which has a tridiagonal coefficient matrix. It is not used in scipy.solve_banded().

Also the example pseudo code you provided in your answer has a mistake. solve_banded(), also requires the number of upper and lower non-zero diagonals.

It is not mandatory to answer.