Use rank-1 update blas kernel for quasi-Newton updates
Opened this issue · 7 comments
We can get some nice speedup if we use LinearAlgebra.BLAS.ger! to do the rank-1 updates of the approximate Hessian or inverse Hessian in quasi-Newton methods instead of H = H + (w * w') / θ which is much slower for large matrices.
We can use fallback implementation though for non-BLAS types.
A bit more advanced would be performing rank-1 updates over the decomposition of H directly, for Direct approximations, to avoid re-factorizing it every time we call \. We have lowrankupdate for Cholesky.
Yes, thanks for opening this issue. You are touching things that I have planned for, but pushed forward in favor of some other functionality, but yes, those things would be cool.
- Standard implementation that's closer to the linear algebra/math to work as fallback
- Fast implementations for
Array's of supported types (say Float32 and Float64, whatever is available) with dispatch - A factorized version of H that is updated. I originally thought I was going to prioritize this, but then I read some not so great benchmarks and pushed it down the list. But if we already have lowrankupdate (I didn't know that) for Cholesky, that seems like a no-brainer.
We also need L-BFGS and maybe even variants of it. We also need L-SR1 which seems nice is TR contexts (there are even closed forms for the spectral decomposition in the memoryless versions).
Yes I can perhaps tackle those issues at some point once the package is more stable. One thing is that for small matrices, the BLAS kernel is actually slower than naive implementation, so we need to keep that in mind.
I think such a check for size should be negligable?
Yes for all reasonably sized matrices it should be.
Your objective and gradient calculation would have to be overwhelmingly simple for such a branch to have any effect I'd say. Either way, we can just have it as an option. Auto will check, but specific choices could in principle be inferred if necessary.