Higher order SVD
Higher-order SVD (in sense of Tucker decomposition) of the matrix $M$ with dimensions $d_1\times d_2\times\cdots\times d_n$ is
$$ M_{i_1,i_2,\dots,i_N} = \sum_{j_1} \sum_{j_2}\cdots \sum_{j_N} s_{j_1,j_2,\dots,j_N} u^{(1)}_{i_1,j_1} u^{(2)}_{i_2,j_2} \dots u^{(N)}_{i_N,j_N}, $$ where $s$ is the core tensor and $u^{(i)}$ is the orthogonal matrix.
The matrix $u^{(i)}$ is left singular vectors of regular SVD of matrix $T$ with dimensions $d_i \times\prod_{j\neq i}d_j$ ($T$ is Flatten of $M$ over all dimensions except $i$-th). Therefore, high-order SVD is just a composition of regular SVDs over each dimension of $M$
hoSVD[m_] := With[{d = ArrayDepth[m]},
{Fold[Transpose[#1, RotateRight@Range[d]].#2 &, m, #], ConjugateTranspose /@ #} &@
Table[Conjugate@First@SingularValueDecomposition@Flatten[m, {{i}, Delete[Range[d], i]}], {i, d}]]
It returns {s, u} where s is the core tensor and u is the list of orthogonal matrices
If $n=2$ it is regular SVD
m = RandomReal[NormalDistribution[], {3, 4}];
{s, u} = hoSVD[m];
Max@Abs[m - Transpose@u[[1]].s.u[[2]]]
1.11022*10^-15
To generalize to the multidimensional case we need a generalization of the Dot
dot[s_, u_, n_] := Transpose[Transpose[s, Ordering[#]].u, #] &@Append[#, n] &@
Delete[#, n] &@Range@ArrayDepth[s];
It is Dot over n-th dimension of s.
Now we can write
m2 = dot[dot[s, u[[1]], 1], u[[2]], 2];
Max@Abs[m - m2]
1.11022*10^-15
Multidimensional case
m = RandomReal[NormalDistribution[], {3, 4, 5}];
{s, u} = hoSVD[m];
m2 = dot[dot[dot[s, u[[1]], 1], u[[2]], 2], u[[3]], 3];
Max@Abs[m - m2]
3.9968*10^-15
m = RandomReal[NormalDistribution[], {3, 4, 5, 6}];
{s, u} = hoSVD[m];
Max@Abs[Fold[dot[#, #2[[1]], #2[[2]]] &, s, Transpose@{u, Range@Length[u]}] - m]
6.21725*10^-15
I worked on this topic today as I'm interested in it for an interpolation problem.
I followed the interlaced computation from the Wikipedia link in the other answers (corresponding paper is here), and improved ybeltukov's answer. This allows the decomposition to be faster when the number of points is high and ranks are low, which can be a common case.
This paper allowed me to understand how to use the decomposition for interpolation. Interpolating a tensor is basically the product of 1D-interpolations of singular vectors across indices (ie. of different columns of the us below) weighted by the core tensor (s). Given that the number of non zero values of the core tensor can be very small compared to the number of dimensions (if the data is not random), the interpolation algorithm can be made very fast.
We have for example
$$ M(x, y, z) = \sum_{j_1} \sum_{j_2} \sum_{j_3} s_{j_1,j_2,j_3} u^{(1)}_{j_1}(x) u^{(2)}_{j_2}(y) u^{(3)}_{j_3}(z) $$
This can actually be useful for integrating a function as well, by integrating splines for example. So if the core tensor is indeed small it's possible to do high dimensional integrations very fast.
Here is the code
defaultTol=1.*^-10;
rankFunction[mat_,tol_:defaultTol]:=Length@Select[Diagonal@mat,Abs@#>tol&];
hoSVD2[m_,maxRank_:Infinity,rankF_:rankFunction]:=
Module[{d,us,Ak,r,u2,u,w,v,rotatePermutation,flateningDimensions},
d=ArrayDepth[m];
flateningDimensions={{1},Delete[Range[d],1]};
rotatePermutation=RotateRight@Range[d];
Ak=m;
us=
Table[
(*we flatten the dimensions other than the first one to have a matrix,and run an SVD on it*)
{u,w,v}=SingularValueDecomposition[Flatten[Ak,flateningDimensions],UpTo[maxRank]];
(*we keep the singular vectors corresponding to non zero singular values*)
r=rankF@w;
u2=u[[All,;;r]];
(*we rotate the kth dimension with u2,and transpose Ak for the next iteration
in order to have the (k+1)th dimension as the leftmost one,
multiplying with u2 can decrease the amount of data of Ak,thus accelerating next iterations*)
Ak=Transpose[(ConjugateTranspose@u2).Ak,rotatePermutation];
u2
,
{k,d}
];
(*s=Ak;*)
{Ak,us}
];
(*u(n).s, multiplying s on the left across index n*)
dot2[s_,u_,n_]:=
Module[{perm,permInv},
permInv=Range@ArrayDepth[s]//Delete[#,n]&//Prepend[#,n]&;
perm=Ordering[permInv];
Transpose[u.Transpose[s,perm],permInv]
];
(*Interpolation function*)
interpolateVector[v_,x_]:=Interpolation[v][x];
interpU[u_,i_,j_,x_]:=interpolateVector[u[[i,All,j]],x[[i]]];
interpS[s_,spos_,u_,x_]:=s[[Sequence@@spos]]*Times@@Table[interpU[u,i,spos[[i]],x],{i,Length@spos}];
interpM[s_,positiveS_,u_,x_]:=interpS[s,#,u,x]&/@positiveS//Total;
And examples
(*example data*)
f[x_,y_,z_]:=Cos[x+y+z] Exp[-0.03(z+y)];
nPoints=6;
m=Table[f[x,y,z],{x,Range[nPoints]},{y,Range[nPoints]},{z,Range[nPoints]}]*1.;
(*reconstruction of the original tensor, works with any number of indices*)
{s,u}=hoSVD2[m];
m2 = Fold[dot2[#, #2[[1]], #2[[2]]]&, s, Transpose@{u, Range@Length[u]}];
Max@Abs[m - m2]
(*interpolation example, works with any number of indices*)
x={1.8,4.2,2.2};
f@@x
positiveS=Position[s,x_/;Abs@x>0.00001];
interpM[s,positiveS,u,x]
(*hoSVD on a big tensor*)
nPoints2=100;
mBig=Table[f[x,y,z],{x,Range[nPoints2]},{y,Range[nPoints2]},{z,Range[nPoints2]}]*1.;
{s,u}=hoSVD2[mBig];//Timing(*finishes in a reasonable amount of time*)
Wikipedia has a perhaps more accessible description.
In multilinear algebra, there does not exist a general decomposition method for multi-way arrays (also known as N-arrays, higher-order arrays, or data-tensors) with all the properties of a matrix singular value decomposition (SVD).
A matrix SVD simultaneously computes
(a) a rank-R decomposition and
(b) the orthonormal row/column matrices.
These two properties can be captured separately by two different decompositions for multi-way arrays. Property (a) is extended to higher order by a class of closely related constructions known collectively as CP decomposition (named after the two most popular and general variants, CANDECOMP and PARAFAC). Such decompositions represent a tensor as the sum of the n-fold outter products of rank-1 tensors, where n is the dimension of the tensor indices.
Property (b) is extended to higher order by a class of methods known variably as Tucker3, N-mode SVD, and N-mode principal component analysis (PCA). (This article will use the general term "Tucker decomposition".) These methods compute the othonormal spaces associated with the different axes (or modes) of a tensor. The Tucker decomposition is also used in multilinear subspace learning as multilinear principal component analysis. This terminology was coined by P. Kroonenberg in the 1980s, but it was later called multilinear SVD and HOSVD (higher-order SVD) by L. De Lathauwer. Historically, much of the interest in higher-order SVDs was driven by the need to analyze empirical data, especial in psychometrics and chemometrics. As such, many of the methods have been independently invented several times, often with subtle variations, leading to a confusing literature. Abstract and general mathematical theorems are rare (though see Kruskal1 with regard to the CP decomposition); instead, the methods are often designed for analyzing specific data types. The 2008 review article by Kolda and Bader2 provides a compact summary of the history of these decompositions, and many references for further reading.
Also, see this package for MathTensor. I think it got updated for v8.
"The creation of MathTensor will open a new chapter in the development of all fields of science and engineering that use serious tensor analysis." - Stephan Wolfram
All this said, it's hard to answer such a general question. Try to provide a specific example and some code of what you've tried.