# Computing Tucker via the HOSVD

## Contents

- Higher-order Singular Value Decomposition (HOSVD) and Sequentially-truncased HOSVD (ST-HOSVD)
- Simple example of usage
- Generate a core with different accuracies for different sizes
- Generate data tensor with core described above
- Compute (full) HOSVD
- Compute low-rank HOSVD approximation
- Verbosity - Getting more or less information.
- Specify the ranks
- Specify the mode order
- Generate bigger data tensor with core described above
- ST-HOSVD compared to HOSVD

## Higher-order Singular Value Decomposition (HOSVD) and Sequentially-truncased HOSVD (ST-HOSVD)

The HOSVD computes a Tucker decomposition of a tensor via a simple process. For each mode k, it computes the r_k leading left singular values of the matrix unfolding and stores those as factor matrix U_k. Then it computes a `ttm` of the original tensor and all the factor matrices to yield the core of size r_1 x r_2 x ... x r_d. The core and factor matrices are used to form the `ttensor`. The values of r_k that lead to a good approximation can be computed automatically to yield a specified error tolerance; this is recommended and the default in our code. The ST-HOSVD is an improvement on the HOSVD that does a TTM in *each* mode before moving on to the next mode. This has the advantage of shrinking the tensor at each step and reducing subsequent computations. ST-HOSVD is the default in the `hosvd` code.

- L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31:279-311, 1966, http://dx.doi.org/10.1007/BF02289464
- L. D. Lathauwer, B. D. Moor and J. Vandewalle, A multilinear singular value decomposition, SIAM J. Matrix Analysis and Applications, 21(4):1253-1278, 2000, http://dx.doi.org/10.1137/S0895479896305696
- N. Vannieuwenhoven, R. Vandebril and K. Meerbergen, A New Truncation Strategy for the Higher-Order Singular Value Decomposition, SIAM J. Scientific Computing, 34(2):A1027-A1052, 2012, http://dx.doi.org/10.1137/110836067

## Simple example of usage

% Create random 50 x 40 x 30 tensor with 5 x 4 x 3 core info = create_problem('Type','Tucker','Num_Factors',[5 4 3],'Size',[50 40 30],'Noise',0.01); X = info.Data; % Compute HOSVD with desired relative error = 0.1 T = hosvd(X,0.1); % Check size of core coresize = size(T.core) % Check relative error relerr = norm(X-full(T))/norm(X)

Computing HOSVD... Size of core: 5 x 4 x 3 ||X-T||/||X|| = 0.00995852 <=0.100000 (tol) coresize = 5 4 3 relerr = 0.0100

## Generate a core with different accuracies for different sizes

We will create a core tensor that has is nearly block diagonal. The blocks are expontentially decreasing in norm, with the idea that we can pick off one block at a time as we increate the prescribed accuracy of the HOSVD. To do this, we use `tenrandblk`.

% Block sizes (need not be cubic). Number of rows is the number % of levels and number of columns is the order of the tensor. bsz = [3 2 1; 2 2 2; 2 3 4]; % Squared norm of each block. Must be length L and sum to <= 1 bsn = [.9 .09 .009]'; % Create core tensor with given block structure and norm 1 G = tenrandblk(bsz,bsn,true);

Created block of size 3 x 2 x 1 with norm (0.948683)^2=0.900000 Created block of size 2 x 2 x 2 with norm (0.300000)^2=0.090000 Created block of size 2 x 3 x 4 with norm (0.094868)^2=0.009000 Created tensor of size 7 x 7 x 7 with off-block-diaognal norm (0.031623)^2=0.001000

```
fprintf('Size of G: %s\n', tt_size2str(size(G)));
```

Size of G: 7 x 7 x 7

## Generate data tensor with core described above

We take the core G and embed into into a larger tensor X by using orthogonal transformations. The true rank of this tensor is equal to the size of G.

