How to efficiently perform matrix inversion more than once $\left(A^TA + \mu I \right)^{-1}$ if $\mu$ is changing but $A$ is fixed? 
How to efficiently perform matrix inversion more than once $\left(A^TA + \mu I \right)^{-1}$ if $\mu \in \mathbb{R}$ is changing but $A \in \mathbb{R}^{n \times m}$ is fixed?


Or do I need to invert the whole matrix in every update of $\mu$?

 A: If the change in $\mu$ is small, you might use the series
$$ (B+\lambda I)^{-1} = (B+\mu I)^{-1} - (\lambda - \mu) (B + \mu I)^{-2} + (\lambda - \mu)^2 (B + \mu I)^{-3} - \ldots $$
A: Do you need to have the matrix in explicit form, or are you happy for it to be factored?  If the latter, then note that adding a constant to the diagonal of a matrix leaves the eigenvectors unchanged and shifts all the eigenvalues by that constant.  That's easy to see since if Ax=$\lambda$x, then (A+$\mu$I)x = ($\lambda+\mu$)x.  If you've done the eigendecomposition, then the new inverse is just V $(\Lambda + \mu I)^{-1}$ V$^T$.  Since the thing being inverted is diagonal, this is very fast.
A: To invert any matrix $\bf X$ you can reformulate it into
"Find the matrix which if multiplied by $\bf X$ gets closest to $\bf I$".
To solve that problem, we can set up the following equation system.
$$\min_{\bf v}\|{\bf M_{R(X)}} {\bf v} - \text{vec}({\bf I})\|_2^2$$
where $\text{vec}({\bf I})$ is vectorization of the identity matrix, $\bf v$ is vectorization of ${\bf X}^{-1}$ which we are solving for  and $\bf M_{R(X)}$ represents multiplication (from the right) by matrix $\bf X$. If no inverse exists, then this should still find a closest approximation to an inverse in the 2-norm sense.
This will be a linear least squares problem that should converge in worst case same number of iterations as matrix size.
But, if we run some iterative solver and have a good initial guess, then it can go much faster than that.
Information on how to construct $\bf M_{R(X)}$ should exist at the wikipedia entry for Kronecker Product.
If you look at this:
$$({\bf B^T\otimes A})\text{vec}({\bf X}) = \text{vec}({\bf AXB}) = \text{vec}({\bf C})$$
We can rewrite it: $$({\bf B^T\otimes A})\underset{\bf v}{\underbrace{\text{vec}({\bf X})}}-\underset{\text{vec}({\bf I})}{\underbrace{\text{vec}({\bf C})}} = {\bf 0}$$
Maybe now it becomes clearer what $\bf M_{R(X)}$ should be.

Edit:
If we combine this with an iterative Krylov subspace solver, we are allowed to choose an initial "guess" for the solution to the equation system.
So let us assume we have found one solution $X_1 = (A^TA + \mu_1 I)^{-1}$.
We can now use $X_1$ as initial guess when solving for $X_2 = (A^TA+\mu_2 I)^{-1}$
