HOSVD in three steps

The way to compute a Tucker decomposition is called HOSVD, and it uses only tools you already have:

1. Unfold the tensor along each axis (Part I).
2. Run SVD on each unfolding; keep the top components. These are the factor matrices.
3. Contract the original tensor against all factor matrices to get the core (section 06).

Us = [np.linalg.svd(unfold(T, ax), full_matrices=False)[0] for ax in range(3)]
r = (2, 2, 3)
Us = [Us[i][:, :r[i]] for i in range(3)]
core  = np.einsum('ijk,ia,jb,kc->abc', T, Us[0], Us[1], Us[2])
recon = np.einsum('abc,ia,jb,kc->ijk', core, Us[0], Us[1], Us[2])