Gkb Golub-Kahan

DFT Discrete Fourier Transformation FFT Fast Fourier Transform DCT Discrete Cosine Transformation

Keywords: HOSVD, decomposition CP, LSQR, T-product, produit cosine, traitement des images et des vidéos, reconnaissance faciale HOSVD, CP decomposition, color image restoration, video restoration, LSQR, face recognition, t-product, cosine product PSF Point Spread Function ALS Alternating Least Squares KPA Kronecker Product Approximation PCA Principal Component Analysis SVD Singular Value Decomposition CP CANDECOMP decomposition

Les tenseurs sont des tableaux multidimensionnels. Ils sont une généralisation des matrices et des vecteurs. Ils fournissent un moyen naturel de représenter les données dans différents domaines, ce qui fait le domaine tensoriel un cadre idéal pour formuler et résoudre de nombreux problèmes dans différents domaines. L'une des applications les plus importantes des tenseurs se trouve dans le domaine du traitement des images, comme la restauration d'images et la reconnaissance des visages, où les images couleur (images RGB) sont présentées comme des tenseurs d'ordre 3 et une vidéo composée d'images couleur est un tenseur d'ordre 4. Mais lorsque nous travaillons dans des espaces de dimension supérieure pour résoudre des problèmes d'ordre supérieur, un ensemble de défis se pose, comme un problème connu sous le nom de "the curse of dimensionality". En effet, lorsque la dimension augmente, les problèmes d'ordre supérieur deviennent plus difficiles (besoin en calcul et en mémoire) car la taille des données d'un tenseur augmente exponentiellement avec l'augmentation de la dimensionnalité du tenseur. Par conséquent, les calculs tensoriels deviennent très coûteux. Dans cette thèse, nous nous sommes concentrés sur la résolution de certains problèmes tensoriels. Les algorithmes proposés sont obtenus en combinant des méthodes itératives telles que la méthode LSQR et des décompositions d'ordre supérieur pour surmonter le problème de la dimensionnalité. Les approches que nous proposons sont appliquées à la restauration d'images et de vidéos. De plus, cette thèse étudie les méthodes de reconnaissance des visages basées sur le format tensoriel, où des outils d'algèbre multi-linéaire tels que la HOSVD (Higher-Order Singular Value Decomposition) ont été utilisés. Nous proposons un nouvel algorithme qui peut être appliqué à une base de données d'images représentées par un tenseur d'ordre 3 ou 4.

Acknowledgements

First of all, I would like to thank Mr. Hassane SADOK, and Mr. Abdeslem Hafid BENTBIB, who have supervised me throughout this thesis and who have shared with me their brilliant intuitions. I would like to thank them for their kindness, their availability despite their numerous duties and for the numerous encouragements they have given me.

I would like to thank all the members of my jury: Mr. Nour Eddine ALAA, Mr. Abdellah BNOUHACHEM, Mr. Hassan SAFOUHI, and Mr. Otmane SOUHAR, for the honor they have given me by agreeing to be the referees of this thesis.

I address all my gratitude to all my friends and to all the people who helped me in the realization of this work. My sincere thanks also goes to Abdelilah HAKIM the director of the laboratory LAMAI-FSTG, and Carole ROSIER the director of the laboratory LMPA-UCLO as well as Khalide JBILOU and Abderrahman BOUHAMIDI to have welcomed me in the research unit and to have allowed me to work under good conditions.

My thanks would be incomplete if I did not mention my dear professors and educators who have accompanied me from kindergarten to my PhD.

Last but not the least, I would like to thank my family: my parents, my sister, and my brother for supporting me throughout writing this thesis and my life in general.

Publications and talks

Publications and submitted papers

• Bentbib, Abdeslem H., Asmaa Khouia, and Hassane Sadok. "The LSQR method for solving tensor least-squares problems." Electronic Transactions on Numerical Analysis 55 (2022): 92-111.

• Bentbib, Abdeslem H., Asmaa Khouia, and Hassane Sadok. "Color image and video restoration using tensor CP decomposition." BIT Numerical Mathematics (2022): 1-22.

• Asmaa khouia, "Higher-Order Tensor Decomposition Applied To Face Recognition" (submited).

List of talks

We presented our works in:

• Numerical Methods for Large Scale Problems which took place: 6-10 June 2022 in Belgrade • Francophone Computer Algebra Days, national days of formal calculation which took place: February 28 -March 4, 2022 in Marseille.

• AHI EVRAN International Conference on Scientific Research, which took place: November 30 -December 1-2, 2021.

Chapter 1 Introduction

Tensors are multidimensional array, and they generalize matrices to higher dimensions.

The order of a tensor is the number of dimensions, also known as ways or modes.

Scalars can therefore be interpreted as zeroth-order tensors, vectors as first-order tensors, and matrices as second-order tensors. The tensors of order three or higher are refereed as higher-order tensors. The following figure shows how we can move from scalars to tensors. Tensors and their decompositions originally appeared in 1927 [START_REF] Frank | The expression of a tensor or a polyadic as a sum of products[END_REF], but have remained untouched by the computer science community until the late 20th century. They have since then spread to numerous other disciplines, including machine learning. In the last decade, the field of tensors has gained a huge interest in many scientific areas.

Examples include computer vision [START_REF] Alex | Multilinear image analysis for facial recognition[END_REF][START_REF] Alex | Multilinear subspace analysis of image ensembles[END_REF], signal processing [START_REF] De | From matrix to tensor: Multilinear algebra and signal processing[END_REF][START_REF] Comon | Tensor decompositions[END_REF], numerical analysis [START_REF] Beylkin | Numerical operator calculus in higher dimensions[END_REF][START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF], neuroscience [START_REF] Christian | Tensorial extensions of independent component analysis for multisubject fmri analysis[END_REF][START_REF] Martınez-Montes | Concurrent eeg/fmri analysis by multiway partial least squares[END_REF] and other more fields. Just as matrices are used to represent linear transformations, tensors can be used to represent more general types of transformations. Tensors are also a natural way of representing multidimensional data, such as images, where grayscale images can be considered as second-order tensors, color images (RGB images) are presented as third-order tensors, and a video composed of color images is a fourth-order tensor. Many of the datasets we deal with today come in a tensor format, such as in medical imaging, where different modalities of medical images of a given organ of a patient are captured in order to make a medical decision. Those images are often structured into a tensor. One of the reasons for the importance of tensors comes from the fact that they can be used to efficiently represent very high-dimensional data. For example, a tensor with 10 6 elements can represent a million-dimensional vector. Thus, tensors represent an essential tool for handling large-scale datasets. In order to handle large size data stored as tensors, higher-order factorizations (originated with Hitchcock in 1927 [START_REF] Frank | The expression of a tensor or a polyadic as a sum of products[END_REF][START_REF] Frank | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF]) are proved as an important tool to extract useful information out of data. Rather than flattening multi-dimensional arrays into matrices and using matrix factorization techniques to preserve the multi-way nature of the data and extract the underlying factors in each dimension, unlike matrix factorizations, the uniqueness of some tensor factorizations can be guaranteed under mild conditions.

