High Order Singular Value Decomposition for Plant Biodiversity Estimation

We propose a new method to estimate plant biodiversity with R{\'e}nyi and Rao indexes through the so called High Order Singular Value Decomposition (HOSVD) of tensors. Starting from NASA multispectral images we evaluate biodiversity and we compare original biodiversity estimates with those realised via the HOSVD compression methods for big data. Our strategy turns out to be extremely powerful in terms of storage memory and precision of the outcome. The obtained results are so promising that we can support the efficiency of our method in the ecological framework.


Introduction
In order to face Earth changes in space and time, satellite imageries are nowadays being used to provide timely information over the whole globe.
From this point of view, starting from the early '70s, technological improvements in remote sensors led to a growth in the number and in the size of available data. For example to store a band of the entire Earth's surface from the MODIS sensor with a low spectral resolution, 5600 m, we need around 99 MB. Therefore for rasters with higher spectral, spatial or radiometric resolutions the memory request increases significantly. One may point out that 99 MB are not much, if compared with the features of modern machines. However the memory needed to process thousands of images might represent a crucial issue. Moreover, in order to perform any index computation over these images, they should be loaded into the computer RAM, which usually has lower capacity and is occupied also by the system application and by our computing script. Since the advent of big data, which include multispectral and hyperspectral images, improving storage and analytical techniques became fundamental. Mathematicians are facing this challenge together with computer scientists, physicists and engineers (cf. e.g. [DLDMV00, VVM12, ADKM11, LC10, STG + 19, BP94, BBK + 16, BDHM17, GSR + 14, SSV + 05, BBCM11, BBC + 19]). The usual structure for storing multispectral images are tensors. For matrices approximation there exists an optimal technique, called Singular Values Decomposition, SVD (cf. [Sch07,Ste93,Mir59,EY36]). One application of SVD is image compression (cf. e.g. [WBWJ15,JXP12,GVL96]). Since we need to store data in tensors, the generalisation of matrices, we will present High Order Singular Value Decomposition (HOSVD) firstly introduced in [DLDMV00, VVM12], a generalisation to tensors of SVD to tensors. There are various techniques to approximate tensors by taking advantage of the HOSVD, those that will be implemented in this paper are the so called Truncated HOSVD (T-HOSVD) and Sequentially Truncated HOSVD (ST-HOSVD). Even if only for some special cases HOSVD provides optimal results ([VNVM14, BDHR17, DOT18]), it is possible to present an estimate of the tensor approximation errors. Indeed the core of the present paper will be application of T-HOSVD and ST-HOSVD techniques, their modern versions and the error made by using them. Indeed in the last section we will apply some possible HOSVD implementations to tensors in which we stored RED and NIR bands. Next from the compressed tensors we get NDVI rasters and over them we compute biodiversity indexes. So we will be able to compare our biodiversity estimations from compressed data with those realised over original data. As far as we know, this is the first attempt to estimate biodiversity with the presented indexes from HOSVD compressed data. Moreover the obtained results are extremely promising, letting us support its efficiency in ecological framework.

Ecological background
One of the main component of Erath's biosphere is biodiversity, which is in strict relationship with the planet status in space and time. Measuring biodiversity over wide spatial scales is difficult concerning time, costs and logistical issues, e.g.: (i) the number of sampling units to be investigated, (ii) the choice of the sampling design, (iii) the need to clearly define the statistical population, (iv) the need for an operational definition of a species community, etc. (cf. [Chi07]). Hence, ecological proxies of species diversity are important for developing effective management strategies and conservation plans worldwide (cf. [RBC + 10]). From this point of view, environmental heterogeneity is considered to be one of the main factors associated to a high degree of biological diversity given that areas with higher environmental heterogeneity can host more species due to the greater number of available niches within them. Therefore, measuring heterogeneity from satellite images might represent a powerful approach to gather information of diversity changes in space and time in an effective manner (cf. [ZFL + 19]). Nowadays, the advent of satellites has made possible having real images of a territory, even with a remarkable quality. This approach is behind the the remote sensing discipline. Many definitions have been proposed during the years for this subject. We refer to the one presented in [CW11,p.6].
Definition 2.1. Remote sensing is the practice of deriving information about Earth's land and water surfaces using images acquired from an overhead perspective, using electromagnetic radiation in one or more regions of the electromagnetic spectrum, reflected or emitted from the Earth's surface.
This field of study is based on a well known physical phenomenon: different materials and different organisms absorb and reflect electromagnetic radiations differently. Most of the modern satellite sensors can acquire multiple images, divided in the so called bands.
Definition 2.2. Bands or channels are the recorded spectral measurements.
Depending on the sensor features, we can acquire images dived into the different bands.
Definition 2.3. Multi-spectral sensor can acquire from 4 up to 10 bands. Hyper-spectral sensor can acquire more than 100 bands.
Data used in the present paper come from the MODerate-Resolution Imaging Spectroradiometer sensor, or MODIS, built on both the satellites Terra and Aqua, cf. [CJ]. MODIS measures 36 bands in the visible and infrared spectrum, at different resolutions. Firstly there are RED and NIR bands with pixel size of 250km and next 5 bands, still with RED and NIR, at 500m of spatial resolution. These are extremely useful for land observation. The remaining bands with 1km of resolution consist of monitoring images from visible spectrum, MIR, NIR and TIR. The data are registered four times a day: twice daytime and twice night-time. Images usually include the entire Earth's surface. From the 60s the scientific community highlighted two important relations between spectral measurements and biomass: • a direct relation between NIR region and biomass, i.e. greater the biomass greater the NIR reflected radiation measured and vice versa; • an inverse relation between RED spectral region and biomass, i.e. greater the biomass, lower the RED visible reflected spectrum measured and vice versa. Therefore the relation between NIR and RED reflectance is central for the vegetal biomass estimation. The aim of vegetation indexes is measuring biomass or vegetative vigour on the base of digital brightness values. The one used in the present article is the Normalized Difference Vegetation Index, or NDVI, presented in 1974 by Rouse and others, cf. [RHSD74].
Definition 2.4. Given a region R, let RED, N IR ∈ M m×n (R) be respectively the RED and the NIR raster band of R imagery. The normalized difference vegetation index of region R is N DV I ∈ M m×n (R) such that . . , m} and for every j ∈ {1, . . . , n}, when it is defined.
The innovative idea of [RDB + 13] was the application of information theory studies and ecological indexes to remote sensed images. The ecological indexes used in the present paper are the Rao and Rényi ones, which we are going to define.
Definition 2.5. Given a spectral image of a sample area, let N be the image radiometric resolution and let p i be the relative abundances of the i-th value for every i ∈ {1, . . . , N }. Fixed a distance function d, we build up a pairwise spectral difference matrix D ∈ M N (R) such that for every i, j ∈ {1, . . . , N }. The Rao's Q index for the sample area is Definition 2.6. Given a sample area with N species and defined p i the relative abundances for every i ∈ {1, . . . , N }, in decreasing order, the Rényi index is The main difference between them is that Rao index takes into account both the frequencies and the numerical values of each pixel. On the other side Rényi index considers only pixels frequencies.
Lastly notice that rasters usually are split into small chunks, called windows, over which the biodiversity chosen index is computed. We need to generate Rényi's and Rao's computation codes to measure how much biodiversity information is lost using approximated tensor data. Since we keep the ecologists moving window approach, our implementation presents also a very basic multi-core modality. The multi-cores Rény's and Rao's implementations crux is presented in the Appendix, section 6.1.

