Generic decomposition for A \ x

I have an algorithm where I have to perform A * x and A \ x a lot of times for the same A::AbstractMatrix. Importantly, A can be dense, triangular, sparse, etc.

I would like to save a generic decomposition object for A which makes this fast and numerically stable.

First I tried lu(A), but that does not work for eg LowerTriangular. So now I am special-casing for these things. Is there a generic solution?

Maybe add a specialization for lu(::LowerTriangular) in the stdlib? It’s hard to imagine a better generic solution than lu.

1 Like

Except maybe QR?

For unequal sizes, sure.

A is square (sorry, forgot to say).

In fact for A::Diagonal or A::AbstractTriangular, I can just leave it alone. Would an lu method that does this make sense? Or complement with an L = I. I am not sure.

Isn’t that the point of factorize?

1 Like

That’s a good point, unfortunately it is not type stable. A version that would be less eager to find the best solution could work better for me.

If you’re factoring a matrix I’m guessing the dynamic dispatch overhead is the least of your problems. Can’t you get by with a function barrier or a parametric type?

Thanks, I will benchmarks and see how it goes. This is all for CI code by the way, where I am also checking inference, but I think putting the function barrier in the right place could work around this.