Consider a linear, differentiable, and invertible function $f_\theta:\mathbb{R}^n\rightarrow \mathbb{R}^n$, with parameters $\theta$, with exactly one fixed point $z_0 = f_\theta(z_0)$. At that fixed point, I'd like to calculate the set of derivatives $\frac{\partial z_0}{\partial \theta}$. Given the function's properties and the Jacobian inverse identity, $J_{f^{-1}}=J_f^{-1}$, is there a computationally efficient way to calculate $\frac{\partial z_0}{\partial \theta}$?
This problem has come up recently in my work. It's very similar to the problem approached in [1] (Supplement A, Theorem 1) and [2] (Section 2.1), which end up with the equation $$\frac{\partial z}{\partial \theta}=\frac{\partial f}{\partial \theta} [I-\frac{\partial f}{\partial z}]^{-1}$$ The main bottleneck here is solving $[I-\frac{\partial f}{\partial z}]^{-1}$. Because the full matrix inverse is $O(n^3)$, and n is large for practical problems, these references solve for this inverse with iterative methods like a Neumann Series[2]. In my problem, however, f is both linear and invertible, which I'm hoping could simplify or speed up the computational process of finding $\frac{\partial z_0}{\partial \theta}$.
For instance, $J_f^{-1}$ can be quickly calculated with the Jacobian inverse identity by simply differentiating $f^{-1}$, which for this function is computationally easy. Using this with the Woodbury identity [3] to find $[I-\frac{\partial f}{\partial z}]^{-1}$ doesn't appear to help with computational complexity, however, since the identity matrix is of course full rank and we thus must perform another full inverse. Is there another approach to the problem that I'm missing, or are the function's properties simply not helpful in this case?
Thanks for your help!
[1] https://arxiv.org/pdf/1909.01377.pdf