The update for the Hessian using BFGS is given by: $$H_{k+1}=(I-\rho_ks_ky_k^T)H_k(I-\rho_k y_ks_k^T)+\rho_ks_k s_k^T$$
where $\rho_k=\dfrac{1}{y_k^Ts_k}$.
Nocedal and Wright, Numerical Optimization say that that update runs in $O(n^2)$, but how is that possible if the update has matrix multiplication? Namely, $(I-\rho_ks_ky_k^T)$ is a matrix that has to be multiplied with $H_k$. Matrix multiplication is $O(n^3)$ (It can be efficiently computed and it is less than $O(n^3)$, but certainly no $O(n^2)$). So, how is it possible to have runtime $O(n^2)$? I know I must be misunderstanding something. I'd appreciate if someone can clarify.