Let $A$ be a $n\times m$ matrix, $X$ be a $d \times m$ matrix and $b_i\in \mathbb R$, $i=1\dots n$. $$ f(x) = 4 XA^\top (I+\Delta_X)^{-3}(B+B\Delta_X-I)A $$ where $B=\sum_ie_ib_ie_i^\top$ and $\Delta_X=\sum_i e_ie_i^\top AX^\top XA^\top e_ie_i^\top$ are diagonal matrices.
I have found that $f(x)$ is the first derivative of $$ \sum_i \left(\frac{1}{1+\|Xa_i\|^2}-b_i\right)^2 $$
What is the derivative $\partial f / \partial X$ ?
I'm trying to get it from this: $$\tag{1}\label{eq1} \frac{\partial^2}{\partial x^2}\left(\frac{1}{1+g_i(x)}-b_i\right)^2 = 2\left(\frac{3-b_i-b_ig_i(x)}{(1+g_i(x))^4}g'(x)^2 + \frac{b_i+b_ig_i(x)-1}{(1+g_i(x))^3}g''(x)\right) $$ where $g_i(x)=\|xa_i\|^2$.
But I'm facing some difficulties with $g'(x)^2$, considering $g'(x)$ is a matrix. Besides, I'm not really sure \eqref{eq1} can help.
Before proceeding to the second derivative, let's check on that first derivative.
Define the diagonal matrix $C$ using the Hadamard product as follows: $$ \eqalign { C &= I\circ (AX^TXA^T) + I \cr } $$ And let $u$ denote the vector of all ones. Then your objective function can be reformulated as $$ \eqalign { F(X) &= \|C^{-1}u- b\|^2 \cr } $$ which is easier to work with. For the first derivative, I obtained $$ \eqalign { f(X) &= \frac {\partial F} {\partial X} \cr &= 4 XA^T \bigg[I\circ(C^{-1}bu^TC^{-1} - C^{-2}uu^TC^{-1})\bigg] A \cr } $$ I'm not sure if this expression is equal to yours or not.
Update with details
Defining $M\!=\!(C^{-1}u\!-\!b)$ for later convenience, let's begin by stating a couple of well-known differentials $$\eqalign{ dC^{-1} &= -C^{-1}\,dC\,C^{-1} \cr d\,\|M\|^2 &= d\,(M:M) = 2\,M:dM \cr } $$ and employ them to find $dF$.
$$\eqalign{ dM &= dC^{-1}u \cr dF &= 2\,M:dM \cr &= 2\,Mu^T:dC^{-1} \cr &= -2\,Mu^T : C^{-1}\,dC\,C^{-1} \cr &= -2\,C^{-1}Mu^TC^{-1} : dC \cr } $$ The next thing we need is the differential of $C$: $$\eqalign{ C &= I\circ\big(AX^TXA^T\big)+I \cr dC &= I\circ\big(A\,\,d(X^TX)\,\,A^T\big) \cr &= I\circ\big(A\,\,(dX^TX + X^TdX)\,\,A^T\big) \cr &= 2\,I\circ\big(A\,\,{\rm sym}(X^TdX)\,\,A^T\big) \cr } $$ I should pause and define the Frobenius product and the sym() operator $$\eqalign{ A:B &= {\rm tr}(A^TB) \cr {\rm sym}(A) &= \frac {1} {2} (A+A^T) \cr } $$ and a few of their more useful properties, especially with respect to the Hadamard product $$\eqalign{ (A\circ B):C &= A:(B\circ C) \cr A:{\rm sym}(B) &= {\rm sym}(A):B \cr {\rm sym}(I\circ A) &= I\circ A = A\circ I \cr A:BC &= B^TA:C \cr &= AC^T:B \cr } $$ Now substitute $dC$ back into the expansion of $dF$.
$$\eqalign{ dF &= -2\,C^{-1}Mu^TC^{-1} : \big(2\,I\circ(A\,\,{\rm sym}(X^TdX)\,\,A^T)\big) \cr &= -4(\,C^{-1}Mu^TC^{-1})\circ I : \big(A\,\,{\rm sym}(X^TdX)\,\,A^T\big) \cr &= -4\,A^T((\,C^{-1}Mu^TC^{-1})\circ I)A : {\rm sym}(X^TdX) \cr &= -4\,A^T((\,C^{-1}Mu^TC^{-1})\circ I)A : (X^TdX) \cr &= -4\,XA^T((\,C^{-1}Mu^TC^{-1})\circ I)A : dX \cr &= \bigg(\frac {\partial F} {\partial X}\bigg) : dX \cr } $$
Update #2
It occurs to me that since $C$ is diagonal and $u$ is a vector of ones, we can get rid of the Hadamard product in the sub-expression $$ \eqalign { I\circ(C^{-1}Mu^TC^{-1}) &= C^{-1}(I\circ(Mu^T))C^{-1} \cr &= C^{-1}{\rm diag}(M)C^{-1} \cr &= C^{-1}(C^{-1}-B)C^{-1} \cr &= C^{-3}(I - BC) \cr } $$ where $B = {\rm diag}(b)$, and all of the matrices commute since they're diagonal.
Plugging this back into the derivative expression $$ \eqalign { \frac {\partial F} {\partial X} &= 4\,XA^TC^{-3}(BC-I)A \cr } $$
I now see that setting $C=(I+\Delta_X)$ recovers your original expression for the derivative. $$ $$ Update #3: The Second Derivative
What we have so far is $$ \eqalign { f &= 4\,XA^TC^{-3}(BC-I)A \cr dC &= I\circ(A(dX^TX+X^TdX)A^T) \cr } $$ Since $f$ is a matrix, I prefer to denote it by $Y$.
Now the idea is to take the differential $$ \eqalign { dY &= 4(dX)A^TC^{-3}(BC-I)A + 4XA^TC^{-3}B(dC)A \cr &+\, 4XA^T\bigg(dC^{-3}\bigg)(BC-I)A \cr \cr &= 4(dX)A^TC^{-3}(BC-I)A + 4XA^TC^{-3}B(dC)A \cr &+\, 4XA^T\bigg(C^{-3}(dC)C^{-1} + C^{-2}(dC)C^{-2} + C^{-1}(dC)C^{-3}\bigg)(I-BC)A \cr } $$ And apply the ${\rm vec}()$ operation $$ \eqalign { {\rm vec}(AXB) &= (B^T\otimes A){\rm vec}(X) \cr } $$ to isolate all of the $dC, dX$ terms as vectors. $$ \eqalign { {\rm vec}(dY) &= 4((A^TC^{-3}(BC-I)A)^T\otimes I){\rm vec}(dX) \cr &+ 4\big(A^T(I-BC)C^{-1}\otimes C^{-3}AX^T\big){\rm vec}(dC) \cr &+ 4\big(A^T(I-BC)C^{-2}\otimes C^{-2}AX^T\big){\rm vec}(dC) \cr &+ 4\big(A^T(I-BC)C^{-3}\otimes C^{-1}AX^T\big){\rm vec}(dC) \cr &+ 4(A^T\otimes XA^TC^{-3}B){\rm vec}(dC) \cr } $$ For notational convenience, collect all the matrix coefficients and re-write this mess as a simple vector equation $$ \eqalign { dy &= G\,dx + H\,dc \cr } $$ But we're still not done, we need to express $dc$ in terms of $dx$. However, we must be careful with notation since ${\rm vec}(dX) = dx$, but ${\rm vec}(dX^T) = P\,dx \neq dx^T$.
The $P$ matrix is a permutation which depends on the dimension of $dx$. $$ \eqalign { dc &= {\rm vec}(I)\circ{\rm vec}(A\,d\!X^TXA^T+AX^Td\!XA^T) \cr &= {\rm vec}(I)\circ(AX^T\otimes A)Pdx + {\rm vec}(I)\circ(A\otimes AX^T)\,dx \cr &= {\rm diag}({\rm vec}(I))\,\bigg((AX^T\otimes A)P + (A\otimes AX^T)\bigg)\,dx \cr &= K\,dx \cr } $$
Finally, we're ready to write down an expression for the second derivative $$ \eqalign { \frac {\partial y} {\partial x} &= G + HK \cr &= \frac {\partial\,{\rm vec}(Y)} {\partial\,{\rm vec}(X)} \cr } $$