a is (d,t), B is (k,d,t)

def _principal_angle(a, B):
  '''a is (d,t), B is (k,d,t)'''
  # TODO: check case for t = d-1
  if a.shape[1] == 1:
    return a.T.dot(B)[0,:,0]

  # find normals that maximize distance when projected
  x1 = np.einsum('abc,adc->abd', B, B).dot(a) - a   # b.dot(b.T).dot(a) - a
  x2 = np.einsum('ab,cad->cbd', a.dot(a.T), B) - B  # a.dot(a.T).dot(b) - b
  xx = np.vstack((x1, x2))

  # batch PCA (1st comp. only)
  xx -= xx.mean(axis=1)[:,None]
  c = np.einsum('abc,abd->acd', xx, xx)
  _, vecs = np.linalg.eigh(c)
  fpc = vecs[:,:,-1]
  fpc1 = fpc[:len(x1)]
  fpc2 = fpc[len(x1):]

  # a.dot(fpc1).dot(b.dot(fpc2))
  lhs = a.dot(fpc1.T).T
  rhs = np.einsum('abc,ac->ab', B, fpc2)
  return np.einsum('ij,ij->i', lhs, rhs)