Fix: SciPy Not Working — Sparse Matrix Errors, Optimization Convergence, and Import Issues

The Error

You import a SciPy function and the symbol isn’t there:

ImportError: cannot import name 'lena' from 'scipy.misc'
ImportError: cannot import name 'imread' from 'scipy.misc'

Or you build a sparse matrix and get a deprecation warning:

SparseEfficiencyWarning: The matrix subclass is not the recommended way to
represent matrices. Please use scipy.sparse.csr_array instead.

Or your optimization runs to the iteration limit without converging:

scipy.optimize.OptimizeResult:
    success: False
    message: 'Maximum number of function evaluations has been exceeded.'
    nit: 1000

Or scipy.integrate.quad returns garbage:

IntegrationWarning: The integral is probably divergent, or slowly convergent.
result: (42.0, 1.5e+08)   # Huge error bar — result is meaningless

Or you install SciPy and pip builds from source for an hour, then fails:

error: Fortran compiler not found

SciPy is a massive library — optimization, integration, statistics, signal processing, sparse matrices, linear algebra, interpolation, and more. Each submodule has its own evolution and deprecations. The migration from scipy.misc (mostly deprecated), the rise of _array types alongside _matrix types, and the 1.11+ API cleanup have changed a lot of commonly-copied StackOverflow answers.

Why This Happens

SciPy builds on NumPy and adds scientific algorithms. Many parts of SciPy predate NumPy best practices — scipy.misc held helper functions that have moved to better homes (scipy.datasets, scipy.signal, imageio). The scipy.sparse module had csr_matrix (subclass of np.matrix) as its main type for decades, but SciPy 1.8+ is transitioning to csr_array (subclass of np.ndarray) which behaves consistently with the rest of NumPy.

Optimization and integration problems often fail because the default algorithms assume well-behaved inputs. Real-world problems have local minima, steep gradients, or singular points that break defaults. Understanding which solver fits your problem is the difference between convergence and garbage.

Platform and Environment Differences

SciPy’s performance and even its API correctness depend heavily on which BLAS/LAPACK backend it’s linked against, which architecture it was compiled for, and whether NumPy was upgraded to 2.0. The same pip install scipy produces measurably different runtime behavior on different machines.

SciPy 1.11+ and NumPy 2 compatibility. SciPy 1.13+ supports NumPy 2.0 fully; SciPy 1.10 and earlier do not. If you upgrade NumPy to 2.x but keep an old SciPy, importing SciPy raises AttributeError: module 'numpy' has no attribute 'float_' because NumPy 2 removed deprecated dtype aliases. Pin both together: pip install "numpy>=2.0" "scipy>=1.13". The reverse — NumPy 1.x with SciPy 1.13+ — works because SciPy 1.13 has backward-compatibility shims. For dtype changes in NumPy 2 that ripple into SciPy code, see NumPy not working.

M1 native BLAS vs Rosetta. Apple Silicon Macs running native arm64 Python get SciPy linked against Apple’s Accelerate framework — fast, native, no extra install. Running x86_64 Python under Rosetta gets OpenBLAS via Rosetta translation, which is 3–5x slower for linear algebra. Check with np.show_config() and scipy.show_config() — Accelerate appears as accelerate_info on native arm64. The fix: install arm64 Python from python.org, then pip install --no-cache-dir scipy to force a fresh native wheel.

Conda vs pip wheel performance. conda install -c conda-forge scipy links against Intel MKL by default. pip install scipy links against OpenBLAS. On Intel CPUs, MKL is 1.5–3x faster for linalg.solve, linalg.svd, and matrix multiplication. On AMD Zen CPUs, the difference is smaller (5–20%) because MKL historically downgraded to slow paths on non-Intel chips — set MKL_DEBUG_CPU_TYPE=5 to force the AVX2 path. ARM machines get the same speed from both because MKL doesn’t run on ARM.

MKL vs OpenBLAS thread management. Both libraries spawn worker threads for matrix operations. By default they spawn N_CPU threads — when you call a SciPy function from within a multiprocessing pool, you get N_CPU * N_CPU total threads and the machine thrashes. Set OMP_NUM_THREADS=1 and MKL_NUM_THREADS=1 before importing SciPy when using multiprocessing, or use threadpoolctl to manage the count dynamically.

