CODEWITHSUNDEEP
pythonscipylinear-algebra
Y. Ryu
Y. Ryu

Reputation: 61

How to use scipy.linalg.solve_banded to solve Ax = b when A is a striped matrix?

To solve a linear equation system

Ax=b

in which A is a striped matrix like

image

who has 5 non-zero diagonals.

But unlike a banded matrix, A has three non-zero diagonals whose offsets are 0, -1, and 1, and two non-zeros diagonals with offsets of -m and m.

I tried directly solving it with

diagonals = [Ap, As, An, Aw, Ae]
A = diags(diagonals, [0, -1, 1, -m, m]).toarray()
x = linalg.solve(A, b)

This method created the whole A. But A is spares so this method wasted too many memories to save zero elements.

So I tried using solve_banded

A = np.zeros((2 * m + 1, len(initial)))
A[0] = Ae
A[m - 1] = An
A[m] = Ap
A[m + 1] = As
A[2 * m] = Aw
x = linalg.solve_banded((m, m), A, b)

This method costs less memories than the previous ones, but it still wasted some on (2m-4) zero vectors. Is there any smarter methods which uses only the five non-zero vectors?

Upvotes: 6

Views: 3691

Answers (1)

Wunderbar
Wunderbar

Reputation: 577

Hey I can partially answer this question. To reduce the memory storage, you can maintain the sparse matrix form, without turning it to large matrix by dropping to.array() command:

A = diags(diagonals, [0, -1, 1, -m, m])

Now sparse matrix has its own spsolve method and hence scipy.sparse.linalg.spsolve(A,b) will work.

To boost performance further, you can sparse.linalg.splu(A).solve(b) which converts A by LU-decomposition and then the SuperLU object again have a solve method. (As I tested, this method is slightly faster than directly spsolve, which also uses a LU-like decomposition. See here for instance.)

The only problem here is in midst of LU-decomposition the memory usage will also be huge due to decomposition.

I was also wondering if there is any way to synthesise the solve_banded method with sparse matrix, as there is no inherent method.

Upvotes: 4

Related Questions