TT-operator

[rkuoyemnbeek]

Hi,

In the original paper of prof Oseledets ( [https://epubs.siam.org/doi/abs/10.1137/090752286?journalCode=sjoce3] ) a TT-operator is defined as (see pg 2312) " A matrix M is said to be in TT-format if its elements are defined as M(i_1,..., i_d, j_1, ..., j_d) = M1(i_1, j_1) ... Md(i_d,j_d) (1) "

Let d = 3, we take random matrices M1 \in \R^{n1 x n1}, M2 \in \R^{n2 x n2}, M3 \in \R^{n3 x n3} and we define the operator M as in (1).

function M_ttm = M_maker( M1, M2, M3)

d = 3;
r = ones(4,1);
n = zeros(d,1);
n(1) = size(M1, 1)^2; % we assume square matrices
n(2) = size(M2, 1)^2;
n(3) = size(M3, 1)^2;

M_tt = tt_tensor();
M_tt.d = d;
M_tt.n = n;
M_tt.r = r;
ps = cumsum( [1; n.*r(2:end).*r(1:end-1)]);
M_tt.ps = ps;
core_M = zeros(ps(end)-1,1);

core_M(1:ps(2)-1) = M1(:);
core_M(ps(2):ps(3)-1) = M2(:);
core_M(ps(3):ps(4)-1) = M3(:);

M_tt.core = core_M;
M_ttm = tt_matrix( M_tt, [size(M1, 1), size(M2, 1), size(M3, 1)], [size(M1, 1), size(M2, 1), size(M3, 1)]);
end

According to the paper, this should be equal this TT-operator back in matrix format, is equal to M1 \otimes M2 \otimes M3

But according to my scripts, this is equal to M3 \otimes M2 \otimes M1

n1 = 3; n2 = 4; n3 = 5;
M1 = randn(n1); M2 = randn(n2); M3 = randn(n3);
M_ttm = M_maker( M1, M2, M3);

% according to the paper, M_ttm should be equal to 
M_test = kron( M1, kron( M2, M3));
norm( full( M_ttm) - M_test) % not equal to 0
% but it seems not right (I assume that the full(.) operator is correctly
% implemented)

% It seems that it is equal to
M_test1 = kron(M3, kron( M2, M1));
norm( full( M_ttm) - M_test1) % equal to 0

For me, it seems that or the full(.) operator is not correct, or the paper is not correct or I am missing something. Can you clarify how it works?

[dolgov]

Hello, that's correct, if you use the built-in kron, you need to permute the arguments. This is a long story how a linear algebra notation ended up incompatible with efficient coding. Arrays in Matlab and Fortran are listed starting from the first dimension, e.g. A(i,j) = a(i+(j-1)n), where a = A(:) is the raw contiguous storage. In contrast, the (standard) Kronecker product (and hence Matlab kron) lists the second argument first, i.e. c((i-1)m + j) = a(i)*b(j), if c = kron(a, b). Whereas 2 arguments are easy to swap, in high dimensions this becomes a mess. It's more convenient to use the left Kronecker product (https://epubs.siam.org/doi/pdf/10.1137/1031127 or http://www.cs.cornell.edu/cv/ResearchPDF/KPHist.pdf) which is compatible with the native array storage, returned in particular by full. The left kron is implemented in tkron function for tt_tensor and tt_matrix, and it should be easier to use than building up the class properties from scratch.

[rkuoyemnbeek]

Hi, Thank you very much for the quick responses to my questions. I was not aware of the different definition of kronecker products in the papers about TT-tensors. I was for some weeks very confused about it, so thanks for clarifying this.

BTW: I am a Phd-Student at KU Leuven in the numerical linear algebra group. I use Tensor-Trains in my research.

regards, Koen