cornellius-gp / linear_operator

A LinearOperator implementation to wrap the numerical nuts and bolts of GPyTorch
MIT License
95 stars 28 forks source link

[Docs] Example in readme is giving `NotPSDError` #99

Open jakee417 opened 3 months ago

jakee417 commented 3 months ago

📚 Documentation/Examples

Is documentation wrong?

In readme example for LinearOperator: https://github.com/cornellius-gp/linear_operator?tab=readme-ov-file#with-linearoperator

from linear_operator.operators import DiagLinearOperator, LowRankRootLinearOperator
# C = torch.randn(1000, 20)
# d = torch.randn(1000)
# b = torch.randn(1000)
A = LowRankRootLinearOperator(C) + DiagLinearOperator(d)  # represents C C^T + diag(d)
torch.linalg.solve(A, b)  # computes A^{-1} b efficiently!

will give

NotPSDError: Matrix not positive definite after repeatedly adding jitter up to 1.0e-06.

Think you know how to fix the docs? (If so, we'd love a pull request from you!)

from linear_operator.operators import DiagLinearOperator, LowRankRootLinearOperator
C = torch.randn(1000, 20)
d = torch.ones(1000) * 1e-9
b = torch.randn(1000)
A = LowRankRootLinearOperator(C) + DiagLinearOperator(d)
torch.linalg.solve(A, b)

should solve the issue

(Link to LinearOperator documentation)

Balandat commented 2 months ago

Hmm yeah d should just not be sampled from a normal distribution - you'll end up adding negative values to the diagonal. Something like d = 1e-9 + torch.rand(1000) should be safe.