Search code examples
pythonmatlabscipylinear-algebrasparse-matrix

Solve Over-determined sparse matrix in Scipy (from Matlab to Python)

Given a large sparse matrix A which are banded or tridiagonals (however it is called) and a vector f, I would like to solve for Z, where AZ = f.

There are 6 diagonals, not clearly shown here. basically there are 6 diagonals, not clearly shown here

A has more M rows than N columns (just by 1, M ~= N), hence it is over-determined. Here is the source Matlab code, and I would like to convert it to its Scipy equivalent.

Matlab

A = A(:,2:end); #less one column
f = f(:);

Z = A\f;
Z = [0;-Z];
Z = reshape(Z,H,W);
Z = Z - min(Z(:));

My attempt on Scipy gives me this, but solving Z with scipy.sparse.linalg lsqr & lsmr is a lot slower than Matlab \ as well as not giving a good enough solution. A is created as a csr_matrix.

Python

A = A[:,1:]
f = f.flatten(1)

Z = la.lsqr(A, f, atol=1e-6, btol=1e-6)
#Z = la.lsmr(A, f)   # the other method i used
Z = Z[0]
Z = np.append([0], np.negative(Z))
Z = np.reshape(Z, (height, width), order='F').copy()
Z = Z - Z.flatten(1).min()

Could anyone recommend a better alternative to solve for Z, that is as effective and fast as Matlab \ ?

Solution

  • This looks like a good candidate for solve_banded.

    Unfortunately, the interface for providing the banded matrix is a little complex. You could start by converting your sparse matrix to DIA format, and work from there.