Mathematical background
The celebrated Schmidt [Sch07,Ste93], Mirsky [Mir59], Eckhart, Young [EY36] theorem states that for every rank-r matrix A ∈ M m×n (C) there exists the best rank s approximation for every s < r in the Frobenius norm, actually the result holds for arbitrary (O m × O n )-invariant norms [Mir60]. The famous generalization of this result [DLDMV00] to higher order tensors fails in general to produce an output that is the "best" approximation of a given tensor, cf. [VVM12]. It succeeds for the so called orthogonally decomposable tensors [BDHR17]. Another type of generalization of the concept of "best rank-r approximation of a tensor" proposed in [DOT18] works for general tensors. Despite the possible non-existence of the best approximation of a tensor under the construction of [DLDMV00], this technique turns out to be very convenient from the computational point of view and it is possible to prove that the outcome is a "quasi-optimal solution. In this paper we will mainly use the technique popularized by L. De Lathauwer and J. Vandewalle in [DLDMV00], the so called High Order Singular Value Decomposition (HOSVD), we will study two types of errors made by the specific approximation, we will make a comparison among them and we will show that they will be small enough to have very precise results.
Here some basic mathematical preliminaries. The number field will always be the complex field of numbers C.
Definition 3.1. Let V be a linear subspace of the tensor space C n 1 ⊗ · · · ⊗ C n d . If for every i ∈ {1, . . . , d}, there exist a subspace V i ⊆ C n i such that Remark that not every subspace of C n 1 ⊗ · · · ⊗ C n d is separable. The structure of a separable tensor space has some consequences on its elements.
Definition 3.2. The multilinear rank of a tensor A ∈ C n 1 ⊗ · · · ⊗ C n d is the d-uple (r 1 , . . . , r d ) with the property that r i is the minimal dimension of a subspace V i ⊂ C n i such that A ∈ V 1 ⊗ · · · ⊗ V d for every i ∈ {1, . . . , d}.
Let A ∈ C n 1 ⊗ · · · ⊗ C n d be a tensor and let (r 1 , . . . , r d ) ∈ N d . We will discuss if there exists, and in case how to determine, a tensor M of multilinear rank lower or equal component-wise than (r 1 , . . . , r d ) which minimizes the Frobenius norm of the tensor difference, i.e.
This problem is also known as Low MultiLinear Rank Approximation (LMLRA). M. Ishteva and L. De Lathauwer firstly stated this problem and introduced this acronym, cf. [IAVHDL11]. Looking for a tensor M ∈ C n 1 ⊗ · · · ⊗ C n d satisfying the stated rank properties means searching a subspace V i ⊂ C n i of dimension r i for every i ∈ {1, . . . , d} such that, if the approximation tensor M exists, it belongs to V 1 ⊗ · · · ⊗ V d . The techniques of T-HOSVD and ST-HOSVD are described by L. De Lathauwer, B. De Moor and J. Vandewalle, in [DLDMV00] and by N. Vannieuwenhoven, R. Vandebril and K. Meerbergen, in [VVM12] respectively. Here we will recall the fundamental aspects for the reader convenience.
More generally, let p q = [1, d] be a partition of d with s = |p| and t = |q|. Then, we can associate to A even more multilinear maps, namely This interpretation of A is called a flattening of A.
The punchline of the HOSVD is to apply the SVD to the flattenings of a given tensor and to use the More-Penrose transform, Definition 3.4, to build an approximated tensor. As already outlined, this technique is optimal for the so called orthogonally decomposable (ODeCo) tensors, cf. [BDHR17], while for other cases there are no evidences that this procedure would lead to "the best multilinear-rank approximation" of a given tensor. It's worth noting that, unlike matrices, tensors of higher order can fail to have best rank-r approximations, cf. [dSL08]. Anyway the various algorithms of the HOSVD are extremely explicit and it is possible to estimate the measure of the error made by thanking this approximation. We will show that for the applied purpose of this paper the approximation is very good in the considered problem.
HOSVD provides a sparse representation of the given tensor, whose costs in terms of storage use can be easily computed. Indeed given A ∈ C n 1 ⊗ · · · ⊗ C n d whose multilinear rank is known to be (r 1 , . . . , r d ), from HOSVD we get that A = (U 1 , . . . , U d ) · C with U i rank-r i orthogonal (n i × r i )matrices and C ∈ C r 1 ⊗ · · · ⊗ C r d the so called core tensor. The storage of each matrix U i costs n i r i memory units for every i ∈ {1, . . . , d} and storing the core tensor C needs d i=1 r i memory units. In conclusion the sparse representation of A costs r i memory units which is extremely better than d i=1 n i . The following proposition reveals the main strategy of the HOSVD.
Proposition 3.6 (HOSVD [VVM12]). Let V = V 1 ⊗ · · · ⊗ V d be a separable tensor subspace of be an orthogonal basis of V for the standard product, and let (U i ) = (u i j i ) i=1,...,d;j i =1,...,r i be the corresponding orthogonal bases for the V i 's, i=1, . . . , d. The projector π i : C n 1 ⊗ · · · ⊗ C n d → C n 1 ⊗ · · · ⊗ C n d is such that for every A ∈ C n 1 ⊗ · · · ⊗ C n d for every i ∈ {1, . . . , d}. Define P V (A) := π 1 . . . π d (A). Then for every A ∈ C n 1 ⊗ · · · ⊗ C n d and for every ρ permutation of d elements, we get that: (1) Corollary 3.7. Under the hypothesis of Proposition 3.6, we have that for every tensor A ∈ C n 1 ⊗ · · · ⊗ C n d Usually, in applications, the multilinear rank of the given tensor is not known. Consequently, working on the exact error is not convenient. Therefore, on the basis of Proposition 3.6 and Corollary 3.7, N. Vannieuwenhoven K. Meerbergen and R. Vandebril developed two new strategies to approximate tensors, cf. [VVM12], when the multilinear rank is not known. The first is based on the approximation of the upper bound of the error, stated in Corollary 3.7. As a matter of fact reducing the upper bound implies reducing the the exact error.
The Truncated Higher Order Singular Value Decomposition (T-HOSVD) has the SVD as key concept. Fixed A ∈ C n 1 ⊗ · · · ⊗ C n d , consider the i-th term of the upper bound summation, i.e.
By the positivity of each term of the previous expression, minimizing the upper bound means minimizing each term of the upper bound summation, i.e. taking the minimum of the norm of the difference between A and its projection in just one direction each time. Because the tensor is approximated only with respect to one direction each time, the minimization problem is mathematically arg min i.e. looking for for the best approximation at rank r i of the i-th flattening of A for every i ∈ {1, . . . , d}. However thanks to Schmidt, Mirsky, Eckhart, Young theorem the problem for matrices has a close solution which is obtained through SVD of the i-th flattening of A, truncated at the r i singular values for every i ∈ {1, . . . , d}. It is clear that the order of projectors application is not significant for the T-HOSVD. This won't be true anymore for the Sequentially Truncated High Order Singular Value (ST-HOSVD). The idea of the ST-HOSVD is minimizing each term of the summation of Proposition 3.6. Let A ∈ C n 1 ⊗ · · · ⊗ C n d be a tensor and let (r 1 , . . . , r d ) be the target multilinear rank of the approximation. The first step is looking for the projector which minimises the first error term of Equation (1): arg min i.e. arg min The last formulation of this first step has a close solution, thanks to Schmidt, Mirsky, Eckhart, Young theorem. We are looking for the best rank r 1 approximation of the 1-st flattening. So we can conclude that the matrix U 1 obtained from the SVD of the first flattening of A truncated at the r 1 -th column is such thatÛ 1 = arg min U 1 ∈O(n 1 ×r 1 ) The core tensor of the first step, Fixed theπ 1 =Û 1Û H 1 , the next step is looking for the projector which minimises the second error term of Equation (1): But thanks to multilinearity, the last equation becomes arg min arg min SinceÛ 1 is fixed, the second step problem becomes arg min 2 whose solution is the matrixÛ 2 obtained through the SVD of the 2-nd flattening of C (1) truncated at the r 2 -th column. Defined the new core tensor we proceed similarly with the next ones. At the (d − 1)-th step the (d − 1)-th core tensor is defined as For the last d-th step, we look for the projectors which minimises the last error term of Equation (1): which using the multilinearity and the projectors properties becomes arg min arg min The solution of the last problem is the matrixÛ d obtained through the SVD of the d-th flattening of the core tensor C (d−1) truncated at the r d -th column.
We can now state both the T-HOSVD and the ST-HOSVD algorithms, whose implementations are in the Appendix, Section 6.2 and Section 6.3 respectively.

