BUG: sparse.linalg.lsqr: variance returned appears systematically incorrect
Describe your issue.
I have encountered unusual behaviour with scipy.sparse.linalg.lsqr.
I am using the variance reported by LSQR (the estimated diagonals of (A.T@A)^(-1)) but the estimated variance is too low. The level that it is too low seems to depend on the number of data points relative to the number of model parameters. For ~1,000 data points and ~500 model parameters the variance reported by SciPy's implementation of LSQR is low by a factor of 5.
In the figure below I am computing the diagonals of (A.T@A)^(-1) directly and comparing to what scipy.sparse.linalg.lsqr returns. If everything was working correctly, I should expect to see this line hover around y=1. Less than 1 indicates LSQR is under-estimating the uncertainties. I set extremely strict tolerances for convergence, and it does not appear that lsqr is stopping early. For small numbers of parameters, LSQR provides the correct variance estimates. As the number of parameters grows, LSQR underestimates the variance by between a factor of 2x or 5x in these examples. However, as the number of model parameters increases further and gets close to the number of data points, the variance from LSQR becomes closer to the truth.
I plotted the summary statistics returned from LSQR (second figure) but they did not obviously explain the behaviour I see.
I could not find any tests in scipy's implementation that verifies whether the estimated variance is realistic or not.
If anyone knows something about this, I would appreciate your input.
Reproducing Code Example
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import linalg
from scipy.linalg import cho_factor, cho_solve
np.random.seed(0)
fig2, ax_metas = plt.subplots(2, 4)
fig, ax = plt.subplots()
for n in (250, 501, 1001):
ms = []
ratios = []
metas = []
pylops_ratios = []
for m in np.unique(np.logspace(0, np.log10(n), 100).astype(int)):
if m > n:
continue
A = np.random.randn(n, m)
x_true = np.random.randn(m)
#noise, var = (0, np.ones(n))
var = np.random.uniform(0.5, 1.5, n)
noise = np.random.randn(n) * 0.1
y = A @ x_true + noise
C = np.diag(var)
Cinv = np.diag(1/var)
# Method 1: Direct solve
direct_cov = np.linalg.solve(A.T @ Cinv @ A, np.eye(m))
direct_var = np.diag(direct_cov)
# Method 2: LSQR with Cholesky
L = np.linalg.cholesky(C)
Linv = np.linalg.inv(L)
A_white = Linv @ A
y_white = Linv @ y
x_chol, *meta, var_chol = linalg.lsqr(A_white, y_white, calc_var=True, show=False, atol=1e-12, btol=1e-12, damp=1e-6, iter_lim=10000000)
ms.append(m)
ratios.append(np.mean(var_chol/direct_var))
metas.append(meta)
ms = np.array(ms)
metas = np.array(metas)
for j, ax_meta in enumerate(ax_metas.flat):
ax_meta.scatter(ms/n, metas[:, j], label=f"n={n}")
ax.scatter(ms/n, ratios, label=f"n={n}")
labels = ("istop", "itn", "r1norm", "r2norm", "anorm", "acond", "arnorm", "xnorm")
for j, ax_meta in enumerate(ax_metas.flat):
ax_meta.set_ylabel(labels[j])
ax_meta.set_xlabel(f"number of parameters / number of data points")
ax.legend()
ax.set_xlabel(f"number of parameters / number of data points")
ax.set_ylabel(f"mean of [(diagonal from scipy.sparse.linalg.lsqr) / (diagonal from direct solve)]")
fig.suptitle("Computing diagonals of (A.T @ C^(-1) @ A)^(-1) from LSQR")Error message
N/ASciPy/NumPy/Python version and system information
1.14.1 2.1.2 sys.version_info(major=3, minor=12, micro=5, releaselevel='final', serial=0)
Build Dependencies:
blas:
detection method: extraframeworks
found: true
include directory: unknown
lib directory: unknown
name: accelerate
openblas configuration: unknown
pc file directory: unknown
version: unknown
lapack:
detection method: extraframeworks
found: true
include directory: unknown
lib directory: unknown
name: accelerate
openblas configuration: unknown
pc file directory: unknown
version: unknown
pybind11:
detection method: config-tool
include directory: unknown
name: pybind11
version: 2.12.0
Compilers:
c:
commands: cc
linker: ld64
name: clang
version: 15.0.0
c++:
commands: c++
linker: ld64
name: clang
version: 15.0.0
cython:
commands: cython
linker: cython
name: cython
version: 3.0.11
fortran:
commands: gfortran
linker: ld64
name: gcc
version: 13.3.0
pythran:
include directory: ../../../../../../private/var/folders/zq/1mkv5ktx0_5cz1rlynrmdmlh0000gn/T/pip-build-env-pw6uu9w1/overlay/lib/python3.12/site-packages/pythran
version: 0.16.1
Machine Information:
build:
cpu: x86_64
endian: little
family: x86_64
system: darwin
cross-compiled: false
host:
cpu: x86_64
endian: little
family: x86_64
system: darwin
Python Information:
path: /private/var/folders/zq/1mkv5ktx0_5cz1rlynrmdmlh0000gn/T/cibw-run-mpfy6fi_/cp312-macosx_x86_64/build/venv/bin/python
version: '3.12'Source: scipy/scipy