import numpy as np
A = np.random.rand(4,4)
k = np.empty((4,4,4))
for ix in range(4):
    for iy in range(4):
        x = A[ix,]
        y = A[iy,]
        sx = np.power(x - x[:,np.newaxis],2)
        sy = np.power(y - y[:,np.newaxis],2)
        k[ix,iy] = (sx + sy).sum().T

k[ix,iy] = (sx + sy).sum(axis=-1)

xij = (A[:, None, :] - A[..., None])**2
k = np.sum(xij[:, None, :, :] + xij, axis=-1)