Results
In this section we will first describe the data used: RED, NIR and NDVI bands of different territories, provided by NASA, [Did18b,Did18a]. The first step in our analysis will be computing the biodiversity index over the NASA NDVI imageries. Then we will generate a new NDVI from RED and NIR using Definition 2.4, since we do not know how NASA generates NDVI from the other two bands. Over these "relative" NDVI we will compute the biodiversity indexes. Next we will generate 3-order tensors with just RED and NIR bands. Then we will approximate these tensors with T-HOSVD and ST-HOSVD. Lastly from the approximated tensors we will extract RED and NIR imageries to get "approximated" NDVI. Over these last NDVI we compute the biodiversity indexes. In conclusion we will measure the error made in estimating biodiversity when we use approximated NDVI instead of NASA NDVI or relative NDVI. 4.1. Data. From MOD13C2v006 and MOD13A3v006, both sensed by MODIS, cf. [CJ], and available at [Did18b,Did18a], we select the RED, NIR rasters and the NASA computed NDVI. MOD13C2v006 is a product characterized by 13 layers, each of which stores an imagery of the Earth with different properties. Each raster is a matrix of 3600 rows and 7200 columns. The side of each pixel corresponds to 5600 m, which is the spatial resolution. The imageries are monthly: they are obtained from the daily data through NASA's algorithms. The chosen data from January 2010 until December 2018 are in hdf format. MOD13A3v006 is a similar product with a higher spatial resolution. While each image from the previous dataset covered the entire Earth's surface, each one from this second dataset covers just 1200 × 1200 m 2 . We select the 20 components, called granules, to compose an Europe's map. We download in GEOTiff the granules from RED, NIR bands and from the NDVI computed by the NASA, from June 2011 until December 2018. Also in this case they are obtained by the NASA scientists from daily data. The final dimension of most of the rasters is 4800 rows and 6000 columns.
We do not talk about every used test elements from MOD13A3v006 dataset, since those of December 2016 and of December 2017 do not include all the 20 granules. Moreover we remark that rasters of December 2012 and December 2015 have the correct dimension, but they have respectively one and two missing areas. Lastly NASA stores the data in 16-bit signed integers. The next step is creating 3-order tensors. Taking advantage of the GDAL dependencies for python, we simply convert rasters into matrix, removing the metadata useless for our aims. Then we store NDVI, RED and NIR into 3-order arrays.
Definition 4.1. Let T E be the set of all the tensor A ∈ R n 1 ⊗ R n 2 ⊗ R 3 with n 1 = 4800 and n 2 = 6000 such that • A ·,·,1 is the NDVI raster, • A ·,·,2 is the RED raster band, • A ·,·,3 is the NIR raster band of Europe dataset for every month and for every year. The cardinality of T E is n E = 91. Similarly let T W be the set of all the tensor A ∈ R m 1 ⊗ R m 2 ⊗ R 3 with m 1 = 3600 and m 2 = 7200 such that • A ·,·,1 is the NDVI raster, • A ·,·,2 is the RED raster band, • A ·,·,3 is the NIR raster band of Earth dataset for every month and for every year. The cardinality of T W is n W = 108.
Since obtaining biodiversity indexes takes long time and needs many resources, we perform all the computation over HPC@UniTrento, the university of Trento cluster. The indexes measured are Rao and Rényi both with window side equal to 11. We maintain the same window side on Europe's and on Earth's images, since the different spatial resolutions lead to similar rasters dimensions. 4.2. NASA and relative NDVI. The first step is computing both Rao and Rényi indexes over the NDVI raster, extracted from the loaded tensor, i.e. over (A k h ) ·,·,1 for every k h ∈ {1, . . . , n h } for every h ∈ {E, W }. Definition 4.3. Let R h be the set of Rao index computed over (A k h ) ·,·,1 for every k h ∈ {1, . . . , n h } for every h ∈ {E, W }. Similarly let I h be the set of Rényi index computed over (A k h ) ·,·,1 for every We use the obtained images as comparison term. Then we compute also an NDVI from the RED and NIR raster, using Definition 2.4, since the algorithm used by NASA for NDVI creation is not public.
Let f h = l • g h be the function that associates to each tensor A ∈ T h the NDVI raster obtained from (A k h ) ·,·,2 and (A k h ) ·,·,3 for every k h ∈ {1, . . . , n h } for every h ∈ {E, W }. Theñ Remark 4.5. Since rasters have only non negative elements for every NIR and RED band, the second case of function l in the previous definition is verified when both elements of NIR and RED rasters are zero. In that case we assign to NDVI the default missing value, M = −3000.
Consequently we perform again index estimation over N h elements for every h ∈ {E, W }.
Remark 4.7. Henceforth we assume that elements of the same set pairs (R h ,R h ) and (I h ,Ĩ h ) are ordered equally for every h ∈ {E, W }.
Even now, we can present some comparison between these two types of estimates. We compute the error A j − B j F for every A j ∈ R i and B j ∈R i for every j ∈ {1, . . . , n i } for every i ∈ {E, W }. Next we also compute the error per pixel dividing the error by the number of pixels. Since we are working with huge matrices, this type of distributed error is useful to understand how big is on average the error made pointwise. In the following table we present some statistics, where e is the error vector and ep is the error per pixel vector. Besides for every v ∈ C n , we set We make on average a 0.6% error per pixel for Europe Rao estimation using self-made NDVI instead of NASA NDVI, while the error is on average of 6% per pixel when Rao is computed over Earth NDVIs. Both the unbiased variance are very small, which means that errors are not very different from the mean. Notice that also the minimum approximation error is greater for Earth than for Europe data. We suppose that this is linked with the water surface greater in World rasters than in Europe's ones. Similarly we compute the error A j − B j F for every A j ∈ I i and B j ∈Ĩ i for every j ∈ {1, . . . , n i } for every i ∈ {E, W }. With the same notation already introduced, we present some statistics for Rényi index. In this case the mean error per pixel for both the dataset is smaller than the previous mean. One possible explication could be that Rényi index takes into account only frequencies while Rao index considers both values and their frequencies. However in this case variance in World error is smaller than in Europe, while in the Rao case there is the opposite situation. Moreover we observe that the minimum and the maximum approximation error in the World dataset for both the indexes is realised by the same element, i.e. April 2018 and May 2014 respectively. This consideration holds also for Europe dataset: the minimum error comes from April 2018 and the maximum from December 2012. 4.3. Approximated NDVI. Before applying the approximation tensors codes, described in the Appendix, Sections 6.2 and 6.3, we have to highlight one limit of the python function svds. It takes as additional parameter k, which is the rank of the wanted approximation and which has to be strictly lower than both the dimensions of the given matrix. So if we had passed just a 3-order tensor such that n 3 = 2, for the third flattening we would have fixed k equal to 1, getting a vector: this is a too low order tensor for our aims. Therefore we decide to increase n 3 up to 3, adding another matrix to our tensor: in the first case we take twice RED band raster and once NIR one, in the second case we take twice NIR band and once RED one.
Definition 4.9. Let R = {10, 50, 100, 500, 1000} be a set with the given order fixed, then the target multilinear rank we choose are r j = (i j , i j , 2) for every i j ∈ R for every j ∈ {1, . . . , 5}. Let T k,h,j be the set of T-HOSVD approximation at multilinear rank r j of tensors from the set T k,h for every h ∈ {E, W }, for every k ∈ {N, R} and for every j ∈ {1, . . . , 5}. Similarly let S k,h,j be the set of ST-HOSVD approximation.
Before presenting results related to indexes computation, we have some data about storage memory use. Since it depends on the core tensor and on the projectors dimensions, which are equal for T-HOSVD and ST-HOSVD, we report only a table for each dataset. For each tensor A ∈ T k,W the ratio between the memory used for storing the core tensor and the projectors over the memory used for storing A is the same, for every k ∈ {N, R}. For tensors of T k,E it holds the same, except for those elements composed by a lower number of granules, for every k ∈ {N, R}. Since they are 2 over 91, we neglect them and in the table we present the ratios of memory usage for each rank approximation. We call these ratios absolute compression ratios, because they have as denominator the memory used to store two time RED band and once NIR band, or vice-versa. In the following table they are present as Abs. R, when RED band raster is repeated and as Abs. N, in the other case. Beside we list also a relative compression ratio, where the denominator is the amount of memory necessary to store once RED and once NIR band. In the table it is Rel. Moreover here and all along the paper for simplicity we write as rank only the significant components of the multilinear rank, i.e. i j for every i j ∈ R for every j ∈ {1, . . . , 5}.