ARM64 sparse matrix support. scipy.sparse works identically on ARM and x86, but some BLAS-accelerated paths (sparse-dense matrix product) are slower on ARM because of less optimized BLAS kernels. For very large sparse problems on ARM servers, consider using scipy.sparse.linalg.cg (iterative solver) instead of spsolve (direct factorization) — the iterative method is less dependent on BLAS performance.

Windows wheel sizes. Windows SciPy wheels are ~50 MB because they bundle a complete LAPACK implementation. Linux wheels are smaller because they assume system BLAS exists. On a fresh Windows install with no Visual C++ Redistributable, the import succeeds (everything is statically linked) — this is the one case where Windows is easier than Linux for SciPy installation. For pip wheel build failures on unusual platforms, see pip could not build wheels.

Cross-compile failures. When deploying to ARM64 Lambda or Graviton EC2 from an x86 dev machine, pip install scipy inside a Linux x86 Docker image produces an x86 wheel that crashes on the ARM target. Use a multi-platform Docker build (docker buildx build --platform linux/arm64) or build directly on the target architecture. The crash is silent — Lambda just kills the function after init timeout.

Fix 1: scipy.misc Import Errors

ImportError: cannot import name 'lena' from 'scipy.misc'
ImportError: cannot import name 'imread' from 'scipy.misc'
ImportError: cannot import name 'imresize' from 'scipy.misc'

Most of scipy.misc was deprecated in 1.0 and removed in 1.10+.

Migration table:

Example — replacing scipy.misc.imread:

# OLD — removed in SciPy 1.3
from scipy.misc import imread
img = imread("photo.jpg")

# NEW — use Pillow
from PIL import Image
import numpy as np
img = np.array(Image.open("photo.jpg"))

# Or imageio
import imageio.v3 as iio
img = iio.imread("photo.jpg")

For test images that used to come from scipy.misc:

from scipy.datasets import face, ascent

img = face()      # RGB raccoon face, shape (768, 1024, 3)
img = ascent()    # Grayscale image, shape (512, 512)

scipy.datasets requires pooch which installs automatically with recent SciPy, but may need pip install pooch on older versions.

Fix 2: Sparse Matrix _matrix vs _array

SparseEfficiencyWarning: The matrix subclass is not the recommended way
to represent matrices. Please use csr_array instead.

SciPy is migrating from scipy.sparse.csr_matrix (which inherits from np.matrix) to scipy.sparse.csr_array (which inherits from np.ndarray). The _array versions behave like regular NumPy arrays — * is element-wise, @ is matrix multiplication.

Old pattern:

from scipy.sparse import csr_matrix
import numpy as np

# csr_matrix uses np.matrix conventions — * is matrix multiply
A = csr_matrix(np.eye(3))
B = csr_matrix(np.eye(3))

A * B   # Matrix multiplication (returns csr_matrix)
A.multiply(B)   # Element-wise

New pattern (SciPy 1.8+):

from scipy.sparse import csr_array
import numpy as np

# csr_array uses ndarray conventions — * is element-wise, @ is matrix multiply
A = csr_array(np.eye(3))
B = csr_array(np.eye(3))

A * B   # Element-wise (same as NumPy)
A @ B   # Matrix multiplication (same as NumPy)

Convert between types:

from scipy.sparse import csr_matrix, csr_array

old_matrix = csr_matrix(some_array)

# Convert to _array form
new_array = csr_array(old_matrix)

# Or go back
back_to_matrix = csr_matrix(new_array)

Sparse format reference:

Convert sparse to dense (only for small matrices):

sparse_mat = csr_array(np.eye(100))

dense = sparse_mat.toarray()   # Convert to numpy ndarray
# WARNING: materializes all zeros too — can OOM on large sparse matrices

Pro Tip: Always use csr_array for new code. SciPy will eventually deprecate and remove csr_matrix. Converting a codebase mid-project is painful because the multiplication semantics differ — fix it early or maintain compatibility by using .multiply() and @ consistently.

Fix 3: optimize.minimize Did Not Converge

from scipy.optimize import minimize

result = minimize(f, x0=[1, 1])
# message: 'Maximum number of function evaluations has been exceeded.'
# success: False

