# Higher order SVD

Does anyone know how to do a higher order SVD in Mathematica ? A good reference seems to be here http://www.sandia.gov/~tgkolda/pubs/pubfiles/SAND2007-6702.pdf but I don't understand their formalism yet, maybe someone has already some experience with this in Mathematica and this could speed up the understanding. Thanks

Another interesting reference is here ftp://ftp.cs.sjtu.edu.cn:990/zhang-lq/Papers/TensorFactorization/2004icpr_tensorRepresentaion_Wang.pdf

-
The paper you reference has great illustrations, that make thinking about tensor decompositions quite intuitive. Thanks for the reference. – Jagra May 22 '12 at 2:27

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

-
Nice, thank you very much for your time, your implementation is quite compact which is good. This can have applications to many fields. – faysou Oct 23 '13 at 22:28

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.

-
This is more of a long comment than an answer I suppose... – J. M. May 22 '12 at 2:16
@J.M. -- Acknowledged, just trying to help the question down the road towards an answer. – Jagra May 22 '12 at 2:29