# How lazy is too lazy? An exercise in pseudoinverses
**URL:** https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020
**Category:** Numerics
**Tags:** linearalgebra, lazy-evaluation
**Created:** [August 26, 2021, 7:33am UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020 "2021-08-26T07:33:16Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@artkuo](https://discourse.julialang.org/u/artkuo)
#### Post date: [August 26, 2021, 7:33am UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020/1 "2021-08-26T07:33:16Z")
Inspired by [this discussion](https://discourse.julialang.org/t/ols-difference-between-and-normal-equation-and-inv-or-pinv/66739), 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](https://github.com/artkuo/LazyPseudoinverse.jl):
```julia
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?
---
### Author:  [@RoyiAvital](https://discourse.julialang.org/u/RoyiAvital)
#### Post date: [August 26, 2021, 2:39pm UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020/2 "2021-08-26T14:39:04Z")
Pay attention that `A' * A` might not be a Positive Definite matrix and hence the Cholesky factorization might not work.
---
### Author:  [@artkuo](https://discourse.julialang.org/u/artkuo)
#### Post date: [August 26, 2021, 3:41pm UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020/3 "2021-08-26T15:41:19Z")
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.
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [August 27, 2021, 12:54pm UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020/4 "2021-08-27T12:54:06Z")
> [@artkuo](#):
>
> 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
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
```
---
### Author:  [@artkuo](https://discourse.julialang.org/u/artkuo)
#### Post date: [August 27, 2021, 11:16pm UTC](https://discourse.julialang.org/t/how-lazy-is-too-lazy-an-exercise-in-pseudoinverses/67020/5 "2021-08-27T23:16:07Z")
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](https://github.com/artkuo/LazyPseudoinverse.jl/tree/singularexception).) 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.