# Orthogonalization schemes in tensor train format

Martina Iannacito, Olivier Coulaud, Luc Giraud

► To cite this version:

Martina Iannacito, Olivier Coulaud, Luc Giraud. Orthogonalization schemes in tensor train format. SIAM Conference on Applied Linear Algebra (LA24), May 2024, Paris, France. hal-04587668

HAL Id: hal-04587668

https://inria.hal.science/hal-04587668

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inria

# Orthogonalization schemes in tensor train format

SIAM Conference on Applied Linear Algebra, May 14, 2024, Paris

Martina Iannacito*

joint work with Olivier Coulaud† and Luc Giraud†

* Università di Bologna

† Centre Inria at the University of Bordeaux

<figure><img src="image_1.png" /></figure>

Finanziato dall'Unione europea
NextGenerationEU

<figure><img src="image_2.png" /></figure>

Ministero
dell'Università
e della Ricerca

<figure><img src="image_3.png" /></figure>

Italiadomani
DIPARTIMENTO DI SCIENZE E TECNOLOGIE

<figure><img src="image_4.png" /></figure>

ALMA MATER STUDIUM
UNIVERSITÀ DI Bologna
Inria

# Orthogonalization schemes in tensor train format

SIAM Conference on Applied Linear Algebra,
May 14, 2024, Paris

Martina Iannacito*

joint work with Olivier Coulaud† and Luc Giraud†

* Università di Bologna

† Centre Inria at the University of Bordeaux

<figure><img src="image_2.png" /></figure>

Finanziato
dall'Unione europea
NextGenerationEU

<figure><img src="image_3.png" /></figure>

Ministero
dell'Università
e della Ricerca

<figure><img src="image_4.png" /></figure>

Italiadomani
DIPARTIMENTO DI SCIENZE E TECNOLOGIE

<figure><img src="image_5.png" /></figure>

ALMA MATER STUDIUM
UNIVERSITÀ DI Bologna
Computing an orthogonal basis is a fundamental sub-problem for numerical simulations
Computing an orthogonal basis is a fundamental sub-problem for numerical simulations
Computing an orthogonal basis is a fundamental sub-problem for numerical simulations

**High-dimensional** simulations arise in

* Stochastic equations

* Uncertainty quantification problems

* Quantum and vibration chemistry

* Optimization

* Machine learning

<figure><img src="image_5.png" /></figure>
Computing an orthogonal basis is a fundamental sub-problem for numerical simulations

**High-dimensional** simulations arise in

* Stochastic equations

* Uncertainty quantification problems

* Quantum and vibration chemistry

* Optimization

* Machine learning

<figure><img src="image_6.png" /></figure>
Example

