In a departure from Guyon & Males (2017)[1]}, we solve this with a least squares pseudo-inversion[2]}, with a regularization constant \(\alpha \)  (which we set to \(\alpha =1\)  for simulations.)
\(\mathbf {F} &=& ((\mathbf {D}^T)^\dagger \mathbf {P}^T)^T \\\mathbf {F} &=& \mathbf {P}\mathbf {D}^T(\mathbf {D}\mathbf {D}^T + \alpha \mathbf {I})^{-1}\) 