Rank
Europe Earth  Table 3. Rate of compression.

Rel
Remark that even with the greater component wise target multilinear rank, we need just a small percentage of memory with respect the one used for storing the entire tensor. In addiction to this, even the relative ratio at the highest multilinear rank present a significant saving in memory use. In order to not get lost during the presentation with the numerous indexes used, we will give a general idea. After having generated new tensors and having approximated them, we compute new NDVIs through function f h of Definition 4.4 applied on T k,h,j ∪ S k,h,j for every h ∈ {E, W }, for every k ∈ {N, R} and for every j ∈ {1, . . . , 5}. We compute biodiversity index over them. Lastly we measure the difference in estimating biodiversity from approximated NDVI and NASA or self-made NDVI. Definition 4.11. Let I R,h,k,j be the set of Rényi index computed over elements of N R,h,t,j for every h ∈ {E, W }, for every t ∈ {T, S} and for every j ∈ {1, . . . , 5}, and call them i j -approximated estimates for every i j ∈ R and for every j ∈ {1, . . . , 5}.
Remark 4.12. Notice that we have 4 indexes for approximated estimates set: 1 st index: indicates the repeated matrix in the starting tensor R for RED and N for NIR; 2 nd index: the belonging dataset, E for Europe and W for Earth; 3 rd index: the approximation algorithm, T for T-HOSVD and S for ST-HOSVD; 4 th index: is associated with the target multilinear rank.
The last step is computing the error with respect to the original estimates, i.e.
A k − C k,j for every A k ∈ I h and for every C k,j ∈ I R,h,T,j ∪ I R,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. Moreover we compute the error with respect to relative estimates, i.e. B k − C k,j for every B k ∈Ĩ h and for every C k,j ∈ I R,h,T,j ∪ I R,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. Now we present some statistics for the errors per pixel with respect to original estimates, stored in vector epO and relative estimates, in vector epR for Europe dataset, for both the decomposition techniques for each target multilinear rank. In the following table we report as rank only the significant component of the multilinear rank, for simplicity.
Remark 4.13. To simplify the discussion henceforth and all along the article the i j -original error will be the error between original estimate and i j -approximated estimate, while the i j -relative error will be the error between relative estimate and i j -approximated estimate for every i j ∈ R.  Table 4. Statistics for Rényi index over N ∈ N R,E,t,j .
Notice that ST-HOSVD original and relative error per pixel is lower than T-HOSVD original and relative error per pixel, when the first two components of target multilinear rank are greater or equal than 500. The variance of both errors is quite low, even if it increases in the last three multilinear ranks. Moreover we underline that the minimum relative error and the minimum original error are equal up to the forth decimal digit. They also decrease, when the components of multilinear rank increase. On the other hand the maximum of relative errors and the maximum of original errors do not coincide. Beside, they increase significantly in the forth and fifth approximation. Finally we notice that the average relative error per pixel is frequently slightly greater than the original one. For T-HOSVD this inequality between original and relative error average happens from the third approximation, while for the ST-HOSVD from the second one. The difference seems to grow for increasing multilinear rank components. This could appear a bit strange, since we expect the contrary. However we have to remark that Rényi index takes into account only raster values frequencies, neglecting the values themselves. In addiction from the complete data, we observe that the relative error of elements with missing granules, tensors of December 2012 and December 2015, is more than 3 times the original error.
We can now list statistics about the Earth dataset. Similarly in vector epO there are the errors per pixel with respect to original Rényi estimates, while in epR with respect to relative estimates for each target multilinear rank.  Table 5. Statistics for Rényi index over N ∈ N R,W,t,j . We remark that even in this case on average ST-HOSVD technique leads to lower original and relative errors than T-HOSVD, for the last two and for the last one target multilinear rank respectively. Moreover we observe that relative error mean is slightly lower than original one, as we expected. Minimum and maximum of both errors decrease for increasing multilinear rank components. Certainly the most stunning value is the variance of relative error per pixel at the last approximation. For both T-HOSVD and ST-HOSVD it is lower than 10 −4 . We include the images associated to Rényi index in the five approximations for the Europe worst case and the Earth best case, both from ST-HOSVD approximation technique.
Example 4.14. Looking at the Rényi index computed over Europe NDVI of February 2013 as it is in Figure 1, we immediately notice that biodiversity seems to be quite high, near 4.5 almost everywhere in Europe. However as remarked in [RMR17], Rényi index tends to overestimate biodiversity. Besides we underline that Rényi index computed over self-made NDVI does not differ much from its computation over NASA NDVI. Next we have in Figure 2A the same index computed over self-made NDVI and the approximated Rényi estimates at different multilinear ranks. We can notice that when the first two components of the multilinear rank grow, in the Rényi estimation some new noising elements appear, leading to high errors. We believe that this type of phenomenon deserves further analysis. However in the internal land of Europe, the biodiversity estimation is quite close to the relative and original one, for multilinear rank components strictly greater than 100.
In the next example we will consider the element of Earth dataset, which realises the minimum error.
Example 4.15. Firstly we show in Figure 3 the Rényi estimate over NASA NDVI of October 2017. As we said in the Example 4.14 Rényi index provides quite high biodiversity values. Indeed also in this case there are many Earth areas with a biodiversity value close to 4.5. Then we present the same index over self-made and approximated NDVI. Even if the small printing dimensions reduce the detail precision, our eyes do not perceive at first glance much difference between the last three index approximations and the index of Figure 3. However with a closer observation, for example, we notice that islands of the Pacific ocean disappear in first approximation and they partially reappear in the last two images. We may conclude that the original and relative error is in this case linked with these missing territories, but also with the different evaluation of biodiversity in Amazon, for example.
Next step in our discussions is computing Rényi index over approximated NDVI of N N,h,t,j for every h ∈ {E, W }, for every t ∈ {T, S} and j ∈ {1, . . . , 5}.
Definition 4.16. Let I N,h,t,j be the set of Rényi index computed over elements of N N,h,t,j for every h ∈ {E, W }, for every t ∈ {T, S} and for every j ∈ {1, . . . , 5}. We decide to call these i j -approximated estimates for every i j ∈ R and for every j ∈ {1, . . . , 5}.  As before, we compute the error with respect to the original estimates, i.e.
A k − C k,j for every A k ∈ I h and for every C k,j ∈ I N,h,T,j ∪ I N,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. Moreover we compute the error with respect to relative estimates, i.e. B k − C k,j for every B k ∈Ĩ h and for every C k,j ∈ I N,h,T,j ∪ I N,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }.
Lastly we report some statistics about i j -original and i j -relative errors per pixel, stored in vector epO and epR for each i j ∈ R. Firstly a table for Europe related data.  Table 6. Statistics for Rényi index over N ∈ N N,E,t,j .
Almost every consideration for statistics of approximated estimates of I R,E,t,j for every t ∈ {T, S} and for every j ∈ {1, . . . , 5} holds also in this case. However we can notice that on average the 100 relative and original error are lower that 500 one, but this is not true anymore for 1000 relative and original. Moreover when the first two multilinear rank components are strictly smaller that 500, T-HOSVD relative and original errors are lower than ST-HOSVD ones. For the statistics over I R,E,t,j elements, this consideration holds only for relative error at the third multilinear rank. The remarks about not full granules elements are true also in this case. We underline that at each multilinear rank both the original and the relative errors on average are smaller in this second case, i.e. applying the procedure to tensors where the NIR band raster is repeated. Lastly some statistical aspects about Earth approximated estimates. As previously, we list a table with mean, variance, min and max for i j -original and i j -relative errors per pixel, stored respectively in vector epO and epR for each i j ∈ R.  Table 7. Statistics for Rényi index over N ∈ N N,W,t,j . Also in this case the considerations presented for I R,W,t,j error statistics hold. Notice that again variance is lower than 10 −4 for 1000-relative error of ST-HOSVD. Moreover also in this case ST-HOSVD is convenient only when the first two multilinear rank components are greater than 500. Lastly original and relative error on average are again smaller in this second case, i.e. when we apply our method to tensors where NIR band is repeated. We will present briefly the Rényi index image over Europe NDVI of February 2013 and over Earth NDVI of October 2017, obtained starting from the correspondent tensors of T N,h for h ∈ {E, W }.
Example 4.17. In Example 4.14 we presented the Rényi index image, which realises the highest original and relative error, starting from RED repeated band. Here we have the image of Rényi index for the same element, with approximation obtained from NIR band repeated. We remark that also starting from twice NIR and once RED band tensor element of February 2013 realises the highest original and relative error. Again we observe an increasing presence of noise in the north Europe area for growing multilinear rank components. Example 4.18. Similarly we display approximation of Rényi index for Earth element of October 2017, obtained from tensors where is repeated twice the NIR band. As in Example 4.14 we underline that the number of detected territories in the Pacific ocean grows when the first two multilinear rank components grow. Moreover comparing Figure 4F and Figure 6E, we notice that the amount of Pacific island present is greater in the second one. Indeed this is the element which minimises the original error also in this second proceeding way. Therefore even in this case we can affirm that choosing a starting tensor with repeated NIR band is more convenient that starting with repeated RED band.

