a = np.random.rand(50, 40, 30, 6)
b = np.empty(a.shape[:-1] + (3, 3), dtype=a.dtype)
b[..., (0), :] = a[..., :3]
b[..., 1:, (0)] = a[..., 1:3]
b[..., (1), 1:] = a[..., 3:5]
b[..., (2), 1:] = a[..., 4:]
b[1, 2, 3, 2] = b[1, 2, 3, 0] + b[1, 2, 3, 1]
pi = gu_pinv(b)
Python, simultaneous pseudo-inversion of many 3x3, singular, symmetric, matrices