如何在Python中比numpy.linalg.lstsq更快求解大规模稀疏超定线性方程组？

Great question! Since your matrix is sparse and the system is overdetermined (way more rows than columns), using numpy.linalg.lstsq is definitely not the most efficient choice—it’s optimized for dense matrices, so it’s wasting resources processing all those zero entries. Here are several targeted approaches to speed up your least squares solve:

1. Use Sparse Matrix Iterative Solvers (Best for General Cases) #

Scipy has dedicated solvers for large sparse least squares problems that avoid converting your matrix to a dense format. lsqr and lsmr are the go-to options here—they’re iterative methods designed specifically for overdetermined sparse systems, and they’ll handle your 150k×140 matrix with ease.

Example code:

import scipy.sparse as sp
from scipy.sparse.linalg import lsqr

# Convert your matrix to CSR format (efficient for row-based operations)
A_sparse = sp.csr_matrix(your_dense_or_raw_sparse_matrix)
# Solve Ax = b
x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var = lsqr(A_sparse, b)

lsqr converges quickly for systems with a small number of columns (like your 140 columns) and uses minimal memory compared to dense solvers.

2. Leverage the Overdetermined Structure with Normal Equations (Fastest if Stable) #

Since your system has way more rows than columns, you can compute the normal equations: (A^T A x = A^T b). The magic here is that (A^T A) is only a 140×140 dense matrix—solving this tiny system is nearly instantaneous, even with basic methods.

Note: This works best if your matrix has a reasonable condition number (not too ill-conditioned). If (A) is very ill-conditioned, the normal equations might amplify numerical errors, but for most real-world problems with 140 columns, this is the fastest approach by far.

Example code:

import numpy as np
import scipy.sparse as sp

A_sparse = sp.csr_matrix(your_matrix)
# Compute small dense normal equation components
AtA = A_sparse.T.dot(A_sparse).toarray()
Atb = A_sparse.T.dot(b)

# Solve with numpy's lstsq (or use Cholesky for even faster results)
x = np.linalg.lstsq(AtA, Atb, rcond=None)[0]

# Alternative: Cholesky decomposition (faster if AtA is positive definite)
L = np.linalg.cholesky(AtA)
x = np.linalg.solve(L.T, np.linalg.solve(L, Atb))

3. Add a Preconditioner (For Ill-Conditioned Systems) #

If your matrix is ill-conditioned (high acond value from lsqr), iterative solvers can converge slowly. A preconditioner can fix this by transforming the system into one with better conditioning. For sparse least squares, an incomplete LU (ILU) preconditioner on (A^T A) works well.

Example code:

from scipy.sparse.linalg import spilu, LinearOperator

# Build ILU preconditioner for A^T A
ilu = spilu(A_sparse.T.dot(A_sparse))
M = LinearOperator((AtA.shape[0], AtA.shape[0]), matvec=ilu.solve)

# Solve with preconditioned lsqr
x, istop, itn, ... = lsqr(A_sparse, b, M=M)

4. GPU Acceleration (If You Have Hardware) #

If you access to an NVIDIA GPU, using GPU-accelerated libraries like CuPy can drastically speed up the solve. CuPy mirrors SciPy’s API, so the code change is minimal:

import cupy as cp
from cupyx.scipy.sparse.linalg import lsqr

# Move data to GPU
A_cupy = cp.sparse.csr_matrix(A_sparse.get())
b_cupy = cp.array(b)

# Solve on GPU
x_cupy = lsqr(A_cupy, b_cupy)[0]

# Move result back to CPU if needed
x = cp.asnumpy(x_cupy)

Quick Recommendation #

Start with the normal equations approach—it’s the fastest for your column count, and you can quickly check if the solution is stable by computing the condition number of (AtA) with np.linalg.cond(AtA). If the condition number is larger than ~1e8, switch to lsqr with a preconditioner.

内容的提问来源于stack exchange，提问作者Ruan Putka