I'm trying to improve the speed of the following iteration to calculate $s_k$:
$$B_k^{-1} = \Bigg( I + \frac{s_{k}s_{k-1}^T}{||s_{k-1}||^2}\Bigg)...\Bigg(I+ \frac{s_1s_0^T}{||s_0||^2}\Bigg) B_0^{-1}\\s_{k+1} = -\frac{B_k^{-1}r_{k+1}}{1+\frac{s_k^TB_k^{-1}r_{k+1}}{||s_k||^2}}$$
If we calculate $B_k^{-1}$ explicitly then each step is roughly $O(n^4)$ to calculate $B_k^{-1}$. If instead we compute $B_k^{-1}r_{k+1}$ by expanding the product from the right, we do $k$ matrix-vector multiplications which is $O(n^3)$, a lot better.
Is there a way to find $s_k$ in the form of a low rank update using the value of previous $s_i$?