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}]]

m = RandomReal[NormalDistribution[], {3, 4}];
{s, u} = hoSVD[m];
Max@Abs[m - Transpose@u[[1]].s.u[[2]]]

dot[s_, u_, n_] := Transpose[Transpose[s, Ordering[#]].u, #] &@Append[#, n] &@
     Delete[#, n] &@Range@ArrayDepth[s];

m2 = dot[dot[s, u[[1]], 1], u[[2]], 2];
Max@Abs[m - m2]

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]

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]

defaultTol=0.00000001;
rankFunction[mat_,tol_:defaultTol]:=Length@Select[Diagonal@mat,Abs@#>tol&];

hoSVD2[m_,rankF_:rankFunction]:=
Module[{d=ArrayDepth[m],us,Ak,r,u2,u,w,v,rotatePermutation,flateningDimensions},
    Ak=m;
    rotatePermutation=RotateRight@Range[d];
    flateningDimensions={{1},Delete[Range[d],1]};

    us=
        Table[
            {u,w,v}=SingularValueDecomposition@Flatten[Ak,flateningDimensions];
            r=rankF@w;
            u2=u[[All,;;r]];
            Ak=Transpose[(ConjugateTranspose@u2).Ak,rotatePermutation];

            u2
            ,
            {i,d}
        ];

    (*s=Ak;*)
    {Ak,us}
];

(*u(n).s*)
dot2[s_,u_,n_]:=
Module[{perm,permInv},
    permInv=Range@ArrayDepth[s]//Delete[#,n]&//Prepend[#,n]&;
    perm=Ordering[permInv];
    Transpose[u.Transpose[s,perm],permInv]
];

interpolateVector[v_,x_]:=Interpolation[v][x];
interpU[u_,i_,j_,x_]:=interpolateVector[u[[i,All,j]],x[[i]]];
interpS[s_,spos_,u_,x_]:=Part[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;

f[x_,y_,z_]:=Cos[x+y]Sin[z]
nPoints=6;
m=Table[f[x,y,z],{x,Range[nPoints]},{y,Range[nPoints]},{z,Range[nPoints]}]*1.;

{s,u}=hoSVD2[m];
m2=dot2[dot2[dot2[s,u[[1]],1],u[[2]],2],u[[3]],3];
Max@Abs[m-m2]

x={2.,4.1,3.1};
f@@x

positiveS=Position[s,x_/;Abs@x>0.00001];
interpM[s,positiveS,u,x]

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