# 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

• 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

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*)

• You might want to give a second argument à la UpTo[maxrank] to SingularValueDecomposition in order to compute a truncated SVD. For example, in your bigM case, a maximal rank for 10 or so should reduce the time for SVD drastically without loosing anything of the approximation quality. – Henrik Schumacher May 19 '19 at 19:02
• Thank you for your comment, I've updated the function accordingly, even putting UpTo[Infinity] speeds up the computation massively. This is quite an impressive algorithm, in terms of what it can do compared to the time it takes. – faysou May 20 '19 at 7:34
• Indeed, UpTo[Infinity] leads to a massive speed up. Interesting... Thanks for letting be know! – Henrik Schumacher May 20 '19 at 10:09
• The algorithm must be iterative and stop when singular values are too small. – faysou May 20 '19 at 10:32
• – faysou May 21 '19 at 8:14

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.'s ennui 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