One of the most important tensor factorizations are the CANDECOMP/PARAFAC (CP) decomposition and the Tucker decomposition. Those decompositions can be considered to be higher-order extensions of the matrix singular value decomposition. The CP decomposition of multi-dimensional arrays was first introduced by Hitchcock in 1927 [START_REF] Frank | The expression of a tensor or a polyadic as a sum of products[END_REF][START_REF] Frank | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF]. It decomposes a tensor as a sum of rank-one order tensors. It has received a lot of attention in the different areas of science and engineering, such as signal processing, data science, and machine learning. One of the famous methods for computing the CP decomposition is the alternating least squares (ALS). It was first proposed by Carroll and Chang [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "eckart-young" decomposition[END_REF]. There are multiple approaches to obtaining the Tucker decomposition, but the most widely used is the one known as the higher-order singular value decomposition (HOSVD) [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF]. The applications of Tucker decompositions are varied, such as face recognition, human motion, data compression, etc.

Image Restoration

Images are produced in order to record or display useful information. They play a very important role in many aspects of our lives. A gray-scale image is stored and represented numerically as a matrix. Each element of the matrix is called a pixel, which means that the dimension of an image is the number of pixels across the height and the width of the image. Each of these pixels contains a numerical value called pixel values. These pixel values represent the intensity of the pixel and they vary from 0 to 255. The numbers closer to zero represent the darker shades, while the numbers closer to 255 represent the lighter or the white shade. The matrix of pixel values is known as the channel, and in the case of a gray-scale image, there is only one channel. In the case of color images, they can be generated from three primary colors: red, green, and blue, which means a color image is composed of three channels: red channel, green channel, and blue channel.

Due to imperfections in the image formation process, the recorded image often represents a degraded version of the original scene. We can distinguish two sources of degradation: the process of image formation and the process of image recording. The degradation due to the process of image formation is usually denoted by blurring It can be caused by relative motion between the camera and the original scene, or by the optical system. The degradation introduced by the recording process is usually denoted by noise and is due to measurement errors. Even with the most recent technologies, however, there are many situations where the gradations of the original scene are too important. This can be, for example, due to the difficult acquisition conditions encountered and sometimes for reasons of cost.

The field of image restoration is concerned with the problem of undoing the effects of imperfections in the image formation process. More specifically, the goal of image restoration is to estimate the proprieties of the imperfect imaging system (blur) from the observed degraded image. In many fields, images represent a key source of data, making high quality imaging systems very essential. As a consequence, the applications of image restoration are various, such as medical imaging [START_REF] Persons | 3d iterative restoration of tomosynthetic images[END_REF][START_REF] Schweiger | Computational aspects of diffuse optical tomography[END_REF], astronomical imaging [START_REF] Mark | Digital image restoration[END_REF][START_REF] Berry | The handbook of astronomical image processing. The handbook of astronomical image processing[END_REF], surveillance [START_REF] Nguyen | Efficient generalized crossvalidation with applications to parametric image restoration and resolution enhancement[END_REF] etc.

In medical imaging, a number of imaging techniques and devices are being invented.

As in any other imaging system, medical image acquisition devices also introduce degradation to the images. Image restoration methods play an important role in improving the quality of images obtained from medical imaging devices.

One of the most common and important applications of image restoration is in the field of astronomy, where the images obtained from space telescopes are subject to many degradations. They were a result of atmospheric turbulence, aberrations of the optical system, relative motion between the camera and the object, and other reasons due to the enormous expense required to obtain such images. The loss of information due to the degradation of astronomical images could be devastating, which makes astronomical imaging an important application of image restoration.

The linear model of the image restoration problem is described by the Fredholm integral equation of the first kind:

g(x, y) = R 2
k(x, y; s, t) f (s, t)dsdt + η(x, y)

(1.1)
where f is the true image, g is the observed image, and η is additive noise. The kernel function k models the blurring operation and is called the point spread function (PSF). In many situations, the blur is assumed to be spatially invariant, which means the kernel operator k satisfies k(x, y; s, t) = k(xs, yt) . In this case, the degradation model given in (1.1) is formulated as a convolution operation:

g(x, y) = (k f )(x, y) + η(x, y) Unfortunately, the model (1.1) is not very useful for image restoration because of the complexity implied by the possibility of having a different PSF k(x, y; s, t) at each coordinate (s, t) of the image, and it is unrealistic to assume that one might be able to estimate a different PSF for each location in the image. In addition, we do not have a precise function definition for g because the observed image is recorded digitally, and thus is known only at discrete values. Moreover, in many cases, it is necessary to estimate k from measured data. Therefore, it is natural to consider the digital image restoration problem

g = K f + η (1.2)
which is obtained from the equation (1.1) by discretizing the functions and approximating integration with a quadrature rule. K is a matrix that represents the blurring operation, and it can be constructed via the point spread function (PSF) and the boundary conditions [START_REF] Per | Deblurring images: Matrices, spectra and filtering[END_REF] since edges are important structures of the true image and they should be preserved during image restoration. In the case of spatially invariant blur, the structure of the matrix Kdepends on the imposed boundary conditions. For example, if we impose zero boundary conditions by assuming that outside the borders of the image is black everywhere, the matrix K is a block Toeplitz with Toeplitz blocks (BTTB). This assumption is useful when dealing with astronomical images, since most of the time it is possible to assume that the outside borders of the image are black.

If the images are assumed to be matrices of size m × n. Then the blurring matrix K is of size N × N , with N = mn the number of pixels of the image, and g, f , and η are vectors of size N. The problem of image restoration has now been reduced to the problem of solving the equation (1.2). There are a large number of approaches providing solutions to the image restoration problem [START_REF] Harry | Digital image restoration[END_REF][START_REF] Angel | Restoration of images degraded by spatially varying pointspread functions by a conjugate gradient method[END_REF][START_REF] Hafid Bentbib | Global golub-kahan bidiagonalization applied to large discrete ill-posed problems[END_REF]. A classical approach for solving Eq (1.

2) is to calculate its least squares solution

min f K f -g (1.3)
But there are many aspects that make solving the problem (1.3) very challenging. In fact, the dimensions of the matrix K can be extremely large since if the observed image is of size 512 × 512, then the matrix K is of size 262144 × 262144. Thus, the problem (1.3) is large-scale, which makes solving it computationally expensive since these computations usually involve matrix vector products. A successful technique for overcoming this problem is to exploit the structure of the matrix K. In particular, when the blur is separable, which means the horizontal and vertical components of the blur can be separated, in this case, the blur kernel satisfies k(x, y; s, t) = k 1 (x, s)k 2 (y, t) , and so the matrix K can be represented as Kronecker product of two matrices K r and K c ,

K = K r ⊗ K c
If the observed images have m × n pixels, then K r and K c are matrices of size n × n and m × m respectively. This Kronecker decomposition of the matrix K reduces the dimension of the problem (1.2) from mn × mn to m × n. In this case, the blurring model can be formulated in this form,

K c FK T r = G,
where G, F ∈ R m×n , and K T r is the transpose of K r . In the non-separable case, one can approximate the matrix K by solving the Kronecker product approximation (KPA)

problem [START_REF] Loan | Approximation with kronecker products[END_REF].

( Kr , Kc ) = arg min

K r ,K c K -K r ⊗ K c . (1.4)
The second challenging problem of the image restoration problem is that the matrix K is severely ill-conditioned, with singular values decaying to and clustering at 0. This means that the problem (1.3) is sensitive to any perturbation. Regularization is needed to avoid computing solutions that are corrupted by noise. There are several regularization techniques that can be used to regularize the problem (1.2). The most popular regularization approach in the field of image restoration research is the Tikhonov regularization. The Tikhonov seeks to determine a useful approximation of f by replacing the minimization problem (1.3) by the problem of the form:

min f K f -g 2 2 + λ L f 2 2 (1.5)
where λ > 0 is the regularization parameter that need to be chosen, and L is the regularization matrix. A common choice of L is the identity matrix. In this case, we obtain Tikhonov regularization in standard form:

min f K f -g 2 2 + λ f 2 2 (1.6)
The problem (1.6) is equivalent to the problem:

min f K f -ĝ 2 where ĝ = Ö g 0 è , K = Ö K √ λI è
, which satisfy the normal equation:

(K T K + λI) f = K T g.
For the choice of the regularization parameter. There exist numerous methods [START_REF] Calvetti | L-curve curvature bounds via lanczos bidiagonalization[END_REF][START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF][START_REF] Fenu | Gcv for tikhonov regularization via global golub-kahan decomposition[END_REF][START_REF] Kindermann | Convergence analysis of minimization-based noise level-free parameter choice rules for linear ill-posed problems[END_REF], such as generalized cross validation (GCV) [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. The GCV parameter λ GCV is computed to minimize the GCV function,

GCV(λ) = ||K f λ -g|| 2 2 (trace(I -KK -1 λ K T )) 2
, where f λ is solution of K λ f = K T g, and

K λ = K T K + λL T L.

Face recognition

In the past several decades, biometrics has gained a lot of attention, and it has been growing rapidly. A biometric system is a pattern recognition system that recognizes a person based on their physiological characteristics (fingerprints, face, hand contour, ... etc.), or behavioral characteristics (signature, ... etc.). These characteristics are called biometric modalities. The following figure shows the different types of biometric modalities, There are different approaches and algorithms that have been developed for face recognition. Principal Component Analysis (PCA) [START_REF] Sirovich | Low-dimensional procedure for the characterization of human faces[END_REF][START_REF] Turk | Eigenfaces for recognition[END_REF] and Linear Discriminant Analysis (LDA) [START_REF] Peter | Eigenfaces vs. fisherfaces: Recognition using class specific linear projection[END_REF] are well known techniques commonly exploited in the field of face recognition. The PCA approach, known as the eigenface method, is one of the popular methods for feature selection and dimension reduction. It was introduced by Matthew Turk and Alex Pentland in 1991. Numerous extensions to the standard PCA method have been developed, such as multi-linear PCA [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF]. Next, we will give a brief description of the eigenfaces algorithm.

Eigenfaces (PCA)

Eigenface is a set of features obtained by principal component analysis (PCA) building on singular value decomposition (SVD) to project the higher-dimensional face-image space to a lower dimension. Let's consider a set of l images (X 1 , . . . , X l ) of dimension m × n. The first step is to convert these images into vectors (x 1 , . . . , x l ) of size N, with

N = mn.
Then, we calculate the mean of all the face vectors x = 1 l l ∑ i=1

x i , and subtract it from each vector x i , i = 1, . . . , l

a i = x i -x.
We define the matrix A of size N × l as follows,

A = [a 1 a 2 • • • a l ]
The covariance matrix is given by

C = 1 l -1 AA T
where A T is the transpose of A, C is symmetric semi-definite and is orthogonally diagonalizable, which can be written as C = UDU T , where D is a diagonal matrix that contain the eigenvalues of C, U is an orthogonal matrix, where each column u i of U is an eigenvector of C. The eigenvectors of the matrix C are referred to as the principal components. After computing the eigenvectors and the eigenvalues of the covariance matrix C, we chose the k eigenvectors of C corresponding to the k largest eigenvalues, with k < l. Using the first k eigenvectors, we can approximate each of the normalized face vectors a i as a linear combination of (u 1 , . . . , u k ).

a i = k ∑ j=1 α i j u j (1.7)
These u j for j = 1, . . . , k are called eigenfaces, since each eigenface can be viewed as a feature. Using (1.7), we can represent each vector a i by its coordinates

α i = (α 1 , . . . , α k )
with respect to the k principal components. We point out that, for numerical reasons, the covariance matrix C and its eigenvectors are not computed. Instead, given the SVD of the matrix A = USV T , we have with α is the coordinate of ȳ with respect to the k principal components. Recently, multiple extensions of the eigenfaces method have been developed, such as the t-SVD approach [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF] defined via the t-product. This approach differ from the traditional eigenfaces method is that the data is represented as third-order tensor of size m × l × n, which mean in this approach the training images are not vectorized, and similarly 1.3. Chapter-by-chapter overview to the traditional method the covariance tensor and its associated are not computed.

C = 1 l -1 US 2 U T , with D = 1 l-1 S 2 ,
Instead, a generalization of the matrix SVD via the t-product is used.

The performance of face recognition algorithms such as eigenfaces and Fisherfaces (LDA based approach) is good when the only variable that counts for image formation is the identity of the person, which means that the faces are captured under controlled conditions. In reality, there are several factors that can affect the image of a given person, such as illumination, that have a great influence on the appearance of the face in the image. Illumination conditions are unavoidable in the real world, especially when views are collected at different times. Other factors that can affect a face image are view angle, expression... Another problem with using face recognition algorithms such as eigenfaces is that in different fields, a natural representation of images is a thirdorder tensor rather than a simple matrix of vectorized images. Multiple approaches have been developed to overcome these problems. Among the solutions proposed is using a tensor representation of the data. For example, we can represent the face database of subjects photographed in different poses under different illuminations and different face expressions as a fifth-order tensor when the images are represented as matrices, or we can represent the face database as a sixth-order tensor when the images are represented as third-order tensors. In order to manipulate the data stored as a tensor, different multi-linear algebra tools are used. For instance, in [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF][START_REF] Alex | Multilinear image analysis for facial recognition[END_REF] the High-Order SVD (HOSVD, see [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF] ) is used to classify the image of an unknown person, and in [START_REF] Brandoni | Tensor-train decomposition for image recognition[END_REF] the authors explore the use of the Tensor-Train decomposition (for TTdecomposition, see [START_REF] Ivan | Tensor-train decomposition[END_REF]) for multi-feature recognition strategies.

Chapter-by-chapter overview

The thesis can be divided into chapters, a brief summary of every chapter is given as follows:

Chapter 3 is the preliminaries chapter in which we present the basics of multi-linear algebra that will be useful in the development of the techniques presented in this thesis and to help the readers understand multi-linear concepts. We introduce some basic tensor and matrix operations together with their properties, and also give a brief introduction to some of the basic tensor decompositions used in later chapters.

In chapter 4, we propose a new approach to image and video restoration. This approach constructs a degradation model based on a tensor representation, where a color image is represented by a third-order tensor, and a video composed of color images is a fourthorder tensor. Applying tensor CP decomposition to our original problem leads to three subproblems. To solve those subproblems, we apply the global LSQR algorithm and a new algorithm based on Golub Kahan bidiagonalization.

In chapter 5, we are interested in finding an approximate solution X of the tensor least squares minimization problem min

X X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) -G where G ∈ R J 1 ×J 2 ו••×J N and A (i) ∈ R J i ×I i (i = 1, . . . , N) are known, and X ∈ R I 1 ×I 2 ו••×I N is the
unknown tensor to be approximated. Our approach is based on two steps. Firstly, we apply the CP or HOSVD decomposition to the right-hand side tensor G. Secondly, we perform the well-known Golub-Kahan bidiagonalization on each coefficient matrix

A (i) (i = 1, . . . , N)
to obtain a reduced tensor least squares minimization problem. This type of equation may appear in color image and video restorations.

In chapter 6, we explore the use of the t-product, the cosine product, and the outer product, applied to face recognition. The proposed approaches are based on using tensor decompositions of an arrangement of images in a database when we add a factor such as illumination, view angle, or expression. Our algorithms can be applied to a database of images represented by a third or fourth-order tensor. intérêt considérable dans de nombreux domaines scientifiques. Citons par exemple la vision par ordinateur [START_REF] Alex | Multilinear image analysis for facial recognition[END_REF][START_REF] Alex | Multilinear subspace analysis of image ensembles[END_REF], le traitement du signal [START_REF] De | From matrix to tensor: Multilinear algebra and signal processing[END_REF][START_REF] Comon | Tensor decompositions[END_REF], l'analyse numérique [START_REF] Beylkin | Numerical operator calculus in higher dimensions[END_REF][START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF], les neurosciences [START_REF] Christian | Tensorial extensions of independent component analysis for multisubject fmri analysis[END_REF][START_REF] Martınez-Montes | Concurrent eeg/fmri analysis by multiway partial least squares[END_REF] les différents domaines de la science et de l'ingénierie, tels que le traitement du signal, la science des données et l'apprentissage automatique. L'une des méthodes les plus connues pour calculer la décomposition CP est la méthode des moindres carrés alternatifs (ALS). Elle a été proposée pour la première fois par Carroll et Chang [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "eckart-young" decomposition[END_REF]. Il existe de multiples approches pour obtenir la décomposition de Tucker, mais la plus utilisée est celle connue sous le nom de décomposition en valeurs singulières d'ordre supérieur (HOSVD) [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF]. Les applications des décompositions de Tucker sont variées, telles que la reconnaissance des visages, le mouvement humain, la compression de données, etc.

Chapter 2

Introduction en français

Restauration d'images

Les images sont produites afin d'enregistrer ou d'afficher des informations utiles. Elles jouent un rôle très important dans de nombreux aspects de notre vie. Une image en niveaux de gris est stockée et représentée numériquement sous forme de matrice. En raison des imperfections du processus de formation de l'image, l'image enregistrée représente souvent une version dégradée de la scène originale. On peut distinguer deux sources de dégradation : le processus de formation de l'image et le processus d'enregistrement de l'image. La dégradation due au processus de formation de l'image est généralement désignée par le terme "flou". Elle peut être causée par le mouvement relatif entre la caméra et la scène originale, ou par le système optique. La dégradation introduite par le processus d'enregistrement est généralement appelée bruit et est due à des erreurs de mesure. Cependant, même avec les technologies les plus récentes, il existe de nombreuses situations où les gradations de la scène originale sont trop importantes. Cela peut être dû, par exemple, aux conditions d'acquisition difficiles rencontrées et parfois pour des raisons de coût.

Le domaine de la restauration d'images s'intéresse au problème de l'annulation des effets des imperfections dans le processus de formation de l'image. Plus précisément, le but de la restauration d'image est d'estimer les propriétés du système d'imagerie imparfait (flou) à partir de l'image dégradée observée. Dans de nombreux domaines, les images représentent une source de données essentielle, ce qui rend les systèmes d'imagerie de haute qualité très indispensables. Par conséquent, les applications de la restauration d'images sont variées, telles que l'imagerie médicale [START_REF] Persons | 3d iterative restoration of tomosynthetic images[END_REF][START_REF] Schweiger | Computational aspects of diffuse optical tomography[END_REF], l'imagerie astronomique [START_REF] Mark | Digital image restoration[END_REF][START_REF] Berry | The handbook of astronomical image processing. The handbook of astronomical image processing[END_REF], la surveillance [START_REF] Nguyen | Efficient generalized crossvalidation with applications to parametric image restoration and resolution enhancement[END_REF] L'une des applications les plus courantes et les plus importantes de la restauration d'images se trouve dans le domaine de l'astronomie, où les images obtenues à partir de télescopes spatiaux sont soumises à de nombreuses dégradations. Elles sont le résultat de turbulences atmosphériques, d'aberrations du système optique, de mouvements relatifs entre la caméra et l'objet, et d'autres raisons dues à l'énorme dépense nécessaire pour obtenir de telles images. La perte d'informations due à la dégradation des images astronomiques pourrait être dévastatrice, ce qui fait de l'imagerie astronomique une application importante de la restauration d'images.

Le modèle linéaire du problème de restauration d'images est décrit par l'équation intégrale de Fredholm: ( Kr , Kc ) = arg min La Tikhonov cherche à déterminer une approximation utile de f en remplaçant le problème de minimisation (2.3) par le problème de la forme :

g(x, y) = R 2 k(x
K r ,K c K -K r ⊗ K c . ( 2
min f K f -g 2 2 + λ L f 2 2 (2.5)
où λ > 0 est le paramètre de régularisation qui doit être choisi, et L est la matrice de régularisation. Un choix courant de L est la matrice d'identité. Dans ce cas, nous obtenons une régularisation de Tikhonov sous forme standard :

min f K f -g 2 2 + λ f 2 2 (2.6)
Le problème (2.6) est équivalent au problème :

min f K f -ĝ 2 où ĝ = Ö g 0 è , K = Ö K √ λI è
, qui satisfont l'équation normale :

(K T K + λI) f = K T g.
Pour le choix du paramètre de régularisation. Il existe de nombreuses méthodes [START_REF] Calvetti | L-curve curvature bounds via lanczos bidiagonalization[END_REF][START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF][START_REF] Fenu | Gcv for tikhonov regularization via global golub-kahan decomposition[END_REF][START_REF] Kindermann | Convergence analysis of minimization-based noise level-free parameter choice rules for linear ill-posed problems[END_REF], comme la validation croisée généralisée (GCV) [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. Le paramètre GCV λ GCV est calculé en minimisant la fonction GCV, Ensuite, nous donnerons une brève description de l'algorithme des faces propres.

GCV(λ) = ||K f λ -g|| 2 2 (trace(I -KK -1 λ K T )) 2 , où f λ est la solution de K λ f = K T g, et K λ = K T K + λL T L.

Reconnaissance des visages

Faces propres (PCA)

Eigenface est un ensemble de caractéristiques obtenues par l'analyse en composantes principales ( PCA ) en s'appuyant sur la décomposition en valeurs singulières ( SVD ) pour projeter l'espace de l'image du visage de dimension supérieure à une dimension inférieure. Considérons un ensemble de l images (X 1 , . . . , X l ) de dimension m × n. La première étape consiste à convertir ces images en vecteurs (x 1 , . . . , x l ) de taille N, avec

N = mn.
Ensuite, nous calculons la moyenne de tous les vecteurs de visage x = 1 l l ∑ i=1

x i , et on la soustrait de chaque vecteur x i , i = 1, . . . , l.

a i = x i -x.
On définit la matrice A de taille N × l comme suit,

A = [a 1 a 2 • • • a l ]
La matrice de covariance est donnée par

C = 1 l -1 AA T
a i = k ∑ j=1 α i j u j (2.7)
Ces u j pour j = 1, . . . , k sont appelés faces propres, puisque chaque face propre peut être considérée comme une caractéristique. En utilisant (2.7), nous pouvons représenter chaque vecteur a i par ses coordonnées α i = (α 1 , . . . , α k ) par rapport aux k composantes principales. Nous soulignons que, pour des raisons numériques, la matrice de covariance C et ses vecteurs propres ne sont pas calculés. Au lieu de cela, étant donné la SVD de la matrice A = USV T , nous avons [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF][START_REF] Alex | Multilinear image analysis for facial recognition[END_REF] le High-Order SVD (HOSVD, voir [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF] ) est utilisé pour classer l'image d'une personne inconnue, et dans [START_REF] Brandoni | Tensor-train decomposition for image recognition[END_REF], les auteurs explorent l'utilisation de la décomposition Tensor-Train (pour la décomposition TT, voir [START_REF] Ivan | Tensor-train decomposition[END_REF]) pour les stratégies de reconnaissance multi-fonctions.

C = 1 l -1 US 2 U T , avec D = 1 l-1 S 2 ,

Aperçu chapitre par chapitre

La thèse peut être divisée en chapitres, un bref résumé de chaque chapitre est donné L'application de la décomposition CP tensorielle à notre problème original conduit à trois sous-problèmes. Pour résoudre ces sous-problèmes, nous appliquons l'algorithme global LSQR et un nouvel algorithme basé sur la bidiagonalisation de Golub Kahan.

Dans le chapitre 5, nous nous intéressons à la recherche d'une solution approximative X du problème de minimisation des moindres carrés du tenseur Chapter 3

min X X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) -G où G ∈ R J 1 ×J 2 ו••×J N et A (i) ∈ R J i ×I i (i = 1, . . . , N) sont connus, et X ∈ R I 1 ×I 2 ו

Preliminaries

In this chapter we present the basics of multi-linear algebra that will be useful in the development of the techniques presented in this thesis and to help the readers understand multi-linear concepts. We introduce some basic tensor and matrix operations together with their properties, and also give a brief introduction to some of the basic tensor decompositions used in later chapters.

Notation

Throughout this thesis we use the following notation, real numbers or scalars are denoted as lowercase letters x, y ∈ R, vectors as boldface lowercase letters x = (x i ) ∈ R n and matrices as boldface capital letters A = a ij ∈ R m×n . Tensors are denoted by boldface Euler script letters X =

x i 1 i 2 •••i N ∈ R I 1 ×I 2 ו••×I N .
The columns and rows of a matrix are written, respectively, as

A = [a 1 • • • a n ] and A =       a T 1 . . . a T m      

Matrix analysis Definition 3.1.

Let A and B be two matrices of the same size, the inner product of A, B is given by

A, B = tr(A T B) (3.1)
where tr(A) denotes the trace of A, the sum of its diagonal elements, and the corresponding norm of (3.1) is the Frobenius norm defined as

A = » A, A Properties 3.2.1.
Let A and B two matrices of appropriate size, then we have

• A T = A • AB ≤ A B • A -B 2 = A 2 + B 2 -A, B

Matrix products

Apart from the standard matrix product, there are multiple other matrix products that play an important role when working with tensors.

Definition 3.2. (Kronecker product)

The Kronecker product of two matrices A ∈ R I×J and B ∈ R T×R is a matrix denoted as For matrices and vectors of appropriate size, the following properties hold:

A ⊗ B ∈ R IT×JR
1. (A ⊗ B) T = A T ⊗ B T , 2. (A ⊗ B)(C ⊗ D) = AC ⊗ BD, 3. A and B orthogonal ⇒ A ⊗ B orthogonal, 4. (A ⊗ B)v = vec Ä BVA T ä , v = vec(V),
For a matrix of size m × n, vec(A) is a vector of size mn, obtained by stacking the columns of the matrix on top of one another.

Proof.

1. Fallows immediately from Definition 3.2.

For

A ∈ R m×n , B ∈ R p×q , C ∈ R n×r , D ∈ R q×s , (A ⊗ B mp×nq )(C ⊗ D nq×rs ) =       a 11 B • • • a 1n B . . . . . . . . . a m1 B • • • a mn B             c 11 D • • • c 1r D . . . . . . . . . c n1 D • • • c nr D       =       ∑ n k=1 a 1k c k1 BD • • • ∑ n k=1 a 1k c kr BD . . . . . . . . . ∑ n k=1 a mk c k1 BD • • • ∑ n k=1 a mk c kr BD       = AC ⊗ BD ∈ R mp×rs .
3. Assuming that A ∈ R m×m and B ∈ R n×n are orthogonal, and using the properties

(1) and ( 2) of Theorem 3.3, we have

(A ⊗ B) T (A ⊗ B) = Ä A T ⊗ B T ä (A ⊗ B) = A T A ⊗ B T B = I m ⊗ I n = AA T ⊗ BB T = (A ⊗ B) Ä A T ⊗ B T ä = (A ⊗ B)(A ⊗ B) T ,
where I m and I n denote identity matrices of order m and n, respectively. From the definition of Kronecker product, obviously

I m ⊗ I n = I mn . 4. For A ∈ R m×n , B ∈ R p×q and v ∈ R nq , (A ⊗ B)v =       a 11 B • • • a 1n B . . . . . . . . . a m1 B • • • a mn B             v 1 . . . v nq       =       a 11 Bv 1:q + a 12 Bv q+1:2q + • • • + a 1n Bv (n-1)q:nq . . . a m1 Bv 1:q + a m2 Bv q+1:2q + • • • + a mn Bv (n-1)q:nq       = vec à ï Bv 1:q Bv q+1:2q • • • Bv (n-1)q:nq ò       a 11 a 21 • • • a m1 . . . . . . . . . a 1n a 2n • • • a mn       í = vec(B ï v 1:q v q+1:2q • • • v (n-1)q:nq ò V A T ) Definition 3.4. ( Khatri-Rao product)
The Khatri-Rao product (KRP), denoted by , of two matrices A ∈ R m×r and B ∈ R p×r is defined as the "column-wise Kronecker product" given by

A B = a 1 ⊗ b 1 |a 2 ⊗ b 2 | . . . | a r ⊗ b r ,
where a j , b j are the j th columns of A and B respectively. By definition A and B have the same number of columns, and A B ∈ R mp×r .

Theorem 3.5.

For matrices and vectors of appropriate size, the following properties hold:

1. (A ⊗ B)(C D) = AC BD, 2. (A B)v = vec B diag(v)A T , If v ∈ R n , diag (v) denotes diagonal n × n matrix with elements v 1 , v 2 .
. . , v n on the diagonal.

Proof.

1. For A ∈ R m×n , B ∈ R p×q , C ∈ R n×r , D ∈ R q×r , using property (2) of Theorem 3.3,
we have

(A ⊗ B mp×nq )(C D nq×r ) = Å A ⊗ B) ï c 1 ⊗ d 1 • • • c r ⊗ d r ò = ï (A ⊗ B) (c 1 ⊗ d 1 ) • • • (A ⊗ B) (c r ⊗ d r ) ò = ï Ac 1 ⊗ Bd 1 • • • Ac r ⊗ Bd r ò = AC BD.
2. For A ∈ R m×n , B ∈ R p×n and v ∈ R n , using property (4) of Theorem 3.3, we have

(A B)v = ï a 1 ⊗ b 1 a 2 ⊗ b 2 • • • a n ⊗ b n ò       v 1 . . . v n       = (a 1 ⊗ b 1 ) v 1 + (a 2 ⊗ b 2 ) v 2 + • • • + (a n ⊗ b n ) v n = vec Ä b 1 v 1 a T 1 ä + vec Ä b 2 v 2 a T 2 ä + • • • + vec Ä b n v n a T n ä = vec Ä b 1 v 1 a T 1 + b 2 v 2 a T 2 + • • • + b n v n a T n ä = vec           ï b 1 b 2 • • • b n ò           v 1 v 2 . . . v n                     a T 1 a T 2 . . . a T n                     . Definition 3.6. (Hadamard product)
The Hadamard product of two equal-size matrices A, B ∈ R m×n is the element-wise product denoted by * and defined as

A * B =           a 11 b 11 a 12 b 12 • • • a 1n b 1n a 21 b 21 a 22 b 22 • • • a 2n b 2n . . . . . . . . . . . . a m1 b m1 a m2 b m2 • • • a mn b mn           ∈ R m×n Proposition 3.7.
For matrices of appropriate size, the following properties hold:

1. A * B = B * A 2. (A * B) T = A T * B T 3. A * (B * C) = (A * B) * C 4. A * (B + C) = A * B + A * C 5. (A B) T (A B) = A T A * B T B

Singular-Value Decomposition

The most widely known and widely used matrix decomposition is the Singular-Value Decomposition (SVD). It has numerous applications in machine learning and statistics.

The SVD can be used in the calculation of matrix operations, such as matrix inverse and data reduction. It can be used in least-squares linear regression, image compression, restoration, and denoising data.

Theorem 3.8. (SVD) Any m × n matrix A, with m ≥ n, can be factorized

A = U Ö Σ 0 è V T
where U ∈ R m×m and V ∈ R n×n are orthogonal, and Σ ∈ R n×n is diagonal,

Σ = diag (σ 1 , σ 2 , . . . , σ n ) σ 1 ≥ σ 2 ≥ • • • ≥ σ n ≥ 0
Proof. See [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF].

The columns u 1 , . . . , u n of U, which form an orthonormal set, are called left singular vectors.. The columns v 1 , . . . , v n of V , which also form an orthonormal set, are called right singular vectors. The SVD yields a decomposition of A as a sum of n rank-one matrices:

A = UΣV = n ∑ i=1 σ i u i v i (3.2)
The factorization (3.2) can equivalently be expressed by the equations

Av i = σ i u i for i = 1, . . . , n.
In a general a matrix may have many different SVDs. However the following proposition shows that all SVDs involve the same singular values. Thus we may speak of the singular values of a matrix A.

Proposition 3.9. Given any SVD of A, the singular values are the square roots of the nonzero eigenvalues of A A or AA (these matrices have the same eigenvalues).

Proof. We show the result for A A. Given a SVD A = UΣV , we have

A AV = Ä UΣV ä Ä UΣV ä V = V Σ U UΣV V = V Σ ΣV V = V Σ 2
It follows that for i = 1, . . . , r each right singular vector v i of A is an eigenvector of A A with non-zero eigenvalue σ 2 i . The remaining columns of V span the eigenspace of A A corresponding to the eigenvalue zero.

Similarly we can show that the left singular values of A are eigenvectors of AA .

QR Decomposition

Another matrix decomposition is QR decomposition, which is a factorization of a matrix into a product of an orthogonal matrix and a triangular matrix. Proof. See [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF].

The QR decomposition can be symbolically illustrated as follows

Tensor Computation

Tensors are multi-dimensional arrays, which are higher-order generalizations of vectors (first-order tensors) and matrices (second-order tensors). The order of a tensor is the number of dimensions, also known as ways or modes. Tensor can be formally defined as

A = (a i 1 i 2 •••i N ) ∈ R I 1 ×I 2 ו••×I N . If I 1 = I 2 = • • • = I N , then the tensor A is called cubic (cubical). An N-th order cubical tensor Y ∈ R l 1 ×I 2 ו••×I N is diagonal if its elements y i 1 ,i 2 ,...,i N = 0 only if i 1 = i 2 = • • • = i N .
We use I to denote the cubical identity tensor with ones on the superdiagonal and zeros elsewhere. A tensor fiber is a one-dimensional fragment of a tensor, obtained by fixing all indices except for one. A matrix column is a mode-1 fiber and a matrix row is a mode-2 fiber.

Third-order tensors have column, row, and tube fibers, denoted by A(:, j, k), A(i, :, k)

and A(i, j, :) respectively (using Matlab notations). A tensor slice is a two-dimensional section of a tensor, obtained by fixing all indices except for two indices. For example in case of a third-order tensor, we use A(i, :, :), A(:, j, :) and A(:, :, k) to denote horizontal, lateral and frontal slices respectively. A(:, :, k)

is more often denoted as A (k) .

A, B ∈ R I 1 ×I 2 ו••×I N is given by A, B = I 1 ∑ i 1 =1 I 2 ∑ i 2 =1 • • • I N ∑ i N =1 A i 1 •••i N B i 1 •••i N .
It follows immediately that

A, A = A 2 = I 1 ∑ i 1 =1 I 2 ∑ i 2 =1
. . .

I N ∑ i N =1 A 2 i 1 ...i N .
This is analogous to the matrix Frobenius norm.

Matricization

Matricization, also known as unfolding or flattening, is the process of reordering the elements of an N-way array into a matrix. For instance, a 2 × 3 × 4 tensor can be arranged as a 6 × 4 matrix or a 2 × 8 matrix, etc. A special case of matricization is n-mode matrix defined as follow, Definition 3.12 ([20, 56]).

The n-mode matrix of a tensor

A ∈ R I 1 ×I 2 ו••×I N is denoted by A (n) ∈ R I n ×(I 1 •••I n-1 I n+1 •••I N )
and arranges the mode-n fibers into columns of a matrix. More specifically, we have

A (n) (i n , j) = A(i 1 , i 2 , ..., i N ), where j = 1 + N ∑ k=1,k =n (i k -1)J k , and J k = k-1 ∏ m=1,m =n I m .
Example 3.3.1. Let X ∈ R 4×3×2 then the three n-mode matrix are

X (1) =          
x

y = vec(Y) = î Y(:, 1) T , Y(:, 2) T , . . . , Y(:, T) T ó T ∈ R IT .
Analogously, the vectorization of a tensor Y is defined as the vectorization of the associated 1-mode unfolded matrix Y (1) :

vec(Y) = vec(Y (1) ).

Tensor products

Definition 3.14. (n-mode product)

Let X ∈ R I 1 ×I 2 •••×I N be an N th order tensor and U ∈ R J×I n be a matrix. Then the n- mode product of X by U, denoted by X × n U, is a tensor of size

I 1 × I 2 × • • • × I n-1 × J × I n+1 × • • • × I N
whose entries are given by:

(X × n U) i 1 •••i n-1 ji n+1 •••i N = I n ∑ i n =1 X i 1 •••i n-1 ji n+1 •••i N U ji n
The idea can also be expressed in terms of the mode-n matrix:

B = X × n U ⇐⇒ B (n) = UX (n) Example 3.3.2.
Let X ∈ R 4×3×2 be defined by its frontal slices

X 1 =           1 5 9 2 6 10 3 7 11 4 8 12           , X 2 =           13
          , And let U =    1 3 5 2 4 6   . Then the 2-mode product of X and U is y = X × 2 U ∈ R 4×2×2
with frontal slices

Y 1 =           61
          , Y 2 =           169 220 178 232 187 244 196 256           Proposition 3.15.
For tensors and matrices of appropriate size, the following properties hold:

1. X = X (n) 2. X × n ( U + V) = X × n U + X × n V 3.
For distinct modes in a series of multiplication, the order of the multiplication is irrele-

vant, i.e X × m U × n V = X × n v × m U
If the modes are the same, then:

X × n U × n V = X × n VU 4. X , Y × n A = ¨X × n A T , Y ∂ , 5. If U is an orthonormal matrix, then (a) Y = X × n U ⇒ X = Y × n U T , (b) X = X × n U , 6. Y = X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) ⇔ Y (n) = A (n) X (n) Ä A (N) ⊗ • • • ⊗ A (n+1) ⊗ A (n-1) ⊗ • • • ⊗ A (1) ä T . 7. X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) ≤ X A (1) A (2) • • • A (N) .
Proof.

1. Follows immediately from Definition 3.11.

Let

Y = X × n ( U + V), we have Y = X × n ( U + V) ⇐⇒ Y (n) = ( U + V)X (n) = UX (n) + VX (n) , which lead to Y = X × n U + X × n V.
3. For distinct n-mode product, we have

X × m U × n V = ∑ i m ∑ i n X i 1 ...i n ...i m ...i N U ji n V ki m = = ∑ i n i m X i 1 ...i n ...i m ...i N U ji n V ki m = ∑ i m i n A i 1 ...i n ...i m ...i N V ki m U ji n = = ∑ i n ∑ i m X i 1 ...i n ...i m ...i N V ki m U ji n = X × n V × m U
If the modes are the same, we have

X × n U × n V = ∑ i n ∑ i n A i 1 ...i n ...i N U i n i n V ki n n = = ∑ i n A i 1 ...i n ...i N ∑ i n U i n i n V ki n = ∑ i n A i 1 ...i n ...i N ∑ i n n V ki n U i n i n = = ∑ i n A i 1 ...i n ...i N W ki n = B i 1 ...k....i N = X × n VU
4. Follows immediately from Definition 3.11.

5. Let X be a N-order tensor of size

I 1 × • • • × I N and U is an orthonormal matrix of size J × I n . Multiplying Y = X × n U by U T , we have Y × n U T = X × n U × n U T = X × n U T U U is orthonormal (U T U = I)
, which shows (a). On the other hand using orthogonal invariance of matrix Frobenius, we have

X × n U F = UX (n) F = X (n) F = X F 6. Let X be I 1 × • • • × I N tensor and A (n) J n × I n matrices
. By Definition 3.12, each

y j 1 •••j N element of J 1 × • • • × J N tensor y is mapped to y j n k element of matrix Y (n) , with k = 1 + N ∑ l=1 l =n j l -1 l-1 ∏ m=1 m =n J m
We will prove the statement by showing that every element

Ä X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) ä j 1 •••j N maps to element A (n) X (n) Ä A (N) ⊗ • • • ⊗ A (n+1) ⊗ A (n-1) ⊗ • • • ⊗ A (1) ä T j n k
, with k as stated. From Definition 3.14, we have

Ä X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) ä j 1 •••j N = I 1 ∑ i 1 =1 • • • I N ∑ i N =1 x i 1 •••i N a (1)
j 1 i 1 • • • a (N) j N i N
On the other hand, by denoting

M n = A (N) ⊗ • • • ⊗ A (n+1) ⊗ A (n-1) ⊗ • • • ⊗ A (1) T , we have Ä A (n) X (n) M n ä j n k = A (n) j n , : X (n) M n [:, k] = I n ∑ i n =1 a (n) j n i n X (n) M n [i n , k] = I n ∑ i n =1 a (n) j n i n În ∑ i=1 X (n) [i n , i] M n [i, k] with În = I 1 • • • I n-1 I n+1 • • • I N . Now, X (n) [i n , i] = x i 1 •••i N , with i = 1 + N ∑ l=1 l =n (i l -1) l-1 ∏ m=1 m =n I m .
From the definition of Kronecker product follows that the same i stands in

M n [i, k] = ã(N) i N j N • • • ã(n+1) i n+1 j n+1 ã(n+1) i n-1 j n-1 • • • ã(1) i 1 j 1 ,
with ã(m) i m j m denoting an element of A (m) T . Using these conclusions, we can rewrite (3.6) as

I 1 ∑ i 1 =1 • • • I N ∑ i N =1 x i 1 •••i N a (1)
j 1 i 1 • • • a (N) j N i N ,
which completes the proof.

7. Follows directly from the previous property.

Definition 3.16. (n-mode tensor-vector product)

The n-mode multiplication of a tensor

Y ∈ R I 1 ×I 2 ו••×I N by a vector a ∈ R I n is denoted by Y × n a and has dimension I 1 × • • • × I n-1 × I n+1 × • • • × I N , that is, Z = Y × n a ∈ R I 1 ו••×I n-1 ×I n+1 ו••×I N ,
Element-wise, we have

z i 1 ,i 2 ,...,i n-1 ,i n+1 ,...,i N = I n ∑ i n =1 y i 1 ,i 2 ,...,i N a i n with Z k 1 ,...,k N = Y i 1 ,...,i N X j 1 ,...,j N , k n = j n + (i n -1)J n , n = 1, . . . , N.

Notion of rank for tensors

The rank of a matrix is the number of linearly independent column vectors, or, equivalently, the number of non-zero singular values. However, this definition of the matrix rank is not directly extensible to tensors.

Definition 3.19. (Rank-one tensor) An N th order tensor A ∈ R I 1 ×I 2 ו••×I N is rank 1 if it
can be written as the outer product of N vectors:

A = v 1 • v 2 • • • • • v N with v k ∈ R i k
This means that:

A (i 1 , i 2 , . . . , i N ) = N ∑ k=1 v k (i k ) for all 1 ≤ i k ≤ I k v k (i k ) denotes the i th k element of vector v k Definition 3.20. (Rank of a tensor)
The rank of a tensor is the minimum number of rank one tensors we need to add up to get it.

Definition 3.21. (Tensor n-rank)

The n-rank is defined as the number of linearly independent n-mode fibres of a tensor

A : rank n (A) = rank A (n) Definition 3.22. (Multi-linear rank)
The multi-linear rank of a tensor A ∈ R I 1 ×I 2 ו••×I N is a vector, noticed rank A and is given by:

rank A = rank A (1) , . . . , rank A (N)
Where rank A (n) , for 1 ≤ n ≤ N, is the rank of the unfolding matrix A (n)

Tensor decompositions

In this section, we give a brief introduction to higher order decompositions. In particular, we focus on two tensor decompositions, which are CP decomposition that approximate a tensor as sum of rank one tensors, and higher order SVD (HOSVD) decomposition.

CP decomposition

