I've rewritten the question in math notation, and I've left the old version in physics bra-ket notation here.
Background
A simple consequence of the singular value decomposition is that any vector $v$ in a vector space $V$ formed by the tensor product of two smaller spaces ("subsystems") $U$ and $W$ of dimension $d_U$ and $d_W$,
$v \in V = U \otimes W$,
has a special decomposition in terms of rank-one tensors (aka product states)
$v = \sum_{i=1}^d \lambda_i u_i \otimes w_i $, $\qquad u_i \in U$, $\qquad w_i \in W$, $\qquad d = \mathrm{min}(d_U,d_W)$
built from the fixed orthonormal bases $\{ u_i \}$ and $\{ w_i \}$.
This decomposition has many nice properties, and it's simple to see that there's no way to generalize it to the case of 3 or more subsystems while keeping all of them. For instance, a generic state in $V = \bigotimes_{n=1}^N V_n$ cannot be expressed as a sum of fewer than $\tilde{d}$ simple tensors, with $\tilde{d} \gg d_n = \mathrm{dim}V_n$ for all $n=1,\ldots,N$.
Vague Question
Is there a useful/natural/canonical decomposition of a vector into a "small" number of orthogonal simple tensors for the case of 3 or more subsystems? At the very least, the GHZ state should take it's canonical form in such a decompostion:
$v_{GHZ} = a_1 \otimes \cdots \otimes a_N + b_1 \otimes \cdots \otimes b_N$, $\qquad a_n, b_n \in V_n$, $\qquad \langle a_n ; b_n \rangle = 0$ for all $n$.
Specific Question
If we guess that the entropy function
$H[\{ p_i \} ] = -\sum_i p_i \ln p_i$
is appropriate, we can define the "minimum-entropy product-state decomposition" to be the decomposition
$v = \sum_{i=1}^{\tilde{d}} \lambda_i v_i$, $\qquad v_i = \bigotimes_{n=1}^N \psi_i^n$, $\qquad \psi_i^n \in V_n$
with the minimum value for the entropy $H[\{ p_i = \lambda_i^2 \} ]$, under the condition that the $\{ v_i \}$ are orthonormal. Note that we allow $\langle \psi_i^n ; \psi_j^n \rangle \neq 0$ for $i \neq j$.
The natural questions to ask are: Is this decomposition generically unique (i.e. unique except for some set of measure zero in the global vector space)? Is it continuous? (What other properties should this decomposition have to satisfy to be useful?)
I am 99% sure that in the case of $N=2$, this reduces to the Schmidt decomposition and that the answers to both questions is "yes".
Is any of this sensitive to our choice of the entropy function $H$, as opposed to some other permutation-invariant and majoritization-preserving function of the spectrum?
Best Answer
Vague question:
The hierarchical higher order SVD would lead you to the canonical form for the GHZ state. When you have $H_1 \otimes H_2 \otimes \ldots \otimes H_n$ you basically first Schmidt decompose your state $\psi$ relative to $H_1 \otimes (H_2 \otimes \ldots \otimes H_N)$ and then continue like this until you finish.
In more generality, the Tucker decomposition is the most general multilinear generalization of the SVD. It is not, however, unique. The hierarchical higher order SVD I mentioned above is a special case of the Tucker decomposition.
Specific question:
I don't know the answer to your question, but I did notice this: If you start with $\psi = \sum_{jk} \gamma_{jk} \alpha_j \otimes \beta_k$ then your entropy is essentially the entropy of the post measurement mixed state with measurement basis being $\alpha_j \otimes \beta_k$. So what you want to look at is which local measurement scheme gives you the minimum post-measurement entropy. So in the bipartite case all you need to do is take a partial trace and whatever the entropy of that matrix, that is the minimum post-measurement entropy. I.e. you get entanglement, which confirms your suspicion that the Schmidt decomposition is what you get in bipartite case. A similar strategy might work in the multipartite case.