% Size of X xsz = [20 20 20]; % Create orthogonal matrices U = cell(3,1); for k = 1:3 V = matrandorth(xsz(k)); U{k} = V(:,1:size(G,k)); end % Create X X = full(ttensor(G,U)); % The norm should be unchanged fprintf('||X||=%f\n',norm(X));

||X||=1.000000

## Compute (full) HOSVD

We compute the ST-HOSVD using the `hosvd` method. We specify the tolerance to close to machine precision. Ideally, it finds a core that is the same size as G.

```
fprintf('ST-HOSVD...\n');
T = hosvd(X,2*sqrt(eps));
```

ST-HOSVD... Computing HOSVD... Size of core: 7 x 7 x 7 ||X-T||/||X|| = 8.3359e-15 <=0.000000 (tol)

## Compute low-rank HOSVD approximation

The norm squared of the first two blocks of G is 0.99, so specifying an error of 1e-2 should yield a core of size 4 x 4 x 3. However, the conservative nature of the algorithm means that it may pick something larger. We can compensate by specifying a larger tolerance.

% Using 1e-2 exactly is potentially too conservative... fprintf('Result with tol = sqrt(1e-2):\n'); T = hosvd(X, sqrt(1e-2),'verbosity',1); % But a small multiple (i.e., |ndims(X)|) usually works... fprintf('Result with tol = sqrt(3e-2):\n'); T = hosvd(X, sqrt(3e-2),'verbosity',1);

Result with tol = sqrt(1e-2): Computing HOSVD... Size of core: 5 x 5 x 5 ||X-T||/||X|| = 0.0517313 <=0.100000 (tol) Result with tol = sqrt(3e-2): Computing HOSVD... Size of core: 3 x 3 x 3 ||X-T||/||X|| = 0.100034 <=0.173205 (tol)

Similarly, lhe norm squared of the first block of G is 0.9, so specifying an error of 1e-1 should result in a core of size 3 x 2 x 1.

% Using 1e-1 exactly is potentially too conservative... fprintf('Result with tol = sqrt(1e-1):\n'); T = hosvd(X, sqrt(1e-1),'verbosity',1); % But a small multiple (i.e., |ndims(X)|) usually works... fprintf('Result with tol = sqrt(3e-1):\n'); T = hosvd(X, sqrt(3e-1),'verbosity',1);

Result with tol = sqrt(1e-1): Computing HOSVD... Size of core: 2 x 2 x 2 ||X-T||/||X|| = 0.20788 <=0.316228 (tol) Result with tol = sqrt(3e-1): Computing HOSVD... Size of core: 1 x 1 x 1 ||X-T||/||X|| = 0.316437 <=0.547723 (tol)

## Verbosity - Getting more or less information.

Setting the verbosity to zero suppresses all output. Cranking up the verbosity gives some insight into the decision-making process...

% Example 1 T = hosvd(X, sqrt(3e-1),'verbosity',10);

Computing HOSVD... ||X||^2 = 1 tol = 0.547723 eigenvalue sum threshold = tol^2 ||X||^2 / d = 0.1 Reverse cummulative sum of evals of Gram matrix: 1: 1.0000 <-- Cutoff 2: 0.1000 3: 0.0331 4: 0.0097 5: 0.0034 6: 0.0004 7: 0.0001 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Reverse cummulative sum of evals of Gram matrix: 1: 0.9000 <-- Cutoff 2: 0.0002 3: 0.0001 4: 0.0000 5: 0.0000 6: 0.0000 7: 0.0000 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Reverse cummulative sum of evals of Gram matrix: 1: 0.8999 <-- Cutoff 2: -0.0000 3: -0.0000 4: -0.0000 5: -0.0000 6: -0.0000 7: -0.0000 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Size of core: 1 x 1 x 1 ||X-T||/||X|| = 0.316437 <=0.547723 (tol)

Example 2

```
T = hosvd(X, sqrt(3*eps),'verbosity',10);
```

