% hierarchical Tucker with 4 indices, h-rank r
r1234=r; r12=r;r34=r;r1=r;r2=r;r3=r;r4=r;
for i1=0:1;for i2=0:1;for i3=0:1;for i4=0:1;
                A(i1+1,i2+1,i3+1,i4+1)=i1+2*i2+4*i3+8*i4;
            end;end;end;end
A=rand(size(A));
T=tensor(A);
a12=tenmat(T,[1,2]);
a34=tenmat(T,[3,4]);
a1=tenmat(T,[1]); a2=tenmat(T,[2]);
a3=tenmat(T,[3]); a4=tenmat(T,[4]);
a1234=tenmat(T,[ ]);  
[u1,s1,v1]=svd(double(a1)); rs1=size(s1);
r1=min(r1,rs1(1));  u1=u1(:,1:r1);% u1,pause
[u2,s2,v2]=svd(double(a2)); rs2=size(s2);
r2=min(r2,rs2(1)); u2=u2(:,1:r2);
[u3,s3,v3]=svd(double(a3)); rs3=size(s3);
r3=min(r3,rs3(1)); u3=u3(:,1:r3);
[u4,s4,v4]=svd(double(a4)); rs4=size(s4);
r4=min(r4,rs4(1)); u4=u4(:,1:r4);
[u12,s12,v12]=svd(double(a12)); rs12=size(s12);
r12=min(r12,rs12(1)); u12=u12(:,1:r12);
[u34,s34,v34]=svd(double(a34)); rs34=size(s34);
r34=min(r34,rs34(1)); u34=u34(:,1:r34);
% s12=s12(1:r12,1:r2); 
% ss=size(s12); ls=ss(1)*ss(2);
% s1234=reshape(s12,[ls,1])
b12=kron(u2,u1)'*u12;
b34=kron(u4,u3)'*u34;
b1234=kron(u34,u12)'*a1234';
va=kron(kron(u4,u3),kron(u2,u1))*kron(b34,b12)*b1234;
norm(va-a1234')