spsolve vs lsmr

About

When solving large-scale linear systems, computational efficiency and memory usage are critical considerations. Sparse matrices, which contain a significant number of zero elements, are commonly used in scientific computing, engineering simulations, and optimization problems to reduce memory footprint and improve performance.

Two widely used solvers for sparse linear systems in SciPy are spsolve and lsmr. While both are designed to handle sparse matrices, they operate in fundamentally different ways. spsolve is a direct solver that leverages LU decomposition to obtain an exact solution, making it highly accurate but potentially memory-intensive. On the other hand, lsmr is an iterative solver based on the least-squares method, which is often more efficient for extremely large systems but may require tuning to achieve convergence.

Experiment spsolve vs lsmr (N=50000)

Code

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import time
import matplotlib.pyplot as plt

N = 50000

np.random.seed(42)
B_sparse = np.random.rand(N)

# Sparse Matrix Construction
diagonals = [
    np.random.rand(N),
    np.random.rand(N-1),
    np.random.rand(N-1)
]
A_sparse = sp.diags(diagonals, offsets=[0, 1, -1], format='csr')

# Solve using spsolve
start_time = time.time()
X_sparse_spsolve = spla.spsolve(A_sparse, B_sparse)
sparse_spsolve_time = time.time() - start_time

# Solve using lsmr
start_time = time.time()
X_sparse_lsmr = spla.lsmr(A_sparse, B_sparse)[0]
sparse_lsmr_time = time.time() - start_time

# Plot results
labels = ["Sparse (spsolve)", "Sparse (lsmr)"]
times = [sparse_spsolve_time, sparse_lsmr_time]

plt.figure(figsize=(8, 5))
plt.bar(labels, times)
plt.ylabel("Time (seconds)")
plt.title("Comparison of Computational Time: spsolve vs lsmr")
plt.savefig("test.jpg")

Result

Experiment Sparse vs Dense (N=1000)

Code

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import time
import matplotlib.pyplot as plt

N = 1000

np.random.seed(42)
B_sparse = np.random.rand(N)

# Sparse Matrix Construction
diagonals = [
    np.random.rand(N),
    np.random.rand(N-1),
    np.random.rand(N-1)
]
A_sparse = sp.diags(diagonals, offsets=[0, 1, -1], format='csr')

# Dense Matrix Construction (from sparse matrix)
A_dense = A_sparse.toarray()
B_dense = B_sparse.copy()

# Solve using spsolve (Sparse)
start_time = time.time()
X_sparse_spsolve = spla.spsolve(A_sparse, B_sparse)
sparse_spsolve_time = time.time() - start_time

# Solve using lsmr (Sparse)
start_time = time.time()
X_sparse_lsmr = spla.lsmr(A_sparse, B_sparse)[0]
sparse_lsmr_time = time.time() - start_time

# Solve using np.linalg.solve (Dense)
start_time = time.time()
X_dense_solve = np.linalg.solve(A_dense, B_dense)
dense_solve_time = time.time() - start_time

# Solve using np.linalg.lstsq (Dense least squares)
start_time = time.time()
X_dense_lstsq = np.linalg.lstsq(A_dense, B_dense, rcond=None)[0]
dense_lstsq_time = time.time() - start_time

# Plot results
labels = ["Sparse (spsolve)", "Sparse (lsmr)", "Dense (solve)", "Dense (lstsq)"]
times = [sparse_spsolve_time, sparse_lsmr_time, dense_solve_time, dense_lstsq_time]

plt.figure(figsize=(10, 5))
plt.bar(labels, times, color=["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728"])
plt.ylabel("Time (seconds)")
plt.title("Comparison of Computational Time: Sparse vs Dense Solvers")
plt.grid(axis="y", linestyle="--", alpha=0.7)
plt.savefig("test.jpg")

Result