# How lazy is too lazy? An exercise in pseudoinverses

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:

1. Iâ€™d appreciate any pointers on smarter ways to do things, if anyone cares to look at the code (linked above).
2. 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?
3. 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`.
4. 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?
3 Likes

Pay attention that `A' * A` might not be a Positive Definite matrix and hence the Cholesky factorization might not work.

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.

Use `cholesky(x, check=false)` and `issucess` to manually check whether it succeeded.

``````julia> using LinearAlgebra

julia> S = rand(4,3) |> x -> x * x'
4Ă—4 Matrix{Float64}:
1.64844   1.5099    0.852409  0.949308
1.5099    1.51767   1.07208   0.850572
0.852409  1.07208   1.16213   0.555333
0.949308  0.850572  0.555333  0.671274

julia> cholesky(S)
ERROR: PosDefException: matrix is not positive definite; Cholesky factorization failed.
Stacktrace:
[1] checkpositivedefinite
@ ~/Documents/languages/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/factorization.jl:18 [inlined]
[2] cholesky!(A::Hermitian{Float64, Matrix{Float64}}, ::Val{false}; check::Bool)
@ LinearAlgebra ~/Documents/languages/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:266
[3] cholesky!(A::Matrix{Float64}, ::Val{false}; check::Bool)
@ LinearAlgebra ~/Documents/languages/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:298
[4] cholesky(A::Matrix{Float64}, ::Val{false}; check::Bool)
@ LinearAlgebra ~/Documents/languages/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:394
[5] cholesky (repeats 2 times)
@ ~/Documents/languages/julia/usr/share/julia/stdlib/v1.8/LinearAlgebra/src/cholesky.jl:394 [inlined]
[6] top-level scope
@ REPL[30]:1

julia> c = cholesky(S, check=false)
Failed factorization of type Cholesky{Float64, Matrix{Float64}}

julia> issuccess(c)
false
``````

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.