Rao index. In this second part Rao index is computed over approximated NDVI rasters.
Definition 4.19. Let R R,h,k,j be the set of Rao index computed over elements of N R,h,t,j for every h ∈ {E, W }, for every t ∈ {T, S} and for every j ∈ {1, . . . , 5}. We call these i j -approximated estimates for every i j ∈ R and for every j ∈ {1, . . . , 5}.
As we did previously, we compute the error with respect to the original estimates, i.e.
A k − C k,j for every A k ∈ R h and for every C k,j ∈ R R,h,T,j ∪ R R,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. Moreover we compute the error with respect to relative estimates, i.e.
B k − C k,j for every B k ∈R h and for every C k,j ∈ R R,h,T,j ∪ R R,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. To describe our results we report some statistics about i j -original and i j -relative error per pixel, stored in vector epO and epR for each i j ∈ R. Firstly we list information for Europe related data. The most evident aspect is the high mean of both original and relative error made, even when the  Table 8. Statistics for Rao index over N ∈ N R,E,t,j .
components of multilinear rank grow. Indeed this average is close to a 20% of error per pixel. If we compare it with the average error per pixel made for Rényi index, we could think that in this case HOSVD is not performant. However we want to remark two critical aspects: firstly Rao index takes into account also the values of NDVI, not only their frequencies. Next we remind that at multilinear rank r 5 we are keeping nearly 15% of the total information, as in Table 3. Therefore even if results on average are not as good as in Rényi case, we do not exclude the power of this method for Rao index computation. Besides in the best case we have both original and relative error per pixel near to 11%, which is appreciable, as we will see. Lastly we remark an interesting and unexpected twist. If in the previous analysis T-HOSVD performed better when multilinear rank components were small with respect to ST-HOSVD, in the present case T-HOSVD provides better result than the other algorithm when the first two multilinear rank components are greater that 500. Notice also that in this case the relative error is on average greater than the original one. This phenomenon is at least in part linked to the two incomplete elements, December 2012 and 2015. Indeed for these elements the relative error is slightly less than twice the original error.
Now we present some statistics for the original and relative errors of Rao estimates for Earth dataset. As before in the following table epO is the vector of the errors per pixel computed with respect to original estimates, while in epR we store the errors per pixel with respect to relative estimates.  Table 9. Statistics for Rao index over N ∈ N R,W,t,j .
We immediately notice that even with the greatest multilinear rank components, the average original and relative error per pixel is quite high if compared with the ones in Rényi case. We again remark that 16% error per pixel on average has been obtained using only 16.4% of the total information. However even in the best case the relative error per pixel is 14.5%, which is quite enough if compared with the minimum error per pixel in the Europe previous case. Besides we notice that the relative error per pixel is slightly lower that the original ones as in the previous Earth case. As before we have that when the first two components of the multilinear rank are lower than 500, ST-HOSVD techniques leads to better results than T-HOSVD. Moreover as in the previous case the variance decreases when the first two multilinear rank components grow.
Before moving to the analysis of error related to Rao index over N N,h,p,j for every h ∈ {E, W }, for every p ∈ {T, S} and for every j ∈ {1, . . . , 5}, we show the best case for Rao index for Europe and the worst case for Earth dataset in the following examples. In both the case we use as approximation technique T-HOSVD.
Example 4.20. Looking at Rao index computed over Europe NDVI of December 2014, in Figure  7, we immediately notice that biodiversity is much lower than in the Rényi case, 4.14. Besides we underline that even in this case Rao index computed over self-made NDVIs is enough close to its computation over NASA NDVI. Lastly we highlight that in this particular NDVI element the north Europe area is entirely set to missing values, which are shown with white colour. Therefore we can assume that the error is the minimum since this peculiarity. Next we have in Figure 8A the same index computed over self-made NDVI and the approximated Rao estimated at different multilinear rank. It is significant that it is hard to distinguish the Rao index computed over NASA NDVI from the ones computed over the forth and the fifth approximated NDVI.
Next we analyse the worst case of Rao index for the Earth dataset.
Example 4.21. Firstly in Figure 9 we present Rao index over NASA NDVI of May 2014. As in the Europe case, biodiversity estimated by Rao index has lower values. In the tropical regions for example we have the highest values, which are close to 2. Moreover there are some areas as the Sahara desert which present extremely low biodiversity values.
Next we show in Figure 10 the Rao index computed over self-made and approximated NDVIs. It is clear observing the approximated Rao pictures and index of 9 that the high error is linked with the overestimation of biodiversity in the North pole area. Indeed it is a remarkable phenomenon in the first four approximations. In Figure 10F it tends to decrease, but it is still appreciable. Besides also the majority of the Pacific ocean islands are missing in all the approximations. These regions were probably related to error also in Example 4.15.  Then we compute the error with respect to the original estimates, i.e.
A k − C k,j for every A k ∈ R h and for every C k,j ∈ R N,h,T,j ∪ R N,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }. Moreover we compute the error with respect to relative estimates, i.e.
B k − C k,j for every B k ∈R h and for every C k,j ∈ R N,h,T,j ∪ R N,h,S,j , for every k ∈ {1, . . . , n h }, for every j ∈ {1, . . . , 5} and for every h ∈ {E, W }.
In conclusion we report some statistics about the original and the relative errors per pixel, stored respectively in vector epO and epR for the Europe dataset. Comparing Table 8 with Table 10, we notice that the average original and relative error per pixel are greater in the second case. In other words, starting from a tensor with repeated NIR band is not convenient for computing Rao index over Europe elements. Besides we underline that the minimum relative and original error are realised when the starting tensor has RED band repeated. We can not affirm nothing about the maximum original and relative error, since for i j ∈ {10, 100} it is lower when the repeated band is the NIR, while for i j ∈ {50, 500, 1000} when it is repeated RED band. Lastly we highlight as always for the Europe dataset that the average original error is lower than the relative one. As in the repeated RED case, we remark that ST-HOSVD is convenient when the first two components of multilinear are strictly lower than 500, otherwise T-HOSVD provides better results. Then we Relative approximation