Computing HOSVD... ||X||^2 = 1 tol = 2.58096e-08 eigenvalue sum threshold = tol^2 ||X||^2 / d = 2.22045e-16 Reverse cummulative sum of evals of Gram matrix: 1: 1.0000 2: 0.1000 3: 0.0331 4: 0.0097 5: 0.0034 6: 0.0004 7: 0.0001 <-- Cutoff 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Reverse cummulative sum of evals of Gram matrix: 1: 1.0000 2: 0.1000 3: 0.0336 4: 0.0097 5: 0.0044 6: 0.0013 7: 0.0003 <-- Cutoff 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Reverse cummulative sum of evals of Gram matrix: 1: 1.0000 2: 0.0998 3: 0.0331 4: 0.0095 5: 0.0049 6: 0.0017 7: 0.0001 <-- Cutoff 8: -0.0000 9: -0.0000 10: -0.0000 11: -0.0000 12: -0.0000 13: -0.0000 14: -0.0000 15: -0.0000 16: -0.0000 17: -0.0000 18: -0.0000 19: -0.0000 20: -0.0000 Size of core: 7 x 7 x 7 ||X-T||/||X|| = 8.3359e-15 <=0.000000 (tol)

## Specify the ranks

If you know the rank you want, you can specify it. But there's no guarantee that it will satisfy the specified tolerance. In such cases, the method will throw a warning.

% Rank is okay T = hosvd(X,sqrt(3e-1),'ranks',bsz(1,:)); % Rank is too small for the specified error T = hosvd(X,sqrt(3e-1),'ranks',[1 1 1]); % But you can set the error to the tensor norm to make the warning go away T = hosvd(X,norm(X),'ranks',[1 1 1]);

Computing HOSVD... Size of core: 3 x 2 x 1 ||X-T||/||X|| = 0.316418 <=0.547723 (tol) Computing HOSVD... Size of core: 1 x 1 x 1 ||X-T||/||X|| = 0.316437 <=0.547723 (tol) Computing HOSVD... Size of core: 1 x 1 x 1 ||X-T||/||X|| = 0.316437 <=1.000000 (tol)

## Specify the mode order

It's also possible to specify the order of the modes. The default is 1:ndims(X).

```
T = hosvd(X,sqrt(3e-1),'dimorder',ndims(X):-1:1);
```

Computing HOSVD... Size of core: 1 x 1 x 1 ||X-T||/||X|| = 0.316437 <=0.547723 (tol)

## Generate bigger data tensor with core described above

Uses the same procedure as before, but now the size is bigger.

% Size of Y ysz = [100 100 100]; % Create orthogonal matrices U = cell(3,1); for k = 1:3 V = matrandorth(ysz(k)); U{k} = V(:,1:size(G,k)); end % Create Y Y = full(ttensor(G,U));

## ST-HOSVD compared to HOSVD

The answers are essentially the same for the sequentially-truncated HOSVD and the HOSVD...

fprintf('ST-HOSVD...\n'); T = hosvd(Y,2*sqrt(eps)); fprintf('HOSVD...\n'); T = hosvd(Y,2*sqrt(eps),'sequential',false);

ST-HOSVD... Computing HOSVD... Size of core: 7 x 7 x 7 ||X-T||/||X|| = 5.80269e-15 <=0.000000 (tol) HOSVD... Computing HOSVD... Size of core: 7 x 7 x 7 ||X-T||/||X|| = 4.24284e-15 <=0.000000 (tol)

But ST-HOSVD may be slightly faster than HOSVD for larger tensors.

fprintf('Time for 10 runs of ST-HOSVD:\n'); tic, for i =1:10, T = hosvd(Y,2*sqrt(eps),'verbosity',0); end; toc fprintf('Time for 10 runs of HOSVD:\n'); tic, for i =1:10, T = hosvd(Y,2*sqrt(eps),'verbosity',0,'sequential',false); end; toc

Time for 10 runs of ST-HOSVD: Elapsed time is 0.359940 seconds. Time for 10 runs of HOSVD: Elapsed time is 0.547478 seconds.