`lssolve` fails where `linsolve` succeeds on `I x = [1, 0]`

[gdalle]

Am I just holding it wrong?

julia> using KrylovKit, LinearAlgebra

julia> linsolve(I(2), [1.0, 0.0])
([1.0, 0.0], ConvergenceInfo: one converged value after 1 iterations and 2 applications of the linear map;
norms of residuals are given by 0.00e+00.)

julia> lssolve(I(2), [1.0, 0.0])
┌ Warning: LSMR lssolve finished without converging after 100 iterations:
│ * ‖ b - A x ‖ = NaN
│ * ‖ x ‖ = NaN
│ * ‖ Aᴴ(b - A x) - λ^2 x ‖ = NaN
│ * number of operations = 201
└ @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/lssolve/lsmr.jl:134
([NaN, NaN], ConvergenceInfo: no converged values after 100 iterations and 201 applications of the linear map;
norms of residuals are given by NaN.)

[Jutho]

Thanks for reporting, this is indeed a bug where a vector becomes exactly zero and is than divided by its norm (which is not very likely in realistic examples, but should not happen nonetheless). I've already added a quick fix, but I still also want to add a test.

[gdalle]

I encountered it during implicit differentiation, when backpropagating basis vectors through a diagonal Jacobian. It didn't feel like such an outlandish use ^^

[Jutho]

No definitely not, that's not what I wanted to say (though I guess with a diagonal matrix, one wouldn't typically use a Krylov method? 😄).

[gdalle]

Thanks for the fix!
In implicit differentiation we don't know the structure of the Jacobian, and the main reason to use Krylov methods is because we want to deal with operators instead of matrices. Diagonal Jacobians are of course a trivial case, but they're good for testing!