Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LSMR fails with NaN output for trivial problems #20

Open
tvercaut opened this issue Dec 2, 2022 · 2 comments
Open

LSMR fails with NaN output for trivial problems #20

tvercaut opened this issue Dec 2, 2022 · 2 comments

Comments

@tvercaut
Copy link

tvercaut commented Dec 2, 2022

Thanks for the nice library. I wanted to try LSMR but my first attempt to use it with a trivial problem failed with NaN output.

Steps to reproduce:

import torch
import torchmin

A = torch.eye(10)
xtrue = torch.zeros((10, 1))
b = A @ xtrue

x = torchmin.lstsq.lsmr.lsmr(A,b)[0]
print(x)

which resulted in

tensor([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])

instead of 0s

@tvercaut
Copy link
Author

tvercaut commented Dec 2, 2022

Note that the same holds with a slightly less trivial but still trivial case:

import torch
import torchmin

A = torch.eye(10)
xtrue = torch.ones((10, 1))
b = A @ xtrue

x = torchmin.lstsq.lsmr.lsmr(A,b)[0]
print(x)

Comparing the source code with that of scipy points to the following issues.

normr and normar are initialiased differently and do not include the "convergence" tests at the very beginning

normar = b.new_tensor(0)
normr = b.new_tensor(0)

instead of
https://github.com/scipy/scipy/blob/2347d9309fbbabb1d3f89d35b7a42d0d53f002b2/scipy/sparse/linalg/_isolve/lsmr.py#L290-L307

Within the main iteration loop, the convergence test is only done every 10 iterations

# ------- Test for convergence --------
if itn % 10 == 0:

instead of every time
https://github.com/scipy/scipy/blob/2347d9309fbbabb1d3f89d35b7a42d0d53f002b2/scipy/sparse/linalg/_isolve/lsmr.py#L409

@tvercaut
Copy link
Author

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant