# 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. The 29th Biennial Conference on Numerical Analysis 2023, Jun 2023, Glasgow, United Kingdom. hal-04394115

HAL Id: hal-04394115

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

Submitted on 15 Jan 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
# 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. The 29th Biennial Conference on Numerical Analysis 2023, Jun 2023, Glasgow, United Kingdom. hal-04394115

HAL Id: hal-04394115

https://inria.hal.science/hal-04394115v1

Submitted on 15 Jan 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>
<figure><img src="image_2.png" /></figure>

Orthogonalization schemes
in tensor train format

29th Biennial Numerical Analysis Conference,
June 28th 2023, Glasgow

Martina Iannacito*
joint work with Olivier Coulaud† and Luc Giraud†

* KU Leuven

† Centre Inria at the University of Bordeaux
<figure><img src="image_2.png" /></figure>

Orthogonalization schemes
in tensor train format

29th Biennial Numerical Analysis Conference,
June 28th 2023, Glasgow

Martina Iannacito*
joint work with Olivier Coulaud† and Luc Giraud†

* KU Leuven

† Centre Inria at the University of Bordeaux
* Stochastic equations

* Uncertainty quantification problems

* Quantum and vibration chemistry

* Optimization

* Machine learning

reduce their problems to

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

$$ \min_{x \in S} ||b - Ax|| $$

$$ Ax = \lambda x $$

$$ Ax = b $$
* Stochastic equations

* Uncertainty quantification problems

* Quantum and vibration chemistry

* Optimization

* Machine learning

reduce their problems to

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

$$ \text{Least-squares} \\ \min_{x \in S} ||b - Ax|| $$

$$ \text{Eigenpairs} \\ Ax = \lambda x $$

$$ \text{Linear systems} \\ Ax = b $$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$q_1 = a_1,$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$q_1 = a_1,$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1,
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1,
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}_h linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}
$$
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{k+1} \rangle}{\|q_j\|^2} q_j
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}
$$

Are these versions equivalent in finite precision arithmetic?
Gram-Schmidt process [Leon et al. 2013]

Given {$a_h$}$_h$ linearly independent vectors, we construct an
orthogonal basis {$q_h$} such that

$$
q_1 = a_1, \quad q_2 = a_2 - \frac{\langle a_2, q_1 \rangle}{\|q_1\|^2} q_1, \dots \quad q_{k+1} = a_{k+1} - \sum_{j=1}^{k} \frac{\langle a_j, q_{j-1} \rangle}{\|q_{j-1}\|^2} q_{j-1}
$$

How to make an algorithm out of it?

let $Q_k = [q_1, \dots, q_k]$ be a $n \times k$ matrix

All projections, then subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}
$$

Projections alternate subtractions

$$
q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}
$$

Are these versions equivalent in finite precision arithmetic?
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10 \quad \text{in fp64}$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{CGS} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix},$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{CGS} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix},$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix} \quad \text{with} \quad \varepsilon = 1e-10 \quad \text{in fp64}$$

$$Q_{\text{CGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix} \quad \text{with} \quad \varepsilon = 1e-10 \quad \text{in fp64}$$

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{CGS} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{\text{CGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$

$$Q_{\text{MGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon/2 \\ 0 & \varepsilon & -\varepsilon/2 \\ 0 & 0 & \varepsilon \end{bmatrix},$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{\text{CGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$

$$Q_{\text{MGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon/2 \\ 0 & \varepsilon & -\varepsilon/2 \\ 0 & 0 & \varepsilon \end{bmatrix},$$
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{\text{CGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$

$$Q_{\text{MGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon/2 \\ 0 & \varepsilon & -\varepsilon/2 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_1, q_3 \rangle = \varepsilon \sqrt{\frac{2}{3}}$$

This example was taken from [Björck 1996].
Not equivalent in general

Classical Gram-Schmidt

All projections, then subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - Q_k Q_k^\top) a_{k+1}$$

$$A = \begin{bmatrix} 1 & 1 & 1 \\ \varepsilon & 0 & 0 \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}$$

with

$$\varepsilon = 1e-10$$

in fp64

Modified Gram-Schmidt

Projections alternate subtractions

$$q_{k+1} \leftarrow (\mathbb{I}_n - q_k \otimes q_k) \dots (\mathbb{I}_n - q_1 \otimes q_1) a_{k+1}$$

$$Q_{\text{CGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon \\ 0 & \varepsilon & 0 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_2, q_3 \rangle = \frac{1}{2}$$

$$Q_{\text{MGS}} = \begin{bmatrix} 1 & 0 & 0 \\ \varepsilon & -\varepsilon & -\varepsilon/2 \\ 0 & \varepsilon & -\varepsilon/2 \\ 0 & 0 & \varepsilon \end{bmatrix}, \langle q_1, q_3 \rangle = \varepsilon \sqrt{\frac{2}{3}}$$

This example was taken from [Björck 1996].
Let $Q_k = [q_1, \dots, q_k]$ be the computed orthogonal basis, then its
**loss of orthogonality** is

$$||\mathbb{I}_k - Q_k^\top Q_k||.$$

It measures the quality in terms of orthogonality of the computed
basis. It is linked with the computational precision $u$ and the
linearly dependency of the input vectors $A_k = [a_1, \dots, a_k]$,
estimated through $\kappa(A_k)$.

<table><thead><tr><th colspan="3">Matrix</th></tr><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 et al. 2002]</td><td><b>Gram</b></td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Giraud et al. 2005]</td><td><b>CGS</b></td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Björck 1967]</td><td><b>MGS</b></td><td>O(uκ(A<sub>k</sub>))</td></tr><tr><td>[Giraud et al. 2005]</td><td><b>CGS2</b></td><td>O(u)</td></tr><tr><td>[Giraud et al. 2005]</td><td><b>MGS2</b></td><td>O(u)</td></tr><tr><td>[Wilkinson 1965]</td><td><b>Householder</b></td><td>O(u)</td></tr></tbody></table>
Let $Q_k = [q_1, \dots, q_k]$ be the computed orthogonal basis, then its
**loss of orthogonality** is

$$||\mathbb{I}_k - Q_k^\top Q_k||.$$

It measures the quality in terms of orthogonality of the computed
basis. It is linked with the computational precision $u$ and the
linearly dependency of the input vectors $A_k = [a_1, \dots, a_k]$,
estimated through $\kappa(A_k)$.

<table><thead><tr><th colspan="3">Matrix</th></tr><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 et al. 2002]</td><td><strong>Gram</strong></td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Giraud et al. 2005]</td><td><strong>CGS</strong></td><td>O(uκ<sup>2</sup>(A<sub>k</sub>))</td></tr><tr><td>[Björck 1967]</td><td><strong>MGS</strong></td><td>O(uκ(A<sub>k</sub>))</td></tr><tr><td>[Giraud et al. 2005]</td><td><strong>CGS2</strong></td><td>O(u)</td></tr><tr><td>[Giraud et al. 2005]</td><td><strong>MGS2</strong></td><td>O(u)</td></tr><tr><td>[Wilkinson 1965]</td><td><strong>Householder</strong></td><td>O(u)</td></tr></tbody></table>
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 3 \\ 7 & 1 & 2 \\ 6 & 8 & 7 \end{bmatrix} $$
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 3 \\ 7 & 1 & 2 \\ 6 & 8 & 7 \end{bmatrix} $$
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

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

+ Better representation

- Curse of dimensionality
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

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

+ Better representation

- Curse of dimensionality
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

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

+ Better representation

- Curse of dimensionality

Approximation techniques were proposed

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

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

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

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

+ Better representation

- Curse of dimensionality

Approximation techniques were proposed

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

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

[ 1 8 4
  9 2 2
  7 1 6 ]

since the last century

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

tensor

[ 1 8 4
  9 2 3 1
  7 1 2 4 9
  6 8 7 ]

+ Better representation

- Curse of dimensionality

Approximation techniques were proposed

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

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

Approximation techniques introduce **norm-wise** compression errors
matrix

$$ \begin{bmatrix} 1 & 8 & 4 \\ 9 & 2 & 2 \\ 7 & 1 & 6 \end{bmatrix} $$

since the last century

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

tensor

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

+ Better representation

- Curse of dimensionality

Approximation techniques were proposed

* Canonical Polyadic

* Tucker

* Hierarchical Tucker

* Tensor-Train

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

Approximation techniques introduce **norm-wise** compression errors
<figure><img src="image_10.png" /></figure>

scalar

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

vector

$r_i$

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

matrix

$r_i$

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

order-3 tensor

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

scalar

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

vector

$r_i$

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

matrix

$r_i$

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

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

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is
<figure><img src="image_14.png" /></figure>

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is
<figure><img src="image_18.png" /></figure>

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

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

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

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

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

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

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

<figure><img src="image_27.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_23.png" /></figure>

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

<figure><img src="image_27.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_28.png" /></figure>

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

<figure><img src="image_32.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 $\mathbf{z}$ is an order-$d$ tensor in TT-format and $\delta \in (0, 1)$, then
$\bar{\mathbf{z}} = \text{TT-rounding}(\mathbf{z}, \delta)$ such that

$$ \|\mathbf{z} - \bar{\mathbf{z}}\| \le \delta \|\mathbf{z}\| $$
<figure><img src="image_28.png" /></figure>

scalar

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

vector

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

matrix

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

order-3 tensor

Let $\mathbf{x} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ be order-$d$ tensor, its TT-representation is

<figure><img src="image_32.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 $\mathbf{z}$ is an order-$d$ tensor in TT-format and $\delta \in (0, 1)$, then
$\bar{\mathbf{z}} = \text{TT-rounding}(\mathbf{z}, \delta)$ such that

$$ \|\mathbf{z} - \bar{\mathbf{z}}\| \le \delta \|\mathbf{z}\| $$
In the Tensor-Train framework

compression precision $\delta$

* norm-wise perturbation

* TT-model [Oseledets 2011]

computational precision $u$

* component-wise perturbation

* IEEE model [Higham 2002]
In the Tensor-Train framework

compression precision $\delta$

* norm-wise perturbation

* TT-model [Oseledets 2011]

computational precision $u$

* component-wise perturbation

* IEEE model [Higham 2002]
In the Tensor-Train framework

compression precision $\delta$

* norm-wise perturbation

* TT-model [Oseledets 2011]

computational precision $u$

* component-wise perturbation

* IEEE model [Higham 2002]

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

* How to design orthogonalization kernels for tensor subspace?

* Does compression affect the loss of orthogonality?

* Are tensor results related with the known linear algebra ones?
In the Tensor-Train framework

compression precision $\delta$

* norm-wise perturbation

* TT-model [Oseledets 2011]

computational precision $u$

* component-wise perturbation

* IEEE model [Higham 2002]

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

* How to design orthogonalization kernels for tensor subspace?

* Does compression affect the loss of orthogonality?

* Are tensor results related with the known linear algebra ones?
$$Q, R = \text{CGS}(A, \delta)$$

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

$$1 \quad \text{for } i = 1, \dots, m \text{ do}$$

$$2 \quad \mathbf{p} = \mathbf{a}_i$$

$$3 \quad \text{for } j = 1, \dots, i - 1 \text{ do}$$

$$4 \quad R(i,j) = \langle \mathbf{a}_i, \mathbf{q}_j \rangle$$

$$5 \quad \mathbf{p} = \mathbf{p} - R(i,j)\mathbf{q}_j$$

$$6 \quad \text{end}$$

$$7 \quad \mathbf{p} = \text{TT-rounding}(\mathbf{p}, \delta)$$

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

$$9 \quad \mathbf{q}_i = 1/R(i, i) \mathbf{p}$$

$$10 \quad \text{end}$$

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

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

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

$$1 \quad \text{for } i = 1, \dots, m \text{ do}$$

$$2 \quad \mathbf{p} = \mathbf{a}_i$$

$$3 \quad \text{for } j = 1, \dots, i - 1 \text{ do}$$

$$4 \quad R(i,j) = \langle \mathbf{p}, \mathbf{q}_j \rangle$$

$$5 \quad \mathbf{p} = \mathbf{p} - R(i,j)\mathbf{q}_j$$

$$6 \quad \text{end}$$

$$7 \quad \mathbf{p} = \text{TT-rounding}(\mathbf{p}, \delta)$$

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

$$9 \quad \mathbf{q}_i = 1/R(i, i) \mathbf{p}$$

$$10 \quad \text{end}$$

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

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

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

$$1 \quad \text{for } i = 1, \dots, m \text{ do}$$

$$2 \quad \mathbf{p} = \mathbf{a}_i$$

$$3 \quad \text{for } j = 1, \dots, i - 1 \text{ do}$$

$$4 \quad R(i,j) = \langle \mathbf{a}_i, \mathbf{q}_j \rangle$$

$$5 \quad \mathbf{p} = \mathbf{p} - R(i,j)\mathbf{q}_j$$

$$6 \quad \text{end}$$

$$7 \quad \mathbf{p} = \text{TT-rounding}(\mathbf{p}, \delta)$$

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

$$9 \quad \mathbf{q}_i = 1/R(i, i) \mathbf{p}$$

$$10 \quad \text{end}$$

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

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

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

$$1 \quad \text{for } i = 1, \dots, m \text{ do}$$

$$2 \quad \mathbf{p} = \mathbf{a}_i$$

$$3 \quad \text{for } j = 1, \dots, i - 1 \text{ do}$$

$$4 \quad R(i,j) = \langle \mathbf{p}, \mathbf{q}_j \rangle$$

$$5 \quad \mathbf{p} = \mathbf{p} - R(i,j)\mathbf{q}_j$$

$$6 \quad \text{end}$$

$$7 \quad \mathbf{p} = \text{TT-rounding}(\mathbf{p}, \delta)$$

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

$$9 \quad \mathbf{q}_i = 1/R(i, i) \mathbf{p}$$

$$10 \quad \text{end}$$

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

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

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

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

2 $\mathbf{p}_k = \mathbf{a}_i$

3 for $k = 1, 2$ do

4 $\mathbf{p}_k = \mathbf{p}_{k-1}$
5 for $j = 1, \dots, i-1$ do

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

7 $\mathbf{p}_k = \mathbf{p}_k - R(i, j)\mathbf{q}_j$

8 end

9 $\mathbf{p}_k = \text{TT-rounding}(\mathbf{p}_k, \delta)$

10 end

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

12 $\mathbf{q}_i = 1/R(i, i) \mathbf{p}_2$

13 end

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

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

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

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

2 $\mathbf{p}_k = \mathbf{a}_i$

3 for $k = 1, 2$ do

4 $\mathbf{p}_k = \mathbf{p}_{k-1}$
5 for $j = 1, \dots, i-1$ do

6 $R(i, j) = \langle \mathbf{p}_k, \mathbf{q}_j \rangle$

7 $\mathbf{p}_k = \mathbf{p}_k - R(i, j)\mathbf{q}_j$

8 end

9 $\mathbf{p}_k = \text{TT-rounding}(\mathbf{p}_k, \delta)$

10 end

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

12 $\mathbf{q}_i = 1/R(i, i) \mathbf{p}_2$

13 end

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

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

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

2 $\mathbf{p}_k = \mathbf{a}_i$

3 for $k = 1, 2$ do

4 $\mathbf{p}_k = \mathbf{p}_{k-1}$
5 for $j = 1, \dots, i-1$ do

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

7 $\mathbf{p}_k = \mathbf{p}_k - R(i, j)\mathbf{q}_j$

8 end

9 $\mathbf{p}_k = \text{TT-rounding}(\mathbf{p}_k, \delta)$

10 end

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

12 $\mathbf{q}_i = 1/R(i, i) \mathbf{p}_2$

13 end

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

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

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

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

2 $\mathbf{p}_k = \mathbf{a}_i$

3 for $k = 1, 2$ do

4 $\mathbf{p}_k = \mathbf{p}_{k-1}$
5 for $j = 1, \dots, i-1$ do

6 $R(i, j) = \langle \mathbf{p}_k, \mathbf{q}_j \rangle$

7 $\mathbf{p}_k = \mathbf{p}_k - R(i, j)\mathbf{q}_j$

8 end

9 $\mathbf{p}_k = \text{TT-rounding}(\mathbf{p}_k, \delta)$

10 end

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

12 $\mathbf{q}_i = 1/R(i, i) \mathbf{p}_2$

13 end

Output: $Q = \{\mathbf{q}_1, \dots, \mathbf{q}_m\}, R$
Let $A = [a_1, \dots, a_m]$, then we look for $A = QR$ with $Q^T Q = I_m$ compute the Gram matrix

$$A^T A = (R^T Q^T) QR = R^T R$$

this is (almost) the **Cholesky** factorization of $A^T A$ that can be written as

$$A^T A = R^T R = LL^T$$

with the Cholesky factor $L = R^T$. Thus, it follows

$$A = QR = QL^T$$

from which we obtain $Q$ by solving a linear system.
Let $A = [a_1, \dots, a_m]$, then we look for $A = QR$ with $Q^T Q = I_m$ compute the Gram matrix

$$A^T A = (R^T Q^T) QR = R^T R$$

this is (almost) the **Cholesky** factorization of $A^T A$ that can be written as

$$A^T A = R^T R = LL^T$$

with the Cholesky factor $L = R^T$. Thus, it follows

$$A = QR = QL^T$$

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

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

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

2. $L = \text{cholesky}(G)$

3. $\{\mathbf{p}_1, \dots, \mathbf{p}_m\}$ from $\mathcal{A}$ and $(L^\top)^{-1}$

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

5. $q_i = \text{TT-rounding}(\mathbf{p}_i, \delta)$

6. end

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

In TT-format the following modifications occur

* $G(i, j)$ is the scalar product of $\mathbf{a}_i$ and $\mathbf{a}_j$

* The inverse of $L^\top$ is explicitly computed

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

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

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

2. $L = \text{cholesky}(G)$

3. $\{\mathbf{p}_1, \dots, \mathbf{p}_m\}$ from $\mathcal{A}$ and $(L^\top)^{-1}$

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

5. $q_i = \text{TT-rounding}(\mathbf{p}_i, \delta)$

6. end

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

In TT-format the following modifications occur

* $G(i, j)$ is the scalar product of $\mathbf{a}_i$ and $\mathbf{a}_j$

* The inverse of $L^\top$ is explicitly computed

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

Given a vector $x \in \mathbb{R}^n$ and a direction $y \in \mathbb{R}^n$, the Householder reflector $H$ reflects $x$ along $y$, i.e.,

$$Hx = ||x||y \quad \text{with} \quad ||y|| = 1.$$

Thanks to its properties, $H$ writes as

$$H = \mathbb{I}_n - \frac{2}{||z||^2} z \otimes z \quad \text{with} \quad z = (x - ||x||y).$$

The practical implementation of the Householder transformation kernel uses the components of the input vectors.

## Remark

The Householder algorithm does **not** readily apply to tensor in TT-formats, because of the compressed nature of this format. It is described in [Iannacito 2022]
## Householder transformation

Given a vector $x \in \mathbb{R}^n$ and a direction $y \in \mathbb{R}^n$, the Householder reflector $H$ reflects $x$ along $y$, i.e.,

$$Hx = ||x||y \quad \text{with} \quad ||y|| = 1.$$

Thanks to its properties, $H$ writes as

$$H = \mathbb{I}_n - \frac{2}{||z||^2} z \otimes z \quad \text{with} \quad z = (x - ||x||y).$$

The practical implementation of the Householder transformation kernel uses the components of the input vectors.

### Remark

The Householder algorithm does **not** readily apply to tensor in TT-formats, because of the compressed nature of this format. It is described in [Iannacito 2022]
Let $\Delta_d$ is the TT-Laplacian, then $\{a_k\}$ are 'Krylov tensors', i.e.,
$$a_{k+1} = \text{TT-rounding}(\Delta_d x_k, \max\_\text{rank} = 1) \text{ with } x_{k+1} = \frac{1}{\|a_{k+1}\|} a_{k+1}$$
Let $\Delta_d$ is the TT-Laplacian, then $\{a_k\}$ are 'Krylov tensors', i.e.,
$$a_{k+1} = \text{TT-rounding}(\Delta_d x_k, \text{max\_rank} = 1) \text{ with } x_{k+1} = \frac{1}{\|a_{k+1}\|} a_{k+1}$$
Let $\Delta_d$ is the TT-Laplacian, then $\{a_k\}$ are 'Krylov tensors', i.e.,
$$a_{k+1} = \text{TT-rounding}(\Delta_d x_k, \max\_\text{rank} = 1) \quad \text{with} \quad x_{k+1} = \frac{1}{\|a_{k+1}\|} a_{k+1}$$

<figure><img src="image_34.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* Gram approach

* CGS

* MGS

* CGS2

* MGS2

* Householder transformation
Let $\Delta_d$ is the TT-Laplacian, then $\{a_k\}$ are 'Krylov tensors', i.e.,
$$a_{k+1} = \text{TT-rounding}(\Delta_d x_k, \max\_\text{rank} = 1) \text{ with } x_{k+1} = \frac{1}{\|a_{k+1}\|} a_{k+1}$$

<figure><img src="image_34.png" /><figcaption>Figure: 20 tensors of order $d=3$ and mode $n=15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* Gram approach

* CGS

* MGS

* CGS2

* MGS2

* Householder transformation
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

<figure><img src="image_35.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* Gram approach

* CGS

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

$$A_k = [\mathrm{vec}(\mathbf{a}_1), \dots, \mathrm{vec}(\mathbf{a}_k)]$$

then we conjecture

$$\mathcal{O}(\delta\kappa^2(A_k))$$
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

<figure><img src="image_35.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* Gram approach

* CGS

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

$$A_k = [\mathrm{vec}(\mathbf{a}_1), \dots, \mathrm{vec}(\mathbf{a}_k)]$$

then we conjecture

$$\mathcal{O}(\delta\kappa^2(A_k))$$
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

<figure><img src="image_36.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* MGS

* $\delta\kappa(A_k)$ with

$$A_k = [\mathrm{vec}(\mathbf{a}_1), \dots, \mathrm{vec}(\mathbf{a}_k)]$$

then we conjecture

$$\mathcal{O}(\delta\kappa(A_k))$$
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

<figure><img src="image_36.png" /><figcaption>Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$</figcaption></figure>

* MGS

* $\delta\kappa(A_k)$ with

$$A_k = [\mathrm{vec}(\mathbf{a}_1), \dots, \mathrm{vec}(\mathbf{a}_k)]$$

then we conjecture

$$\mathcal{O}(\delta\kappa(A_k))$$
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

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

Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$

* CGS2

* MGS2

* Householder transformation
then we conjecture

$$ \mathcal{O}(\delta) $$
Loss of orthogonality for {$\mathbf{a}_k$} 'Krylov tensors'

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

Figure: 20 tensors of order $d = 3$ and mode size $n = 15$, compression precision $\delta = 10^{-5}$, computational precision $u = \mathcal{O}(10^{-16})$

* CGS2

* MGS2

* Householder transformation
then we conjecture

$$ \mathcal{O}(\delta) $$
<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(μκ<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(μκ<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(μκ(A<sub>k</sub>))</td><td>O(δκ(A<sub>k</sub>))</td></tr><tr><td>CGS2</td><td>O(μ)</td><td>O(δ)</td></tr><tr><td>MGS2</td><td>O(μ)</td><td>O(δ)</td></tr><tr><td>Householder</td><td>O(μ)</td><td>O(δ)</td></tr></tbody></table>

with *u* the computational precision, *δ* the compression precision
<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(μκ<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(μκ<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(μκ(A<sub>k</sub>))</td><td>O(δκ(A<sub>k</sub>))</td></tr><tr><td>CGS2</td><td>O(μ)</td><td>O(δ)</td></tr><tr><td>MGS2</td><td>O(μ)</td><td>O(δ)</td></tr><tr><td>Householder</td><td>O(μ)</td><td>O(δ)</td></tr></tbody></table>

with *u* the computational precision, *δ* the compression precision
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

> the compression precision $\delta$ independent from the machine architecture

> the compression precision $\delta$ replacing the computational one $u$

> the compression acting norm-wise rather than component-wise

More detailed results can be found at [Coulaud et al. 2022]
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

> the compression precision $\delta$ independent from the machine architecture

> the compression precision $\delta$ replacing the computational one $u$

> the compression acting norm-wise rather than component-wise

More detailed results can be found at [Coulaud et al. 2022]
* All the 6 kernels can be generalized to the TT-framework

* Loss of orthogonality bounds appears to hold true with

> the compression precision $\delta$ independent from the machine architecture

> the compression precision $\delta$ replacing the computational one $u$

> the compression acting norm-wise rather than component-wise

More detailed results can be found at [Coulaud et al. 2022]

What is left out?

* Theoretical proof of the loss of orthogonality bounds

* 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 precision $\delta$ independent from the machine architecture

> the compression precision $\delta$ replacing the computational one $u$

> the compression acting norm-wise rather than component-wise

More detailed results can be found at [Coulaud et al. 2022]

What is left out?

* Theoretical proof of the loss of orthogonality bounds

* Investigating the quality of the tensor subspace spanned by the orthogonal basis
<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>
<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>
[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: 10.1016/B978-0-12-824447-0.00007-8.

[2] Å. Björck. *Numerical Methods for Least Squares Problems*. Society for Industrial and Applied Mathematics, 1996. DOI: 10.1137/1.9781611971484.

[3] Å. 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.

[4] O. Coulaud, L. Giraud, and M. Iannacito. *On some orthogonalization schemes in Tensor Train format*. Tech. rep. RR-9491. Inria Bordeaux - Sud-Ouest, Nov. 2022.

[5] L. Giraud, J. Langou, and M. Rozložnik. "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: 10.1016/j.camwa.2005.08.009.

[6] N. J. Higham. *Accuracy and Stability of Numerical Algorithms*. Second. Society for Industrial and Applied Mathematics, 2002. DOI: 10.1137/1.9780898718027.

[7] M. Iannacito. "Numerical linear algebra and data analysis in large dimensions using tensor format". Theses. Université de Bordeaux, Dec. 2022.

[8] S. J. Leon, Å. Björck, and W. Gander. "Gram-Schmidt orthogonalization: 100 years and more". In: *Numerical Linear Algebra with Applications* 20.3 (2013), pp. 492-532. DOI: 10.1002/nla.1839.

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

[10] 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/S1046827500370883.

[11] 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: 10.1016/B978-0-12-824447-0.00007-8.

[2] Å. Björck. *Numerical Methods for Least Squares Problems*. Society for Industrial and Applied Mathematics, 1996. DOI: 10.1137/1.9781611971484.

[3] Å. 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.

[4] O. Coulaud, L. Giraud, and M. Iannacito. *On some orthogonalization schemes in Tensor Train format*. Tech. rep. RR-9491. Inria Bordeaux - Sud-Ouest, Nov. 2022.

[5] L. Giraud, J. Langou, and M. Rozložnik. "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: 10.1016/j.camwa.2005.08.009.

[6] N. J. Higham. *Accuracy and Stability of Numerical Algorithms*. Second. Society for Industrial and Applied Mathematics, 2002. DOI: 10.1137/1.9780898718027.

[7] M. Iannacito. "Numerical linear algebra and data analysis in large dimensions using tensor format". Theses. Université de Bordeaux, Dec. 2022.

[8] S. J. Leon, Å. Björck, and W. Gander. "Gram-Schmidt orthogonalization: 100 years and more". In: *Numerical Linear Algebra with Applications* 20.3 (2013), pp. 492-532. DOI: 10.1002/nla.1839.

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

[10] 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/S1046827500370883.

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

Questions?
Thanks for the attention.

Questions?