def gu_pinv(a, rcond=1e-15):
    a = np.asarray(a)
swap = np.arange(a.ndim)
swap[[-2, -1]] = swap[[-1, -2]]
u, s, v = np.linalg.svd(a)
cutoff = np.maximum.reduce(s, axis=-1, keepdims=True) * rcond
mask = s > cutoff
s[mask] = 1.0 / s[mask]
s[~mask] = 0
return np.einsum('...uv,...vw->...uw', np.transpose(v, swap) * s[..., (None
    ), :], np.transpose(u, swap))
Python, одновременная псевдоинверсия многих 3x3, сингулярных, симметричных матриц