The default algorithm (BFGS) assumes a smooth, unconstrained function with well-behaved gradients. Many real problems break this assumption.

Choose the right algorithm:

from scipy.optimize import minimize

# For smooth unconstrained problems — BFGS (default), L-BFGS-B (memory-efficient)
result = minimize(f, x0, method='L-BFGS-B')

# For non-smooth functions — Nelder-Mead (gradient-free)
result = minimize(f, x0, method='Nelder-Mead')

# For problems with bounds — L-BFGS-B with bounds
result = minimize(f, x0, method='L-BFGS-B', bounds=[(0, 1), (0, 10)])

# For problems with constraints — SLSQP or trust-constr
from scipy.optimize import LinearConstraint, NonlinearConstraint

linear_con = LinearConstraint([1, 1], -np.inf, 1)   # x0 + x1 <= 1
result = minimize(f, x0, method='trust-constr', constraints=[linear_con])

# For global optima — differential evolution (no gradient required)
from scipy.optimize import differential_evolution
result = differential_evolution(f, bounds=[(0, 10), (0, 10)])

Provide analytical gradients for 10–100x speedup when possible:

from scipy.optimize import minimize
import numpy as np

def f(x):
    return x[0]**2 + x[1]**2

def grad_f(x):
    return np.array([2*x[0], 2*x[1]])

result = minimize(f, x0=[1, 1], jac=grad_f, method='BFGS')   # Uses analytical gradient

Increase iteration limits:

result = minimize(
    f, x0,
    method='L-BFGS-B',
    options={
        'maxiter': 10000,      # Maximum iterations
        'maxfun': 100000,       # Maximum function evaluations
        'gtol': 1e-8,           # Gradient tolerance (convergence criterion)
        'ftol': 1e-10,          # Function value tolerance
    },
)

Scale your variables if they span different orders of magnitude:

# Variables: price in $1000s, count in integers 0-1
# WRONG — BFGS struggles with mismatched scales
def f(x):
    price, count = x
    return (price - 50)**2 + (count - 0.5)**2

# CORRECT — normalize to similar scales
def f_scaled(x):
    price_scaled, count = x
    price = price_scaled * 100
    return (price - 50)**2 + (count - 0.5)**2

Common Mistake: Checking only result.x without checking result.success. A failed optimization still returns values — they’re just wrong. Always check:

result = minimize(f, x0)
if not result.success:
    raise RuntimeError(f"Optimization failed: {result.message}")
print(f"Converged at x = {result.x}")

Fix 4: scipy.integrate.quad Warnings

IntegrationWarning: The integral is probably divergent, or slowly convergent.

quad is adaptive quadrature — it works great for smooth integrands but struggles with singularities, oscillations, and unbounded intervals.

For integrals with singularities, split the interval:

from scipy.integrate import quad
import numpy as np

# Singular at x=0 — poor convergence
def f(x):
    return 1 / np.sqrt(x)

# WRONG — quad struggles near the singularity
result, err = quad(f, 0, 1)   # Warning

# CORRECT — use the `points` argument to mark singularities
result, err = quad(f, 0, 1, points=[0])   # Treats 0 as a breakpoint

# Or compute analytically near the singular point

For oscillatory integrals, use a specialized method:

from scipy.integrate import quad

# Oscillating integrand
def g(x):
    return np.cos(1000 * x) * np.exp(-x)

# WRONG — default quad misses oscillations
result, err = quad(g, 0, 10)

# CORRECT — weight by oscillatory factor
result, err = quad(lambda x: np.exp(-x), 0, 10, weight='cos', wvar=1000)

For infinite intervals:

# Use np.inf and -np.inf explicitly
result, err = quad(lambda x: np.exp(-x**2), -np.inf, np.inf)
# Should return ~sqrt(pi) = 1.7724...

Vector-valued integration:

from scipy.integrate import quad_vec

def f_vec(x):
    return np.array([np.sin(x), np.cos(x), x**2])

result, err = quad_vec(f_vec, 0, 1)

Fix 5: Interpolation and Extrapolation

from scipy.interpolate import interp1d

x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 4, 1, 3, 5])
f = interp1d(x, y)
f(10)   # ValueError: A value in x_new is above the interpolation range