(B)
Component rank 10 compression Component rank 50 compression Component rank 100 compression Component rank 500 compression Component rank 1000 compression  Table 11. Statistics for Rao index over N ∈ N N,W,t,j .
In this case comparing Table 9 with Table 11, we notice that starting tensors with repeated NIR band provides on average better results. Indeed for both the decomposition techniques at the fifth approximation the average original and relative error per pixel is around 17.5% in the first case and 16.5% in the second one. This 1% difference between the RED and NIR case is present also in the minimum and maximum relative and original error per pixel for both the decomposition techniques. Besides we underline that also in this case ST-HOSVD is convenient when the first two multilinear components are lower than 100, otherwise T-HOSVD is preferable. In conclusion to these analysis we want to empathise that an average error of 16.5% per pixel is not much, if we remind that we use only 16.4% of the total information available. Lastly as before, we display the Rao index computed for December 2014 and May 2014, starting from tensors with repeated NIR bands.
Example 4.23. Firstly we remark that December 2014 realises the minimum original and relative error, also starting from tensors with repeated NIR band. In Figure 11 there are approximated Rao estimates at different multilinear rank, for this NIR repeated case. As before, our eyes hardly perceive differences between Rao computed over NASA or self-made NDVI and over the last three approximated NDVIs.
Example 4.24. Firstly we remark that May 2014 do not realise the minimum original and relative error starting from tensors with repeated NIR band. Therefore these are not the worst approximations in the Earth dataset, in this second case. In Figure 12 there are approximated Rao estimates at different multilinear rank, for this NIR repeated case. As in Figure 10 the problematic areas are the North pole and the Pacific ocean. Indeed we notice that there is a biodiversity overestimation in the Arctic regions, which decreases for increasing multilinear rank components. Moreover even in the fifth approximation there are missing island in the Pacific ocean. However we can not easily notice much difference with respect to Figure 10.

