patrick-kidger
Bug: NormalCG jvp can do CG on singular systems
NormalCG applies CG to the normal equations A* A x = A* b, and CG's documentation requires A* A to be nonsingular. This is the case if A is full rank and tall, but is not if A is wide, even if A is full rank. The problem is the jvp requires a normalcg solve with A replaced by A*, which is wide if A is tall.
The fix is simple: the issue is only that the wrong normal equations are used for the wide case, here they are A A^* y = b, where x = A^* y. In the wide A case, x is the minimum norm solution, and in the tall case, x is the solution in the least squares sense. In both cases, x is the pseudoinverse solution. Fix should be trivial
See also #156
As, good spot. This is what we document... but it's also very much a footgun. 😄 In the same spirit as #156 I'd be happy to take a PR on this if you have the appetite?
One possibility to obviate both this issue and #156 is to make a lx.Normal wrapper taking a solver as input and returning another solver, and setting lx.NormalCholesky = lx.Normal(lx.Cholesky()) and lx.NormalCG() = lx.Normal(lx.CG()) as the exact same logic is needed there looks to little benefit of inlining the logic. This could in principle be useful for any solver, and in general upgrades one working on nonsingular positive definite matrices to one working on full rank matrices. I can submit a PR if you'd like to see this approach.
Ah, agreed. I like the sound of that approach. Very happy to take a PR on that!