The efficient solution of large-scale multiterm linear matrix equations is a challenging task in numerical linear algebra, and it is a largely open problem. We propose a new iterative scheme for symmetric and positive definite operators, significantly advancing methods such as truncated matrix-oriented Conjugate Gradients (CG). The new algorithm capitalizes on the low-rank matrix format of its iterates by fully exploiting the subspace information of the factors as iterations proceed. The approach implicitly relies on orthogonality conditions imposed over much larger subspaces than in CG, unveiling insightful connections with subspace projection methods. The new method is also equipped with memory-saving strategies. In particular, we show that for a given matrix Y, the action L(Y) in low rank format may not be evaluated exactly due to memory constraints. This problem is often underestimated, though it will eventually produce Out-of-Memory breakdowns for a sufficiently large number of terms. We propose an ad-hoc randomized range-finding strategy that appears to fully resolve this shortcoming. Experimental results with typical application problems illustrate the potential of our approach over various methods developed in the recent literature.
@techreport{palitta2025SsCG,title={A subspace-conjugate gradient method for linear matrix equations},author={Palitta, Davide and Iannacito, Martina and Simoncini, Valeria},year={2025},month=jan,eprint={2501.02938},archiveprefix={arXiv},primaryclass={math.NA},keywords={report},}
2022
HAL
On some orthogonalization schemes in Tensor Train format
Olivier Coulaud, Luc Giraud, and Martina Iannacito
In the framework of tensor spaces, we consider orthogonalization kernels to generate an orthogonal basis of a tensor subspace from a set of linearly independent tensors. In particular, we investigate numerically the loss of orthogonality of six orthogonalization methods, namely Classical and Modified Gram-Schmidt with (CGS2, MGS2) and without (CGS, MGS) re-orthogonalization, the Gram approach, and the Householder transformation. To tackle the curse of dimensionality, we represent tensor with low rank approximation using the Tensor Train (TT) formalism, and we introduce recompression steps in the standard algorithm outline through the TT-rounding method at a prescribed accuracy. After describing the algorithm structure and properties, we illustrate numerically that the theoretical bounds for the loss of orthogonality in the classical matrix computation round-off analysis results are maintained, with the unit round-off replaced by the TT-rounding accuracy. The computational analysis for each orthogonalization kernel in terms of the memory requirement and the computational complexity measured as a function of the number of TT-rounding, which happens to be the computational most expensive operation, completes the study.
@techreport{Coulaud2022b,title={{On some orthogonalization schemes in Tensor Train format}},author={Coulaud, Olivier and Giraud, Luc and Iannacito, Martina},url={hal.archives-ouvertes.fr/hal-03850387},type={Research Report},number={{RR-9491}},institution={Inria center at the University of Bordeaux},year={2022},month=nov,keywords={report},hal_id={hal-03850387},hal_version={v2},}
HAL
The backward stable variants of GMRES in variable accuracy
Emmanuel Agullo, Olivier Coulaud, Luc Giraud, and
3 more authors
In the context where the representation of the data is decoupled from the arithmetic used to process them, we investigate the backward stability of two backward-stable implementations of the GMRES method, namely the so-called Modified Gram-Schmidt (MGS) and the Householder variants. Considering data may be compressed to alleviate the memory footprint, we are interested in the situation where the leading part of the rounding error is related to the data representation. When the data representation of vectors introduces componentwise perturbations, we show that the existing backward stability analyses of MGS-GMRES and Householder-GMRES still apply. We illustrate this backward stability property in a practical context where an agnostic lossy compressor is employed and enables the reduction of the memory requirement to store the orthonormal Arnoldi basis or the Householder reflectors. Although technical arguments of the theoretical backward stability proofs do not readily apply to the situation where only the normwise relative perturbations of the vector storage can be controlled, we show experimentally that the backward stability is maintained; that is, the attainable normwise backward error is of the same order as the normwise perturbations induced by the data storage. We illustrate it with numerical experiments in two practical different contexts. The first one corresponds to the use of an agnostic compressor where vector compression is controlled normwise. The second one arises in the solution of tensor linear systems, where low-rank tensor approximations based on Tensor-Train is considered to tackle the curse of dimensionality.
@techreport{Agullo2022,title={{The backward stable variants of GMRES in variable accuracy}},author={Agullo, Emmanuel and Coulaud, Olivier and Giraud, Luc and Iannacito, Martina and Marait, Gilles and Schenkels, Nick},url={https://hal.science/hal-03776837},type={Research Report},number={RR-9483},pages={1-77},institution={Inria center at the University of Bordeaux},year={2022},month=sep,keywords={report},hal_id={hal-03776837},hal_version={v3},}
HAL
A robust GMRES algorithm in Tensor Train format
Olivier Coulaud, Luc Giraud, and Martina Iannacito
We consider the solution of linear systems with tensor product structure using a GMRES algorithm. In order to cope with the computational complexity in large dimension both in terms of floating point operations and memory requirement, our algorithm is based on low-rank tensor representation, namely the Tensor Train format. In a backward error analysis framework, we show how the tensor approximation affects the accuracy of the computed solution. With the bacwkward perspective, we investigate the situations where the (𝑑+1)-dimensional problem to be solved results from the concatenation of a sequence of 𝑑-dimensional problems (like parametric linear operator or parametric right-hand side problems), we provide backward error bounds to relate the accuracy of the (𝑑+1)-dimensional computed solution with the numerical quality of the sequence of 𝑑-dimensional solutions that can be extracted form it. This enables to prescribe convergence threshold when solving the (𝑑+1)-dimensional problem that ensures the numerical quality of the 𝑑-dimensional solutions that will be extracted from the (𝑑+1)-dimensional computed solution once the solver has converged. The above mentioned features are illustrated on a set of academic examples of varying dimensions and sizes.
@techreport{Coulaud2022a,title={{A robust GMRES algorithm in Tensor Train format}},author={Coulaud, Olivier and Giraud, Luc and Iannacito, Martina},url={https://hal.science/hal-03776529v4},type={Research Report},number={RR-9484},pages={1-48},institution={Inria center at the University of Bordeaux},year={2022},month=sep,hal_id={hal-03776529},hal_version={v4},keywords={report},}
2021
HAL
Extension of Correspondence Analysis to multiway data-sets through High Order SVD: a geometric framework
Olivier Coulaud, Alain Franc, and Martina Iannacito
This paper presents an extension of Correspondence Analysis (CA) to tensors through High Order Singular Value Decomposition (HOSVD) from a geometric viewpoint. Correspondence analysis is a well-known tool, developed from principal component analysis, for studying contingency tables. Different algebraic extensions of CA to multi-way tables have been proposed over the years, nevertheless neglecting its geometric meaning. Relying on the Tucker model and the HOSVD, we propose a direct way to associate with each tensor mode a point cloud. We prove that the point clouds are related to each other. Specifically using the CA metrics we show that the barycentric relation is still true in the tensor framework. Finally two data sets are used to underline the advantages and the drawbacks of our strategy with respect to the classical matrix approaches.
@techreport{Coulaud2021,title={{Extension of Correspondence Analysis to multiway data-sets through High Order SVD: a geometric framework}},author={Coulaud, Olivier and Franc, Alain and Iannacito, Martina},url={https://hal.science/hal-03418404},type={Research Report},number={RR-9429},institution={Inria center at the University of Bordeaux; INRAE},year={2021},month=nov,keywords={report},hal_id={hal-03418404},hal_version={v1},}
Peer-reviewed journal
2025
ETNA
A note on TT-GMRES for the solution of parametric linear systems
Olivier Coulaud, Luc Giraud, and Martina Iannacito
Electronic Transactions on Numerical Analysis, Jan 2025
We study the solution of linear systems with tensor product structure using the Generalized Minimal RESidual (GMRES) algorithm. To manage the computational complexity of high-dimensional problems, our approach relies on low-rank tensor representation, focusing specifically on the Tensor Train format. We implement and experimentally study the TT-GMRES algorithm. Our analysis bridges the heuristic methods proposed for TT-GMRES by Dolgov [Russian J. Numer. Anal. Math. Modelling, 28 (2013), pp. 149-172] and the theoretical framework of inexact GMRES by Simoncini and Szyld [SIAM J. Sci. Comput. 25 (2003), pp. 454-477]. This approach is particularly relevant in a scenario where a (d+1)-dimensional problem arises from concatenating a sequence of d-dimensional problems, as in the case of a parametric linear operator or parametric right-hand-side formulation. Thus, we provide backward error bounds that link the accuracy of the computed (d+1)-dimensional solution to the numerical quality of the extracted d-dimensional solutions. This facilitates the prescription of a convergence threshold ensuring that the d-dimensional solutions extracted from the (d+1)-dimensional result have the desired accuracy once the solver converges. We illustrate these results with academic examples across varying dimensions and sizes. Our experiments indicate that TT-GMRES retains the theoretical rounding-error properties observed in matrix-based GMRES.
@article{Coulaud2025Note,title={A note on TT-GMRES for the solution of parametric linear systems},author={Coulaud, Olivier and Giraud, Luc and Iannacito, Martina},journal={Electronic Transactions on Numerical Analysis},publisher={Kent State University Library},volume={62},pages={163-187},year={2025},month=jan,doi={10.1553/etna_vol62s163},keywords={intjour},tag={journal},}
2021
BUMI
High order singular value decomposition for plant diversity estimation
Alessandra Bernardi, Martina Iannacito, and Duccio Rocchini
Bollettino dell’Unione Matematica Italiana, Dec 2021
We propose a new method to estimate plant diversity with Rényi and Rao indexes through the so called High Order Singular Value Decomposition (HOSVD) of tensors. Starting from NASA multi-spectral images we evaluate diversity and we compare original diversity estimates with those realized via the HOSVD compression methods for big data. Our strategy turns out to be extremely powerful in terms of memory storage and precision of the outcome. The obtained results are so promising that we can support the efficiency of our method in the ecological framework.
@article{Bernardi2021,author={Bernardi, Alessandra and Iannacito, Martina and Rocchini, Duccio},title={High order singular value decomposition for plant diversity estimation},journal={Bollettino dell'Unione Matematica Italiana},year={2021},month=dec,day={01},volume={14},number={4},pages={557-591},issn={2198-2759},doi={10.1007/s40574-021-00300-w},keywords={natjour},tag={journal},}
Comm. Eco.
Measuring diversity from space: a global view of the free and open source rasterdiv R package under a coding perspective
Elisa Thouverai, Matteo Marcantonio, Giovanni Bacaro, and
7 more authors
The variation of species diversity over space and time has been widely recognised as a key challenge in ecology. However, measuring species diversity over large areas might be difficult for logistic reasons related to both time and cost savings for sampling, as well as accessibility of remote ecosystems. In this paper, we present a new R package - rasterdiv - to calculate diversity indices based on remotely sensed data, by discussing the theory behind the developed algorithms. Obviously, measures of diversity from space should not be viewed as a replacement of in situ data on biological diversity, but they are rather complementary to existing data and approaches. In practice, they integrate available information of Earth surface properties, including aspects of functional (structural, biophysical and biochemical), taxonomic, phylogenetic and genetic diversity. Making use of the rasterdiv package can result useful in making multiple calculations based on reproducible open source algorithms, robustly rooted in Information Theory.
@article{Thouverai2021,author={Thouverai, Elisa and Marcantonio, Matteo and Bacaro, Giovanni and Re, Daniele Da and Iannacito, Martina and Marchetto, Elisa and Ricotta, Carlo and Tattoni, Clara and Vicario, Saverio and Rocchini, Duccio},title={Measuring diversity from space: a global view of the free and open source rasterdiv R package under a coding perspective},journal={Community Ecology},year={2021},month=apr,day={01},volume={22},number={1},pages={1-11},issn={1588-2756},doi={10.1007/s42974-021-00042-x},keywords={intjour},tag={journal},}
Glob. Ecol.
From zero to infinity: Minimum to maximum diversity of the planet by spatio-parametric Rao’s quadratic entropy
Duccio Rocchini, Matteo Marcantonio, Daniele Da Re, and
20 more authors
The majority of work done to gather information on the Earth’s biodiversity has been carried out using in-situ data, with known issues related to epistemology (e.g., species determination and taxonomy), spatial uncertainty, logistics (time and costs), among others. An alternative way to gather information about spatial ecosystem variability is the use of satellite remote sensing. It works as a powerful tool for attaining rapid and standardized information. Several metrics used to calculate remotely sensed diversity of ecosystems are based on Shannon’s information theory, namely on the differences in relative abundance of pixel reflectances in a certain area. Additional metrics like the Rao’s quadratic entropy allow the use of spectral distance beside abundance, but they are point descriptors of diversity, that is they can account only for a part of the whole diversity continuum. The aim of this paper is thus to generalize the Rao’s quadratic entropy by proposing its parameterization for the first time. Innovation The parametric Rao’s quadratic entropy, coded in R, (a) allows the representation of the whole continuum of potential diversity indices in one formula, and (b) starting from the Rao’s quadratic entropy, allows the explicit use of distances among pixel reflectance values, together with relative abundances. Main conclusions The proposed unifying measure is an integration between abundance- and distance-based algorithms to map the continuum of diversity given a satellite image at any spatial scale. Being part of the rasterdiv R package, the proposed method is expected to ensure high robustness and reproducibility.
@article{Rocchini2021b,author={Rocchini, Duccio and Marcantonio, Matteo and Da Re, Daniele and Bacaro, Giovanni and Feoli, Enrico and Foody, Giles M. and Furrer, Reinhard and Harrigan, Ryan J. and Kleijn, David and Iannacito, Martina and Lenoir, Jonathan and Lin, Meixi and Malavasi, Marco and Marchetto, Elisa and Meyer, Rachel S. and Moudry, Vítězslav and Schneider, Fabian D. and Šímová, Petra and Thornhill, Andrew H. and Thouverai, Elisa and Vicario, Saverio and Wayne, Robert K. and Ricotta, Carlo},title={From zero to infinity: Minimum to maximum diversity of the planet by spatio-parametric Rao’s quadratic entropy},journal={Global Ecology and Biogeography},volume={30},number={5},pages={1153-1162},keywords={intjour},doi={10.1111/geb.13270},url={onlinelibrary.wiley.com/doi/abs/10.1111/geb.13270},year={2021},month=mar,tag={journal},}
MEE
rasterdiv—An Information Theory tailored R package for measuring ecosystem heterogeneity from space: To the origin and back
Duccio Rocchini, Elisa Thouverai, Matteo Marcantonio, and
32 more authors
Ecosystem heterogeneity has been widely recognized as a key ecological indicator of several ecological functions, diversity patterns and change, metapopulation dynamics, population connectivity or gene flow. In this paper, we present a new R package—rasterdiv—to calculate heterogeneity indices based on remotely sensed data. We also provide an ecological application at the landscape scale and demonstrate its power in revealing potentially hidden heterogeneity patterns. The rasterdiv package allows calculating multiple indices, robustly rooted in Information Theory, and based on reproducible open-source algorithms.
@article{Rocchini2021a,author={Rocchini, Duccio and Thouverai, Elisa and Marcantonio, Matteo and Iannacito, Martina and Da Re, Daniele and Torresani, Michele and Bacaro, Giovanni and Bazzichetto, Manuele and Bernardi, Alessandra and Foody, Giles M. and Furrer, Reinhard and Kleijn, David and Larsen, Stefano and Lenoir, Jonathan and Malavasi, Marco and Marchetto, Elisa and Messori, Filippo and Monkeywordhi, Alessandro and Moudrý, Vítězslav and Naimi, Babak and Ricotta, Carlo and Rossini, Micol and Santi, Francesco and Santos, Maria J. and Schaepman, Michael E. and Schneider, Fabian D. and Schuh, Leila and Silvestri, Sonia and Ŝímová, Petra and Skidmore, Andrew K. and Tattoni, Clara and Tordoni, Enrico and Vicario, Saverio and Zannini, Piero and Wegmann, Martin},title={rasterdiv—An Information Theory tailored R package for measuring ecosystem heterogeneity from space: To the origin and back},journal={Methods in Ecology and Evolution},volume={12},number={6},pages={1093-1102},keywords={intjour},doi={10.1111/2041-210X.13583},url={besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.13583},year={2021},month=feb,tag={journal},}
Other works
2022
thesis
Numerical linear algebra and data analysis in large dimensions using tensor format
Martina Iannacito
Inria center at the University of Bordeaux, France, Dec 2022
This work aims to establish which theoretical properties of classical linear algebra techniques, developed in two different contexts, namely numerical linear algebra and data analysis, are maintained and which are lost when extended to low-rank tensors. In the numerical linear algebra part, we experimentally study the effects of rounding errors on an iterative solver (Generalised Minimal RESidual) and several orthogonalization methods (Classical and Modified Gram-Schmidt and Householder transformation, among others) when extended to tensors by the Tensor Train (TT) formalism. In all the considered algorithms, we introduce additional rounding steps, with the TT-rounding compression algorithm, to cope with memory constraints, which are always crucial when dealing with tensors. Our experimental tests suggest that for these algorithms, the classical bounds due to the propagation of rounding errors remain valid, replacing the accuracy of the arithmetic with that of the TT-rounding algorithm. In the data analysis part, we geometrically study the generalization of correspondence analysis to multi-way tables by the Tucker tensor decomposition technique, thus contributing to the understanding of Multi-Way Correspondence Analysis (MWCA). Examples of MWCA applied to datasets complement the theoretical results. In particular, we perform the MWCA on an original ecological dataset provided to us in the framework of the Malabar project of the BioGeCo (INRAE) and Pleiade (Inria) teams.
@phdthesis{Iannacito2022,author={Iannacito, Martina},year={2022},month=dec,title={Numerical linear algebra and data analysis in large dimensions using tensor format},url={theses.fr/s349733},school={Inria center at the University of Bordeaux, France},keywords={thesis},}
2019
thesis
HOSVD FOR MULTISPECTRAL IMAGES. A numerical approach to the plant biodiversity estimate.
Remote sensing has already proved its power in estimating plant biodiversity. Since plants reflect and absorb in a specific way sunlight, ecologists use multispectral photos to estimate plant biodiversity of a certain area through biodiversity indexes, e.g., Rényi’s and Rao’s index. Working with these big data is often expensive in terms of storage memory requirements. A possible solution is compressing multispectral data through tensor techniques, as High Order Singular Value Decomposition. We apply two recent variants of HOSVD to spectral images and we compute biodiversity indexes over the compressed data generated. For the Rényi index compressing data is extremely advantageous, since the information loss in estimating biodiversity is very low. For Rao’s index the results are not as satisfactory, even if they are promising. Using approximatively 15% of the original data information, in the first case we estimate an average error per pixel between 5.5% and 7.5%, while in the second case the average error per pixel is between 17% and 19%.
@mastersthesis{Iannacito2019,author={Iannacito, Martina},year={2019},month=jul,title={{HOSVD FOR MULTISPECTRAL IMAGES. A} numerical approach to the plant biodiversity estimate.},url={webapps.unitn.it/Biblioteca/it/Web/RichiestaConsultazioneTesi/365719},school={University of Trento, Italy},keywords={thesis},type={thesis},}