I have have a material model, defining the deviatoric stress for a nonlinear fluid: $\boldsymbol{\sigma}_{\mathrm{dev}} = f(\dot{\boldsymbol{\varepsilon}}_{\mathrm{dev}})$
Now I wish to find the strain-rate for a given stress (inverse problem solved with Newton iterations). However, I come across a problem with the tangent relation, \begin{equation} \boldsymbol{\mathsf E} := \frac{\partial \boldsymbol{\sigma}_{\mathrm{dev}}}{\partial\dot{\boldsymbol{\varepsilon}}_{\mathrm{dev}}} \end{equation} due to the fact that $\boldsymbol{\mathsf E}$ is singular, e.g. in Voigt notation the matrix would look like this for simple Newtonian flow: \begin{equation} \underline{E} = \mu\begin{pmatrix}\frac23 &-\frac13 &-\frac13& 0& 0& 0 \\ -\frac13 &\frac23 &-\frac13& 0& 0& 0\\ -\frac13 &-\frac13 &\frac23& 0& 0& 0 \\ 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & 1\\ \end{pmatrix} \end{equation} (my actual matrix is not as simple as in the example, but will be deviatoric)
Dealing with symmetry is simple enough, as you can just leave out the symmetric components. But with deviatoric tangents, the "inverse" I would like to find is actually: \begin{equation} \boldsymbol{\mathsf E}_\mathrm{inv} : \boldsymbol{\mathsf E} = \boldsymbol{\mathsf I}_\mathrm{dev} \end{equation} where $\boldsymbol{\mathsf I}_\mathrm{dev} := \boldsymbol{\mathsf I} - \frac13 \boldsymbol{I}\otimes\boldsymbol{I}$.
So my question is, how can I find $\boldsymbol{\mathsf E}_\mathrm{inv}$ ?
Currently I do it by adding a small number to the upper diagonal to make it invertible by normal procedure, but this ruins the convergence speed of my iterations (evaluating my $f$ is costly).
I've realized how increadibly simple the answer is, since \begin{equation} (\boldsymbol{I} \otimes \boldsymbol{I}) : \dot{\boldsymbol{\varepsilon}} = \boldsymbol{0} \end{equation} the volumetric part can be added to the tangent relation \begin{equation} \boldsymbol{\sigma} = \boldsymbol{\mathsf{E}} : \dot{\boldsymbol{\varepsilon}} = \underbrace{(\boldsymbol{\mathsf{E}} + \alpha\;\boldsymbol{I}\otimes\boldsymbol{I})}_{\hat{\boldsymbol{\mathsf E}}} : \dot{\boldsymbol{\varepsilon}} \end{equation} for some nonzero choice of $\alpha$, then construct \begin{equation} \boldsymbol{\mathsf E}_\mathrm{inv} = \hat{\boldsymbol{\mathsf E}}^{-1} \end{equation} which will be invertible and not cause any convergence problems.