Solution (source code)

= Solution

Write the inverse <tensor> as
$$
Y_{ij}=\delta_{ij}+\lambda A_{ij}+\lambda^2C_{ij}+\lambda^3D_{ij}+O(\lambda^4).
$$
Insert this into $(\delta_{ik}+\lambda X_{ik})Y_{kj}=\delta_{ij}$ and equate powers of $\lambda$. One obtains successively
$$
A_{ij}=-X_{ij},\qquad C_{ij}=-X_{ik}A_{kj}=X_{ik}X_{kj},\qquad
D_{ij}=-X_{ik}C_{kj}=-X_{ik}X_{kl}X_{lj}.
$$
Thus the <Neumann series> gives
$$
\boxed{Y_{ij}=\delta_{ij}-\lambda X_{ij}+\lambda^2X_{ik}X_{kj}-\lambda^3X_{ik}X_{kl}X_{lj}+O(\lambda^4).}
$$
Every contraction is an ordinary ordered <matrix product>. No transpose or symmetry assumption is used. Multiplying the truncated series on either side by $I+\lambda X$ leaves $I-\lambda^4X^4$, which verifies the error order. For a fixed finite matrix this expansion converges whenever $|\lambda|\,\|X\|<1$ in a submultiplicative <matrix norm>.