I am trying to compute the partial derivatives of the following discrete time dynamic system with respect to $A$ and $B$ (for linearization purposes):
$$ B_{t+1} = yf(A_tx)^\dagger\\ A_{t+1} = f^{-1}(B_t^\dagger y)x^\dagger $$
Where $f$ is an element-wise invertible function, $x$ and $y$ are static fat matrices, and $A$ and $B$ are matrices of requisite size.
In Derivative of pseudoinverse with respect to original matrix, @user1551 introduces an $e$ and and $E$ that is used to help with the differentiation - but I haven't been able to follow well enough to do it on my own.
Define $F=f(AX)$, then it appears you are trying to solve $$Y = BF$$ using some kind of alternating directions iteration.
But given any $A$ you can solve explicitly for $B$ $$\eqalign{ B &= YF^++M(I-FF^+) \cr }$$ Conversely, given any $B$ you can solve for $A$ $$\eqalign{ F &= B^+Y+L(I-B^+B) \cr G &= f^{-1}(F) \cr A &= GX^++K(I-XX^+) \cr }$$ NB: In the above $(K,L,M)$ are arbitrary matrices. In your iteration it appears you have set them all to zero, i.e. $$\eqalign{ B &= YF^+ \cr F &= B^+Y \cr A &= GX^+\cr }$$
Although you can solve for either matrix given the other, in order to solve for them simultaneously you'll need a constraint equation.