Conclusion
We have shown that the presented approach is extremely convenient for Rényi index estimation to save storage memory. As reported in Tables 4,6, 5 and 7 the average error per pixel is around 5.5% for Earth dataset and is close to 7.6% for the Europe's one when using respectively 16.4% and 14.8% of the total tensor information. Moreover we have seen that starting from tensor with repeated NIR bands is more convenient than with RED repeated, in the Rényi case. For Rao case starting with twice NIR rasters is more convenient only for Earth dataset. Lastly we have remarked Component rank 10 compression Component rank 50 compression Component rank 100 compression that in Rényi case T-HOSVD makes on biodiversity estimation lower error than ST-HOSVD when the first two multilinear rank components are relatively low.
In the Rao case we face greater on average original and relative error per pixel. Indeed in Tables  8 and 9 we have an average error per pixel close to 19.5% and 17% for Europe and for Earth dataset respectively. In this case we have noticed that the roles of T-HOSVD and ST-HOSVD are inverted.
Indeed T-HOSVD makes on biodiversity estimation lower error than ST-HOSVD with the first two multilinear rank components are relatively great.
In conclusion we believe that the presented work might help ecologists in their remote sensed plant biodiversity estimation. Indeed fixed a certain accuracy, they can compress through T-HOSVD and ST-HOSVD the tensors with NIR and RED bands, to save storage memory and at the same time computing with fixed tolerance the biodiversity estimates. Lastly if it was possible to directly download decomposed tensors, ecologists would benefit of a significant computer memory saving. 6. Appendix 6.1. Rao's and Rényi's codes. Remind that for every raster of size (m, n), fixed l the side of the moving window, the biodiversity index is computed l × m × n times. To speed up the entire work we decide to implement a parallel version of the biodiversity index algorithm. When a computer executes a parallel functions its cores perform independently different tasks at the same time. In our case we want that each core of the used machine works on different position of the moving window. Consequently to implement this mechanism we needed a parallel computing library compatible with Python, the chosen programming language. We preferred Joblib, cf. [Var ], since its ease of use. Other two used libraries are itertools, to create iterable elements and spatial from SciPy, to compute the distance element wise for two matrices. The main Rao computation code is #### computation import numpy as np import scipy from scipy import spatial import itertools #### parallelisation import joblib from joblib import Parallel, delayed import multiprocessing def raop(rw): def raout (cl, rw = rw, rasterm = rasterm, missing = missing, w = w, distance_m → = distance_m): tw_labels, tw_values = np.unique(rasterm[(rw-w):(rw+w+1),(cl-w):(cl+w+1)], → return_counts=True) if len(np.where(tw_labels == missing)) != 0: tw_values = np.delete( tw_values, np.where(tw_labels == missing)) tw_labels = np.delete( tw_labels, np.where(tw_labels == missing)) if len(tw_values) > 1: d1 = spatial.distance.cdist(np.diag(tw_labels), np.diag(tw_labels), → distance_m) p = tw_values/np.sum(tw_values) p1 = np.zeros((len(tw_values),len(tw_values))) comb = np.array([x[0]*x[1] for x in list(itertools.combinations(p, 2))]) p1[np.triu_indices(len(tw_values), k=1)] = comb p1[np.tril_indices(len(tw_values), k=-1)] = comb return ((np.sum(np.multiply(p1,d1)))) elif len(tw_values) == 1: return (((0))) else: return ((missing)) Raout = Parallel(n_jobs = NcCores)(delayed(raout)(cl) for cl in range(w,c-w)) return (Raout) out[w:(r-w), w:(c-w)] = (np.asarray(Parallel(n_jobs = NcCores)(delayed(raop)(rw) → for rw in range(w,(r-w)))).reshape(r-2*w,c-2*w)) Since the moving window scrolls over rows and over columns, we define two parallelised functions one inside the others, raop whose variable is just the row index and raout. This second function takes as variable the column index and has other set parameters: rasterm the raster, rw the row index, missing a value used in the raster when data are not present, w the window side and a function distance_m to compute the distance between raster values. The first step is storing the values and their frequencies of the raster area covered by the moving window in the arrays tw_labels and tw_values respectively. Then we check if there is the missing value in the considered area: if it is present, we remove it and its frequency from the storing arrays. Next if there are at least two different elements in the tw_values array, we compute their distance with the function spatial.distance.cdist and save the result in the d1 matrix. In p we put the relative frequencies, obtained from the absolute ones. We define a 0 matrix p1, with the same size of d1. We compute and store in comb all the possible combination between different elements of vector p. We assign these values to the upper and lower triangular part of the matrix p1. Lastly we return the sum of the elements of the product matrix between d1 and p1. Notice that the product matrix will coherently have diagonal null elements, since the distance between a pixel value and itself is 0.
If there is just one element in vector tw_values, we return 0 since the difference between a pixel value and itself is 0. Otherwise, if tw_values is an empty vector, we return missing. Inside function raop we define function raout and we do the first Parallel call. Outside it we do the second Parallel call. Notice that this function is thought to compute Rao's index when the moving window is completely contained into the raster. When this condition is not satisfied, to speed up computation, there is a special implementation of Rao's, available here [Ian19]. With a similar approach, we implemented Rényi's index, whose core is #### computation import numpy as np import scipy #### parallelisation import joblib from joblib import Parallel, delayed import multiprocessing def IndexOP (rw): def IndexOut (cl, rw = rw, rasterm = rasterm, missing = missing, w = w, alpha= → alpha, base = base): tw_labels, tw_values = np.unique(rasterm[(rw-w):(rw+w+1),(cl-w):(cl+w+1)], → return_counts=True) if len(np.where(tw_labels == missing)[0]) != 0: tw_values = np.delete( tw_values, np.where(tw_labels == missing)) if len(tw_values) != 0: p = tw_values/np.sum(tw_values) if np.log(np.sum(p**alpha)) == 0: return(0) else: return((1/(1-alpha)) * np.log(np.sum(p**alpha)) / np.log(base)) else: return (missing) Index_Out = Parallel(n_jobs = NcCores)(delayed(IndexOut)(cl) for cl in range(w, → c-w)) return (Index_Out) out[w:(r-w), w:(c-w)] = np.asarray(Parallel(n_jobs = NcCores)(delayed(IndexOP)(rw → ) for rw in range(w,r-w))).reshape(r-2*w,c-2*w) Similarly we declare function IndexOP, with just row index as variable and inside it function IndexOut. This one takes as parameters some equals to raout's ones. The variable proper of IndexOut are alpha, the Rényi's parameter, and base, the logarithm base. The code structure is similar to raout. We compute and store the values, present inside the raster area covered by the moving window, and their frequencies. We check and in case delete missing value and its frequency. If there is still an element into tw_values, we compute the relative frequencies and return the Rényi index. We test if the index final value is 0 and return 0 to avoid sign problem. If tw_values is an empty vector, we return missing value. We do a first Parallel call of IndexOut inside IndexOP and a second one of IndexOP outside it. Even in this case the presented code will work only when the moving window is entirely contained in the raster. To speed up the computations, we develop a special version of this code which works when the moving window is not fully contained. It is available here [Ian19].
6.2. Approximation codes for T-HOSVD. In order to implement the approximation Algorithm 1 of T-HOSVD we need a python library expressly developed to manage tensors. Moreover we also want that this library to interact properly with NumPy the python most used library for numerical computation, cf. [Oli06]. Therefore we choose TensorLy, cf. [KPAP19]. This recently developed library, firstly presented in 2016, wants to make tensor study and manipulation easy and accessible. Besides their creators projected TensorLy to perfectly match with other famous python libraries, as NumPy, SciPy. They developed most of the main tensor operations and related functions. To compute the singular value decomposition of tensors flattening, the central step in T-HOSVD and ST-HOSVD algorithm, we need a python function able to manage huge arrays. This is not a trivial task, since development of python functions is left to singular initiative. We decide to use svds function from SciPy sparse linear algebra function, cf. [JOP ]. This implementation of SVD takes advantage of matrix sparsity in performing the matrix-vector multiplication. The final implementation of T-HOSVD is import numpy as np import tensorly as tl from tensorly import decomposition as decompose from tensorly import tenalg as Tla import scipy from scipy import spatial from scipy.sparse.linalg import svds This function takes as input variable a tensor T and a list or a tuple, whose values are the target multilinear rank components. Inside the T_hosvd we declare an empty list, L, in which at the i-th step we store matrix U_i from the thin SVD of the i-th flattening for every i ∈ {1, . . . , d}. The dim variable is a tuple containing the size of tensor T. Then for each direction we compute the flattening and its SVD. After the for loop, we get the core tensor with the multilinear product between a list of matrices, our projectors, and the original tensor T. Notice that there is the option projector, with True as default value, to return a list L containing the projectors matrices together with the standard result, i.e. the core tensor.
6.3. Approximation codes for ST-HOSVD. The ST-HOSVD implementation of Algorithm 2 is import numpy as np import tensorly as tl from tensorly import decomposition as decompose from tensorly import tenalg as Tla import scipy from scipy import spatial from scipy.sparse.linalg import svds def ST_hosvd(T, rango, projector = The input arguments of ST-HOSVD and T-HOSVD implementations are the same. The first difference with T_hosvd is the initial declaration of core tensor, set equal to T. Next we highlight that only if the projector variable is True, we initialize an empty list to store the projectors matrices. Another difference is the computation of the core tensor inside the for loop with the component wise product between matrix and tensor, fixed a certain direction. The basic output is still the final core tensor.