Let A ∈ R I 1 ×I 2 ו••×I N be an N th -order tensor. The CP decomposition [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Henk | Towards a standardized notation and terminology in multiway analysis[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF] of

A is given by A = R ∑ r=1 a (1) r • a (2) r • • • • • a (N) r ,
where a

(k)
r are vectors of size I k with 1 ≤ k ≤ N and R is a positive integer. A CP decomposition of a tensor A is called an exact CP decomposition if R = rank(A), with rank(A) [START_REF] Tamara | Tensor decompositions and applications[END_REF] represent the rank of the tensor A defined as the smallest number of rank-one tensors that generate A as their sum. Unlike matrices, who's the best rank-R approximation is given by the leading R factors of the SVD, the rank of a specific given tensor is hard to define [START_REF] Håstad | Tensor rank is np-complete[END_REF]. In practice, the rank of a tensor is determined numerically by fitting various rank-R CP models. But an interesting property associated with CP decomposition for higher-order tensors is uniqueness under some conditions [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Joseph | Rank, decomposition, and uniqueness for 3-way and n-way arrays[END_REF].

If we define A n = î a (n) 1 a (n) 2 • • • a (n) R ó
for n ∈ {1, . . . , N}, the CP decomposition can be symbolically written as

A = A 1 • A 2 • • • • • A N ,
the matrices A n ∈ R I n ×R are called factor matrices. Often, the vectors a (n) r are chosen such that a (n) r = 1. In this case, the CP decomposition is written as

A = R ∑ r=1 λ r a (1) r • a (2) r • • • • • a (N) r ,
where λ r is a scalar that compensates for the magnitudes of vectors a

(n)

r . Using the n-mode multiplication of a tensor by a matrix, we obtain the following representations:

A = Λ × 1 A 1 × 2 • • • × N A N ,
where

Λ ∈ R R×Rו••×R , with entries Λ i 1 ,••• ,i N =      λ r for i 1 = i 2 = • • • = i N = r, 0 otherwise.
For a given integer R, there are many algorithms to compute CP decomposition. The most popular approach is to apply the alternating least squares method (ALS) [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of "eckart-young" decomposition[END_REF][START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF].

Tucker decomposition Definition 3.23. A tensor

A ∈ R I 1 ×I 2 ו••×I N is said to be in Tucker format if it can be represented as A = S × 1 U (1) × 2 U (2) × 3 • • • × N U (N) were S ∈ R R 1 ×R 2 ו••×R N is
a tensor called the core tensor and U (n) ∈ R I n ×R n are matrices called the factor matrices, with R n = rank n (A), for n = 1, 2, . . . , N. Element-wise, There are multiple approaches to compute the Tucker decomposition, the well-known approach is the method known as the higher-order singular value decomposition (HOSVD), that generalize the matrix SVD.

a i 1 i 2 •••i N = R 1 ∑ r 1 =1 R 2 ∑ r 2 =1 • • • R N ∑ r N =1 s r 1 r 2 •••r N u (1) i 1 r 1 u (2) i 2 r 2 • • • u (N) i N r N .
Theorem 3.24.

Any tensor A ∈ R I 1 ×I 2 ×...×I N can be expressed as the product:

A = S × 1 U (1) × 2 U (2) . . . × × N U (N)
With the properties:

1. U (n) = Ä U (n) 1 , U (n) 2 , . . . , U (n) I n ä ∈ R I n ×I n are orthogonal matrices.
2. The sub-tensors S i n =α of S ∈ R I 1 ×I 2 ×...×I N have the following properties:

(a) All-orthogonality: two sub-tensors < S i n =α , S i n =β >= 0 for all possible values of n, α and β subject to α = β.

(b) Ordering: S i n =1 ≥ S i n =2 ≥ • • • ≥ S i n =I n ≥ 0 for all possible values of n.
The matrix representation of HOSVD can be obtained by unfolding A and S :

A (n) = U (n) • S (n) • Ä U (n+1) ⊗ U (n+2) ⊗ . . . ⊗ U (N) ⊗ U (1) ⊗ U (2) . . . ⊗ U (n-1) ä
The Frobenius-norms S i n =α , symbolized by σ (n) α , are n-mode singular values of A and the vectors U

(n)

i are i th n-mode singular vectors.

The idea of the HOSVD method is to unfold the tensor onto each of its modes and perform SVD on each n-mode matrix of the tensor. The left singular vectors will be the factor matrices.

Algorithm 1 HOSVD for computing Tucker decomposition 1: procedure HOSVD(A)

2:

for n = 1, 2, 3, • • • , N do 3: U (n) ←-left singular vectors of A (n) 4:
end for 5:

S = A × 1 U T 1 × 2 U T 2 × 3 • • • × N U T N . 6:
return S, U (1) ,• • • , U (N) .

Introduction

The field of image restoration is concerned with the problem of undoing the effects of imperfections in the image formation process. More specifically, the goal of image restoration is to remove blur and noise from a degraded image to recover an approximation of the original image. The well-known mathematical model associated with image restoration is formulated as follows

g = Kx. (4.1)
If the images are assumed to have m × n pixels, then g ∈ R mn and x ∈ R mn are vectors that denote the corrupted and the true images respectively, and K ∈ R mn×mn is a matrix that denotes the blurring operator. This model is obtained from discretization of Fredholm integral equations

g(u, v) = R 2 k(u, v; s, t)x(s, t)dsdt + η(u, v), (4.2)
where k models the blurring operation, called the point spread function (PSF) and η is additive noise. The solution of (4.1) is very sensitive to perturbations, since the matrix K is severely ill-conditioned, with singular values decaying and clustering at 0. Due to the presence of the noise, we cannot solve the problem (4.1) directly. Instead, we use Tikhonov regularization [START_REF] Bouhamidi | Sylvester Tikhonov-regularization methods in image restoration[END_REF][START_REF] Heinz W Engl | Regularization Methods for The Stable Solution of Inverse Problems[END_REF][START_REF] Werner Engl | Regularization of Inverse Problems[END_REF] which replaces the minimization problem min

x Kxg 2 2 by a problem of the form:

min x Kx -g 2 2 + λ 2 Lx 2 2 , (4.3)
where λ > 0 is the regularization parameter, and L is the regularization operator. To recover multichannel images from their corrupted observations, we use the RGB representation, by treating each channel; see [START_REF] Per | Deblurring Images: Matrices, Spectra, and Filtering[END_REF][START_REF] Wen | domaines, ce qui fait le domaine tensoriel un cadre idéal pour formuler et résoudre de nombreux problèmes dans différents domaines. L'une des applications les plus importantes des tenseurs se trouve dans le domaine du traitement des images, comme la restauration d'images et la reconnaissance des visages[END_REF], under the assumption that we have the same within-channel blurring(i.e., the same PSF) in all three channels. The problem of recovering multi-channels images was discussed in several works, such as in [START_REF] Abdeslem H Bentbib | Solution methods for linear discrete ill-posed problems for color image restoration[END_REF], where they present approaches to determine a low-rank approximation of the linear system of equations AX = B, based on Golub-Kahan bidiagonalization or block Golub-Kahan bidiagonalization.

In this chapter, we discuss the use of the tensor representation of multichannel images seen as third or fourth-order tensors, and the approximation of the blur matrix as a Kronecker product, to construct the degradation model expressed by the following tensor equation

X × 1 K c × 2 K r = G, (4.4)
where X , G ∈ R m×n×p denote the original image and the degraded image respectively, and × i for i = 1, 2 denote the i-mode matrix product [START_REF] Tamara | Tensor decompositions and applications[END_REF]. K r ∈ R n×n and K c ∈ R m×m define within-channel blurring and they model the horizontal within blurring and the vertical within blurring matrices, respectively. We define an approximate solution X of the tensor equation (4.4) as sum of s rank-one tensors (CP decomposition [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF]) which can be written in case of third-order tensor as follows

X = s ∑ l=1 a l • b l • c l ≡ A, B, C , (4.5)
where the symbol • denotes the outer product [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF], s is a positive integer and a l ∈ R m , b l ∈ R n , c l ∈ R p , for l = 1, . . . , s, are respectively the constituent vectors of the corresponding factor matrices A, B, and C. The advantage of using the CP model is that it is unique under mild assumptions [START_REF] Joseph | Rank, decomposition, and uniqueness for 3-way and n-way arrays[END_REF], and the number of terms in the right-hand side of the decomposition (4.5) can be decided freely without any restriction. Besides, the interpretation of the matrices A, B, and C is easier than other tensor decomposition such as HOSVD [START_REF] Tamara | Tensor decompositions and applications[END_REF] or tensor train decomposition [START_REF] Ivan | Tensor-train decomposition[END_REF] since we deal with matrices instead of tensors. In practice, the best value of s is determined numerically by fitting various rank-s CP models. Based on this approximation, the degradation model becomes

K c A, K r B, C = G. (4.6)
We use the alternating least squares (ALS [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF]) method to solve (4.6), this leads us to three subproblems of the form EXF T = H. We then use a method based on Golub Kahan bidiagonalization to solve them. Attractive results are obtained for numerical applications.

Degradation model

Tensors are multi-dimensional arrays, which are higher-order generalizations of matrices and vectors [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF]. Tensors provide a natural way to represent multidimensional data whose entries are indexed by several continuous or discrete variables, they can be used in several applications, such as images restoration. For instance, a color image is represented as a third-order tensor and a video comprised of color images is seen as a fourth-order tensor. We assume the image to be represented by an array of m × n pixels in each one of the p channels (third-order tensor of size m × n × p). Let g (i) ∈ R mn and

x (i) ∈ R mn for i = 1, 2, . . . , p represent respectively the i-th channel of the corrupted and the original image. The degradation model is of the form

Kx (i) = g (i) , i = 1, 2, . . . , p, (4.7)
where K ∈ R mn×mn represents within channel blurring, which is assumed to be the same in all channels. Using the Kronecker product approximation of the blurring matrix in [START_REF] Kamm | Kronecker product and SVD approximations in image restoration[END_REF][START_REF] James G Nagy | Kronecker product approximations for image restoration with reflexive boundary conditions[END_REF][START_REF] James | Kronecker product approximation for preconditioning in three-dimensional imaging applications[END_REF][START_REF] Loan | Approximation with kronecker products[END_REF] or by assuming that the PSF is separable (the horizontal and vertical components of the blur are separated) [START_REF] Per | Deblurring Images: Matrices, Spectra, and Filtering[END_REF], the matrix K given in (4.7) can be decomposed into a Kronecker product of two matrices K r = (k r i,j ) 1≤i,j≤n and K c .

K = K r ⊗ K c = k r 1,1 K c . . . . . . k r 1,n K c . . . k r n,1 K c . . . . . . k r n,n K c . ( 4.8)
In the non-separable case, one can approximate the matrix K by solving the Kronecker product approximation (KPA) problem [START_REF] Loan | Approximation with kronecker products[END_REF].

( Kr , Kc ) = arg min

K r ,K c K -K r ⊗ K c . ( 4.9)
The degradation model using tensor representation is described in the following proposition Proposition 4.1. The blur model associated with a p-channel image is given by the tensor equation

X × 1 K c × 2 K r = G, (4.10)
where X , G ∈ R m×n×p denote the original image and the degraded image respectively.

Proof. Using Kronecker product notation, (4.7) can be written as

(I p ⊗ K)x = g, (4.11)
where x and g are defined by

x =          
x (1) x (2) . . .

x (p)           , g =           g (1)
g (2) . . .

g (p)          
, Using (4.8), equation (4.11) can be expressed as

(I p ⊗ K r ⊗ K c )x = g, (4.12)
which is equivalent to K c [X (:, :, 1), X (:, :, 2), . . . , X (:, :, p)](I p ⊗ K T r ) = [G(:, :, 1), G(:, :, 2), . . . , G(:, :, p)], (4.13) where G ∈ R m×n×p and X ∈ R m×n×p are third-order tensors defined by G (i) , and X (i) for i = 1, . . . , p, respectively as its frontal slices, with g (i) = vec(G (i) ) and x (i) = vec(X (i) )

for i = 1, . . . , p.

From (4.13) and based on the transpose of each frontal slice, we obtain the following equation K r [X (:, :, 1) T , X (:, :, 2) T , . . . , X (:, :, p) T ](I p ⊗ K T c ) = [G(:, :, 1) T , G(:, :, 2) T , . . . , G(:, :, p) T ].

(4.14)

Using the relations [X (:, :, 1) T , X (:, :, 2) T , . . . , X (:, :, p) T ] = X (2) , and [G(:, :, 1) T , G(:, :, 2) T , . . . , G(:, :, p) T ] = G (2) , equation (4.14) can be expressed as

K r X (2) (I p ⊗ K T c ) = G (2) . (4.15)
Applying tensor properties (see [START_REF] Tamara | Tensor decompositions and applications[END_REF] p.8), equation (4.15) can be written as

(X × 1 K c × 2 K r × 3 I p ) (2) = G (2) , (4.16)
which leads to the equation

X × 1 K c × 2 K r = G. (4.17)
Let us define the linear operator

M : R m×n×p → R m×n×p X → X × 1 K c × 2 K r .
Then, the tensor minimization problem becomes min

X M(X ) -G . (4.18)
The general least squares problem min

X X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) -G was treated
in the next chapter, where we work on the coefficient matrices A (i) (i = 1, . . . , N).

Degradation model using the CP decomposition

The CP decomposition [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Henk | Towards a standardized notation and terminology in multiway analysis[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF] factorizes a tensor into a sum of component rankone tensors. For example, given a third-order tensor X ∈ R I×J×P , we wish to write it as

X = s ∑ l=1 a l • b l • c l ≡ A, B, C , (4.19)
where the symbol • denotes the outer product [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF], s is a positive integer and a l ∈ R I , b l ∈ R J , c l ∈ R P , for l = 1, . . . , s are respectively the constituent vectors of the corresponding factor matrices A, B, and C. The CP decomposition can be represented in matricized forms by applying unfolding representations of the tensor X :

X (1) ∼ = A(C B) T , X (2) ∼ = B(C A) T , X (3) ∼ = C(B A) T .
It is often convenient to assume that all vectors have unit length so that we can use the modified Harshman's PARAFAC model given by

X = s ∑ l=1 λ l a l • b l • c l ≡ λ, A, B, C . Proposition 4.2. Let X = a 1 • a 2 • • • • • a N be a rank-one tensor and V a set of N matrices V 1 , V 2 , . . . , V N . Then we have X × 1 V 1 × 2 V 2 • • • × N V N = V 1 a 1 • V 2 a 2 • • • • • V N a N . (4.20)
Proof. The proof follows from properties of the n-mode matrix product.

Using Proposition 4.2, the degradation model given in (4.10) is formulated by

K c A, K r B, C = G.
Then, the tensor minimization problem defined in (4.18) becomes min

A,B,C G -K c A, K r B, C . (4.21)
In the presence of noise, the minimization problem defined in (4.21) is replaced by a regularized problem min

A,B,C G -K c A, K r B, C + α A A + α B B + α C C , (4.22)
where α A , α B , and α C are nonnegative regularization parameters.

Alternating Least Squares (ALS)

In order to estimate the factor matrices A, B and C defined in (4.21), we use ALS (Alternating Least Squares) algorithm [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF][START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF]. Let T be a third-order tensor defined by the following canonical decomposition

T = K c A, K r B, C , we have T (1) ∼ = K c A(C K r B) T , T (2) ∼ = K r B(C K c A) T , T (3) ∼ = C(K r B K c A) T .
Applying the unfolding representations of the tensors T and G. The minimization problem (4.21) is equivalent to the following three expressions min

A,B,C G (1) -K c A(C K r B) T , min A,B,C G (2) -K r B(C K c A) T , min A,B,C G (3) -C(K r B K c A) T . (4.23)
As a result, instead of solving (4.21) for the three variables one time, we can use the three equations given in (4.23) by fixing all factor matrices but one each time. Thus, given three initial factor matrices A 0 , B 0 and C 0 , the ALS method solves the three leastsquares subproblems in (4.23) to obtain the factor matrices A, B, and C. Starting from the initial guesses A 0 , B 0 , and C 0 , the ALS approach fixes B and C to solve for A, then fixes A and C to solve for B, and then fixes A and B to solve for C. This process continues iteratively until some convergence criterion is satisfied. Therefore, this method translates the original generalistic minimization problem to a three subproblems where each one is just a least-squares problem. The approach is summarized in the following algorithm.

Algorithm 2 ALS-Algorithm

1: procedure ALS-ALGORITHM(G, K r , K C , A 0 ,B 0 ,C 0 ) 2: for k = 0, 2, 3, • • • , M -1 do 3: A k+1 = argmin A G (1) -K c A(C k K r B k ) T 4: B k+1 = argmin B G (2) -K r B(C k K c A k+1 ) T 5: C k+1 = argmin C G (3) -C(K r B k+1 K c A k+1 ) T 6:
end for

where E ∈ R l 1 ×l 2 , F ∈ R l 3 ×l 4 and H ∈ R l 1 ×l 3 . For example, in the first equation defined in (4.23), we have:

E = K c , F = (C k K r B k ), and H = G (1) .
The Global LSQR is an iterative regularized method that can be used to solve the problem (4.24), this method is a generalization of the classical LSQR [START_REF] Christopher | Lsqr: An algorithm for sparse linear equations and sparse least squares[END_REF], which is based on the bidiagonalization procedure of Golub and Kahan. After applying m steps of the Golub-Kahan procedure, the problem is reduced to solving the minimization problem min

y m ||β 1 e 1 -T m y m || 2 , (4.25)
where β 1 is equal to H , e 1 ∈ R m represents the first unit vector, and T m is a lower bidiagonal matrix. The problem (4.25) is solved by computing the QR factorization of the matrix T m (e.g using m-Givens rotation). More details concerning Global LSQR can be found in [START_REF] Jbilou | Global FOM and GMRES algorithms for matrix equations[END_REF][START_REF] Peng | A matrix LSQR iterative method to solve matrix equation AXB= C[END_REF][START_REF] Toutounian | Global least squares method (gl-lsqr) for solving general linear systems with several right-hand sides[END_REF].

In the presence of noise, we introduce Tikhonov regularization method by considering the following minimization problem:

min X H -EXF T + λ LX , (4.26)
where L is the regularization operator, and λ represents the regularization parameter.

In case of LSQR method, the Tikhonov regularization is included at each iteration, by replacing the minimization problem given in (4.25) by the following minimization min

y m ||β 1 e 1 -T m y m || 2 + λ Ly m 2 . (4.27)
This goes through a class of methods called hybrid projection methods that combine iterative projection methods and variational regularization methods (see [START_REF] Chung | Computational methods for large-scale inverse problems: a survey on hybrid projection methods[END_REF]). The suitable parameter λ can be found using different techniques such as generalized cross validation (GCV) [START_REF] Fenu | GCV for Tikhonov regularization by partial SVD[END_REF][START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF][START_REF] Howard | Tikhonov Regularization for Large Scale Problems[END_REF], and L-curve criterion [START_REF] Calvetti | Estimation of the Lcurve via Lanczos bidiagonalization[END_REF][START_REF] Hanke | Regularization methods for large-scale problems[END_REF]. Note that both Chung and Gazzola have done a lot of work on choosing regularization parameters for hybrid methods, e.g., [START_REF] Chung | A weighted GCV method for Lanczos hybrid regularization[END_REF][START_REF] Gazzola | Embedded techniques for choosing the parameter in Tikhonov regularization[END_REF]. If the regularization operator L is equal to the identity matrix, the minimization problem given in (4.27) can be written as

   T m λI    y m -    β 1 e 1 0    2 . (4.28)
Then, the minimization problem (4.27) is equivalent to

|| ē1 -Tm y m || 2 , (4.29)
where Tm =

   T m λI   , ē1 =    β 1 e 1 0    .
The problem (4.29) is solved by applying the QR factorization of the matrix Tm (2m-Givens rotation). To compute a suitable regularization parameter, we will use the discrepancy principle, discussed in [START_REF] Werner Engl | Regularization of Inverse Problems[END_REF].

In the next section, we propose a new approach for solving the minimization problem (4.24) based on the Golub-Kahan bidiagonalization for the matrices E and F, instead of using the Golub Kahan bidiagonalization associated with the linear operator L(X) = EXF T .

Simpler case

: l = rank(H) = 1
In this section, we assume that the matrix H in equation (4.24) is of rank-one, i.e., it can be written as follows

H = h 1 h T 2 .
Applying the Golub-Kahan bidiagonalization algorithm (described in Algorithm 3 below) of the matrix E, taking as initial vector h 1 , leads to the following relations

U k+1 (β 1 e 1 ) = h 1 , EZ k = U k+1 B k , (4.30)
where

U k+1 , Z k = [z 1 z 2 • • • z k ] are orthonormal matrices, β 1 = ||h 1 ||
B k =               α 1 β 2 α 2 β 3 α 3 . . . . . . β k α k               .
In the same way, applying the Golub-Kahan bidiagonalization algorithm of the matrix F, taking as initial vector h 2 leads to the following relations

Q k+1 ( β1 e 1 ) = h 2 , FV k = Q k+1 C k , (4.31)
where 1: input: matrix A ∈ R m×n , initial vector g

Q k+1 , V k = [v 1 v 2 • • • v k ] are
2: β 1 = ||g|| 2 3: u 1 = g/β 1 4: v1 = A T u 1 5: α 1 = || v1 || 2 6: v 1 = v1 /α 1 7: for k = 2, 3, • • • do 8: ûk = Av k-1 -α k-1 u k-1 9
:

β k = || û|| 2 10: u k = ûk /β k 11: vk = A T u k -β k v k-1 12
:

α k = || vk || 2 13: v k = vk /α k 14: end for
The method consists of searching for a low-rank approximation of the form

X k = x k (1) x k (2) T = Z k Y k V T k , (4.32)
which may be written as

X k = Z k y k (1) y k (2) T V k T = x k (1) x k (2) T , (4.33)
where x

(1)

k = Z k y (1)
k , and x

k = V k y (2) k . Theorem 4.3. Let R k = H -EX k F T . Minimizing the residual norm ||R k || is equivalent to minimizing ||(β 1 e 1 )( β1 e T 1 ) -B k y k (1) y k (2) T C k T ||. (2)
Proof. Consider the approximate solution X k = Z k y k (1) y k (2) T V T k . Then the corresponding residual is given by

R k = H -EX k F T , = h 1 h T 2 -EZ k y k (1) y k (2) T V k T F T , = U k+1 (β 1 e 1 )( β1 e T 1 )Q k+1 T -U k+1 B k y k (1) y k (2) T C k T Q k+1 T , = U k+1 (β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T Q T k+1 .
Since U k+1 and Q k+1 are orthonormal, we get

||R k || = ||(β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T ||.
The following minimization problem is accomplished by using the QR factorization of the matrices B k and C k .

||(β 1 e 1 )( β1 e 1 T ) -B k y k (1)

y k (2) T C k T ||. (4.34) Thus Q (1) k B k =    R (1) k 0    , Q (2)
k C k =    R (2) k 0    , and
Q (1) k (β 1 e 1 ) = Ö f (1) k φ(1) k+1 è , Q (2)
k ( β1 e 1 ) = Ö f (2) k φ(2) k+1 è , where Q (1)
k and Q

k are the product of k-Givens rotation to eliminate the subdiagonals elements of the matrices B k and C k respectively. The minimizer y

k and y

k of (4.34)

can then be obtained from

f k (1) f k (2) T = R k (1) y k (1) y k (2) T R k (2) T , (4.35)
which can be split into two subproblems R

(1)

k = f (1)
k , and R

k y

(2)

k = f (2) k .
Therefore an approximate solution is formed as 1) x k (2) T , where x

x k = x k (
k and x

k are obtained from the following relations

     x (1) k = x (1) k-1 + φ (1) k d (1) k x (1) 0 = 0 ,      x (2) k = x (2) k-1 + φ (2) k d (2) k x (2) 0 = 0 , d (1)
k and d

k can be updated by

       d (1) k = 1 ρ (1) k (z k -d (1) k-1 θ (1) k ) d (1) 0 = 0 ,        d (2) k = 1 ρ (2) k (v k -d (2) k-1 θ (2) k ) d (2)
R k 2 = φ(2) k+1 f (1) k 2 + φ(1) k+1 f (2) k 2 + ( φ(1) k+1 φ(2) k+1 ) 2 .
Proof. Using Theorem 4.3 and since Q

k and Q

k are orthonormal. We have

R k 2 = (β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T 2 , = Q k (1) (β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T Q k (2) T 2 , =    f (1) k φ(1) k+1    ï f k (2) T φ(2) k+1 ò -    R (1) k 0    y k (1) y k (2) T ï R k (2) T 0 ò 2 . Since R (1)
k y

(1)

k = f (1)
k , and R

k y

(2)

k = f (2) k . We have R k 2 =    0 φ(2) k+1 f (1) k φ(1) k+1 f k (2) T φ(1) k+1 φ(2) k+1    2 , = φ(2) k+1 f (1) k 2 + φ(1) k+1 f (2) k 2 + ( φ(1) k+1 φ(2) k+1 ) 2 .

General case: l = rank(H) > 1

In this section, we assume that the matrix H in equation (4.24) is of rank l, i.e, it can be written as follows

H = l ∑ i=1 h 1 (i) h 2 (i) T .
We set

H 1 = î h (1) 1 , h (2)
1 , . . . , h

, and

H 2 = î h (1) 2 , h (2)
2 , . . . , h (l) 2 ó . In the numerical examples, we compute H 1 and H 2 by applying the SVD to the matrix H.

Applying the Global Golub-Kahan bidiagonalization to the matrices E and F, starting with H 1 and H 2 , leads to the following relations

U k+1 ((β 1 e 1 ) ⊗ I l ) = H 1 , EZ k = U k+1 (B k ⊗ I l ), (4.36)
Q k+1 (( β1 e 1 ) ⊗ I l ) = H 2 , FV k = Q k+1 (C k ⊗ I l ), (4.37)
where U k+1 , Z k , Q k+1 and V k are F-orthonormal matrices. An approximate solution is given by

X k = l ∑ i=1 x 1 (i) x 2 (i) T = (Z k (y k (1) ⊗ I l ))(V k (y k (2) ⊗ I l )) T = X k (1) X k (2) T , (4.38)
where X

(1)

k = Z k (y (1) k ⊗ I l ), and
X (2) k = V k (y (2) k ⊗ I l ). Lemma 4.5. Let S k = S (1) S (2) • • • S (k)
be an F-orthonormal matrix with S (i) ∈ R m×l for i=1,. . .,k, and let X ∈ R l.k×n and Y ∈ R k×k . We have the following relations:

S k X ≤ X , (4.39)
S k (Y ⊗ I l ) = Y . (4.40)
Proof. We have

S k X 2 = m ∑ i=1 n ∑ j=1 (S k X) 2 ij , = n ∑ j=1 (S k X) :j 2 , with (S k X) ij = k.n ∑ p=1 (S k ) ip X pj = k ∑ r=1 S (r)
i: x rj , and x rj is defined using Matlab notation as follows x rj = X((r -1)l + 1 : rl, j). As consequence r) x rj , which leads to

(S k X) :j = k ∑ r=1 S (
S k X 2 = n ∑ j=1 k ∑ r=1 S (r) x rj 2 ≤ n ∑ j=1 k ∑ r=1 x rj 2 ≤ X 2 .
Therefore (4.39) is achieved. On the other hand, we have:

(Y ⊗ I l ) :p = Y :i ⊗ (I l ) :j with p = j + (i -1)l f or i = 1, . . . , k, and j = 1, . . . , l.

Then

S k (Y ⊗ I l ) 2 = l.k ∑ p=1 (S k (Y ⊗ I l )) :p 2 = l ∑ j=1 k ∑ i=1 k ∑ r=1 S (r) Y ri ⊗ (I l ) :j 2 = l ∑ j=1 k ∑ i=1 k ∑ r=1 Y ri S (r) ⊗ (I l ) :j 2 .
Since S k is F-orthonormal matrix. we obtain

S k (Y ⊗ I l ) 2 = k ∑ i=1 k ∑ r=1 Y 2 ri = Y .
This proves (4.40). Proof. Consider the approximate solution

X k = (Z k (y k (1) ⊗ I l ))(V k (y k (2) ⊗ I l )) T .
Then the associated residual is given by

R k = H -EX k F T , = H 1 H 2 T -E(Z k (y k (1) ⊗ I l ))(V k (y k (2) ⊗ I l )) T F T , = U k+1 ((β 1 e 1 )( β1 e T 1 ) ⊗ I l ) -(B k y k (1) y k (2) T C k T ⊗ I l ) Q k+1 T .
Since U k+1 and Q k+1 are F-orthonormal, by applying the relations in Lemma 4.5, we get

R k ≤ (β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T .
The following minimization problem is accomplished by using the QR factorization of the matrices B k and C k (see section 4.4.1).

(β 1 e 1 )( β1 e 1 T ) -B k y k (1)

y k (2) T C k T . (4.41)
The algorithm can be summarized as follows Algorithm 4

1: input: matrix E ∈ R l 1 ×l 2 , F ∈ R l 3 ×l 4 , H ∈ R l 1 ×l 3 , r 2: Decompose H as H = H 1 H T 2 where H 1 ∈ R l 1 ×r , H 2 ∈ R l 3 ×r 3: Compute β 1 (1) = ||H 1 ||, U 1 (1) = H 1 /β 1 (1) , α 1 (1) = ||E T U 1 (1) || and V 1 (1) = E T U 1 (1) /α 1 (1) 4: Compute β 1 (2) = ||H 2 ||, U 1 (2) = H 2 /β 1 (2) , α 1 (2) = ||F T U 1 (2) || and V 1 (2) = F T U 1 (2) /α 1 (2) 5: Set W 1 (l) = V 1 (l) , φ(l) 1 = β 1 (l) , ρ(l) 1 = α 1 (l) , n (l) 0 = 0, for l = 1, 2 6: for i = 1, 2, . . . , itermax do 7: W(1) i = EV i (1) -α (1) i U (1) i , W (2)
i = FV i (2) -α i (2) U i (2) 8:
for l = 1, 2 do 9:

β i+1 (l) = || W(l) i || 10: U i+1 (l) = W(l) i /β i+1 (l) 11:
S(l

) i = E T U i+1 (l) -β (l) i+1 V i (l) 12: α (l)
i+1 = || S(l) i || 13: V i+1 (l) = S(l) i /α i+1 (l)
14:

ρ i (l) = … Ä ρ(l) i ä 2 + Ä β i+1 (l) ä 2 15: c i (l) = ρ(l) i /ρ i (l)
16:

s i (l) = β i (l) /ρ i (l) 17: θ i+1 (l) = s i (l) α i+1 (l) , 18: ρ(l) i+1 = c i (l) α i+1 (l)
19:

φ i (l) = c i (l) φ(l) i 20: φ(l) i+1 = s (l) i φ(l) i 21: X (l) i = X (l) i-1 + φ (l) i ρ (l) i W (l) i 22: W (l) i+1 = V (l) i+1 - φ (l) i+1 ρ (l) i W (l) i 23:
end for 24: respectively.

X i = X i (1) X i (2) T 25: n i (l) = n i-1 (l) + φ i (l) 2 , for l = 1, 2 26: If » φ(1)2 i+1 n i (2) + φ(2)2 i+1 n i (1) + φ(1)2 i+1 φ (
In the presence of noise, the minimization problem defined in (4.34) is replaced by a regularized problem min y k (1) ,y k (2) (β 1 e 1 )( β1

e 1 T ) -B k y k (1) y k (2) T C k T + λ L 1 y k (1) y k (2) T L 2 T , (4.42)
where L 1 and L 2 are the regularization operators, and λ is the regularization parameter.

Which is equivalent to the following minimization problem min

y (1) k ,y (2) k ||(β 1 β1 )e 1 ⊗ e 1 -(C k ⊗ B k )(y (2) k ⊗ y (1) k )|| 2 + λ||(L 2 ⊗ L 1 )y (2) k ⊗ y (1) k || 2 , (4.43)
which can be written as min

y (1) k ,y (2) k    C k ⊗ B k λL 2 ⊗ L 1    y (2) k ⊗ y (1) k -    (β 1 β1 )e 1 ⊗ e 1 0    2 . ( 4.44)
If we consider the particular case where the matrix L 2 reduces to the identity I k , and

L 1 = B k .
The minimization problem (4.44) can be written as min y

k Ö   C k λI k    ⊗ B k è (y (2) k ⊗ y (1) k ) -    (β 1 β1 )e 1 ⊗ e 1 0    2 , (1) k ,y (2)
which is equivalent in matrix form to min y k (1) ,y k (2) (β 1 e 1 )( β1

e 1 T ) -B k y k (1) y k (2) T Ĉk T , (4.46)
where Ĉk =

   C k λI k   .
The minimization problem (4.46) is accomplished by using the QR factorization of the matrices B k and Ĉk , Q

k B k =    R (1) k 0    , Q (1)
k (β 1 e 1 ) = Ö f (1) k φ(1) k+1 è , (1)
and

Q (2) k Ĉk = Q (2) k    C k λI k    =    R (2) k 0    Q (2) k ( β1 e 1 ) = à f (2) k φ(2) k+1 g k í . ( 4.48)
The QR factorization (4.48) is formed similarly to the case λ = 0, except that two rotations are needed instead of one [START_REF] Christopher | Algorithm 583: LSQR: Sparse linear equations and least squares problems[END_REF].

Proposition 4.7. An upper bound of the residual norm ||R k || 2 is given as follows

R k 2 ≤ ḡk 2 f (1) k 2 + φ(1) k+1 f (2) k 2 + φ(1) k+1 ḡk 2 , where ḡk = Ö φ(2) k+1 g k è .
Proof. Similar to the proof of Proposition 4.4.

A parameter selection method

For an appropriate selection of the regularization parameter λ appearing in (4.43), we suggest to use the GCV method [START_REF] Bouhamidi | Sylvester Tikhonov-regularization methods in image restoration[END_REF][START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF]. The GCV method is a well-known method for computing the regularization parameter λ, the regularization parameter λ GCV is chosen to minimize the GCV function.

GCV(λ) = ||K f λ -g|| 2 2 (trace(I -KK -1 λ K T )) 2
, where

K λ = K T K + λL T L, and f λ is solution of K λ f = K T g. Let K = (C k ⊗ B k ) and L = L 2 ⊗ L 1 ,
the GCV function can be simplified for Tikhonov regularization method using the generalized singular value decompositions (GSVD) of the pairs (C k , L 2 ) and (B k , L 1 ). Thus, there exist orthogonal matrices

U 1 , U 2 , V 1 , V 2 and invertible matrices W 1 , W 2 such that U T 1 C k W 1 = C 1 = diag(c 1,1 , • • • , c k,1 ), U T 2 B k W 2 = C 2 = diag(c 1,2 , • • • , c k,2 ),
and

V T 1 L 1 W 1 = S 1 = diag(s 1,1 , • • • , s k,1 ), V T 2 L 2 W 2 = S 2 = diag(s 1,2 , • • • , s k,2 ).
Then the GSVD of the pair (K, L) is given by

U T KW = C = diag(c 1 , • • • , c N ), V T LW = S = diag(s 1 , • • • , s N ),
where

U = U 2 ⊗ U 1 , V = V 2 ⊗ V 1 , C = C 2 ⊗ C 1 , S = S 2 ⊗ S 1 .

Numerical examples

In this section, we provide some numerical results that illustrate the performance of the two algorithms for solving the problem given in (4. [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF], with application to image and video restoration.

To determine the effectiveness of our method, we evaluate the relative error:

Relative error = ||X (k) -X true || ||X true || ,
where X (k) denotes the approximate restoration. In addition, we evaluate the Signalto-Noise Ratio (SNR) defined by

SNR = 10log 10 ||X true -E(X true )|| 2 ||X (k) -X true || 2 ,
where E(X true ) denotes the mean gray-level of the uncontaminated image X true . All experiments are performed on an AMD Ryzen 7 4700U 2 GHz with 16 Ram using Matlab 2013a.

Example 1

In the first example, we present the numerical results recovered by Algorithm 4 and global LSQR when applied to the restoration of gray-scale and color images that have been contaminated by Defocus blur, whose entries are given by:

k(i, j) =      1 πr 2 , if » i 2 + j 2 ≤ r 0, otherwise , with r = 3.
In the case of gray-scale images, the blurred image has been built by the product K c XK T r . Table 4.1 compares the computing time (in seconds) and the relative error of the computed restorations for the max number of iterations (Iter) showed in the table.

Example 2

This example illustrates the performance of Algorithm 2(Algo 3) applied to the restoration of a 3-channel RGB color image that has been contaminated by blur and noise ( generated by Matlab's randn function) with different noise levels. This noise level is defined as follows ν = E Ĝ where E is a tensor that represents the noise in G, i.e., G := Ĝ + E , and Ĝ is the noise-free image associated with the original image X . In this example, we set s = 180 fixed for noise levels ν = 10 -3 and ν = 5 × 10 -3 , and s = 250 for noise level ν = 10 -2 . The true image of size 384 × 512 × 3 is blurred by Gaussian blur of size 3 × 3 with deviation σ = 2, or by motion blur of length 21 and angle equal to 0. For comparison with the existing approaches in the literature, we compare our approach with the methods proposed in [START_REF] Guide | On tensor GMRES and Golub-Kahan methods via the T-product for color image processing[END_REF]. We refer to these methods as T-global GMRES and T-global Golub Kahan. These algorithms represent the tensor version of GMRES and Golub Kahan bidiagonalization using the T-product, when they are applied for solving the problem arising from recovering multichannel images that have the following form

K c * X * K T r = G, (4.49)
where X , G ∈ R m×n×p denote the original image and the degraded image respectively, K c ∈ R m×m×p , and K r ∈ R n×n×p such that K c (:, :, 1) = K c , K c (:, :, i) = 0, for i = 2, . . . , p and K r (:, :, 1) = K r , K r (:, :, i) = 0, for i = 2, . . . , p. The following table compares the computing CPU-time, the relative errors and the SNR of the computed restorations for a fixed number of iterations. Kahan is faster than the other two algorithms stated in 4.3. We point out that when we increase the number of iterations for the T-global Golub Kahan algorithm, for noise levels ν = 10 -2 and ν = 5 × 10 -3 , the associated algorithm fails to converge.

Example 3

Video processing is not a continuous process but a discrete one. That means each time we deal with videos, we are dealing with the sequence of frames themselves, each frame is just an image, which might be represented as an m × n array of pixels in case of gray-scale videos, and m × n × 3 array in case of color videos. In this example, we evaluate the effectiveness of our algorithm applied to the restoration of a grayscale The following table shows the resulting relative error of the previous example and time of execution

Chapter 5

The LSQR method for solving tensor least squares problem

In this chapter, we are interested in finding an approximate solution X of the tensor least squares minimization problem min

X X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) -G where G ∈ R J 1 ×J 2 ו••×J N and A (i) ∈ R J i ×I i (i = 1, . . . , N) are known, and X ∈ R I 1 ×I 2 ו••×I N is
the unknown tensor to be approximated. Our approach is based on two steps. Firstly, we apply the CP or HOSVD decomposition to the right-hand side tensor G. Secondly, we perform the well-known Golub-Kahan bidiagonalization to each coefficient matrix

A (i) (i = 1, .
. . , N) to obtain a reduced tensor least squares minimization problem. This type of equations may appear in color image and video restorations as we described below. Some numerical tests are performed to show the effectiveness of our proposed method.

Introduction

The LSQR algorithm of Paige and Sanders [START_REF] Christopher | Lsqr: An algorithm for sparse linear equations and sparse least squares[END_REF] is one of the most efficient algorithms for solving the following problem

Ax = b,
or the linear least squares minimization,

min x Ax -b 2 ,
where A is a matrix of size m × n and b is a vector of size m. The LSQR method is analytically equivalent to the conjugate gradient method applied to the associated normal equation. The LSQR algorithm is based on the bidiagonalization procedure of Golub and Kahan [START_REF] Golub | Calculating the singular values and pseudoinverse of a matrix[END_REF]. In the few past years, many researches have generated the LSQR algorithm for solving various equations. For instance, in [START_REF] Toutounian | Global least squares method (gl-lsqr) for solving general linear systems with several right-hand sides[END_REF] the authors proposed a global version of LSQR (GL-LSQR) to obtain an approximate solution of the matrix equation

AX = B with A ∈ R n×n and B ∈ R n×s .
Recently, a Golub-Kahan bidiagonalization based on tensor format is presented in [START_REF] Fatemeh | Golub-kahan bidiagonalization for ill-conditioned tensor equations with applications[END_REF][START_REF] Karimi | Global least squares method based on tensor form to solve linear systems in kronecker format[END_REF] to solve the following tensor equation

A(X ) = B, with X is a tensor of size n d × n d-1 × • • • × n 1 ,
and A is a linear operator defined by

A : R n d ×n d-1 ו••×n 1 -→ R n d ×n d-1 ו••×n 1 X -→ A(X ) = I ∑ i=1 X × 1 A i,d × 2 A i,d-1 × 3 • • • × d A i,1 .
Here, × i for i = 1, . . . , d denote the i-mode product defined below. In those references, the authors used Golub-Kahan bidiagonalization to the linear operator A. For an extensive survey on the subject of higher-order tensors we refer to [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF].

This chapter is concerned with the numerical solution of a tensor least squares problem of the form:

min X L (X ) -S × 1 G (1) × 2 G (2) × 3 • • • × N G (N) , (5.1)
where L is a linear tensor operator, defined by

L : R I 1 ×I 2 ו••×I N -→ R J 1 ×J 2 ו••×J N X -→ X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) ,
where

S ∈ R m 1 ×m 2 ו••×m N , G (i) ∈ R J i ×m i and A (i) ∈ R J i ×I i (i = 1, .
. . , N) are known, and

X ∈ R I 1 ×I 2 ו••×I N
is the unknown tensor to be approximated. We point out that our approach can be generally applied to solve the least square problem:

min X L(X ) -G , (5.2)
for an arbitrary right-hand side tensor G ∈ R J 1 ×J 2 ו••×J N , by decomposing the tensor G

using CP [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF] or Tucker decomposition [START_REF] Ledyard | Implications of factor analysis of three-way matrices for measurement of change[END_REF][START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF] better known as the higher-order SVD (HOSVD) [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF]. It is easy to verify that if the right-hand side tensor G is written in CP or HOSVD format, then the solution X can also be written in CP or HOSVD format. The tensor least squares minimization problem (5.2) is a generalization of the equations arise in color image and video restoration which we discuss in section 5.5. It is not difficult to verify that (5.2) is equivalent to the following minimization problem: 1) , where ⊗ denotes the Kronecker product. The operator 'vec' stacks the columns of a matrix or tensor to form a vector. If I i = J i , for i = 1, . . . , N, the eigenvalues of the matrix A arise as a product of eigenvalues of the matrices A (i) (i = 1, . . . , N). The set of the spectrum of A denoted by λ(A) is given by (see [START_REF] Roger A Horn | Topics in matrix analysis[END_REF]).

min X Avec(X ) -vec(G) , with A = A (N) ⊗ • • • ⊗ A (2) ⊗ A (
λ(A) = {λ 1 λ 2 • • • λ N such that λ i ∈ λ(A (i) ), i = 1, . . . , N},
which leads to the following result Lemma 5.1. If I i = J i , for i = 1, . . . , N, then the solution of the tensor problem (5.2) is unique if and only if

λ 1 λ 2 • • • λ N = 0, for all λ i ∈ λ(A (i) ), i = 1, . . . , N.
When X is an order 2 tensor, that is a matrix X, the tensor least squares minimization (5.2) becomes

min X A (1) XA (2) T -G .
In this work, we are interested in finding an approximate solution of the tensor least squares problem (5.1). Our approach is based on performing the well-known Golub-Kahan bidiagonalization to the pair of matrices A (i) , G (i) for i = 1, . . . , N. More generally, we are interested on solving the tensor problem (5.2), by written the right-hand side tensor G in CP or HOSVD format. Using this approach we can solve the problem where R is the CP rank.

The remainder of the chapter is organized as follows. In the next section, we introduce notations adopted in this chapter and some basic definitions, properties related to tensors. In section 5.3, we construct an approximate solution of the minimization problem (5.2) based on Golub-Kahan bidiagonalization. We work on coefficient matrices A (i) (i = 1, . . . , N) by taking the right-hand side tensor G in rank one format, then we generalized to the case where the right-hand side tensor is approximated using HOSVD decomposition in section 5.4. An example of application to image and video restoration is given in section 5.5. Finally, numerical examples are presented in section 5.6, to show the effectiveness of the proposed approach.

Notations and preliminary concepts

In this section, we recall some of the basic properties about tensors, that will be used in the chapter.

Definition 5.2 ([20, 56]). The n-mode matrix of a tensor

A ∈ R I 1 ×I 2 ו••×I N is denoted by A (n) ∈ R I n ×(I 1 •••I n-1 I n+1 •••I N )
and arranges the mode-n fibers into columns of a matrix.

More specifically, we have

A (n) (i n , j) = A(i 1 , i 2 , ..., i N ), where j = 1 + N ∑ k=1,k =n (i k -1)J k , and J k = k-1 ∏ m=1,m =n I m .
Remark 5.2.1. Let A ∈ R I 1 ×I 2 ×I 3 , we can express the n-mode matrix using slices:

A Proposition 5.3 ([20, 56]). Let X ∈ R I 1 ×I 2 ו••×I N be a N-th order tensor and {U i } 1≤i≤N a set of matrices with U i ∈ R J i ×I i , i = 1, . . . , N.

1.

(X × N i=1 U i ) (n) = U n X (n) (U N ⊗ . . . ⊗ U n+1 ⊗ U n-1 ⊗ . . . ⊗ U 1 ) T . 2. vec(X × N i=1 U i ) = (U N ⊗ U N-1 ⊗ • • • ⊗ U 1 )vec(X ). Proposition 5.4. Let Y ∈ R I 1 ×I 2 ו••×I N and X ∈ R J 1 ×J 2 ו••×J N be two N order tensors, we have the following result Y ⊗ X = Y X .
Proof. Easy to verify using Definition 3.18.

Proposition 5.5 ([20]). Let (a i ) 1≤i≤N be a family of N vectors of sizes I i with i = 1, . . . , N.

We have the following relation

vec(a 1 • a 2 • • • • • a N ) = a N ⊗ a N-1 ⊗ • • • ⊗ a 1 . Proposition 5.6 ([9]). Let Y ∈ R I 1 ×I 2 ו••×I N and X ∈ R J 1 ×J 2 ו••×J N be two tensors, U n ∈ R K×I n and V n ∈ R L×J n . We have (Y ⊗ X ) × n (U n ⊗ V n ) = (Y × n U n ) ⊗ (X × n V n ).

Rank one approximation

In the following section, we assume that the right-hand side tensor G in (5.2) is of rank one, which mean, it can be written in this form

G = g (1) • g (2) • • • • • g (N) ,
where g (i) ∈ R J i for i = 1, . . . , N. Applying Golub-Kahan bidiagonalization algorithm to the pairs (A (i) , g (i) ) for i = 1, . . . , N, leads to the following relations for i = 1, . . . , N, U

1 e 1 ) = g (i) ,

A (i) V (i) k = U (i) k+1 B (i) k , (5.3)
where

U (i) k+1 = [u (i) 1 , . . . , u (i) k+1 ] and V (i) k = [v (i) 1 , . . . , v (i)
k ] are orthonormal basis, e 1 is the first unit vector of R k and B

k is a bidiagonal matrix defined as follow

B (i) k =                  α (i) 1 β (i) 2 α (i) 2 β (i) 3 α (i) 3 . . . . . . β (i) k α (i) k β (i) k+1                  . Proposition 5.7. Let X (k) = Y k × 1 V (1) k × 2 V (2) k × 3 • • • × N V (N) k be an approximate solution of (5.2) with Y k ∈ R k×kו••×k .
The corresponding residual R k can be expressed as

R k = (β (1) 1 e 1 • • • • • β (N) 1 e 1 -Y k × 1 B (1) k × 2 • • • × N B (N) k ) × 1 U (1) k+1 × 2 • • • × N U (N) k+1 . (5.4)
And the residual norm is given by

R k = β (1) 1 e 1 • • • • • β (N) 1 e 1 -Y k × 1 B (1) k × 2 • • • × N B (N) k . ( 5.5)
Proof. Using the relations (5.3), we get

R k = G -X k × 1 A (1) × 2 A (2) × 3 • • • × N A (N) = G -Y k × 1 V (1) k × 2 V (2) k × 3 • • • × N V (N) k × 1 A (1) × 2 A (2) × 3 • • • × N A (N) = G -Y k × 1 A (1) V (1) k × 2 A (2) V (2) k × 3 • • • × N A (N) V (N) k = G -Y k × 1 U (1) k+1 B (1)
k × 2 U (2) k+1 B (2)
k × 3 • • • × N U (N) k+1 B (N) k = G -Y k × 1 B (1) k × 2 B (2) k × 3 • • • × N B (N) k × 1 U (1) k+1 × 2 U (2) k+1 × 3 • • • × N U (N) k+1 .
On the other hand, we have

G = g (1) • g (2) • • • • • g (N) = U (1) k+1 (β (1)
1 e 1 ) • U (2) k+1 (β (2)
1 e 1 ) • • • • • U (N) k+1 (β (N) 1 e 1 ) = (β (1) 1 e 1 • β (2) 1 e 1 • • • • • β (N) 1 e 1 ) × 1 U (1) k+1 × 2 U (2) k+1 × 3 • • • × N U (N) k+1 ,
which shows (5.4). On the other hand U (i)

k+1 for i = 1, . . . , N are orthonormal matrices, which proves (5.5).

The method determines the tensor Y k which minimize R k ,

Y k = argmin Y β (1) 1 e 1 • β (2) 1 e 1 • • • • • β (N) 1 e 1 -Y × 1 B (1) k × 2 B (2) k × 3 • • • × N B (N) k . (5.6)
This minimization problem is accomplished by using the QR decomposition to each matrix B (i)

k for i = 1, . . . , N. Q (i) k B (i) k = Ö R (i) k 0 è , Q (i) k (β (i)
1 e 1 ) = Ö f (i) k φ(i) k+1 è
, where the matrix

Q (i)
k is a product of k Givens rotation chosen to eliminate the subdiagonal elements β (

k+1 and

R (i) k =               ρ (i) 1 θ (i) 2 ρ (i) 2 θ (i) 3 . . . . . . ρ (i) k-1 θ (i) k ρ (i) k              
, and f

(i) k =               φ (i) 1 φ (i) 2 . . . φ (i) k-1 φ (i) k               .
Then, the minimizer Y k of the minimization problem (5.6) can be obtained from the following equation

Y k × 1 R (1) k × 2 R (2) k × 3 • • • × N R (N) k = f (1) k • f (2) k • • • • • f (N) k .
Therefore an approximate solution is given by

X (k) = Y k × 1 V (1) k × 2 V (2) k × 3 • • • × N V (N) k = V k (1) R k (1) -1 f k (1) • V (2) k R (2) k -1 f k (2) • • • • • V (N) k R (N) k -1 f (N) k = x (1) k • x (2) k • • • • • x (N) k , where x k (i) = V (i) k R (i) k -1 f k (i) for i = 1, . . . , N.
Which can be formulated in this form (see [START_REF] Christopher | Lsqr: An algorithm for sparse linear equations and sparse least squares[END_REF]),

     x (i) k = x (i) k-1 + φ (i) k d (i) k , x (i) 0 = 0, where d (i)
k can be updated using the following expression:

       d (i) k = 1 ρ (i) k (v (i) k -θ (i) k d (i) k-1 ), d (i) 0 = 0.
The following lemma is used to prove Theorem 5.9.

Lemma 5.8. Let F and F be two N order tensors of size k

+ 1 × • • • × k + 1 defined by: F =    f (1) k φ(1) k+1    •    f (2) k φ(2) k+1    • • • • •    f (N) k φ(N) k+1    , F =    f (1) k 0    •    f (2) k 0    • • • • •    f (N) k 0    . Then ¨F , F ∂ = F 2 .
Proof.

¨F , F ∂ = k+1 ∑ i 1 =1 k+1 ∑ i 2 =1 • • • k+1 ∑ i N =1 F i 1 ...i N Fi 1 ...i N .
From Definition 3.17, it is easy to verify that Fi

1 ...i N =      F i 1 ...i N for 1 ≤ i 1 , i 2 , . . . , i N ≤ k 0 otherwise. Then ¨F , F ∂ = k ∑ i 1 k ∑ i 2 • • • k ∑ i N F i 1 ...i N F i 1 ...i N = k ∑ i 1 k ∑ i 2 • • • k ∑ i N F 2 i 1 ...i N = F 2 .
In the following proposition, we give an upper bound of the residual norm.

Theorem 5.9. The residual norm ||R k || 2 satisfy the following inequality

R k 2 ≤ P k N ∑ i=1 φ(i)2 k+1 f (i) k 2 ,
where

f (i) k =    f k (i) φ(i) k+1    and P k = N ∏ i=1 f (i) k 2 .
Proof. We have

R k 2 = Y k × 1 B (1) k × 2 • • • × N B k (N) -β 1 (1) e 1 • • • • • β 1 (N) e 1 2 = Ä Y k × 1 B (1) k × 2 • • • × N B k (N) -β (1) 1 e 1 • • • • • β (N) 1 e 1 ä × 1 Q k (1) × 2 • • • × N Q k (N) 2 = Y k × 1    R (1) k 0    × 2 • • • × N    R (N) k 0    -    f (1) k φ(1) k+1    • • • • •    f (N) k φ(N) k+1    2 .
Let Y k be the solution to the following problem:

Y k × 1 R (1) k × 2 R (2) k × 3 • • • × N R (N) k = f (1) k • f (2) k • • • • • f (N) k , then R k 2 =    f (1) k φ(1) k+1 
  •    f (2) k φ(2) k+1    • • • • •    f (N) k φ(N) k+1    -    f (1) k 0    •    f (2) k 0    • • • • •    f (N) k 0    2 .
From Lemma 5.8, the residual norm can be expressed as follow:

R k 2 = F -F 2 = F 2 -2 ¨F , F ∂ + F 2 = F 2 -F 2 = k+1 ∑ i 1 =1 k+1 ∑ i 2 =1 • • • k+1 ∑ i N =1 F 2 i 1 i 2 ...i N - k ∑ i 1 =1 k ∑ i 2 =1 • • • k ∑ i N =1 F 2 i 1 i 2 ...i N ≤ k+1 ∑ i 2 =1 • • • k+1 ∑ i N =1 F 2 k+1i 2 ...i N + k+1 ∑ i 1 =1 • • • k+1 ∑ i N =1 F 2 i 1 k+1...i N + • • • + k+1 ∑ i 1 =1 k+1 ∑ i 2 =1 • • • k+1 ∑ i N-1 =1 F 2 i 1 i 2 ...k+1 .
Using Matlab notation, we have

R k || 2 ≤ F (k + 1, :, . . . , :) 2 + F (:, k + 1, . . . , :) 2 + • • • + F (:, :, . . . , k + 1) 2 .
With F (:, . . . , k + 1, . . . , :)

i-th index f or i = 1, . . . , N, are expressed by

F (:, . . . , k + 1, . . . , :) = φ(i) k+1 f (1) k • • • • • f (i-1) k • f (i+1) k • • • • • f (N) k .
Using Proposition 5.5, we have vec(F (:, . . . , k + 1, . . . , :

)) = φ(i) k+1 f (N) k ⊗ • • • ⊗ f (i+1) k ⊗ f (i-1) k • • • ⊗ f (1) k .
Then

||R k || 2 ≤ N ∑ i=1 φ(i)2 k+1 f (1) k 2 • • • f (i-1) k 2 f (i+1) k 2 • • • f (N) k 2 .

Approximation in HOSVD format

In this section, we assume that the right-hand side tensor G of (5.2) is written in HOSVD format:

G = S × 1 G (1) × 2 • • • × N G (N) ,
where

S ∈ R m 1 ×m 2 ו••×m N and G (i) ∈ R J i ×m i for i = 1, . . . , N. The CP decomposition is a particular case of HOSVD, when m 1 = m 2 = • • • = m N = R and S = I R , where I R
is the tensor identity. Applying global Golub-Kahan bidiagonalization [START_REF] Toutounian | Global least squares method (gl-lsqr) for solving general linear systems with several right-hand sides[END_REF] to the pairs (A (i) , G (i) ) for i = 1, . . . , N, leads to the following relations:

U k+1 (i) (β 1 (i) e 1 ⊗ I m i ) = G (i) , A (i) V (i) k = U (i) k+1 (B (i) k ⊗ I m i ).
(

with

U (i) k+1 = î U (i) 1 U (i) 2 • • • U (i) k+1 ó and V (i) k = î V (i) 1 V (i) 2 • • • V (i) k ó are F-orthonormal.
The method consists in searching an approximate solution of the form

X (k) = Y k × 1 V (1) k × 2 V (2) k × 3 • • • × N V (N) k ,
where Y k is solution to the following minimization problem

min Y S × 1 (β (i) 1 e 1 ⊗ I m 1 ) × 2 • • • × N (β (N) 1 e 1 ⊗ I m N ) -Y × 1 (B (1) k ⊗ I m 1 ) × 2 • • • × N (B (N) k ⊗ I m N ) . (5.8)
In particular, when S is reduce to I R and m 1 = . . . = m N = R, we have the following proposition Proposition 5.10. Let X

(k) = (Y k ⊗ I R ) × 1 V (1) k × 2 V (2) k • • • × N V (N) k with Y k ∈ R k×kו••×k ,
be an approximate solution of (5.2), with the right-hand side tensor G is written in CP decomposition format. Then the corresponding residual R k can be expressed as

R k = (β (1) 1 e 1 • • • • • β (N) 1 e 1 -Y k × 1 B (1) k × 2 • • • × N B (N) k ) ⊗ I R × 1 U (1) k+1 × 2 • • • × N U (N) k+1 .
(5.9)

In this case, The method chooses the tensor Y k which satisfy the minimization problem

Y k = argmin Y β (1) 1 e 1 • β (2) 1 e 1 • • • • • β (N) 1 e 1 -Y × 1 B (1) k × 2 B (2) k × 3 • • • × N B (N) k . (5.10)
Proof. Using the relations (5.7), we get (1) k+1 (B

R k = G -X × 1 A (1) × 2 A (2) × 3 • • • × N A (N) = G -(Y k ⊗ I R ) × 1 V (1) k × 2 V (2) k × 3 • • • × N V (N) k × 1 A (1) × 2 A (2) × 3 • • • × N A (N) = G -(Y k ⊗ I R ) × 1 U
k ⊗ I R ) × 2 U (2) k+1 (B (1)
k ⊗ I R ) × 3 • • • × N U (N) k+1 (B (N) k ⊗ I R ) = G -Y k × 1 B (1) k × 2 B (2) k × 3 • • • × N B (N) k ⊗ I R × 1 U (1) k+1 × 2 U (2) k+1 × 3 • • • × N U (N) k+1 . (2)
On the other hand, we have

G = I R × 1 G (1) × 2 • • • × N G (N) = I R × 1 U (1) k+1 (β (1)
1 e 1 ⊗ I R ) × 2 • • • × N U (N) k+1 (β (N) 1 e 1 ⊗ I R ) = (β (1) 1 e 1 • • • • • β (N) 1 e 1 ⊗ I R ) × 1 U (1) k+1 × 2 • • • × N U (N) k+1 ,
which prove (5.9).

The minimization problem (5.8) is accomplished by using QR decomposition associated to the matrices B (i)

k for i = 1, . . . , N.
Then, the minimizer Y k of the problem (5.8), can be obtained from the following equation:

Y k × 1 (R (1) k ⊗ I m 1 ) × 2 • • • × N (R (N) k ⊗ I m N ) = S × 1 ( f (1) k ⊗ I m 1 ) × 2 • • • × N ( f (N) k ⊗ I m N ).
Therefore an approximate solution is formed as

X (k) = S × 1 V (1) k × 2 • • • × N V (N) k = S × 1 V k (1) (R k (1)-1 f k (1) ⊗ I m 1 ) × 2 • • • × N V (N) k (R k (N)-1 f k (N) ⊗ I m N ) = S × 1 X (1) k × 2 • • • × N X (N) k , (5.11)
where

X k (i) = V k (i) (R (i) k -1 f k (i) ⊗ I m i ) for i = 1, . . . , N. Which can be formulated in this form,      X (i) k = X (i) k-1 + φ (i) k D (i) k X (i) 0 = 0, (5.12) where D (i)
k can be updated using the following expression:

       D (i) k = 1 ρ (i) k (V (i) k -θ (i) k D (i) k-1 )
D (i) 0 = 0. Lemma 5.11. Let V (i) k = î V (i) 1 V (i) 2 • • • V (i) k ó be an F-orthonormal basis with V (i)
l ∈ R I i ×m i for l = 1, . . . , k, and let X be an N order tensors of size km

1 × • • • × km N , and Z ∈ R kו••×k . Then X × i V (i) k X , (5.13)
(Z ⊗ I m i ) × i V (i) k = Z . (5.14)
Proof. See [START_REF] Ah Bentbib | Krylov subspace projection method for sylvester tensor equation with low rank right-hand side[END_REF].

Theorem 5.12. The residual norm ||R k || 2 satisfy the following inequality,

R k 2 ≤ S 2 P k N ∑ i=1 φ(i)2 k+1 f (i) k 2 ,
where

f (i) k =    f k (i) φ(i) k+1    and P k = N ∏ i=1 || f (i) k || 2 .
Proof. Using the first relation in Lemma 5.11, we have

R k 2 Y k × 1 (B (1) k ⊗ I m 1 ) × 2 • • • × N (B (N) k ⊗ I m N ) -S × 1 (β (i) 1 e 1 ⊗ I m 1 ) × 2 • • • × N (β (N) 1 e 1 ⊗ I m N ) 2 . Let Z = Y k × 1 (B (1) k ⊗ I m 1 ) × 2 • • • × N (B (N) k ⊗ I m N ) -S × 1 (β (i) 1 e 1 ⊗ I m 1 ) × 2 • • • × N (β (N) 1 e 1 ⊗ I m N ) 2 .
We have

Z = Y k × N i=1 Ö   R (i) k 0    ⊗ I m i è -S × N i=1 Ö   f (i) k φ(i) k+1    ⊗ I m i è 2 .
Let Y k be the solution to the following equation:

Y × 1 (R (1) k ⊗ I m 1 ) × 2 • • • × N (R (N) k ⊗ I m N ) = S k × 1 ( f (1) k ⊗ I m 1 ) × 2 • • • × N ( f (N) k ⊗ I m N ).
Then

Z = S × N i=1 Ö   f (i) k φ(i) k+1    ⊗ I m i è -S × N i=1 Ö   f (i) k 0    ⊗ I m i è 2 . Let f (i) =    f (i) k φ(i) k+1    , and
f (i) =    f (i) k 0    i = 1, . . . , N.
Using Proposition 5.3, we have

Z ≤ S 2 Ä f (N) ⊗ I m N ä ⊗ • • • ⊗ Ä f (1) ⊗ I m 1 ä -f (N) ⊗ I m N ⊗ • • • ⊗ f (1) ⊗ I m 1 2 2 = S 2 Ä f (N) ⊗ • • • ⊗ f (1) ä ⊗ I m N ⊗ • • • ⊗ I m 1 -f (N) ⊗ • • • ⊗ f (1) ⊗ I m N ⊗ • • • ⊗ I m 1 2 2 = S 2 Ä f (N) ⊗ • • • ⊗ f (1) -f (N) ⊗ • • • ⊗ f (1) ä ⊗ I m N ⊗ • • • ⊗ I m 1 I M 2 2 = S 2 Ä f (N) ⊗ • • • ⊗ f (1) -f (N) ⊗ • • • ⊗ f (1) ä ⊗ I M 2 2 = S 2 f (N) ⊗ • • • ⊗ f (1) -f (N) ⊗ • • • ⊗ f (1) 2 2 = S 2 vec Ä f (1) • • • • • f (N) ä -vec f (1) • • • • • f (N) 2 2 = S 2 f (1) • • • • • f (N) -f (1) • • • • • f (N) 2 .
The result achieved using Theorem 5.9.

The discussed approach can be summarized in Algorithm 5 given below.

Algorithm 5

1: input: Coefficient matrices A (i) for i = 1, . . . , N. The right-hand side tensor G.

2: output: An approximate solution X k to (5.2). for doi = 1, 2, 3 . . . , N 8:

3: Decompose G as S × 1 G (1) × 2 G (2) × 3 • • • × N G (N) 4: Set β 1 (i) = ||G (i) ||, U 1 (i) = G (i) /β 1 , α 1 (i) = ||A (i) T U 1 (i) || and V 1 (i) = A (i) T U (i) 1 /α (i) 1 5: Set W (i) 1 = V (i) 1 , φ(i) 1 = β (i) 1 , ρ(i) 1 = α (i)
Wj = A (i) V (i) j -α (i) j U (i) j , β (i) j+1 = || W(i) j ||, U (i) j+1 = W(i) j /β (i) j+1 9:
S(i)

j = A (i) T U j+1 (i) -β j+1 (i) V j (i) , α (i) j+1 = || S(i) j ||, V j+1 (i) = S(i) j /α j+1 (i) 10: ρ j (i) = ρ(i)2 j + β j+1 (i)2 11: c j (i) = ρ(i) j /ρ j (i) 12: s (i) j = β (i) j /ρ (i) j 13: θ (i) j+1 = s (i) j α (i) j+1 14: ρ(i) j+1 = c (i) j α (i) j+1 15: φ (i)
j = c (i) j φ(i) j 16: n j (i) = n j-1 (i) φ j (i) 2 17: φ(i) j+1 = s (i) j φ(i) j 18: X (i) j = X (i) j-1 + φ (i) j ρ (i) j W (i) j 19: W j+1 (i) = V j+1 (i) -φ i+1 (i) ρ i (i) W j (i) 20: n
(i) j = n j (i) + φ(i)2 j+1 21:
end for 22:

X j = S × 1 X j (1) × 2 • • • × N X j (N) 23: If S Õ N ∏ i=1 n(i) j N ∑ i=1 φ(i)2 j+1 n(i)
(i) j denote f (i) j 2 .
which lead to K c [X (:, :, 1, 1), . . . , X (:, :, 3, k)](I 3 ⊗ K T r ) = [G(:, :, 1, 1), . . . , G(:, :,

which is equivalent to

K c [X (:, :, 1, 1), • • • , X (:, :, 3, k)](I 3 ⊗ I k ⊗ K T r ) = [G(:, :, 1, 1), • • • , X (:, :, 3, k)] (5.20)
The transpose of (5.20) is given by K r [X (:, :, 1, 1) T , . . . , X (:, :, 3, k) T ](

I 3 ⊗ I k ⊗ K T c ) =
î G(:, :, 1, 1) T , . . . , X (:, :, 3, k) T ó (5.21)

As consequence, equation (5.21) can be expressed as

K r X (2) (I 3 ⊗ I k ⊗ K T c ) = G (2) (5.22)
As conclusion the blur model associated to color videos is given by the tensorial equation

X × 1 K c × 2 K r = G (5.
K (i)T c X (:, :, i)K (i) r = G(:, :, i), for i = 1, 2, • • • , p.
By imposing periodic boundary conditions, we assume that the matrices K where F is the Fourier matrix, F is the complex conjugate transpose of FΛ = diag (λ 1 , λ 1 , λ 2 , . . . , λ n ).

where λ k are the eigenvalues of C.

The blurring model associated to the restoration of multi-channel images is described in the following proposition:

Proposition 5.15. Imposing periodic boundary conditions, the blur model associated to pchannel images is given by the tensorial equation

X × 1 F 1 × 2 F * 2 = S,
where F 1 ∈ R m×m , F 2 ∈ R n×n are Fourier matrix and S is defined by S(:,

:, i) = Λ (i)-1 c F 1 G(:, :, i)F * 2 Λ (i)-1 r . with K (i) c = F 1 Λ (i) c F * 1 and K (i) r = F 2 Λ (i) r F * 2 for i = 1, • • • , p.

Numerical examples

In this section, we perform some numerical tests to show the effectiveness of the approach described in this chapter. The first part is devoted to solving the problem (5.2) for given matrices A (i) , i = 1, . . . , N. In the second part, we present some results of application to image and video restoration. In order to solve the problem (5.2), we decompose the right-hand side tensor G ∈ R J 1 ×J 2 ו••×J N using HOSVD or CP decomposition.

In the next tables, we use Algorithm 5-CP to denote Algorithm 5 by decomposing the right-hand side tensor G ∈ R J 1 ×J 2 ו••×J N using CP decomposition as follows,

G = G (1) • G (2) • • • • • G (N) ,
where G (i) ∈ R J i ×R for i = 1, . . . , N. In addition, we use Algorithm 5-HOSVD to denote Algorithm 5 by decomposing G using HOSVD decomposition as follows,

G = S × 1 G (1) × 2 • • • × N G (N) ,
where S ∈ R m 1 ×m 2 ו••×m N , and G (i) ∈ R J i ×m i . All experiments are performed on a 2.7

GHz Intel(R) Core i5 and 8 Ram with Matlab 2016a. In all the tables below, "Iter" stands for the number of iterations.

Part 1:

In this part, we present two numerical examples in order to show the effectiveness of our approach for solving the problem (5.2), for given matrices A (i) , i = 1, . . . , N.

Example 1

In the first example, the coefficient matrices A (i) , i = 1, . . . , N are generated using the following Matlab command:

A (i) = gallery( cycol , [n p], l),
with l = 20. In this case, A (i) are n-by-p matrix with cyclically repeating columns, with rank cannot exceed l. We construct the right-hand side tensor so that all the entries of the exact solution X are equal to one. The following table display the results obtained.

The used stopping criterion is

R k ≤
where is a given tolerance equal to 10 -6 , with number max of iterations equal to 30.

In this example, we decompose the right-hand side tensor using CP decomposition. We point out that the CPU time includes the required time for computing CP decomposition and the construction of the solution X (k) (5.11).

Example 2

In this example, we keep the same data as the previous example, except for the coefficient matrices A (i) , i = 1, . . . , N, they are taken from [START_REF] Lloyd N Trefethen | Numerical linear algebra[END_REF] and they have the same size n:

A (i) = eye(n) + 0.5 sqrt(n) rand(n),
where eye and rand are Matlab functions corresponding to the identity and random matrix, respectively. We compared our approach with GLS-BTF described in [START_REF] Toutounian | Global least squares method (gl-lsqr) for solving general linear systems with several right-hand sides[END_REF]. The numerical results are listed in Table 5.2. and they are obtained by applying Algorithm 5, For this example, the stopping criterion is considered as

R k G < 10 -10 ,
with number max of iterations equal to 160. We point out the CPU-time cover the required time for computing HOSVD decomposition of the right-hand side tensor and the time of construction of the solution X (k) . In the next table 'error' denotes the upper bound of the residual norm described in Theorem 5.12. In this example , we set m i = J i = n for i = 1, . . . , N. The previous table demonstrates the efficiency of our approach especially in term of time of execution. Note that for GLS-BTF method we did not give the exact error and the residual norm, when N = 3, n = 400, and when N = 4, n = 100, due to large CPU-time needed to execute its associated algorithm. In order to show the quality of bound given in Theorem 5.12, we compare between the residual norm and 'error' which represent the upper bound of the residual norm.

We point out that when solving the problem (5.2) for higher dimensions, the approximate solution X (k) (5.11) is not explicitly computed, only the coefficient matrices X (i) k , i = 1, . . . , N are constructed. The numerical results are shown in Table 5.3. We keep the same matrices A (i) , i = 1, . . . , N defined in this example, under the assumption that the right-hand side tensor is written in CP decomposition format:

G = G (1) • G (2) • • • • • G (N) ,
with

G (i) = rand(J i , R), i = 1 . . . , N. The used stopping criterion is R k ≤
where is a given tolerance equal to 10 -10 . with number max of iterations equal to 25 when N = 3, and equal to 30 when N = 4.

Part 2: application to image and video restoration

In this part, we provide some numerical results that illustrate the performance of the approach described in this work applied to the problem of image restoration. To determine the effectiveness of our methods, we evaluate the Relative error.

Relative error = ||X (k) -X true || ||X true || , where X (k) denotes the computed restoration. In addition we evaluate the Signal-to-Noise Ratio (SNR) defined by

SNR = 10log 10 ||X true -E(X true )|| 2 ||X (k) -X true || 2 ,
where E(X true ) denotes the mean gray-level of the uncontaminated image X true . In the following examples, we point out that the CPU time covers both the time of decomposition of the right-hand side tensor G and the construction of the solution. In the next example, we set m = m 1 = m 2 and m i = J i for i = 3, . . . , N, with N = 3 in case of color images and grayscal videos, and N = 4 in case of color videos.

Example 1

This example illustrates the performance of Algorithm 5 applied to the restoration of a 3-channel RGB color image that has been contaminated by Gaussian blur whose function is given by

k(s, t) = 1 2πα 2 exp ß - 1 2α 2 Ä s 2 + t 2 ä ™ ,

Example 2

In this example, we illustrate the effectiveness of our approach applied to the restoration of gray-scale and color videos, seen as third order and fourth order tensors, respectively. The following table shows the results obtained after 20 iterations of Algorithm 5. For completeness, the Figure 5.3 shows the results obtained from the restoration of a video of size 360 × 640 × 30 using Algorithm 5-HOSVD. Note that in the last three experiments displayed in Table 5.5, we present only the results associated to the restoration obtained with algorithm 5-HOSVD, due to the time needed to built the approximate solution tensor.

Conclusion

In this chapter, we proposed a new approach to solve the tensor least squares minimization problem (5.2). We worked under the assumption that the right-hand side tensor is written (or approximated) using CP or higher order singular value decomposition (HOSVD) format. Our goal is to solve the problem (5.2) for higher dimensions, by applying Golub-Kahan bidiagonalization process to each coefficient matrix A (i) for i = 1, . . . , N, and using an LSQR-like method to construct the approximate solution.The presented numerical examples show the effectiveness of the proposed approach.

Chapter 6

Tensor products with application to face recognition

In this chapter, we explore the use of the t-product, the cosine product, and the outer product, applied to face recognition. The proposed approaches are based on using tensor decompositions of an arrangement of images in a database, when we add a factor such as illumination, view angle, or expression. Our algorithms can be applied to a database of images represented by a third or fourth-order tensor. In the numerical results, we compare our approaches with some of the existing methods based on tensor format.

Introduction

Face recognition represents one of the most successful applications of image analysis. It has recently gained increasing interest [START_REF] Brandoni | Tensor-train decomposition for image recognition[END_REF][START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF][START_REF] Anil | Handbook of face recognition[END_REF]. It is due to several applications such as security and video-surveillance. In the past few years, there have been multiple methods being developed to achieve high performance. The PCA (Principal Component Analysis), often goes by the name "eigenfaces" is one of those techniques which is based on matrix decomposition. The eigenfaces form a basis set of all images used to built the covariance matrix. This generates dimension reduction by letting the smaller set of basis images to represent the original training images. Classification can be obtained by comparing how faces are represented by the basis set. This approach was generalized in [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF] using the t-product where the images are considered as twodimensional rather than vectorized objects. However, this method works better when all pictures are taken under similar conditions, and it does not perform well when several factors are varied such as illumination, view angel, and expression. Recently, the face recognition problem was considered using a tensor model. By letting the modes of the tensor represent a different viewing condition, e.g., illumination or facial expression, the precision of the recognition algorithms is improved compared to the PCA method. This tensor representation was recently considered in several works. For instance, in [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF][START_REF] Alex | Multilinear image analysis for facial recognition[END_REF] the High-Order SVD (HOSVD, see [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF] ) is used to classify the image of an unknown person, and in [START_REF] Brandoni | Tensor-train decomposition for image recognition[END_REF] the authors explore the use of the Tensor-Train decomposition (for TT-decomposition, see [START_REF] Ivan | Tensor-train decomposition[END_REF]) for multi-feature recognition strategies.

In this chapter, we are interested in using the tensor decompositions: QR and SVD defined via the t-product or the cosine product, and the CP decomposition defined via the outer product. The approach based on the t-product and the cosine product differs from the one given in [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF], in terms of the data representation. In fact, by adding a factor such as illumination, view angel or expression to construct a dataset of images, we use third-order tensors to represent the dataset of images, when the images are in vectorized format, or as fourth-order tensors, when the images are considered as a two-dimensional array. This approach can also be extended to the case of fifth-order tensors when the images are considered as third-order arrays. Note that in [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF], the representation is limited to the case of third-order tensors.

The remainder of the chapter is organized as follows. In sections 6.2 and 6.3, we recall some properties associated to the t-product and the cosine product between third-order tensors, at the same time, we establish some properties for general higher-order tensors.

In section 6.4, we give a brief introduction to the CP decomposition. In section 6.5, we present new algorithms for tensor face recognition problem, using t-product, cosine product, and CP decomposition. Finally, numerical examples are presented in Section 6.6 that show the effectiveness of the proposed approaches.

t-product, definition and properties

In this part, we recall some definitions and properties related to the t-product in case of third-order tensors [START_REF] Kilmer | Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging[END_REF][START_REF] Kilmer | Factorization strategies for third-order tensors[END_REF]. Let A ∈ R n 1 ×n 2 ×n 3 be a third-order tensor, the operators bcirc, unfold and fold are described by: bcirc

(A) =           A (1) A (n 3 ) A (n 3 -1) • • • A (2) A (2) A (1) A (n 3 ) • • • A (3) . . . . . . . . . . . . . . . A (n 3 ) A (n 3 -1) • • • A (2) A (1)           , unfold(A) =           A (1)
A (2) . . . Before stating some properties associated to the block circulant matrices, it is helpful to look at the simplest form of a block-circulant matrix, that is a circulant matrix [START_REF] Robert | Toeplitz and circulant matrices: A review[END_REF][START_REF] Kra | On circulant matrices[END_REF]. Definition 6.1. A square matrix C of size n × n is a circulant matrix if it has the following form:

A (n 3 )           , bcirc is
C =           r 1 r n r n-1 • • • r 2 r 2 r 1 r n • • • r 3 . . . . . . . . . . . . . . . r n r n-1 • • • r 2 r 1          
, where each row is a cyclic shift of the row above to the right.

Follows from a simple calculations, that Z i e l = e i+l , with all indices are interpreted mod n 3 , which mean if i + l > n 3 , then i + l becomes l with i + l = l + n 3 .

bcirc(A)e k = bcirc(A)e l ⊗ e m = e l ⊗ A (1) e m + e l+1 ⊗ A (2)

e m + • • • + e l-1 ⊗ A (n 3 ) e m = e l ⊗ n 1 ∑ i=1 A (1) (i, m)e i + e l+1 ⊗ n 1 ∑ i=1 A (2) (i, m)e i + • • • + e l-1 ⊗ n 1 ∑ i=1 A (n 3 ) (i, m)e i = n 1 ∑ i=1 Ä A (1) (i, m)e l + A (2) (i, m)e l+1 + • • • + A (n 3 ) (i, m)e l-1 ä ⊗ e i .
On the other hand, we have

n 1 ∑ i=1 n 2 ∑ j=1 circ(A(i, j, :)) ⊗ e i e T j e k = n 1 ∑ i=1 n 2 ∑ j=1 circ(A(i, j, :)) ⊗ e i e T j e l ⊗ e m = n 1 ∑ i=1 circ(A(i, m, :))e l ⊗ e i = n 1 ∑ i=1 Ä A (1) (i, m)I + A (2) (i, m)Z + • • • + A (n 3 ) (i, m)Z n 3 -1 ä e l ⊗ e i = n 1 ∑ i=1 Ä A (1) (i, m)e l + A (2) (i, m)e l+1 + • • • + A (n 3 ) (i, m)e l-1 ä ⊗ e i .
Therefore, the k-th column of bcirc(A) and Theorem 6.3. The matrix bcirc(A) is block diagonalizable. We have:

(F n 3 ⊗ I n 1 )bcirc(A)(F n 3 ⊗ I n 2 ) =           D (1)
D (2) . . .

D (n 3 )           . Proof. It follows from: bcirc(A) = n 1 ∑ i=1 n 2 ∑ j=1 circ(A(i, j, :)) ⊗ e i e T j .
Definition 6.4 ([52, 54]). The t-product * t between two tensors A ∈ R n 1 ×n 2 ×n 3 and B ∈ R n 2 ×m 2 ×n 3 is a tensor given by:

C = A * t B = fold(bcirc(A)unfold(B)), (6.2)
C is of size n 1 × m 2 × n 3 .

Since bcirc(A) is a block diagonalized using the discrete Fourier transform. The tproduct can be computed as

bcirc(A)unfold(B) = (F n 3 ⊗ I n 1 )(F n 3 ⊗ I n 1 )bcirc(A)(F n 3 ⊗ I n 2 )(F n 3 ⊗ I n 2 )unfold(B) = (F n 3 ⊗ I n 1 )           D (1)
D (2) . . .

D (n 3 )           unfold( B), with unfold( B)) = (F n 3 ⊗ I n 2 )unfold(B)
. This can be computed by applying the Fast Fourier Transform (FFT) along the tubes of B. Using the FFT along the third dimension, the t-product can be computed using the following algorithm Algorithm 6 t-Product Computing Using FFT [START_REF] Kilmer | Factorization strategies for third-order tensors[END_REF] 1

: input: A ∈ R n 1 ×n 2 ×n 3 , B ∈ R n 2 ×m 2 ×n 3 2: output : C = A * t B ∈ R n 1 ×m 2 ×n 3 3: Â ←-f f t(A, [] , 3), B ←-f f t(B, [] , 3)
4: for i = 1, 2, . . . , n 3 do 5:

Ĉ(:, :, i) = Â(:, :, i) B(:, :, i)

6: end for

7: C ←-i f f t( Ĉ [] ,
3)

The command f f t(A, [] , 3) is a Matlab command that computes the fast Fourier transform (FFT) on the third dimension of a multiway array. In the next part of this section, we introduce some properties and definitions associated to the t-product, that we are

The t-SVD can be computed using the fast Fourier transform, as described in the following algorithm î Û (:, :, i), Ŝ(:, :, i), V(:, :, i) ó = svd( Â(:, :, i))

6: end for

7: U ←-i f f t( Û [] , 3), S ←-i f f t( Ŝ [] , 3), V ←-i f f t( V [] ,
3) Theorem 6.11 ([54]). (t-QR) Let A be an n 1 × n 2 × n 3 real-valued tensor. Then A can be factored as

A = Q * t R, (6.7)
where Q is an orthogonal tensor of size n 1 × n 1 × n 3 and R is an f-upper triangular tensor of size n 1 × n 2 × n 3 .

The t-QR can be computed using the fast Fourier transform as follows î Q(:, :, i), R(:, :, i) ó = qr( Â(:, :, i))

6: end for

7: Q ←-i f f t( Q [] , 3), R ←-i f f t( R [] , 3)

Generalization of the t-product to higher-order tensors

The t-product can be extended to higher-order tensors, in a recursive manner(see [START_REF] Carla D Martin | An order-p tensor factorization with applications in imaging[END_REF]).

For instance, in case of fourth-order tensors, the t-product is defined as follows Definition 6.12. Let A ∈ R n 1 ×n 2 ×n 3 ×n 4 and B ∈ R n 2 ×m 2 ×n 3 ×n 4 be two fourth-order tensors. The t-product A * t B is a fourth-order tensor of size n 1 × m 2 × n 3 × n 4 given by:

C = A * t B = fold(bcirc(A) * t unfold(B)), (6.8)
where

bcirc(A) =           A (1) A (n 4 ) A (n 4 -1) • • • A (2) A (2) A (1) A (n 4 ) • • • A (3) . . . . . . . . . . . . . . . A (n 4 ) A (n 4 -1) • • • A (2) A (1)           and unfold(A) =           A (1)
A (2) . . .

A (n 4 )          
, with A (i) = A(:, :, :, i), for i = 1, . . . , n 4 . Lemma 6.13. Let A ∈ R n 1 ×n 2 ×n 3 ×n 4 be a fourth-order tensor, we have:

bcirc(A) = n 1 ∑ i=1 n 2 ∑ j=1 n 3 ∑ k=1 circ(A(i, j, k, :)) ⊗ e i • e j • e k . ( 6.9)
Similarly, for higher-order tensors we get:

bcirc(A) = n 1 ∑ i 1 =1 n 2 ∑ i 2 =1 • • • n N-1 ∑ i N-1 =1 circ(A(i 1 , i 2 , • • • , i N-1 , :)) ⊗ e i 1 • e i 2 • • • • • e i N-1 , (6.10) with A is a N-order tensor of size n 1 × n 2 × • • • × n N .
Proof. For definiteness, we focus on proving (6.9), (6.10) is proved in recursive way. In analogue to the case of circulant matrices, we can write bcirc(A) as follows,

bcirc(A) = I ⊗ A (1) + Z ⊗ A (2) + • • • + Z n 3 -1 ⊗ A (n 3 ) = n 3 ∑ i=1 Z i-1 ⊗ A (i) Let e m ∈ R n 3 . Then bcirc(A) × 3 e T m = n 3 ∑ i=1 Z i-1 ⊗ A (i) × 3 e T m = n 3 ∑ i=1 Z i-1 ⊗ A (i) (:, :, m).
Based on Lemma 6.2, we obtain bcirc

(A) × 3 e T m = n 1 ∑ i=1 n 2 ∑ j=1 circ(A(i, j, m, :)) ⊗ e i • e j .
On the other hand, we have

n 1 ∑ i=1 n 2 ∑ j=1 n 3 ∑ k=1 circ(A(i, j, k, :)) ⊗ e i • e j • e k × 3 e T m = n 1 ∑ i=1 n 2 ∑ j=1 circ(A(i, j, m, :)) ⊗ e i • e j ,
which proves that the m-th frontal slice of bcirc(A) and

n 1 ∑ i=1 n 2 ∑ j=1 n 3 ∑ k=1 circ(A(i, j, k, :)) ⊗ e i •
e j • e k are equal. This proves (6.9).

Theorem 6.14. The tensor bcirc(A) is block diagonalizable.

Proof. It follows from:

bcirc(A) = n 1 ∑ i=1 n 2 ∑ j=1 n 3 ∑ k=1 circ(A(i, j, k, :)) ⊗ e i • e j • e k .
To compute the t-product, we apply the FFT along the fourth dimension. The associated algorithm is expressed as follows Algorithm 9 Fourth-order t-Product using FFT

1: input: A ∈ R n 1 ×n 2 ×n 3 ×n 4 , B ∈ R n 2 ×m 2 ×n 3 ×n 4 2: output: C = A * t B ∈ R n 1 ×m 2 ×n 3 ×n 4 3: Â ←-f f t(A, [] , 4), B ←-f f t(B, [] , 4) 4: for i = 1, 2, . . . , n 4 do 5:
Ĉ(:, :, :, i) = Â(:, :, :, i) * t B(:, :, :, i)

6: end for 7: C ←-i f f t( Ĉ [] , 4)
Remark 6.2.1. In the same way, notations, definitions and decompositions associated to the t-product can be extended to higher-order tensors. Proposition 6.9 can be generalized also to higher-order tensors.

Cosine transform product, definition and properties

Let A ∈ R n 1 ×n 2 ×n 3 be a third-order tensor, we define

mat(A) =           A (1) A (2) • • • A (n 3 ) A (2) A (1) • • • A (n 3 -1) . . . . . . . . . . . . A (n 3 ) A (n 3 -1) • • • A (1)           +           A (2) • • • A (n 3 ) 0 . . . . . . . . . A (n 3 ) A (n 3 ) 0 . . . . . . 0 A (n 3 ) • • • A (2)          
, is block Toepltiz-plus-Hankel matrix. Here 0 denotes the zero matrix of size n 1 × n 2 .

The command ten(.) is defined as the inverse of the mat operation:

ten(mat(A)) = A. The cosine product can be computed without forming explicitly the block matrices in the previous definition. Indeed, let y ∈ R n then mat(y) is a n × n Toeplitz-plus Hankel matrix. Let C n denote the n × n orthogonal DCT matrix, we have

C n mat(y)C T n = diag(d),
where d is a vector of eigenvalues can be computed as

d = W -1 (C n mat(y)e 1 ),
where W = diag(C n (:, 1)). Since mat(y)e 1 = (I + Z)vec(y), with Z is the n × n circulant up-shift matrix ( In Matlab, we have: Z = diag(ones(n -1, 1), 1) ). Therefore

d = W -1 (C n (I + Z)vec(y)).
Let L be an operator defined as

L : R n -→ R n y -→ L(y) = W -1 C n (I + Z) M vec(y). Definition 6.16. Let A be an n 1 × n 2 × n 3 tensor. Then L(A) = Â is n 1 × n 2 × n 3 tensor
whose tube fibers âi,j are computed as âi,j = Â(i, j, :) = L(A(i, j, :

)) i = 1, . . . , n 1 , j = 1, . . . , n 2 .
Based on the definition of the n-mode product, we can notice that

L(A) = A × 3 M.
Since mat(A) and mat(B) are block diagonalized using the discrete cosine transform.

Theorem 6.17. The matrix mat(A) is block diagonalizable by the matrix pair C n 3 ⊗ I n 1 and C T n 3 ⊗ I n 2 , with diagonal blocks given by the frontal slices of L(A).

Proof. See [START_REF] Kernfeld | Tensor-tensor products with invertible linear transforms[END_REF] Using the previous theorem, the cosine product can be computed as

mat(A)mat(B) = (C n 3 ⊗ I n 1 )(C n 3 ⊗ I n 1 )mat(A)(C n 3 ⊗ I n 1 )(C n 3 ⊗ I n 1 )mat(B)(C n 3 ⊗ I n 1 )(C n 3 ⊗ I n 1 ) = (C n 3 ⊗ I n 1 )           Â(1) Â (2)
. . .

Â(n 3 )                     B(1) B (2)
. . .

B(n 3 )           (C n 3 ⊗ I n 1 ).
The cosine product can be computed using the following algorithm:

Algorithm 10 Cosine Product: third-order tensor [START_REF] Kernfeld | Tensor-tensor products with invertible linear transforms[END_REF] 1

: input: A ∈ R n 1 ×n 2 ×n 3 , B ∈ R n 2 ×m 2 ×n 3 2: output: C = A * c B ∈ R n 1 ×m 2 ×n 3 3: Â ←-L(A), B ←-L(B)

Generalization of the cosine product to higher-order tensors

In the same way as the t-product, the cosine product can be extended to higher-order tensors. In fact, in case of fourth-order tensors, we have:

Definition 6.21. Let A be an n 1 × n 2 × n 3 × n 4 tensor. Then L(A) = Â is n 1 × n 2 × n 3 ×
n 4 is a tensor whose tube fibers âi,j,k are computed as âi,j,k = Â(i, j, k, :) = L(A(i, j, k, :)) i = 1, . . . , n 1 , j = 1, . . . , n 2 k = 1, . . . , n 3 .

Another way to compute  is to notice that

L(A) = A × 4 M.
The cosine product can be computed using the following algorithm Algorithm 11 Cosine Product: fourth-order tensor

1: input: A ∈ R n 1 ×n 2 ×n 3 ×n 4 , B ∈ R n 2 ×m 2 ×n 3 ×n 4 2: output: C = A * c B ∈ R n 1 ×m 2 ×n 3 ×n 4 3: Â ←-L(A), B ←-L(B)
4: for i = 1, 2, . . . , n 4 do 5:

Ĉ(:, :, :, i) = Â(:, :, :, i) * c B(:, :, :, i)

6: end for 7: C ←-L -1 ( Ĉ) Remark 6.3.1. In the same way as the t-product, the matrix SVD and QR can be generalized using the cosine product for higher-order tensors.

CP decomposition

Let A ∈ R I 1 ×I 2 ו••×I N be an N th -order tensor. The CP decomposition [START_REF] Cichocki | Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation[END_REF][START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Henk | Towards a standardized notation and terminology in multiway analysis[END_REF][START_REF] Tamara | Tensor decompositions and applications[END_REF] of

A is given by A = R ∑ r=1 a (1) r • a (2) r • • • • • a (N) r ,
where a with rank(A) [START_REF] Tamara | Tensor decompositions and applications[END_REF] represent the rank of the tensor A defined as the smallest number of rank-one tensors that generate A as their sum. Unlike matrices, who's the best rank-R approximation is given by the leading R factors of the SVD, the rank of a specific given tensor is hard to define [START_REF] Håstad | Tensor rank is np-complete[END_REF]. In practice, the rank of a tensor is determined numerically by fitting various rank-R CP models. But an interesting property associated with CP decomposition for higher-order tensors is uniqueness under some conditions [START_REF] Richard A Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Joseph | Rank, decomposition, and uniqueness for 3-way and n-way arrays[END_REF].

If we define

A n = î a (n) 1 a (n) 2 • • • a (n) R ó
for n ∈ {1, . . . , N}, the CP decomposition can be symbolically written as

A = A 1 • A 2 • • • • • A N ,
the matrices A n ∈ R I n ×R are called factor matrices. Often, the vectors a (n) r are chosen such that a (n) r = 1. In this case, the CP decomposition is written as

A = R ∑ r=1 λ r a (1) r • a (2) r • • • • • a (N) r ,
where λ r is a scalar that compensates for the magnitudes of vectors a (n) r . Using the n-mode multiplication of a tensor by a matrix, we obtain the following representation:

A = S × 1 A 1 × 2 • • • × N A N ,
where S ∈ R R×Rו••×R , with entries

S i 1 ,••• ,i N =      λ r for i 1 = i 2 = • • • = i N = r, 0 otherwise.
An another representation for the CP decomposition is using slices:

A(:, :, i 3 , • • • , i N ) = A 1 ΛD(A(:, i 3 )) . . . D(A(:, i N ))A T 2 ,

where D(A(:, i)) = diag(A(:, i)), Λ = diag(λ), and λ = [λ 1 , . . . , λ r ].

Application to face recognition

In this section, we are going to present an application to face recognition of the tensor decompositions discussed in this chapter. We first use the CP decomposition and then the QR and SVD decompositions defined via the t-product or the cosine product.

Assuming that we have a collection of images of n p persons, where each image is an m i 1 × m i 2 array with m i 1 m i 2 = n i . Further, we assume that each person has been photographed with n e different facial expressions. Given an image z of an unidentified person in an unknown expression, represented by a vector in R n i or by a matrix of size m i 1 × m i 2 , the aim is to determine which of the n p persons the new image is closest to.

CP decomposition classification 6.5.1.1 Third-order classification

Consider A ∈ R n i ×n e ×n p . The CP decomposition associated to A is given by A = S × i F × e G × p H, where × i , × e , × p are the 1-mode, 2-mode, and 3-mode multiplications, respectively. If we fix the expression e and the person p, the image of a person p in expression e is given by A(:, e, p) = FΛD g e h T p , (6.12)

where g e = G(e, :), h p = H(p, :), and Λ = D(λ). Then (6.12) can be expressed as follows where M e = F 2 F e 1 .

The symbol denotes the Khatri-Rao product [START_REF] Tamara | Tensor decompositions and applications[END_REF]. Then, we classify z ∈ R m 1 ×m 2 as person p in expression ê, where ( ê, p) = argmin

Cosine product and t-product decompositions classification

In the next sections, * denote the t-product or the cosine product. The collection of images is stored as a third-order tensor:

A ∈ R n p ×n e ×n i , or as a fourth-order tensor:

A ∈ R n p ×n e ×m 1 ×m 2 .

Third-order classification

Using the tensor QR factorization, the tensor A ∈ R n p ×n e ×n i can be written as follows

A = Q * R. (6.16)
For a fixed expression e, and using Theorem 6.9, A(:, e, :) can be written as The classification strategy is then analogous to the case of QR factorization, given an new image to be classified z ∈ R n i . The distance of z from the person p is given by α e -U (p, :, :) , (6.24) where α e is solution of the problem min α e Zα e * S * V T (:, e, :) , (6.25) with Z ∈ R 1×1×n i , such that Z(1, 1, :) = z.

Fourth-order classification

In this section, we extend the classification method defined previously to the case of fourth-order tensors. We consider images as matrices of size m 1 × m 2 , instead of vectors. Using the t-QR or the cosine QR factorization, the tensor A ∈ R n p ×n e ×m 1 ×m 2 can be written as follows A = Q * R. end for 10: end for

Numerical results

In this section, we apply the algorithms mentioned above on a database of faces of n p persons in n e different expressions. Each image can be either considered as an m 1 × m 2 matrix, or as a m 1 m 2 vector, so that a database can be represented either by a fourthorder tensor or by a third-order tensor. In this work, we consider the datasets described in the following table we present the success rate of all considered algorithms applied to 'Orl' dataset, for s% = 50%. All percentages correspond to 3 consecutive run, in each run we test 30 images chosen randomly from the test set. We compare the approaches described in this chapter with the approach based on HOSVD described in [START_REF] Eldén | Matrix methods in data mining and pattern recognition[END_REF], and the Principal Component Analysis (PCA-eingfaces) generalized to the tensor case (see [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF]), where the data is represented by a third-order tensor A ∈ R m 1 ×M×m 3 , with M represent the number of images. We point out, that in the next tables, CP 3D and CP 4D refer to three and four dimensional CP decomposition classification respectively. (t-SVD 3D, t-QR 3D

) and (t-SVD 4D, t-QR 4D ) refer to three and four dimensional classification based on t-SVD and t-QR decompositions respectively. In the same way, (dct-SVD 3D, dct-QR 3D ) and (dct-SVD 4D, dct-QR 4D ) refer to three and four dimensional classification based on cosine SVD and cosine QR factorizations respectively. As we can remark, the t-QR 3D and the dct-QR 3D algorithms give a better classification performance than the other tensor based approaches. In term of CPU time, we can observe that the CP 3D algorithm require less CPU time than the other approaches.

Numerical results with compression via truncated HOSVD

In this section, we present the numerical results associated with 3D classifications, with compression via truncated HOSVD. Definition 6.22. Let A ∈ R I 1 ×I 2 ו••×I N be an N th -order tensor. The Tucker decomposition (often referred to as the higher-order SVD (HOSVD) [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF][START_REF] Ledyard | Some mathematical notes on three-mode factor analysis[END_REF] ) of A is defined as:

A = V × 1 U 1 × 2 U 2 × 3 • • • × N U N ,
where V ∈ R R 1 ×R 2 ו••×R N is the core tensor, U n ∈ R I n ×R n are factor matrices, and

(R 1 , . . . , R N ) is the multi-linear rank of the tensor A, where R n = rank(A (n) ). In case of HOSVD decomposition, the factor matrices U n , for n = 1, . . . , N are orthonormal.

For computing the orthonormal factors U n for n = 1, . . . , N, and the tensor V. We compute the SVD associated with each n-mode matrix of the tensor A:

A (n) = U n Σ n V T n ,

Conclusion

In this chapter, we have proposed approaches for face recognition. Based on applying the tensor CP decomposition or QR/SVD decompositions defined using the t-product or the cosine product to a dataset of images. Our approaches are applied to a database represented by third or fourth-order tensors. In the numerical examples, we compare our approaches with PCA method based on tensor format and HOSVD based approach. In 3D classifications, we obtain comparable or better recognition results with the compressed versions of the algorithms.

Figure 1 . 1 :

 11 Figure 1.1: 0D, 1D, 2D and 3D tensors.

Figure 1 . 2 :

 12 Figure 1.2: A color image representation as a third-order tensor .

Figure 1 . 3 :

 13 Figure 1.3: gray-scale and color videos representation as a third-order and fourthorder tensors.

Figure 1 . 4 :

 14 Figure 1.4: A multi-modal MRI dataset of a patient.

Figure 1 . 5 :

 15 Figure 1.5: Tensor Representation of a database.

Figure 1 . 6 :

 16 Figure 1.6: Gray-scale image representation as a matrix.

Figure 1 . 7 :

 17 Figure 1.7: Degradation Model

Figure 1 . 8 :

 18 Figure 1.8: Example of biometric modalities

  and U contains the eigenvectors of C. Given an unknown face Y ∈ R m×n which we want to classify. The first step is to subtract the face from the mean, ȳ = yx, with y = vec(Y). Then we project the normalized vector ȳ into the eigenspace U k ȳ. At the end, we classify the person in the image Y as the person X l in the training set, with l satisfy l = argmin i αα i

  Les tenseurs sont des tableaux multidimensionnels, et ils généralisent les matrices à des dimensions supérieures. L'ordre d'un tenseur est le nombre de dimensions, également appelées voies ou modes. Les scalaires peuvent donc être interprétés comme des tenseurs d'ordre zéro. les vecteurs comme des tenseurs de premier ordre et les matrices comme des tenseurs de second ordre. Les tenseurs d'ordre trois ou plus sont appelés tenseurs d'ordre supérieur. La figure suivante montre comment on peut passer des scalaires aux tenseurs.

Figure 2 . 1 :

 21 Figure 2.1: Tenseurs 0D, 1D, 2D et 3D

  et d'autres domaines encore. Tout comme les matrices sont utilisées pour représenter les transformations linéaires, les tenseurs peuvent être utilisés pour représenter des types de transformations plus généraux. Les tenseurs sont également un moyen naturel de représenter des données multidimensionnelles, comme les images, où les images en niveaux de gris peuvent être considérées comme des tenseurs de deuxième ordre, les images en couleurs (images RVB) sont présentées comme des tenseurs de troisième ordre, et une vidéo composée d'images en couleurs est un tenseur de quatrième ordre.

Figure 2 . 2 :

 22 Figure 2.2: Une représentation d'une image couleur sous forme de tenseur du troisième ordre .

Figure 2 . 3 : 4 .

 234 Figure 2.3: représentation des vidéos en gris et en couleur sous forme de tenseurs d'ordre 3 et 4

Figure 2 . 4 :

 24 Figure 2.4: Ensemble de données IRM multimodales d'un patient .

Figure 2 . 5 :Chapter 2 .

 252 Figure 2.5: Représentation tensorielle d'une base de données .

Figure 2 . 6 :

 26 Figure 2.6: Représentation d'une image à échelle de gris sous forme de matrice .

  etc. En imagerie médicale, un certain nombre de techniques et de dispositifs d'imagerie sont inventés. Comme dans tout autre système d'imagerie, les dispositifs d'acquisition d'images médicales introduisent également des dégradations dans les images. Les méthodes de restauration d'images jouent un rôle important dans l'amélioration de la qualité des images obtenues à partir de dispositifs d'imagerie médicale.

  et U contient les vecteurs propres de C. Étant donné un visage inconnu Y ∈ R m×n que l'on veut classer. La première étape consiste à soustraire le visage de la moyenne, ȳ = yx, Ensuite, nous projetons le vecteur normalisé ȳ dans l'espace propre U k ȳ. Enfin, nous classons la personne dans l'image Y comme la personne X l dans l'ensemble de données, avec l satisfaisant à l = argmin i αα i avec α est la coordonnée de ȳ par rapport aux k composantes principales. Récemment, de multiples extensions de la méthode des faces propres ont été développées, telles que l'approche t-SVD [40] définie via le t-produit. Cette approche diffère de la méthode traditionnelle des faces propres car les données sont représentées sous la forme d'un tenseur d'ordre 3 de taille m × l × n, ce qui signifie que dans cette approche, les images utilisées ne sont pas vectorisées et, comme dans la méthode traditionnelle, le tenseur de covariance et son associé ne sont pas calculés. Au lieu de cela, une généralisation de la SVD matricielle via le t-produit est utilisée. Les performances des algorithmes de reconnaissance de visages tels que les eigenfaces et Fisherfaces (approche basée sur LDA) sont bonnes lorsque la seule variable qui compte dans la formation de l'image est l'identité de la personne, ce qui signifie que les visages sont capturés dans des conditions contrôlées. En réalité, il existe plusieurs facteurs qui peuvent affecter l'image d'une personne donnée, comme l'illumination, qui ont une grande influence sur l'apparence du visage dans l'image. Les conditions d'illumination sont inévitables dans le monde réel, en particulier lorsque les vues sont collectées à différents moments. D'autres facteurs qui peuvent affecter l'image d'un visage sont l'angle de vue, l'expression... Un autre problème lié à l'utilisation d'algorithmes de reconnaissance des visages tels que les visages propres est que dans différents domaines, une représentation naturelle des images est un tenseur du troisième ordre plutôt qu'une simple matrice d'images vectorisées. De multiples approches ont été développées pour surmonter ces problèmes. Parmi les solutions proposées, l'utilisation d'une représentation tensorielle des données. Par exemple, nous pouvons représenter la base de données de visages de sujets photographiés dans différentes poses sous différents éclairages et différentes expressions du visage comme un tenseur d'ordre 5. Par exemple, nous pouvons représenter la base de données de visages de sujets photographiés dans différentes poses sous différents éclairages et différentes expressions du visage comme un tenseur d'ordre 5 lorsque les images sont représentées comme des matrices, ou nous pouvons représenter la base de données de visages comme un tenseur d'ordre 6 lorsque les images sont représentées comme des tenseurs d'ordre 3. Afin de manipuler les données stockées sous forme de tenseur, différents outils d'algèbre multi-linéaire sont utilisés. Par exemple, dans

  comme suit : Le chapitre 3 est le chapitre des préliminaires dans lequel nous présentons les bases de l'algèbre multi-linéaire qui seront utiles dans le développement des techniques présentées dans cette thèse et pour aider les lecteurs à comprendre les concepts multi-linéaires. Nous introduisons certaines opérations tensorielles et matricielles de base ainsi que leurs propriétés, et nous donnons également une brève introduction à certaines des décompositions tensorielles de base utilisées dans les chapitres suivants. Dans le chapitre 4, nous proposons une nouvelle approche de la restauration d'images et de vidéos. Cette approche construit un modèle de dégradation basé sur une représentation tensorielle, où une image couleur est représentée par un tenseur de troisième ordre, et une vidéo composée d'images couleur est un tenseur de quatrième ordre.

and defined as A ⊗ B = a 11 B

 11 . . . . . . 1J B . . . a I1 B . . . . . . a I J B Theorem 3.3.

Theorem 3 .

 3 10. (QR decomposition)Any matrix A ∈ R m×n , m ≥ n, can be transformed to upper triangular form by an orthogonal matrix. The transformation is equivalent to a decompositionA = Q Ö R 0 èwhere Q ∈ R m×m is orthogonal and R ∈ R n×n is upper triangular. If the columns of A are linearly independent, then R is nonsingular.

Figure 3 . 1 :

 31 Figure 3.1: Symbolic illustration of the QR decomposition.

Figure 3 . 2 :

 32 Figure 3.2: Special forms of third-order tensors: identity tensor.

Figure 3 . 3 :

 33 Figure 3.3: Fibers of a third-order tensor.

Figure 3 . 4 :Definition 3 . 11 .

 34311 Figure 3.4: Slices of a third-order tensor.

Figure 3 . 5 :

 35 Figure 3.5: CP decomposition of a third-order tensor.

Figure 3 . 6 :

 36 Figure 3.6: Tucker representation of a third-order tensor.

7 :

 7 end procedureWe point out that a rank (r 1 , • • • , r N ) approximation with r n ≤ R n for n = 1, . . . , N can be obtained simply by restricting the factor matrices U n to the first r n columns (truncated SVD) for n = 1, . . . , N, and by restricting the core tensor S.In this chapter we propose a new approach to image and video restoration. This approach constructs a degradation model based on a tensor representation, where a color image is represented by a third-order tensor, and a video composed of color images is a fourth-order tensor. Applying tensor CP decomposition to our original problem leads to three subproblems. To solve those subproblems, we apply global LSQR algorithm, and a new algorithm based on Golub Kahan bidiagonalization. Some numerical tests are presented to show the effectiveness of the proposed methods.

4

 4 Solving the problem EXF T = HThe minimization problems given in (4.23) have the following form min X H -EXF T , (

  and B k is a lower bidiagonal matrix. The coefficients β 2 , . . . , β k and α 1 , . . . , α k determined by Algorithm 3 define the lower bidiagonal matrix

Algorithm 3

 3 orthonormal matrices, β1 = ||h 2 || and C k is a lower bidiagonal matrix. Golub-Kahan bidiagonalization of the matrix A.

0 = 0 .Proposition 4 . 4 .

 044 The residual norm ||R k || 2 can be computed as follows

Theorem 4 . 6 .

 46 Let R k = H -EX k F T . Then minimizing the residual norm ||R k || is equivalent to minimizing (β 1 e 1 )( β1 e 1 T ) -B k y k (1) y k (2) T C k T .

Figure 4 . 1 :

 41 Figure 4.1: Example 2. (a) Exact image, (b) Degraded image(motion blur + noise ν = 10 -3 ), Restored images by (c) Algorithm 2(Algo 3)

  video of size 200 × 200 × 20 that have been contaminated by both blur (Gaussian blur of size 9 × 9 with deviation σ = 4) and additive noise with noise level equal to 10 -2 . The following figure shows the third frame associated with the true, the blurred, and the restored video.

Figure 4 . 2 :

 42 Figure 4.2: Example 3: Frame no. 3: (a) Original frame, (b) Blurred and noisy frame, (c) Restored frame by Algorithm2(Algo 3)

( 5 . 2 )

 52 for higher orders, since we are dealing with matrices instead of tensors. In fact, when the dimension increases the problem becomes harder (computation and memory requirement), since the data size of a tensor increases exponentially with the increase of the dimensionality of the tensor. As consequence tensor computations become much expensive. For example the n-mode product given in Definition 3.14 has a computational complexity of O(J N ∏ i=1 i =n I i ). In addition, by writing the approximate solution in CP or HOSVD decomposition format, we reduce the required memory. For instance, the CP decomposition transforms the storage complexity of an I N tensors into O(NRI),

1 )

 1 = [A(:, 1, :), A(:, 2, :), . . . , A(:, I 2 , :)] A (2) = A(:, :, 1) T , A(:, :, 2) T , . . . , A(:, :, I 3 ) T A (3) = [A(1, :, :), A(2, :, :), . . . , A(I 1 , :, :)] .

1 6 :

 6 for doj = 1, 2, 3 . . . , k 7:

j is small enough then stop 24 : end for Remark 5 . 4 . 1 .

 24541 In the line 23 of Algorithm 5, we compute the upper bound of the residual norm ||R k || given in Theorem 5.12 where n

23 ) 5 . 5 . 1 .

 23551 RemarkIn case where each channel of the p channels are not affected by the same blur. The blurring model associated is of the form

r for i = 1 , 2 ,Lemma 5 . 14 .

 12514 • • • , p are circulant matrices. Using the following well known propriety associated to the diagonalization of circulant matrices. Let C ∈ R n×n be a circulant matrix, Then C is diagonalized matrix. More precisely, C = F * ΛF,

Figure 5 . 1 :

 51 Figure 5.1: Example 1. (a) Exact image, (b) blurred image.

Figure 5 . 2 :

 52 Figure 5.2: Example 1. restored images for R=m=150 (a) approach based on CP decomposition, (b) approach based on HOSVD decomposition.

Figure 5 . 3 :

 53 Figure 5.3: Example 2. First row: original frames. Second row: Blurred frames. Third row: Restored frames.

  a block circulant matrix of size n 1 n 3 × n 2 n 3 , and the command unfold(A) takes the tensor A into a block n 1 n 3 × n 2 matrix. The operation fold takes unfold(A) back to the tensor form: fold(unfold(A)) = A.

  j, :)) ⊗ e i e T j are equal. This proves (6.1).

Definition 6 .

 6 15 ([50]). The cosine product * c between two tensors A ∈ R n 1 ×n 2 ×n 3 and B ∈ R n 2 ×m 2 ×n 3 is a tensor of size n 1 × m 2 × n 3 given by C = A * c B = ten(mat(A)mat(B)).(6.11)

4 :6: end for 7 :

 47 for i = 1, 2, . . . , n 3 do 5:Ĉ(:, :, i) = Â(:, :, i) B(:, :, i) C ←-L -1 ( Ĉ) Properties 6.3.1.[START_REF] Kernfeld | Tensor-tensor products with invertible linear transforms[END_REF] 1. The cosine product is associative: Let A, B and C be third-order tensors of appropriate size, thenA * c (B * c C) = (A * c B) * c C.

2 .Proposition 6 . 18 .Proposition 6 . 19 .Lemma 6 . 20 .

 2618619620 Let A, B be third-order tensors such that A * c B and B T * c A T are defined, then(A * c B) T = B T * c A T . [50] Let Ĵ ∈ R n 1 ×n 1 ×n 3 be so that Ĵ (:, :, i) = I n 1 ×n 1 . Then L -1 ( Ĵ ) = J is the identity element under L Let A ∈ R n 1 ×n 2 ×n 3 , B ∈ R n 2 ×m 2 ×n 3 . Then (A * c B)(i, j, ; ) = A(i, :, ; ) * c B(:, j, ; ) i = 1, . . . , n 1 , j = 1, . . . , m 2 .Proof. We have J (i, :, :) * c A = A(i, :, :) f or i = 1, . . . , n 1 and B * c J (:, j, :) = B(:, j, :) f or j = 1, . . . , m 2 . Then (A * c B)(i, j, ; ) = J (i, :, :) * c (A * c B) * c J (:, j, :) = (J (i, :, :) * c A) * c (B * c J (:, j, :)) = A(i, :, ; ) * c B(:, j, ; ). Let A ∈ R n 1 ×n 2 ×n 3 , B ∈ R n 2 ×m 2 ×n 3 and C ∈ R m 2 ×p 2 ×n 3 , we have (A * c B * c C)(i, j, :) = A(i, :, :) * c B * c C(:, j, :).

r

  are vectors of size I k with 1 ≤ k ≤ N and R is a positive integer. A CP decomposition of a tensor A is called an exact CP decomposition if R = rank(A),

A

  (:, e, p) = C e h T p ,where C e = FΛD g e . The columns of C e are basis vectors for the expression e and the row p of H, i.e., h p , holds the coordinates of the image of person p in this basis. Furthermore, the same h p holds the coordinates of the images of person p in all expression bases. Then, if the image z that we want classify is an image of person p in expression e, then the coordinates of z in that basis are equal to h p . Thus, we can classify z by computing its coordinates in all the expression bases and testing, for each expression, whether the coordinates of z coincide or nearly coincide with the elements of any row of H. The coordinates of z in the expression basis e can be obtained by solving the least squares problem min α e C e α ez for e = 1, . . . , n e . Then, for each e = 1, . . . , n e and for each p = 1, . . . , n p , we compute d(e, p) = α eh T p . We classify z as person p in expression ê, where ( ê, p) = argmin e,p d(e, p).

6. 5 . 12 1: 3 : 5 : 6 : 9 :F 1

 51235691 Application to face recognition 123 Algorithm input: Collection of images A ∈ R n i ×n e ×n p , z ∈ R n i image to classify 2: Decompose A as A = S × i F × e G × p H for e = 1; 2, . . . , n e do ez with C e = FΛD g e for p = 1; 2, . . . , n p do Compute d(e, p) = α e -Classify z as person p in expression ê, where ( ê, p) = argmin e,p d(e, p)6.5.1.2 Fourth-order classificationNow, we consider images as matrices of size m 1 × m 2 , which mean, the collection of n p persons in n e expressions is stored as fourth-order tensor A ∈ R m 1 ×m 2 ×n e ×n p . The CP decomposition of the tensor A is expressed as followsA = S × i 1 F 1 × i 2 F 2 × e G × p H,where S ∈ R R×R×R×R , × i 1 , × i 2 , × e , × p are the 1-mode, 2-mode, 3-mode, and 4-mode multiplications, respectively. Thus, the image of a person p in expression e in the database is given by A(:, :, e, p) = F 1 ΛD g e D h p F T 2 , (6.13)where g e = G(e, :),h p = H(p, :). If we set C = L × i 1 F 1 × i 2 F 2 × e G, (6.13) becomes A(:, :, e, p) = C e × p h p , where C e = C(:, : e, :). So to classify the person z ∈ R m 1 ×m 2 , we compute for each e = 1, . . . , n e and for each p = 1, . . . , n p d(e, p) = α eh T p , where α e is solution of the minimization problem min α e C e × p α ez for e = 1, . . . , n e . Equivalently, we can set α e as solution of the following minimization problem min α e D g e X α e F T 2z where X α e = D (α e ) . (6.14) If we set F e 1 = F 1 D g e , the minimization problem (6.14) becomes min α e F e 1 X α e F T 2z where X α e = D (α e ) . (6.15) Using the vec operator, the minimization problem (6.15) is equivalent to the following minimization min α e M e α evec(z)

1 : 4 : 5 : 6 : 9 :

 14569 input: Collection of images A ∈ R m 1 ×m 2 ×n e ×n p , z ∈ R m 1 ×m 2 image to classify 2: Decompose A as A = S × i 1 F 1 × i 2 F 2 × e G × pH 3: for e = 1, 2, . . . , n e do Solve min α e M e α evec(z) where M e = F 2 F e 1 , and F e 1 = F 1 D g e for p = 1, 2, . . . , n p do Compute d(e, p) = α e -Classify z as person p in expression ê, where ( ê, p) = argmin e,p d(e, p)

A 14 1: 4 : 6 :

 1446 (:, e, :) = Q * R(:, e, :), f or e = 1, . . . , n e . (6.17) Thus, the image of a person p in expression e is expressed by A(p, e, :) = Q(p, :, :) * R(:, e, :), f or e = 1, . . . , n e , p = 1, . . . , n p . (6.18) Given an image z ∈ R n i of an unknown person in an unknown expression, that we want to classify. The coordinates of z in the expression basis 'e' can be found by solving the following least squares problem min α e Zα e * R(:, e, :) , (6.19) with Z ∈ R 1×1×n i , such that Z(1, 1, :) = z. Then, for each e = 1, . . . , n e and for each p = 1, . . . , n p , we compute: α e -Q(p, :, :). . (6.20) We can then recognize the person in the image by attributing the label corresponding to the closest match. Algorithm input: Collection of images A ∈ R n p ×n e ×n i , z image to classify 2: Decompose A as A = Q * R 3: Z(1, 1, :) = z(:) for e = 1, 2, . . . , n e do e * R(:, e, :) for p = 1, 2, . . . , n p do 7: if then α e -Q(p, :, :) < tol, then classify z as the person p In the same way. Applying the tensor SVD decomposition to the third-order tensor A, that contain the images data leads to the following relation A = U * S * V T . (6.21) Fixing a particular value of the expression parameter, i.e j = e corresponds to using V(:, e, :). In other word, we have A = U * S * V T (:, e, :). (6.22) Thus, the image of a person p in the expression e is defined by A(:, e, p) = U (p, :, :) * S * V T (:, e, :). (6.23)

(6. 26 ) 15 1: 4 : 6 :

 261546 For a fixed expression e, A(:, e, :, :) can be written as A(:, e, :, :) = Q * R(:, e, :, :), f or e = 1, . . . , n e . (6.27) Thus, the image of a person p in the expression e is expressed by A(p, e, :, :) = Q(p, :, :, :) * R(:, e, :, :), f or e = 1, . . . , n e , p = 1, . . . , n p . (6.28) Given an image z ∈ R m 1 ×m 2 of an unknown person in an unknown expression, that we want to classify. The coordinates of z in the expression basis can be found by solving the following least squares problem min α e Zα e * R(:, e, :, :) , (6.29) with Z ∈ R 1×1×m 1 ×m 2 , such that Z(1, 1, :, :) = z. Then, for each e = 1, . . . , n e and for each p = 1, . . . , n p , we compute: α e -Q(p, :, :, :) . (6.30) We can then recognize the person in the image by assigning the label corresponding to the closest match. Algorithm input: Collection of images A ∈ R n p ×n e ×m 1 ×m 2 , z image to classify 2: Decompose A as A = Q * R 3: Z(1, 1, :, :) = z for e = 1, 2, . . . , n e do e * R(:, e, :, :) for p = 1, 2, . . . , n p do 7:if then α e -Q(p, :, :, :) < tol, then classify z as the person p

Figure 6 . 1 :

 61 Figure 6.1: Sample face images of three individuals from the Face96 dataset

Figure 6 . 2 :

 62 Figure 6.2: Sample face images of three individuals from the Orl dataset

Algorithm 16 1 : 7 : 8 :Figure 6 . 3 :

 17863 Figure 6.3: Time of decomposition, and time of recognition of t-QR 3D, dct-QR 3D for different values of the truncation parameter π tested on Orl dataset, with s% = 50%.

Figure 2.7: Modèle de dégradation

  Si les images observées ont m × n pixels, alors K r et K c sont des matrices de taille n × n

	où f est l'image réelle, g est l'image observée, et η est un bruit additif. La fonction l'intégration avec une règle de quadrature. K est une matrice qui représente l'opération peut être représentée comme produit de Kronecker de deux matrices K r et K c ,
	noyau k modélise l'opération de flou et est appelée fonction d'étalement du point (PSF). de flou, et elle peut être construite via la fonction d'étalement du point (PSF) et les
	Dans de nombreuses situations, le flou est supposé être spatialement invariant, ce qui conditions aux limites [39] puisque les bords sont des structures importantes de l'image K = K r ⊗ K c
	signifie que l'opérateur de noyau k satisfait à k(x, y; s, t) = k(x -s, y -t). Dans ce cas, le réelle et qu'ils doivent être préservés pendant la restauration de l'image. Dans le cas
	modèle de dégradation donné dans (2.1) est formulé comme une opération de convo-d'un flou spatialement invariant, la structure de la matrice K dépend des conditions
	lution : aux limites imposées. Par exemple, si nous imposons des conditions aux limites nulles et m × m respectivement. Cette décomposition de Kronecker de la matrice K réduit la
	g(x, y) = (k f )(x, y) + η(x, y) en supposant qu'en dehors des frontières de l'image, tout est noir, la matrice K est un dimension du problème (2.2) de mn × mn à m × n. Dans ce cas, le modèle de dégrada-
	bloc Toeplitz avec des blocs Toeplitz (BTTB). Cette hypothèse est utile lorsqu'on traite tion peut être formulé sous cette forme,
	des images astronomiques, car la plupart du temps, il est possible de supposer que les
	bords extérieurs de l'image sont noirs.	K c FK T r = G,
	Si l'on suppose que les images sont des matrices de taille m × n. Alors la matrice de où G, F ∈ R m×n , et K T r est le transposé de K r . Dans le cas non-séparable, on peut
	flou K est de taille N × N, avec N = mn le nombre de pixels de l'image, et g, f , et η approximer la matrice K en résolvant le problème de l'approximation du produit de
	sont des vecteurs de taille N. Le problème de la restauration d'image est maintenant Kronecker (KPA) [76].
	réduit au problème de la résolution de l'équation (2.2). Il existe un grand nombre
	d'approches fournissant des solutions au problème de la restauration d'images [1, 2, 8].
	Une approche classique pour résoudre l'équation (2.2) consiste à calculer sa solution en
	Ainsi, le problème de restauration d'image est souvent appelé déconvolution. Le sym-résolvant ce problème moindres carrés
	bole " " indique une convolution bidimensionnelle. L'objectif de la restauration est
	d'obtenir une estimation f de l'image originale f . Malheureusement, le modèle (2.1) n'est pas très utile pour la restauration d'images en min f K f -g (2.3)
	raison de la complexité impliquée par la possibilité d'avoir une PSF différente k(x, y; s, t) Mais de nombreux aspects rendent la résolution du problème (2.3) très difficile. En
	à chaque coordonnée (s, t) de l'image, et il n'est pas réaliste de supposer que l'on effet, les dimensions de la matrice K peuvent être extrêmement grandes puisque si
	puisse estimer une PSF différente pour chaque emplacement dans l'image. De plus, l'image observée a une taille de 512 × 512, alors la matrice K a une taille de 262144 ×
	nous ne disposons pas d'une définition précise de la fonction pour g car l'image ob-262144. Ainsi, le problème (2.3) est à grande échelle, ce qui rend sa résolution coû-
	servée est enregistrée numériquement, et n'est donc connue que pour des valeurs dis-teuse en termes de calcul puisque ces calculs impliquent généralement des produits
	crètes. De plus, dans de nombreux cas, il est nécessaire d'estimer k à partir de données vectoriels matriciels. Une technique efficace pour surmonter ce problème consiste à
	mesurées. Par conséquent, il est naturel de considérer le problème de restauration exploiter la structure de la matrice K. En particulier, lorsque le flou est séparable, ce
	d'images numériques suivant qui signifie que les composantes horizontale et verticale du flou peuvent être séparées,
	, y; s, t) f (s, t)dsdt + η(x, y) g = K f + η dans ce cas, le noyau du flou satisfait k(x, y; s, t) = k 1 (x, s)k 2 (y, t), et donc la matrice K (2.1)	(2.2)
	qui est obtenue à partir de l'équation (2.1) en discrétisant les fonctions et en approchant

  où A T est le transposé de A, C est symétrique, semi-définie et diagonalisable orthogonalement, ce qui peut s'écrire comme suit : C = UDU T , où D est une matrice diagonale qui contient les valeurs propres de C, U est une matrice orthogonale, où chaque colonne u i de U est un vecteur propre de C. Les vecteurs propres de la matrice C sont appelés les composantes principales. Après avoir calculé les vecteurs propres et les valeurs propres de la matrice de covariance C, nous choisissons les k vecteurs propres de C correspondant aux k plus grandes valeurs propres, avec k < l. En utilisant les k premiers vecteurs propres, nous pouvons approximer chacun des vecteurs de visage normalisés a i comme une combinaison linéaire de (u 1 , . . . , u k ).

  111 x 121 x 131 x 112 x 122 x 132 x 211 x 221 x 231 x 212 x 222 x 232 x 311 x 321 x 331 x 312 x 322 x 332 x 411 x 421 x 431 x 412 x 422 x 432 x 111 x 211 x 311 x 411 x 112 x 212 x 312 x 412 x 121 x 221 x 321 x 421 x 122 x 222 x 322 x 422 x 131 x 231 x 331 x 431 x 132 x 232 x 332 x 432 x 111 x 211 x 311 x 411 . . . x 131 x 231 x 331 x 431 x 112 x 212 x 312 x 412 . . . x 132 x 232 x 332 x 432

	í
	X (2) =
	and
	X (3) =
	
	        

à à í Definition 3.13. The vectorization of matrix Y ∈ R I×T is defined as

  In line 26 of Algorithm 4, we compute an upper bound of the residual norm ||R k || using Proposition 4.4, where n(1) i and n

	(2) i denote || f (1) i || 2 and || f (2) i || 2

2)2 i+1 is small enough then stop 27: end for Remark 4.4.1.

Table 4 . 1 :

 41 Results for Example 1

	Size	Method	Iter Relative error CPU-time(sec)
		Global LSQR	30	5.5 × 10 -2	0.1
	128 × 128	Algorithm 4	20	2.45 × 10 -2	0.09
		Global LSQR	40	2.74 × 10 -2	0.78
	400 × 318	Algorithm 4	30	1.0 × 10 -2	0.71
	The next table shows the results obtained from the restoration of color images with
	Algorithm 2, using Algorithm 4 to solve the three subproblems, or using global LSQR
	denoted respectively by Algorithm 2(Algo 3), and Algorithm 2(G-LSQR), with s = 250.
	We point out that in the next table, and in all the tables below, "Iter" stands for the
	number of iterations and has the following form Iter1(iter2)", where "Iter1" indicates
	the number of iterations of Algorithm 2, and "iter2" indicates the number of iterations
	of Algorithm 4 or global LSQR.

Table 4 . 2 :

 42 Results for Example 1

	Size	Method	Iter Relative error CPU-time(sec)
		Algorithm 2(G-LSQR) 4(30)	4.32 × 10 -2	4.79
	227 × 303 × 3	Algorithm 2(Algo 3) 5(20)	3.1 × 10 -2	3.64
		Algorithm 2(G-LSQR) 4(30)	6.07 × 10 -2	12.2
	384 × 512 × 3	Algorithm 2(Algo 3) 5(20)	5.71 × 10 -2	9.72

Table 4 . 3 :

 43 Comparison of image restoration on image 'peppers' blurred by Gaussian blurBased on the testes reported in Table4.3 and many unstated tests, we remark that our proposed algorithm works effectively for image restoration problems both in terms of the SNR and the relative error. In terms of the CPU-time, we note that T-global Golub

	Noise level	Method	Iter	SNR Relative error CPU-time(sec)
		T-global GMRES	10(m=10) 21.6	4.91 × 10 -2	20.89
		T-global Golub Kahan	15	20.44	5.55 × 10 -2	8.21
	0.001	Algorithm 2(Algo 3)	3(25)	21.86	4.81 × 10 -2	6.3
		T-global GMRES	10(m=10) 19.86	6.20 × 10 -2	21.89
		T-global Golub Kahan	6	18.17	7.37 × 10 -2	8.05
	0.005	Algorithm 2(Algo 3)	4(20)	20.45	5.66 × 10 -2	7.66
		T-global GMRES	10(m=10) 19.61	6.24 × 10 -2	21.8
		T-global Golub Kahan	8	18.2	7.34 × 10 -2	4
	0.01	Algorithm 2(Algo 3)	4(25)	20.47	5.65 × 10 -2	10.7

Table 4 .4: Results for Example 3

 4

	s	Iter Relative error CPU-time(sec)
	100 5(40)	6.78 × 10 -2	10.9
	200 5(40)	5.99 × 10 -2	14.82
	300 4(40)	5.77 × 10 -2	11.74

Table 5 . 1 :

 51 Numerical results for Example 1.

	N	n	p	R Iter	R k	X k -X	CPU-time(sec)
		200 100 5	21	3.14 × 10 -8 7.66 × 10 -13	0.18
	3	400 300 10	30	1.95 × 10 -6 2.46 × 10 -12	0.78
		500 400 10	30	8.40 × 10 -6 4.37 × 10 -12	5.16
		50	50 10	23	2.53 × 10 -7 8.16 × 10 -12	0.11
	4	100 50 20	30	1.92 × 10 -7 3.31 × 10 -12	0.17

Table 5 . 2 :

 52 Comparison of algorithm 5 and GLS-BTF.

	Method	N	n	Iter	R k	error	X k -X	CPU-time(sec)
			100	27	3.42 × 10 -10 1.34 × 10 -9	3 × 10 -11	0.21
		3	400	30	4.72 × 10 -8 7.39 × 10 -8 1.39 × 10 -9	17.94
	Algorithm5-HOSVD		50	27	6.12 × 10 -10 2.81 × 10 -9 4.60 × 10 -11	2.01
		4	100	27	6.34 × 10 -8 1.33 × 10 -7 6.55 × 10 -10	35
			100 131 4.06 × 10 -6		5.41 × 10 -6	11.71
		3
	GLS-BTF [71]		400	_	_		_	_
			50	160	8.5 × 10 -5		1.18 × 10 -4	106.18
		4
			100	_	_		_	_

Table 5 . 3 :

 53 Numerical results for Example 2 with n = 1000, 10.000.

	N	n	R Iter	error	CPU-time(sec)
		1000	10	25	1.97 × 10 -9	0.22
	3	10,000 10	25	3.45 × 10 -7	10.27
		1000	10	30	2.07 × 10 -10	0.34
	4	10,000 5	30	3.03 × 10 -8	16.13

Table 5 . 4 :

 54 Results for Example 1.

R =

m Method Iter SNR Relative error CPU-time(sec)

		(a)				(b)
		Algorithm 5-CP	40	19.51	6.32 × 10 -2	1.96
	100	Algorithm 5-HOSVD	40	22.44	4.51 × 10 -2	0.63
		Algorithm 5-CP	20	21.77	4.87 × 10 -2	2.75
	150	Algorithm 5-HOSVD	20	24.23	3.67 × 10 -2	0.54

Table 5 . 5 :

 55 Comparison of the performance of the two approaches for videos restoration of different sizes.

	Size	R = m	Method	Relative error CPU-time(sec)
		150	Algorithm 5-CP	2.79 × 10 -2	1.15
	200 × 200 × 30	150	Algorithm 5-HOSVD	3.12 × 10 -2	0.25
		200	Algorithm 5-CP	7.79 × 10 -2	4.56
	360 × 640 × 30	200	Algorithm 5-HOSVD	3.13 × 10 -2	1.19
	360 × 640 × 100	200	Algorithm 5-HOSVD	5.75 × 10 -2	3.49
	200 × 200 × 3 × 20	150	Algorithm 5-HOSVD	3.13 × 10 -2	0.68
	360 × 640 × 3 × 30	200	Algorithm 5-HOSVD	3.17 × 10 -2	6.82

Table 6 . 1 :

 61 Size, number of expressions and persons of each databases.

	Dataset Pixel(matrix format) Pixel(vector format) n p	n e
	Face96	196 × 196	38416	119 19
	Orl	92 × 112	10304	40 10

Table 6 . 2 :

 62 Comparison results -Orl dataset

	Method	Success rate Time of decomposition	Time of recognition
	CP 3D (R=40)	87,78	0.085	0.084
	CP 4D (R=35)	88,89	1.53	0.11
	t-QR 3D	94.46	0.52	1.88
	t-QR 4D	93.33	1.28	3.18
	dct-QR 3D	94.46	0.608	2,32
	dct-QR 4D	91.13	1.68	4,32
	HOSVD 3D	86.66	0.24	0.096
	PCA-eignfaces	85.56	0.35	0.03

Table 6 . 3 :

 63 Success rate associated to 3D classification algorithms tested on Face96 dataset, with s% = 25%.

Method Success rate Time of decomposition Time of recognition

	CP 3D	92.22	0.29	0.35
	HOSVD 3D	92.22	1.762	1.11
	t-QR 3D	92.22	0.13	0.36
	t-SVD 3D	93.33	0.18	0.37
	dct-QR 3D	92.22	0.14	0,38
	dct-SVD 3D	95.56	0.15	0,38
	PCA-eignfaces	90.56	2.82	0.15

LIST OF FIGURES

• 9th International Conference (online) on Applied Analysis and Mathematical Modeling, held: June 11-13, 2021.

• Scientific and Technical Days (JST 2021), which took place: March 22-25, 2021.

• International Hybrid Conference on Numerical Analysis and Optimization Days (JANO '13), which took place: 22-24 February 2021.

• Inter-laboratory research day 2020 JRIL, December 24, 2020.

• Doctoral Days of the Ecole Normale Supérieure of Rabat, December 24 and 25, 2019.

Definition 3.17. (Outer product)

The outer product of the tensors Y ∈ R

Z i 1 ,i 2 ,...,i N ,j 1 ,j 2 ,...,j M = Y i 1 ,i 2 ,...,i N X j 1 ,j 2 ,...,j M .

As special cases, the outer product of two vectors a ∈ R I and b ∈ R J yields a rank-one matrix

and the outer product of three vectors: a ∈ R I , b ∈ R J , and c ∈ R Q yields a third-order rank-one tensor:

Example of application to image and video restoration

Many applications require the solution of the problem of the form (5.2). The problem (5.2) may appear in color image and video restoration. We recall that the blur model associated to p-channel images (color images and gray-scale videos) is given by

where X , G ∈ R m×n×p denote the original image and the degraded image respectively.

For more details about the degradation model associated to p-channel images, see chapter 4 section 4.2. This model is constructed under the assumption that the blur is the same in all channels.

In case of color videos, video restoration is the problem of restoring a sequence of k color images (frames). Each frame is represented by a third-order tensor of size m × n × 3 which mean, a color video can be represented as a fourth-order tensor of size m × n × 3 × k. Under the same assumptions used to construct the blurring model associated to p-channel images, the degradation model associated color videos is given by Proposition 5.13. The blur model associated with color videos is given by the tensor equation

where X , G ∈ R m×n×3×k denote the original and the degraded color video.

Proof. Let G(:, :, i, j) for i = 1, 2, 3 and j = 1, 2, • • • , k denotes the gray scale image that constitute the channel of each frame of the blurred video. Assuming that the blur is spatially invariant. The blurring model is modeled by Kvec(X (:, :, i, j)) = vec(G(:, :, i, j)) for i = 1, 2, 3 and

Using the approximation of the blurring matrix as Kronecker product, we have

Chapter 6. Tensor products with application to face recognition

The matrix C is clearly determined by its first column r = (r 1 , r 2 , . . . , r n ). Therefore, the above circulant matrix is also denoted by circ(r). Let Z be an n × n defined by:

Then it is seen easily that

The most important property of circulant matrices [START_REF] Gene | Matrix computations[END_REF] is that they are diagonalizable by the Fourier matrix F n :

where Λ = diag(F n r) is a diagonal matrix with the eigenvalues of C, F * n is the conjugate transpose of F n . Just as in the case of circulant matrices, a block circulant matrix can be diagonalized by the Discrete Fourier Transform (DFT) [START_REF] Kilmer | A third-order generalization of the matrix svd as a product of third-order tensors[END_REF]. In this case a block circulant matrix can be a block diagonalized. Lemma 6.2. Let A ∈ R n 1 ×n 2 ×n 3 be a third-order tensor, we have:

Proof. Let e k ∈ R n 2 n 3 , then e k can be written as e l ⊗ e m , with e m ∈ R n 2 , e l ∈ R n 3 , and

In analogue to the case of circulant matrices, we have

going to use in this work.

Proposition 6.5. The tensor I n 1 n 1 n 3 , whose first frontal slice is the n 1 × n 1 identity matrix, and whose other frontal slices are all zeros, is the identity element of the t-product:

Proof. Easy to verify from Definition 6.4. Definition 6.6. A third-order tensor is f-diagonal if each frontal slice is a diagonal matrix. Likewise, a tensor is f-upper triangular or f-lower triangular if each frontal slice is upper or lower triangular, respectively. Definition 6.7. The transpose of A ∈ R n 1 ×n 2 ×n 3 is a tensor A T ∈ R n 2 ×n 1 ×n 3 obtained by transposing each of the frontal slices and then reversing the order of transposed frontal slices 2 through n 3 .

Definition 6.8. Q ∈ R n 1 ×n 1 ×n 3 is said to be an orthogonal tensor if

Properties 6.2.1. [START_REF] Kilmer | Factorization strategies for third-order tensors[END_REF] 1. The t-product is associative: Let A, B and C be third-order tensors of appropriate size, then

2. Let A, B be third-order tensors such that A * B and B T * t A T are defined, then

Proposition 6.9. Let A ∈ R n 1 ×n 2 ×n 3 , B ∈ R n 2 ×m 2 ×n 3 be third-order tensors, for i =

and put

We point out that a rank (r 1 , . . . , r N ) approximation with r n ≤ R n for n = 1, . . . , N can be obtained simply by restricting the factor matrices U n to the first r n columns (truncated SVD) for n = 1, . . . , N, and by restricting the core tensor V. We write the HOSVD of the tensor A ∈ R n p ×n e ×n i in the following form:

with V ∈ R n p ×n e ×m i , U ∈ R n i ×m i and m i < n i . In application, we set m i = n e n p . Since in our illustrating examples, the numbers of persons and different expressions are small to number of pixels of images in vector format. Applying the t-QR or the cosine QR to the core tensor V, we obtain

Following from Definition 3.14 and from Propositions 6.9 and 6.19, if we fix an expression 'e' and a person 'p', we have A(p, e, :) = (Q(p, :, :) * R(:, e, :)) × 3 U.

So, for a given image z ∈ R n i , we have to solve the least squares problems

Conclusions and perspectives

In this chapter, we present a brief review of the results established in the previous chapters

Summary of results

In this thesis, we proposed algorithms for solving higher-order problems with applications to image restoration and face recognition.

In chapter 4, we constructed a degradation model associated with color images and videos. This model was constructed based on CP decomposition. Then we used the Alternating Least Squares (ALS) method to split the degradation model into three subproblems. To solve those sub-problems, we used an algorithm based on Golub Kahan bidiagonalization.

In chapter 5, we solved the tensor least squares minimization problem

-G , by decomposing the right-hand tensor using HOSVD or CP decomposition, and by working with the coefficient matrices A (i) , for i= 1, . . . , N. In order to overcome the problem of the curse of dimensionality.

In chapter 6, we introduced approaches to face recognition defined via the t-product, the cosine product, and the outer product. Our approaches are applied to a database represented by third or fourth-order tensors.

Perspectives

During this work, multiple questions were raised, which represent interesting axes for future work:

• Applying the face recognition algorithms to a database contaminated by blur or noise.

• Using tensor decompositions for large data completion.

Abstract

Tensors are multi-dimensional arrays. They are a generalization of matrices and vectors. They provide a natural way to represent the data in different fields, which make tensor field a great framework for formulating and solving many problems in different areas. One of the most important applications of tensors is in the field of image processing such as image restoration and face recognition, where color images (RGB images) are presented as third-order tensors and a video composed of color images is a fourth-order tensor. But when we work in higher dimension spaces to solve higherorder problems, a set of challenges arise such as a problem known as "the curse of dimensionality ". In fact, when the dimension increase the higher-order problems become harder (computation and memory requirement ) since the data size of a tensor increases exponentially with the increase of the dimensionality of the tensor. As consequence tensor computations become much expensive. In this thesis, we have focused on solving some tensor problems. The suggested algorithms are obtained by combining iterative methods such as LSQR method and higher-order decompositions to overcome the problem of dimensionality. Our suggested approaches are applied to the restoration of images and videos.

In addition, this thesis studies the face recognition methods based on the tensor format, where multi-linear algebra tools such as the HOSVD (Higher-Order Singular Value Decomposition) have been used. We suggest a new algorithm that can be applied to a database of images represented by a third or fourth-order tensor.

Keywords: HOSVD, CP decomposition, color image restoration, video restoration, LSQR, face recognition, t-product, cosine product.