Inspired by this discussion, I worked up a self-education project to let a naive user type inv(A'A)A'b expecting to do a pseudoinverse by hand. But the package lazily avoids doing any inversion or matrix multiplication, and does the numerically superior A\b. A short demo of (unregistered) LazyPseudoinverse.jl:
julia> using LazyPseudoinverse
julia> A = [1 2 3; 4 5 6; 7. 8 10; 11 13 19];
julia> b = [5, 8, 10, 11.];
julia> inv(A'A)A'b # looks like normal matrix operation
3-element Vector{Float64}:
-4.679738562091531
8.222222222222245
-2.3333333333333317
julia> inv(A'A)A'b == A\b # it's really a QR solve!
true
A few other features:
inv(A'A) is lazy, and if materialized uses a Cholesky decomposition
inv(A'A)b is equivalent to (A'A)\b using Cholesky instead of LU inverse
inv(A'A)A' is lazy, and if materialized does a pinv
inv(A'A)B' also avoids inverse, and again uses Cholesky unless B === A
A'A\b also uses Cholesky instead of LU factorization
inv(A'A)A'b does none of the above, only A\b with QR
This project is not intended to be useful in any way except as a proof of concept, and for me to learn about parametric types. As such, the package isn’t and won’t be fully functional–I only defined enough methods to demo a few basic operations (in the README). However, a few questions:
I’d appreciate any pointers on smarter ways to do things, if anyone cares to look at the code (linked above).
I committed type piracy by hijacking Base.*(::Adjoint,::AbstractMatrix). This was to make A'A lazy and aware of its positive definiteness (PD). I felt pretty bad about doing this, but are there some examples of virtuous and unavoidable type piracy?
In the world of matrices, patterns like A'A and AA' occur often. Would it be helpful for Julia to recognize them automatically and take advantage of PD? I like how I is smart and automatically handles zeros(3,3)+I.
When to be lazy and when to be eager? There’s a (much better) PD package, PDMats.jl. It actually eagerly computes the Cholesky decomposition upon construction of a PDMat, making it immediately ready for inverting/solving. In the present project, I couldn’t actually decide whether it would be better to emulate that package, or to be lazy. Any thoughts on when to be eager? And how much laziness is too much laziness?
Yeah, I didn’t build in a special check for that. If A'A is singular, it looks like inv says ERROR: SingularException(3), but cholesky gives ERROR: PosDefException. So maybe a bit less informative. It might be possible to override that error for this case, and try to give the same-looking error.
That’s nice, check=false makes it easy to do cholesky when possible, and otherwise to fail with SingularException just like a traditional inv(A'A). (This is in its own branch.) Although I’m not sure it was wrong that my previous version threw PosDefException, since that is equally true as SingularException.
This does highlight an issue with trying to be too clever/lazy. A user could type inv(A'A) and get a SingularException, yet lazy inv(A'A)A' and inv(A'A)A'b both still work, because the respective algorithms pinv and \ are designed to give an answer no matter what. I guess the decision was made long ago for these not to yield error or warning.