$$
\left\{
\begin{array}{l}
\mathcal{L}(u) = f \quad \text{in} \quad \Omega \\
u = f_0 \quad \text{in} \quad \partial\Omega
\end{array}
\right.
\quad \text{for} \quad \Omega \subseteq \mathbb{R}^{n_1 \times \cdots \times n_d}
$$

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

$$
\begin{bmatrix}
1 & 8 & 4 \\
9 & 52 & 3 \\
7 & 26 & 4 \\
6 & 8 & 7
\end{bmatrix}
$$

$$
\mathcal{A}\mathcal{X} = \mathcal{B}
$$

where $\mathcal{A}$ is a multilinear operator and $\mathcal{B}$ is a tensor.
Example

$$
\left\{
\begin{array}{l}
\mathcal{L}(u) = f \quad \text{in } \Omega \\
u = f_0 \quad \text{in } \partial\Omega
\end{array}
\right.
\quad \text{for} \quad \Omega \subseteq \mathbb{R}^{n_1 \times \cdots \times n_d}
$$

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

$$
\begin{bmatrix}
1 & 8 & 4 \\
9 & 52 & 3 \\
7 & 26 & 4 \\
6 & 8 & 7
\end{bmatrix}
$$

$$
\mathcal{A}\mathcal{X} = \mathcal{B}
$$

where $\mathcal{A}$ is a multilinear operator and $\mathcal{B}$ is a tensor.
The storage cost increases **exponentially** with respect to the number of variables.

<figure><img src="image_8.png" /></figure>

Several tensor decomposition techniques exist

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

<figure><img src="image_9.png" /><figcaption>Figure: from [Bi et al. 2022]</figcaption></figure>
The storage cost increases **exponentially** with respect to the number of variables.

<figure><img src="image_9.png" /></figure>

Several tensor decomposition techniques exist

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

<figure><img src="image_10.png" /><figcaption>Figure: from [Bi et al. 2022]</figcaption></figure>
The storage cost increases **exponentially** with respect to the number of variables.

<figure><img src="image_10.png" /></figure>

Several tensor decomposition techniques exist

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

<figure><img src="image_11.png" /><figcaption>Figure: from [Bi et al. 2022]</figcaption></figure>
The storage cost increases **exponentially** with respect to the number of variables.

<figure><img src="image_11.png" /></figure>

Several tensor decomposition techniques exist

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

<figure><img src="image_12.png" /><figcaption>Figure: from [Bi et al. 2022]</figcaption></figure>
<figure><img src="image_12.png" /></figure>

scalar

<figure><img src="image_13.png" /></figure>

vector

<figure><img src="image_14.png" /></figure>

matrix

<figure><img src="image_15.png" /></figure>

order-3 tensor
<figure><img src="image_13.png" /></figure>

scalar

<figure><img src="image_14.png" /></figure>

vector

<figure><img src="image_15.png" /></figure>

matrix

<figure><img src="image_16.png" /></figure>

order-3 tensor
<figure><img src="image_16.png" /></figure>

scalar

<figure><img src="image_17.png" /></figure>

vector

<figure><img src="image_18.png" /></figure>

matrix

<figure><img src="image_19.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is
<figure><img src="image_17.png" /></figure>

scalar

<figure><img src="image_18.png" /></figure>

vector

<figure><img src="image_19.png" /></figure>

matrix

<figure><img src="image_20.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is
<figure><img src="image_20.png" /></figure>

scalar

<figure><img src="image_21.png" /></figure>

vector

<figure><img src="image_22.png" /></figure>

matrix

<figure><img src="image_23.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is

<figure><img src="image_24.png" /></figure>

+ The storage cost is $O(dnr^2)$ with $r = \max\{r_i\}$, said **TT-ranks**
<figure><img src="image_21.png" /></figure>

scalar

<figure><img src="image_22.png" /></figure>

vector

<figure><img src="image_23.png" /></figure>

matrix

<figure><img src="image_24.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is

<figure><img src="image_25.png" /></figure>

+ The storage cost is $O(dnr^2)$ with $r = \max\{r_i\}$, said **TT-ranks**
<figure><img src="image_25.png" /></figure>

scalar

<figure><img src="image_26.png" /></figure>

vector

$r_i$

<figure><img src="image_27.png" /></figure>

matrix

$r_i$

<figure><img src="image_28.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is

<figure><img src="image_29.png" /></figure>

+ The storage cost is $O(dnr^2)$ with $r = \max\{r_i\}$, said **TT-ranks**

- Linear combinations increase the TT-ranks
<figure><img src="image_26.png" /></figure>

scalar

<figure><img src="image_27.png" /></figure>

vector

$r_i$

<figure><img src="image_28.png" /></figure>

matrix

$r_i$

<figure><img src="image_29.png" /></figure>

order-3 tensor

The TT-representation of $\mathcal{X}$ an $(n_1 \times \cdots \times n_d)$ tensor is

<figure><img src="image_30.png" /></figure>

+ The storage cost is $O(dnr^2)$ with $r = \max\{r_i\}$, said **TT-ranks**

- Linear combinations increase the TT-ranks
TT-rounding [Oseledets 2011]

If $\mathcal{X}$ is a tensor in TT-format and $\delta \in (0, 1)$, then $\bar{\mathcal{X}} = \text{TT-rounding}(\mathcal{X}, \delta)$ such that

$$ \|\mathcal{X} - \bar{\mathcal{X}}\| \le \delta \|\mathcal{X}\| $$

where $\max(\text{TT-ranks}(\bar{\mathcal{X}})) < (\max \text{TT-ranks}(\mathcal{X}))$
TT-rounding [Oseledets 2011]

If $\mathcal{X}$ is a tensor in TT-format and $\delta \in (0, 1)$, then $\bar{\mathcal{X}} = \text{TT-rounding}(\mathcal{X}, \delta)$ such that

$$ \|\mathcal{X} - \bar{\mathcal{X}}\| \leq \delta \|\mathcal{X}\| $$

where $\max(\text{TT-ranks}(\bar{\mathcal{X}})) < (\max \text{TT-ranks}(\mathcal{X}))$
Inexactness and orthogonality

TT-rounding [Oseledets 2011]

If $\mathcal{X}$ is a tensor in TT-format and $\delta \in (0, 1)$, then $\bar{\mathcal{X}} = \text{TT-rounding}(\mathcal{X}, \delta)$ such that

$$ \|\mathcal{X} - \bar{\mathcal{X}}\| \le \delta \|\mathcal{X}\| $$

where $\max(\text{TT-ranks}(\bar{\mathcal{X}})) < (\max \text{TT-ranks}(\mathcal{X}))$

<figure><img src="image_30.png" /></figure>

Question

How the introduction of compression steps affects the **loss of orthogonality**

$$ ||\mathbb{I}_k - \mathcal{Q}_k^\top \mathcal{Q}_k|| $$

of computed basis $\mathcal{Q}_k$?
Inexactness and orthogonality

TT-rounding [Oseledets 2011]

If $\mathcal{X}$ is a tensor in TT-format and $\delta \in (0, 1)$, then $\bar{\mathcal{X}} = \text{TT-rounding}(\mathcal{X}, \delta)$ such that

$$ \|\mathcal{X} - \bar{\mathcal{X}}\| \le \delta \|\mathcal{X}\| $$

where $\max(\text{TT-ranks}(\bar{\mathcal{X}})) < (\max \text{TT-ranks}(\mathcal{X}))$

<figure><img src="image_31.png" /></figure>

Question

How the introduction of compression steps affects the **loss of orthogonality**

$$ ||\mathbb{I}_k - \mathcal{Q}_k^\top \mathcal{Q}_k|| $$

of computed basis $\mathcal{Q}_k$?
The **loss of orthogonality** of the computed orthogonal basis $\mathbf{Q}_k$

* depends on the computational precision $u$

* may be affected by the collinearity of the input vectors $\mathbf{A}_k = [\mathbf{a}_1, \dots, \mathbf{a}_k]$, expressed by $\kappa(\mathbf{A}_k)$
The **loss of orthogonality** of the computed orthogonal basis $Q_k$

* depends on the computational precision $u$

* may be affected by the collinearity of the input vectors $A_k = [a_1, \dots, a_k]$, expressed by $\kappa(A_k)$
The loss of orthogonality of the computed orthogonal basis $Q_k$

* depends on the computational precision $u$

* may be affected by the collinearity of the input vectors $A_k = [a_1, \dots, a_k]$, expressed by $\kappa(A_k)$

<table><thead><tr><th>Source</th><th>Algorithm</th><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th></tr></thead><tbody><tr><td>[Stathopoulos and Wu 2002]</td><td>Gram</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>CGS</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Björck 1967]</td><td>MGS</td><td>O(uκ(A<sub>k</sub>))</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>CGS2</td><td>O(u)</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>MGS2</td><td>O(u)</td></tr><tr><td>[Wilkinson 1965]</td><td>Householder</td><td>O(u)</td></tr></tbody></table>
The loss of orthogonality of the computed orthogonal basis $\mathbf{Q}_k$

* depends on the computational precision $u$

* may be affected by the collinearity of the input vectors $\mathbf{A}_k = [\mathbf{a}_1, \dots, \mathbf{a}_k]$, expressed by $\kappa(\mathbf{A}_k)$

<table><thead><tr><th>Source</th><th>Algorithm</th><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th></tr></thead><tbody><tr><td>[Stathopoulos and Wu 2002]</td><td>Gram</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>CGS</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Björck 1967]</td><td>MGS</td><td>O(uκ(A<sub>k</sub>))</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>CGS2</td><td>O(u)</td></tr><tr><td>[Giraud, Langou, and Rozložník 2005]</td><td>MGS2</td><td>O(u)</td></tr><tr><td>[Wilkinson 1965]</td><td>Householder</td><td>O(u)</td></tr></tbody></table>
1