Default interp1d doesn’t extrapolate. Provide fill_value or bounds_error=False:

from scipy.interpolate import interp1d

# Fill outside with NaN (silent, not an error)
f = interp1d(x, y, bounds_error=False, fill_value=np.nan)

# Fill with specific value
f = interp1d(x, y, bounds_error=False, fill_value=0)

# Extrapolate (linearly beyond the data)
f = interp1d(x, y, fill_value='extrapolate')

# Use different values above and below the range
f = interp1d(x, y, bounds_error=False, fill_value=(0, 100))  # below, above

interp1d is deprecated in SciPy 1.14+ — use the newer interfaces:

# OLD — still works but deprecated
from scipy.interpolate import interp1d
f = interp1d(x, y, kind='cubic')

# NEW — recommended replacement
from scipy.interpolate import make_interp_spline, CubicSpline, PchipInterpolator

# For cubic splines
spline = CubicSpline(x, y)

# For monotonic cubic Hermite (overshoot-free)
pchip = PchipInterpolator(x, y)

# For linear interpolation
from scipy.interpolate import RegularGridInterpolator
# (for N-D) or just use np.interp for 1D:
y_new = np.interp(x_new, x, y)   # 1D linear, no scipy needed

2D interpolation — interp2d is also deprecated:

# OLD
from scipy.interpolate import interp2d
f = interp2d(x, y, z)

# NEW
from scipy.interpolate import RegularGridInterpolator

interp = RegularGridInterpolator((x, y), z)
z_new = interp([[x1, y1], [x2, y2]])   # Points to evaluate at

Fix 6: scipy.stats API Changes

AttributeError: 'scipy.stats.norm' object has no attribute 'interval'
# Or
TypeError: interval() missing 1 required positional argument: 'confidence'

SciPy 1.9+ changed the interval() method signature. The keyword argument was renamed:

from scipy.stats import norm

# OLD
ci = norm.interval(alpha=0.95, loc=0, scale=1)

# NEW
ci = norm.interval(confidence=0.95, loc=0, scale=1)

Hypothesis tests in scipy.stats now return named tuples with proper fields:

from scipy.stats import ttest_ind

# Most tests return an object with .statistic and .pvalue
result = ttest_ind(group1, group2)
print(result.statistic)   # t statistic
print(result.pvalue)       # p-value
print(result.df)           # degrees of freedom (SciPy 1.11+)
print(result.confidence_interval(0.95))   # CI method

Common distributions:

from scipy import stats

# Normal distribution
stats.norm.pdf(x, loc=0, scale=1)       # PDF
stats.norm.cdf(x, loc=0, scale=1)       # CDF
stats.norm.ppf(0.95)                     # Inverse CDF (quantile)
stats.norm.rvs(size=100)                 # Random samples

# Student's t
stats.t.pdf(x, df=10)

# Chi-squared
stats.chi2.pdf(x, df=5)

# Fit a distribution to data
mu, sigma = stats.norm.fit(data)

Random number generation — prefer the new Generator API:

# OLD — uses global state
from scipy import stats
np.random.seed(42)
samples = stats.norm.rvs(size=100)

# NEW — explicit random state
rng = np.random.default_rng(seed=42)
samples = stats.norm.rvs(size=100, random_state=rng)

For NumPy random number patterns, see NumPy not working.

Fix 7: Installation — Fortran Compiler Missing

error: Fortran compiler not found.
error: library dfftpack has Fortran sources but no Fortran compiler found

SciPy contains Fortran code (LAPACK, FFTPACK, etc.) that needs compilation. Most users get pre-built wheels, but on unusual platforms pip falls back to source builds.

Install from a wheel — the normal path:

pip install --upgrade pip setuptools wheel
pip install scipy

For conda (usually easiest on Windows and older Linux):

conda install -c conda-forge scipy

If you must build from source, install a Fortran compiler:

# Debian/Ubuntu
sudo apt install gfortran liblapack-dev libblas-dev

# macOS
brew install gcc   # Provides gfortran

# Windows — much harder, prefer conda or wheels

For general pip build failures, see pip could not build wheels.

Fix 8: scipy.signal — Filter Design and FFT

