Another potential approach is to use laziness to gain some similar advantages to Linnea. I had a proof-of-concept (posted previously; apologies for self-promotion) that lazily recognizes inv(A'A)*A'b as a pseudo-inverse, similar to Fig 3 of the Linnea article. Except instead of a graph, it stores info such as positive definiteness and matrix solve in types, and then tries to do the right thing when the operations resolve into an expected pattern. Re-posting from before,
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
If the pattern ended up as inv(A'A)C, it would still recognize positive definiteness and use Cholesky factorization. In principle, any use of inv could be done lazily and automatically do \ when right-multiplied by something. The lazy approach is probably very limited in what it can do compared to Metatheory.jl or Linnea, but for simple things it is quite straightforward and doesn’t need to know the sizes in advance.