# TT-orthogonalization kernels
1

# TT-orthogonalization kernels
$$Q, R = \text{CGS}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{\mathcal{A}_1, \dots, \mathcal{A}_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P = \mathcal{A}_i$

3 for $j = 1, \dots, i - 1$ do

4 $R(i, j) = \langle \mathcal{A}_i, Q_j \rangle$

5 $P = P - R(i, j)Q_j$

6 end

7 $P = \text{TT-rounding}(P, \delta)$

8 $R(i, i) = \|P\|$

9 $Q_i = 1/R(i, i) P$

10 end

Output: $Q = \{Q_1, \dots, Q_m\}, R$

$$Q, R = \text{MGS}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{\mathcal{A}_1, \dots, \mathcal{A}_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P = \mathcal{A}_i$

3 for $j = 1, \dots, i - 1$ do

4 $R(i, j) = \langle P, Q_j \rangle$

5 $P = P - R(i, j)Q_j$

6 end

7 $P = \text{TT-rounding}(P, \delta)$

8 $R(i, i) = \|P\|$

9 $Q_i = 1/R(i, i) P$

10 end

Output: $Q = \{Q_1, \dots, Q_m\}, R$

They readily write in TT-format.
$Q, R = \text{CGS}(\mathcal{A}, \delta)$

Input: $\mathcal{A} = \{\mathcal{A}_1, \dots, \mathcal{A}_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P = \mathcal{A}_i$

3 for $j = 1, \dots, i - 1$ do

4 $\mathbf{R}(i, j) = \langle \mathcal{A}_i, \mathcal{Q}_j \rangle$

5 $P = P - \mathbf{R}(i, j)\mathcal{Q}_j$

6 end

7 $P = \text{TT-rounding}(P, \delta)$

8 $\mathbf{R}(i, i) = \|P\|$

9 $\mathcal{Q}_j = 1/\mathbf{R}(i, i) P$

10 end

Output: $\mathcal{Q} = \{\mathcal{Q}_1, \dots, \mathcal{Q}_m\}, R$

$\mathcal{Q}, R = \text{MGS}(\mathcal{A}, \delta)$

Input: $\mathcal{A} = \{\mathcal{A}_1, \dots, \mathcal{A}_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P = \mathcal{A}_i$

3 for $j = 1, \dots, i - 1$ do

4 $\mathbf{R}(i, j) = \langle P, \mathcal{Q}_j \rangle$

5 $P = P - \mathbf{R}(i, j)\mathcal{Q}_j$

6 end

7 $P = \text{TT-rounding}(P, \delta)$

8 $\mathbf{R}(i, i) = \|P\|$

9 $\mathcal{Q}_i = 1/\mathbf{R}(i, i) P$

10 end

Output: $\mathcal{Q} = \{\mathcal{Q}_1, \dots, \mathcal{Q}_m\}, R$

They readily write in TT-format.
$$ \mathcal{Q}, \mathbf{R} = \text{CGS2}(\mathcal{A}, \delta) $$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $\mathcal{P}_k = A_i$

3 for $k = 1, 2$ do

4 $\mathcal{P}_k = \mathcal{P}_{k-1}$

5 for $j = 1, \dots, i - 1$ do

6 $\mathbf{R}(i, j) = \langle \mathcal{P}_{k-1}, \mathcal{Q}_j \rangle$

7 $\mathcal{P}_k = \mathcal{P}_k - \mathbf{R}(i, j)\mathcal{Q}_j$

8 end

9 $\mathcal{P}_k = \text{TT-rounding}(\mathcal{P}_k, \delta)$

10 end

11 $\mathbf{R}(i, i) = \|\mathcal{P}_2\|$

12 $\mathcal{Q}_i = 1/\mathbf{R}(i, i) \mathcal{P}_2$

13 end

Output: $\mathcal{Q} = \{\mathcal{Q}_1, \dots, \mathcal{Q}_m\}, \mathbf{R}$

$$ \mathcal{Q}, \mathbf{R} = \text{MGS2}(\mathcal{A}, \delta) $$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $\mathcal{P}_k = A_i$

3 for $k = 1, 2$ do

4 $\mathcal{P}_k = \mathcal{P}_{k-1}$

5 for $j = 1, \dots, i - 1$ do

6 $\mathbf{R}(i, j) = \langle \mathcal{P}_k, \mathcal{Q}_j \rangle$

7 $\mathcal{P}_k = \mathcal{P}_k - \mathbf{R}(i, j)\mathcal{Q}_j$

8 end

9 $\mathcal{P}_k = \text{TT-rounding}(\mathcal{P}_k, \delta)$

10 end

11 $\mathbf{R}(i, i) = \|\mathcal{P}_2\|$

12 $\mathcal{Q}_i = 1/\mathbf{R}(i, i) \mathcal{P}_2$

13 end

Output: $\mathcal{Q} = \{\mathcal{Q}_1, \dots, \mathcal{Q}_m\}, \mathbf{R}$
$$Q, R = \text{CGS2}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P_k = A_i$

3 for $k = 1, 2$ do

4 $P_k = P_{k-1}$

5 for $j = 1, \dots, i - 1$ do

6 $\mathbf{R}(i, j) = \langle P_{k-1}, Q_j \rangle$

7 $P_k = P_k - \mathbf{R}(i, j)Q_j$

8 end

9 $P_k = \text{TT-rounding}(P_k, \delta)$

10 end

11 $\mathbf{R}(i, i) = \|\mathcal{P}_2\|$

12 $Q_i = 1/\mathbf{R}(i, i) P_2$

13 end

Output: $\mathcal{Q} = \{Q_1, \dots, Q_m\}, R$

$$\mathcal{Q}, R = \text{MGS2}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1 for $i = 1, \dots, m$ do

2 $P_k = A_i$

3 for $k = 1, 2$ do

4 $P_k = P_{k-1}$

5 for $j = 1, \dots, i - 1$ do

6 $\mathbf{R}(i, j) = \langle P_k, Q_j \rangle$

7 $P_k = P_k - \mathbf{R}(i, j)Q_j$

8 end

9 $P_k = \text{TT-rounding}(P_k, \delta)$

10 end

11 $\mathbf{R}(i, i) = \|\mathcal{P}_2\|$

12 $Q_i = 1/\mathbf{R}(i, i) P_2$

13 end

Output: $\mathcal{Q} = \{Q_1, \dots, Q_m\}, R$
Let $\mathbf{A} = [\mathbf{a}_1, \dots, \mathbf{a}_m]$, then we look for $\mathbf{A} = \mathbf{QR}$ with $\mathbf{Q}^\top\mathbf{Q} = \mathbb{I}_m$ compute the Gram matrix

$$ \mathbf{A}^\top\mathbf{A} = (\mathbf{R}^\top\mathbf{Q}^\top)\mathbf{QR} = \mathbf{R}^\top\mathbf{R} $$

this is (almost) the **Cholesky** factorization of $\mathbf{A}^\top\mathbf{A}$ that can be written as

$$ \mathbf{A}^\top\mathbf{A} = \mathbf{R}^\top\mathbf{R} = \mathbf{LL}^\top $$

with the Cholesky factor $\mathbf{L} = \mathbf{R}^\top$. Thus, it follows

$$ \mathbf{A} = \mathbf{QR} = \mathbf{QL}^\top $$

from which we obtain $\mathbf{Q}$ by solving a linear system.
Let $\mathbf{A} = [\mathbf{a}_1, \dots, \mathbf{a}_m]$, then we look for $\mathbf{A} = \mathbf{QR}$ with $\mathbf{Q}^\top\mathbf{Q} = \mathbb{I}_m$ compute the Gram matrix

$$ \mathbf{A}^\top\mathbf{A} = (\mathbf{R}^\top\mathbf{Q}^\top)\mathbf{QR} = \mathbf{R}^\top\mathbf{R} $$

this is (almost) the **Cholesky** factorization of $\mathbf{A}^\top\mathbf{A}$ that can be written as

$$ \mathbf{A}^\top\mathbf{A} = \mathbf{R}^\top\mathbf{R} = \mathbf{LL}^\top $$

with the Cholesky factor $\mathbf{L} = \mathbf{R}^\top$. Thus, it follows

$$ \mathbf{A} = \mathbf{QR} = \mathbf{QL}^\top $$

from which we obtain $\mathbf{Q}$ by solving a linear system.
$$Q, R = \text{Gram}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1. $G$ is $(m \times m)$ Gram matrix from $\mathcal{A}$

2. $L$ = cholesky($G$)

3. $\{P_1, \dots, P_m\}$ from $\mathcal{A}$ and $(L^\top)^{-1}$

4. for $i = 1, \dots, m$ do

5. $Q_i = \text{TT-rounding}(P_i, \delta)$

6. end

Output: $Q = \{Q_1, \dots, Q_m\}, R$

In TT-format the following modifications occur

* $G(i, j) = \langle A_i, A_j \rangle$

* $(L^\top)^{-1}$ is explicitly computed

* $P_i$ is constructed as a linear combination of $\mathcal{A}$ elements
$$Q, R = \text{Gram}(\mathcal{A}, \delta)$$

Input: $\mathcal{A} = \{A_1, \dots, A_m\}, \delta \in \mathbb{R}_+$

1. $G$ is $(m \times m)$ Gram matrix from $\mathcal{A}$

2. $L$ = cholesky($G$)

3. $\{P_1, \dots, P_m\}$ from $\mathcal{A}$ and $(L^\top)^{-1}$

4. for $i = 1, \dots, m$ do

5. $Q_i = \text{TT-rounding}(P_i, \delta)$

6. end

Output: $Q = \{Q_1, \dots, Q_m\}, R$

In TT-format the following modifications occur

* $G(i, j) = \langle A_i, A_j \rangle$

* $(L^\top)^{-1}$ is explicitly computed

* $P_i$ is constructed as a linear combination of $\mathcal{A}$ elements
# Householder transformation

Let $x$ and $y \in \mathbb{R}^n$ be vectors with $||y|| = 1$, the Householder reflector $H$ reflects $x$ along $y$, i.e., $Hx = ||x||y$. The Householder reflector $H$ is

$$H = I_n - \frac{2}{||z||^2}zz^\top \quad \text{with} \quad z = (x - ||x||y).$$

## Remark

* The Householder algorithm uses some of the components of the input vectors

* The Householder algorithm does **not** readily apply to tensor in TT-formats

* The approximation is applied to

* the tensor defining the Householder reflector

* the reflected tensors

* the tensor of the orthogonal basis
# Householder transformation

Let $x$ and $y \in \mathbb{R}^n$ be vectors with $||y|| = 1$, the Householder reflector $H$ reflects $x$ along $y$, i.e., $Hx = ||x||y$. The Householder reflector $H$ is

$$H = I_n - \frac{2}{||z||^2}zz^\top \quad \text{with} \quad z = (x - ||x||y).$$

## Remark

* The Householder algorithm uses some of the components of the input vectors

* The Householder algorithm does **not** readily apply to tensor in TT-formats

* The approximation is applied to

* the tensor defining the Householder reflector

* the reflected tensors

* the tensor of the orthogonal basis
2

Numerical experiments
2

Numerical experiments
We construct a set of 'Krylov TT-tensors', {$A_k$}, such that

* $A_1$ is randomly generated with $\max(\text{TT-rank}) = 1$

* $\chi_k = 1/||A_k|| A_k$

* $A_{k+1} = \text{TT-rounding}(\Delta_d \chi_k, \max\_\text{rank} = 1)$ where $\Delta_d$ is the TT-Laplacian

The perturbations are due to

* the inexactness from the approximation at accuracy $\delta = 10^{-5}$

* the rounding error from the floating point system at precision $u = O(10^{-16})$
We construct a set of 'Krylov TT-tensors', {$A_k$}, such that

* $A_1$ is randomly generated with $\max(\text{TT-rank}) = 1$

* $\chi_k = 1/||A_k|| A_k$

* $A_{k+1} = \text{TT-rounding}(\Delta_d \chi_k, \max\_\text{rank} = 1)$ where $\Delta_d$ is the TT-Laplacian

The perturbations are due to

* the inexactness from the approximation at accuracy $\delta = 10^{-5}$

* the rounding error from the floating point system at precision $u = O(10^{-16})$
Loss of orthogonality for {$A_k$}

<figure><img src="image_31.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

• Gram approach

• CGS

• MGS

• CGS2

• MGS2

• Householder transformation
Loss of orthogonality for {$A_k$}

<figure><img src="image_32.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

• Gram approach

• CGS

• MGS

• CGS2

• MGS2

• Householder transformation
Loss of orthogonality for {$A_k$}

<figure><img src="image_32.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

* Gram approach

* CGS

* $\delta\kappa^2(A_k)$ with

$$A_k = [\text{vec}(A_1), \dots, \text{vec}(A_k)]$$
Loss of orthogonality for {$A_k$}

<figure><img src="image_33.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

* Gram approach

* CGS

* $\delta\kappa^2(A_k)$ with

$$A_k = [\text{vec}(A_1), \dots, \text{vec}(A_k)]$$
Loss of orthogonality for {$A_k$}

<figure><img src="image_33.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

• MGS

• $\delta\kappa(A_k)$ with

$$A_k = [\mathrm{vec}(A_1), \dots, \mathrm{vec}(A_k)]$$
Loss of orthogonality for {$A_k$}

<figure><img src="image_34.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</figcaption></figure>

• MGS

• $\delta\kappa(A_k)$ with

$$A_k = [\mathrm{vec}(A_1), \dots, \mathrm{vec}(A_k)]$$
Loss of orthogonality for {$A_k$}

<figure><img src="image_34.png" /></figure>

• CGS2

• MGS2

• Householder transformation

<caption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</caption>
Loss of orthogonality for {$A_k$}

<figure><img src="image_35.png" /></figure>

• CGS2

• MGS2

• Householder transformation

<caption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$</caption>
<table><thead><tr><th rowspan="2">Algorithm</th><th>Matrix, theoretical</th><th>TT-format, conjecture</th></tr><tr><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th></tr></thead><tbody><tr><td>Gram</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td><td>O(δκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>CGS</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td><td>O(δκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>MGS</td><td>O(uκ(A<sub>k</sub>))</td><td>O(δκ(A<sub>k</sub>))</td></tr><tr><td>CGS2</td><td>O(u)</td><td>O(δ)</td></tr><tr><td>MGS2</td><td>O(u)</td><td>O(δ)</td></tr><tr><td>Householder</td><td>O(u)</td><td>O(δ)</td></tr></tbody></table>

with *u* the computational precision, *δ* the compression accuracy
<table><thead><tr><th rowspan="2">Algorithm</th><th>Matrix, theoretical</th><th>TT-format, conjecture</th></tr><tr><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th><th>||I<sub>k</sub> - Q<sub>k</sub><sup>T</sup>Q<sub>k</sub>||</th></tr></thead><tbody><tr><td>Gram</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td><td>O(δκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>CGS</td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td><td>O(δκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>MGS</td><td>O(uκ(A<sub>k</sub>))</td><td>O(δκ(A<sub>k</sub>))</td></tr><tr><td>CGS2</td><td>O(u)</td><td>O(δ)</td></tr><tr><td>MGS2</td><td>O(u)</td><td>O(δ)</td></tr><tr><td>Householder</td><td>O(u)</td><td>O(δ)</td></tr></tbody></table>

with *u* the computational precision, *δ* the compression accuracy
<table><thead><tr><th></th><th>cost in fp operations</th><th>cost in TT-rounding</th></tr></thead><tbody><tr><td>Gram</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>CGS</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>MGS</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>CGS2</td><td>O(4nm<sup>2</sup>)</td><td>2m</td></tr><tr><td>MGS2</td><td>O(4nm<sup>2</sup>)</td><td>2m</td></tr><tr><td>Householder</td><td>O(2nm<sup>2</sup> - 2m<sup>3</sup>/3)</td><td>4m</td></tr></tbody></table>

where *m* is the basis size and *n* is the length of each element
<table><thead><tr><th></th><th>cost in fp operations</th><th>cost in TT-rounding</th></tr></thead><tbody><tr><td>Gram</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>CGS</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>MGS</td><td>O(2nm<sup>2</sup>)</td><td>m</td></tr><tr><td>CGS2</td><td>O(4nm<sup>2</sup>)</td><td>2m</td></tr><tr><td>MGS2</td><td>O(4nm<sup>2</sup>)</td><td>2m</td></tr><tr><td>Householder</td><td>O(2nm<sup>2</sup> - 2m<sup>3</sup>/3)</td><td>4m</td></tr></tbody></table>

where *m* is the basis size and *n* is the length of each element
3

Outlook
3

Outlook
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

* the compression accuracy $\delta$ replacing the computational precision $u$

* the compression accuracy $\delta$ can be read as a computational precision with two peculiarities

> it is independent from the machine architecture

> it controls the perturbation norm-wise rather than component-wise

More detailed results can be found at [Coulaud2022-Ortho]
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

* the compression accuracy $\delta$ replacing the computational precision $u$

* the compression accuracy $\delta$ can be read as a computational precision with two peculiarities

> it is independent from the machine architecture

> it controls the perturbation norm-wise rather than component-wise

More detailed results can be found at [Coulaud2022-Ortho]
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

* the compression accuracy $\delta$ replacing the computational precision $u$

* the compression accuracy $\delta$ can be read as a computational precision with two peculiarities

> it is independent from the machine architecture

> it controls the perturbation norm-wise rather than component-wise

More detailed results can be found at [Coulaud2022-Ortho]

What is left out?

* Theoretical proof of the loss of orthogonality bounds in the inexact framework

* Investigating the quality of the tensor subspace spanned by the orthogonal basis
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

* the compression accuracy $\delta$ replacing the computational precision $u$

* the compression accuracy $\delta$ can be read as a computational precision with two peculiarities

> it is independent from the machine architecture

> it controls the perturbation norm-wise rather than component-wise

More detailed results can be found at [Coulaud2022-Ortho]

## What is left out?

* Theoretical proof of the loss of orthogonality bounds in the inexact framework

* Investigating the quality of the tensor subspace spanned by the orthogonal basis
[1] Y. Bi et al. "Chapter 1 - Tensor decompositions: computations, applications, and challenges". In: *Tensors for Data Processing*. Ed. by Y. Liu. Academic Press, 2022, pp. 1-30. ISBN: 978-0-12-824447-0. DOI: https://doi.org/10.1016/B978-0-12-824447-0.00007-8. URL: https://www.sciencedirect.com/science/article/pii/B9780128244470000078.

[2] Á. Björck. "Solving linear least squares problems by Gram-Schmidt orthogonalization". In: *BIT Numerical Mathematics* 7.1 (Mar. 1967), pp. 1-21. DOI: 10.1007/BF01934122.

[3] L. Giraud, J. Langou, and M. Rozložník. "The loss of orthogonality in the Gram-Schmidt orthogonalization process". In: *Computers & Mathematics with Applications* 50.7 (2005). Numerical Methods and Computational Mechanics, pp. 1069-1075. DOI: https://doi.org/10.1016/j.camwa.2005.08.009.

[4] I. V. Oseledets. "Tensor-Train Decomposition". In: *SIAM Journal on Scientific Computing* 33.5 (2011), pp. 2295-2317. DOI: 10.1137/090752286.

[5] A. Stathopoulos and K. Wu. "A Block Orthogonalization Procedure with Constant Synchronization Requirements". In: *SIAM Journal on Scientific Computing* 23.6 (2002), pp. 2165-2182. DOI: 10.1137/S1064827500370883.

[6] J. H. Wilkinson. *The algebraic eigenvalue problem*. en. Numerical Mathematics and Scientific Computation. Oxford, England: Clarendon Press, 1965.
[1] Y. Bi et al. "Chapter 1 - Tensor decompositions: computations, applications, and challenges". In: *Tensors for Data Processing*. Ed. by Y. Liu. Academic Press, 2022, pp. 1-30. ISBN: 978-0-12-824447-0. DOI: https://doi.org/10.1016/B978-0-12-824447-0.00007-8. URL: https://www.sciencedirect.com/science/article/pii/B9780128244470000078.

[2] Á. Björck. "Solving linear least squares problems by Gram-Schmidt orthogonalization". In: *BIT Numerical Mathematics* 7.1 (Mar. 1967), pp. 1-21. DOI: 10.1007/BF01934122.

[3] L. Giraud, J. Langou, and M. Rozložník. "The loss of orthogonality in the Gram-Schmidt orthogonalization process". In: *Computers & Mathematics with Applications* 50.7 (2005). Numerical Methods and Computational Mechanics, pp. 1069-1075. DOI: https://doi.org/10.1016/j.camwa.2005.08.009.

[4] I. V. Oseledets. "Tensor-Train Decomposition". In: *SIAM Journal on Scientific Computing* 33.5 (2011), pp. 2295-2317. DOI: 10.1137/090752286.

[5] A. Stathopoulos and K. Wu. "A Block Orthogonalization Procedure with Constant Synchronization Requirements". In: *SIAM Journal on Scientific Computing* 23.6 (2002), pp. 2165-2182. DOI: 10.1137/S1064827500370883.

[6] J. H. Wilkinson. *The algebraic eigenvalue problem*. en. Numerical Mathematics and Scientific Computation. Oxford, England: Clarendon Press, 1965.