Mathematical proof for why gradient descent algorithm always converges
It looks like you put some emphasis on mathematical rigor. It just happens my motto is currently "from math and larger than math".

Yes, there are many mathematical proofs for conclusions along the line that gradient descent method converges to local minimum when the steps are small enough. We will need some precise statement on the condition of $f$ and the steps taken. Let me quote from Wikipedia, 


  With certain assumptions on the function $F$ (for example, $F$ convex and $\nabla F$ Lipschitz) and particular choices of $\gamma$ (e.g., chosen either via a line search that satisfies the Wolfe conditions or the Barzilai-Borwein method shown as following),
  $$ \gamma_n=\frac{(x_n-x_{n-1})^T[\nabla F(x_n) - \nabla F(x_{n-1})]}
{||\nabla F(x_n) - \nabla F(x_{n-1})||^2}$$
  convergence to a local minimum can be guaranteed. When the function $F$ is convex, all local minima are also global minima, so in this case gradient descent can converge to the global solution.


For some actual statements and proofs in their full clear detail, you can read a course note by Ryan Tibshirani at CMU. Unfortunately, because of the required mathematical rigor,  they are too long to be included here.