from scipy.signal import butter, filtfilt
import numpy as np

# Design a Butterworth lowpass filter
fs = 1000       # Sample rate (Hz)
cutoff = 100    # Cutoff frequency (Hz)
order = 4

# Old API (SciPy <1.12) — returns tuple
# b, a = butter(order, cutoff, btype='low', fs=fs)

# Modern: use second-order sections for numerical stability
sos = butter(order, cutoff, btype='low', fs=fs, output='sos')

# Apply filter
from scipy.signal import sosfiltfilt
filtered = sosfiltfilt(sos, noisy_signal)

sosfiltfilt vs filtfilt — SOS (second-order sections) is more numerically stable for high-order filters. Use it by default.

FFT for signal analysis:

from scipy.fft import fft, fftfreq

# Faster than np.fft for large arrays
y_fft = fft(signal)
freqs = fftfreq(len(signal), d=1/fs)   # d = sample spacing

# Only positive frequencies (real signal)
mask = freqs >= 0
import matplotlib.pyplot as plt
plt.plot(freqs[mask], np.abs(y_fft[mask]))
plt.xlabel("Frequency (Hz)")

For statistical plotting and distribution visualization that pairs naturally with scipy.stats, see Seaborn not working.

Still Not Working?

Performance — Use SciPy Functions, Not Python Loops

SciPy routines are C/Fortran under the hood. A loop in Python calling scipy.stats.norm.pdf for each element is slow; call it once on an array:

import numpy as np
from scipy.stats import norm

x = np.linspace(-3, 3, 100000)

# SLOW
y = np.array([norm.pdf(xi) for xi in x])

# FAST — vectorized
y = norm.pdf(x)

Linear Algebra — scipy.linalg vs numpy.linalg

scipy.linalg has more functionality (eigenvalue problems, LU decomposition options, specialized solvers) but numpy.linalg is enough for most cases:

import numpy as np
from scipy import linalg

# Both work
x = np.linalg.solve(A, b)
x = linalg.solve(A, b)

# scipy.linalg adds options like assume_a='pos' for positive-definite matrices
x = linalg.solve(A, b, assume_a='pos')   # Faster for symmetric positive-definite

# Scipy has specialized decompositions NumPy lacks
L, low = linalg.lu_factor(A)
x = linalg.lu_solve((L, low), b)

Sparse Linear Solvers

For sparse systems, use scipy.sparse.linalg:

from scipy.sparse import csr_array
from scipy.sparse.linalg import spsolve, cg

A_sparse = csr_array(A)   # Convert to sparse if A has many zeros
x = spsolve(A_sparse, b)   # Direct solve (LU)

# Iterative methods for very large systems
x, info = cg(A_sparse, b)   # Conjugate gradient (for symmetric positive definite)

Integration with scikit-learn

Many scikit-learn internals use SciPy sparse matrices. When fit() returns surprising types, check whether you’re getting a sparse or dense result — .toarray() materializes the dense form (and can OOM on large data), while .tocsr() keeps it sparse.

OverflowError in scipy.stats Hypothesis Tests

For very small p-values, scipy.stats returns 0.0 instead of underflowing to a near-zero float. This is intentional but masks the magnitude of the result. Use the log_p versions where available (scipy.stats.norm.logsf, scipy.stats.chi2.logsf) to get the log-probability and reason about it without overflow. For Bayesian work where you need exact tail probabilities, switch to mpmath for arbitrary precision.

LinAlgError: Singular matrix From solve and inv

scipy.linalg.solve and inv fail when the matrix is singular or ill-conditioned. Check the condition number first: np.linalg.cond(A) — anything above 1e10 means the matrix is numerically singular even if mathematically invertible. Use the pseudoinverse (np.linalg.pinv) or a regularized solver (scipy.linalg.lstsq with cond=None) for robust solutions. In optimization, this often shows up because two parameters are perfectly correlated; drop one.

Scientific Plotting Inside Jupyter

When running SciPy inside Jupyter, plots from plt.show() block the kernel by default. Add %matplotlib inline at the top of the notebook for static images, or %matplotlib widget for interactive plots that don’t block. For Matplotlib backend selection and inline rendering issues that affect SciPy visualizations, see Matplotlib not working.