Paper deep dive
Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection
Florian Grivet, Louise Travรฉ-Massuyรจs
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 98%
Last extracted: 3/22/2026, 5:45:22 AM
Summary
This paper evaluates three matrix inversion update methodsโDirect Inversion (DI), Iterative Sherman-Morrison (ISM), and Woodbury Matrix Identity (WMI)โfor streaming outlier detection using the Christoffel function. By deriving theoretical computational costs and validating them with Python simulations, the authors provide a quantitative rule for selecting the optimal method based on matrix size and update rank: ISM is best for rank-1 updates, WMI is efficient for small updates relative to matrix size, and DI is preferred for larger updates.
Entities (5)
Relation Signals (4)
Direct Inversion โ comparedwith โ Iterative Sherman-Morrison
confidence 100% ยท This technical note aims to compare three different updating methods: Direct Inversion (DI), Iterative Sherman-Morrison (ISM), and Woodbury Matrix Identity (WMI)
Woodbury Matrix Identity โ excelsfor โ small updates
confidence 100% ยท WMI excels for small updates relative to matrix size
Iterative Sherman-Morrison โ optimalfor โ rank-1 updates
confidence 100% ยท ISM is optimal for rank-1 updates
DyCF โ uses โ Christoffel function
confidence 100% ยท Ducharlet et al. [2] introduce DyCF, a novel outlier detection algorithm for data streams leveraging the Christoffel function.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Outlier detection identifies data points that deviate significantly from expected patterns, revealing anomalies that may require special attention. Incorporating online learning further improves accuracy by continuously updating the model to reflect the most recent data. When employing the Christoffel function as an outlier score, online learning requires updating the inverse of a matrix following a rank-k update, given the initial inverse. Surprisingly, there is no consensus on the optimal method for this task. This technical note aims to compare three different updating methods: Direct Inversion (DI), Iterative Sherman-Morrison (ISM), and Woodbury Matrix Identity (WMI), to identify the most suitable approach for different scenarios. We first derive the theoretical computational costs of each method and then validate these findings through comprehensive Python simulations run on a CPU. These results allow us to propose a simple, quantitative, and easy-to-remember rule that can be stated qualitatively as follows: ISM is optimal for rank-1 updates, WMI excels for small updates relative to matrix size, and DI is preferable otherwise. This technical note produces a general result for any problem involving a matrix inversion update. In particular, it contributes to the ongoing development of efficient online outlier detection techniques.
Tags
Links
- Source: https://arxiv.org/abs/2603.16697v1
- Canonical: https://arxiv.org/abs/2603.16697v1
Trouble viewing inline? Open PDF directly โ
Full Text
44,712 characters extracted from source content.
Expand or collapse full text
Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection Florian Grivet a,b,โ , Louise Travรฉ-Massuyรจs b,1 a CNES, 18 Av. Edouard Belin, Toulouse, 31400, Toulouse, France b LAAS-CNRS, University of Toulouse, CNRS, 7 Av. du Colonel Roche, 31400, Toulouse, France A R T I C L E I N F O Keywords: Outlier Detection Christoffel function Rank update Matrix Inversion Sherman-Morrison Woodbury Matrix Identity Computational cost Abstract Outlier detection identifies data points that deviate significantly from expected patterns, revealing anomalies that may require special attention. Incorporating online learning further improves accuracy by continuously updating the model to reflect the most recent data. When employing the Christoffel function as an outlier score, online learning requires updating the inverse of a matrix following a rank-ํ update, given the initial inverse. Surprisingly, there is no consensus on the optimal method for this task. This technical note aims to compare three different updating methods: Direct Inversion (DI), Iterative Sherman-Morrison (ISM), and Woodbury Matrix Identity (WMI), to identify the most suitable approach for different scenarios. We first derive the theoretical computational costs of each method and then validate these findings through comprehensive Python simulations run on a CPU. These results allow us to propose a simple, quantitative, and easy-to-remember rule that can be stated qualitatively as follows: ISM is optimal for rank-1 updates, WMI excels for small updates relative to matrix size, and DI is preferable otherwise. This technical note produces a general result for any problem involving a matrix inversion update. In particular, it contributes to the ongoing development of efficient online outlier detection techniques. 1. Introduction The detection of outliers in data streams has become increasingly important in a wide range of applications, from fraud detection to quality control in manufacturing. In such settings, data arrive sequentially and often at high rates, mak- ing online learning approaches particularly attractive. These methods continuously update models as new observations become available, allowing anomaly detection systems to adapt to evolving data distributions while maintaining strong performance over time. Among recent approaches to anomaly detection in data streams [10, 17], Ducharlet et al. [2] introduce an outlier scoring mechanism based on the Christoffel function (CF) [6]. This score is defined in terms of the inverse of a symetric positive definite moment matrix associated with the data. In streaming settings, this matrix is updated sequentially via rank-ํ corrections as new observations arrive. While the resulting CF scores are invariant to the specific in- verse update strategy employed โ up to numerical precision โ the choice of update method has a substantial impact on computational cost, numerical stability, and scalability. These considerations are critical in streaming settings, where efficiency directly constrains real-time applicability. Several strategies are available for updating matrix in- verses after rank-ํ corrections, including Direct Inversion (DI), Iterative Sherman-Morrison (ISM) [11], and the Wood- bury Matrix Identity (WMI) [15]. Despite their widespread โ Corresponding author florian.grivet@cnes.fr (F. Grivet); louise@laas.fr (L. Travรฉ-Massuyรจs) ORCID(s): 0009-0007-7096-3258 (F. Grivet); 0000-0002-5322-8418 (L. Travรฉ-Massuyรจs) 1 Head of ANITI chair ADDX use, there is currently no clear quantitative guidance on which method is preferable under different conditions, such as varying matrix size ํ or update rank ํ. This lack of guidance can lead to inefficient implementations that unnec- essarily limit the practicality of CF-based anomaly detection in streaming environments. This article is presented as a technical note aimed at addressing this gap by comparing inverse update strategies in the context of Christoffel-function-based outlier detection. It does not propose a new scoring model, but instead analyzes how different matrix update methods may affect computa- tional efficiency. The contributions of this technical note are summarized as follows: - This note introduces the Christoffel function, explores its key properties, and presents DyCF, a frugal stream- ing outlier detection method inspired by these founda- tions, that motivates the paperโs work. - This work derives the computational costs of three matrix inverse rank-k update methods, namely Direct Inversion, Iterative Sherman-Morrison, and Wood- bury Matrix Identity. - Summarizing and comparing the three theoretical computational costs yields a unified reference, and the findings are validated through comprehensive Python simulations run on a CPU. - As a key takeaway, this note offers a simple, quantita- tive, and easy-to-remember rule, expressed in terms of the matrix dimension ํ and the update rank ํ, for selecting among the three rank-ํ matrix inverse update strategies implemented in Python on CPU. Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 1 of 10 arXiv:2603.16697v1 [cs.LG] 17 Mar 2026 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection The technical note is organized as follows. Section 2 briefly reviews the Christoffel function, highlighting the properties relevant to anomaly detection. Section 3 then discusses its use in streaming outlier detection, empha- sizing the need for efficient inverse updates under rank-ํ corrections. Section 4 presents the Direct Inversion, Iter- ative ShermanโMorrison, and Woodbury Matrix Identity approaches, detailing their algorithms and theoretical com- putational costs. Section 5 compares the theoretical compu- tational costs of the three methods. Section 6 analyzes theo- retical predictions with empirical results to derive practical implementations. Finally, Section 7 summarizes the main findings, resumes practical guidance, discusses limitations, and points at interesting topics for future work. 2. The Christoffel function The CF originates from the theory of approximation and orthogonal polynomials. Lasserre and Pauwels [6] demonstrated that the CF is related to a sum-of-squares (SOS) polynomial whose sublevel set effectively captures the shape of a dataset. Building on this discovery, Lasserre and Pauwels [7], and Lasserre et al. [8] developed a com- prehensive theoretical framework for data analysis, and in particular, anomaly detection. This section explores key properties of both the theoret- ical CF (referred to as the population Christoffel function) and its empirical counterpart. 2.1. The population Christoffel function Let ํฑ = ( ํฅ 1 ,ํฅ 2 ,โฏ,ํฅ ํ ) โโ ํ . To define polynomials, we adopt the multi-index notation ํผ = ( ํผ ํ ) ํ=1...ํ โโ ํ , such that the monomial ํฑ ํผ of total degree ํํํ(ํฑ ํผ ) = |ํผ| = โ ํ ํ=1 ํผ ํ is given by ํฑ ํผ = ํฅ ํผ 1 1 ํฅ ํผ 2 2 โฏํฅ ํผ ํ ํ . In short form, we denote the set of ํ-variate polynomials byโ[ํฑ]. The dimension ofโ ํ [ํฑ], the space of ํ-variate polynomials of degree at most ํ, is given by ํ ํ (ํ) = ( ํ + ํ ํ ) . Let ํ ํ โถ 1 โค ํ โค ํ ํ (ํ) be a basis ofโ ํ [ํฑ]. We denote ํฃ ํ โถโ ํ โถโ ํ ํ (ํ) ํฑโผ ( ํ 1 (ํฑ),ํ 2 (ํฑ),โฏ,ํ ํ ํ (ํ) (ํฑ) ) ํ The monomials in ํฃ ํ (ํฑ) are graded in the lexicographic order 2 . Let ฮฉ โโ ํ be a compact set, with non-empty interior. Let ํ be a Borel measure supported on ฮฉ and define the associated moment matrix. Definition 2.1 (The moment matrix). The moment matrix of degree ํ โโ, associated with measure ํ, denoted by ํ ํ (ํ) โโ ํ ํ (ํ)รํ ํ (ํ) , is defined as ํ ํ (ํ) = โซ โ ํ ํฃ ํ (ํฑ) ํฃ ํ (ํฑ) ํ ํํ(ํฑ).(1) 2 lexicographic order: monomials are first ordered according to as- cending total degree |ํผ|, and then using lexicographic order on variables considering ํฑ 1 = ํ, ํฑ 2 = ํ, etc. Note that this matrix is symetric positive definite, thus non-singular for all ํ (see [6, Section 2.2] or [13, Remark 2.3] for the proof). The population Christoffel function is defined as follows. Definition 2.2 (The population Christoffel function). The population Christoffel function of degree ํ โโ, associated with the measure ํ, denoted by ฮ ํ ํ (ํฑ), is defined as ฮ ํ ํ (ํฑ) = ํํํ ํโโ ํ [ํฑ] โซ ฮฉ ํ 2 (ํณ) ํํ(ํณ), ํ(ํฑ) = 1 . (2) Now, for any polynomial ํ โโ ํ [ํฑ], there exists some ํฉ โโ ํ ํ (ํ) such that ํ(ํฑ) = ํฉ ํ ํฃ ํ (ํฑ) for any ํฑ โโ ํ . Thus, the objective function becomes โซ ฮฉ ํฉ ํ ํฃ ํ (ํณ) ํฃ ํ (ํณ) ํ ํฉ ํํ(ํณ) = ํฉ ํ ํ ํ (ํ) ํฉ, so that ฮ ํ ํ (ํฑ) = ํํํ ํฉโโ ํ ํ (ํ) ํฉ ํ ํ ํ (ํ) ํฉ,ํฉ ํ ํฃ ํ (ํฑ) = 1 . (3) The Christoffel-Darboux Kernel, which is defined below, is related to the Christoffel function. Definition 2.3 (The Christoffel-Darboux Kernel). The Christoffel-Darboux Kernel (CD-Kernel) associated with the measure ํ, denoted by ํพ ํ ํ (ํฑ, ํฒ), is defined as (ํฑ, ํฒ)โฆ ํพ ํ ํ (ํฑ, ํฒ) = ํฃ ํ (ํฑ) ํ ํ ํ (ํ) ํฃ ํ (ํฒ),(4) while the polynomial ํ ํ,ํ reads ํฑโฆ ํ ํ,ํ (ํฑ) = ํพ ํ ํ (ํฑ, ํฑ) = ํฃ ํ (ํฑ) ํ ํ ํ (ํ) โ1 ํฃ ํ (ํฑ). (5) ํ ํ,ํ is a sum-of-squares (SOS) polynomial of degree2ํ. An interesting property of this SOS polynomial is its behavior inside and outside its support ฮฉ. Lasserre et al. [8, Lemma 4.3.1] quantifies at least the exponential growth with ํ for data points outside the support, while inside, it is at most polynomial [8, Lemma 4.3.2]. The CF of degree ํ โโ defined in equation (2) can be rewritten as ฮ ํ ํ (ํฑ) = 1 ํ ํ,ํ (ํฑ) = 1 ํฃ ํ (ํฑ) ํ ํ ํ (ํ) โ1 ํฃ ํ (ํฑ) (6) and we have ฮ ํ ํ (ํฑ) โ1 = ํ ํ,ํ (ํฑ),(7) so that ฮ ํ ํ (ํฑ) โ1 inherits from the properties of ํ ํ,ํ (ํฑ). 2.2. The empirical Christoffel function In practical applications, the measure ํ is unknown. Let ๎ be a cloud of ํ data points ํฑ โโ ํ sampled from the theoretical measure ํ supported on ฮฉ. We define the discrete measure ํ ํ supported on ๎ such that ํ ํ = 1 ํ โ ํฑโ๎ ํฟ ํฑ where ํฟ ํฑ corresponds to the Dirac measure at ํฑ. The empirical version of the moment matrix can be written as ํ ํ (ํ ํ ) = 1 ํ โ ํฑโ๎ ํฃ ํ (ํฑ) ํฃ ํ (ํฑ) ํ .(8) Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 2 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection Lasserre et al. [8, Corollary 6.3.5] guarantees that the matrix ํ ํ (ํ ํ ) is invertible if the size of ๎: ํ, is greater than ํ ํ (ํ), e.g. the number of samples is greater than the size of the moment matrix. Definition 2.4 (The empirical Christoffel function). Under the condition |๎| = ํ > ํ ํ (ํ), the empirical CF is defined as ฮ ํ ํ ํ (ํฑ) = 1 ํฃ ํ (ํฑ) ํ ํ ํ (ํ ํ ) โ1 ํฃ ํ (ํฑ) .(9) According to Lasserre and Pauwels [7, Theorem 3.13], the empirical CF converges to the population CF as ํ increases: โฮ ํ ํ ํ โฮ ํ ํ โ โ = ํ ํขํ ํฑโโ ํ |ฮ ํ ํ ํ (ํฑ) โ ฮ ํ ํ (ํฑ)| โถ ํโโ 0 ํ.ํ . 3. The Christoffel function for outlier detection in data streams Ducharlet et al. [2] introduce DyCF, a novel outlier de- tection algorithm for data streams leveraging the Christoffel function. This algorithm supports online learning using a rank-ํ update. This section details the anomaly detection method based on the CF as well as the online learning principle of the method. 3.1. The Christoffel function for outlier detection As mentioned in section 2, the polynomial ํ ํ,ํ (ํฑ), and hence ฮ ํ ํ (ํฑ) โ1 , effectively captures the shape of the underlying dataset. Furthermore, as explained above, there exists a dichotomy in the growth behavior of ฮ ํ ํ (ํฑ) โ1 : it exhibits at most polynomial growth when ํฑ is within the support ฮฉ and at least exponential growth when ํฑ is outside ฮฉ. Since ฮ ํ ํ ํ (ํฑ) โ1 converges to ฮ ํ ํ (ํฑ) โ1 , these properties are preserved for finite datasets. Consequently, ฮ ํ ํ ํ (ํฑ) โ1 is well-suited as a scoring function for outlier detection. One can define a level set, or threshold, ํพ ํ,ํ such that all points ํฑ โโ ํ with a value of ฮ ํ ํ ํ (ํฑ) โ1 higher than ํพ ํ,ํ are considered as outliers. This define the scoring function ํ ํ,ํ (ํฑ) = ฮ ํ ํ ํ (ํฑ) โ1 ํพ ํ,ํ ,(10) where a point is detected as an outlier if ํ ํ,ํ โฅ 1. 3.2. Online learning in DyCF Considering ๎ as a dataset, ํ ํ (ํ ํ ) can be seen as a summary or an encoding of this dataset. In a data stream, if a new nominal instance ํฑ arrives, we can improve the performance of the algorithm by integrating this point into the database ๎ and updating the moment matrix ํ ํ (ํ ํ ). This is called online learning. The first update method would be to recalculate the moment matrix ํ ํ (ํ ํ ) using equation (8) and to invert it. However, this last operation is very costly. Nevertheless, we can notice that: ํ ํ (ํ ํ+1 ) = 1 ํ + 1 โ ํณโ๎โชํฑ ํฃ ํ (ํณ)ํฃ ํ (ํณ) ํ = 1 ํ + 1 ( ํ ํ ํ (ํ ํ ) + ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ) . Thus, we can use the Sherman-Morrisonโs formula or the Woodbury Matrix Identity to avoid recalculating the inverse of ํ ํ (ํ ํ+1 ), which is an ํ ํ (ํ) ร ํ ํ (ํ) matrix. When the update involves ํ โโ โ new data points ํฑ โโ ํ , the online learning phase goes as follows in three steps: (i) We denormalize ํ ํ (ํ ํ ) to obtain either ํ = ํ ํ (ํ ํ ) ร ํ, or ํ โ1 = ํ ํ (ํ ํ ) โ1 โํ. (i) The DI method calculates ํ โ1 ํขํํํํกํํ from ํ using ํ โ1 ํขํํํํกํํ = ( ํ + ํ โ ํ=1 ํฃ ํ ( ํฑ (ํ) ) ํฃ ํ ( ํฑ (ํ) ) ํ ) โ1 (11) and the ISM and WMI methods calculate ํ โ1 ํขํํํํกํํ from ํ โ1 using Sherman-Morrison (SM) and WMI formulas, e.g. equations (13) and (14) given below, respectively. (i) Renormalize ํ โ1 ํขํํํํกํํ to obtain the updated inverse moment matrix ํ ํ (ํ ํ+ํ ) โ1 = ํ โ1 ํขํํํํกํํ ร (ํ + ํ). The Sherman-Morrison formula and the Woodbury Matrix Identity โ Suppose ํด โโ ํรํ is an invertible square matrix and ํข,ํฃ โโ ํ are column vectors. Then ํด + ํขํฃ ํ is invertible if and only if 1 + ํฃ ํ ํด โ1 ํข โ 0. In this case, the SM formula [1] states the following (ํด + ํขํฃ ํ ) โ1 = ํด โ1 โ ํด โ1 ํขํฃ ํ ํด โ1 1 + ํฃ ํ ํด โ1 ํข ,(12) which, in our case, becomes: ( ํ + ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ) โ1 = ํ โ1 โ ํ โ1 ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ํ โ1 1+ํฃ ํ (ํฑ) ํ ํ โ1 ํฃ ํ (ํฑ) . (13) Now, let us recall the Woodbury Matrix Identity (WMI): ( ํด + ํํถํ ) โ1 = ํด โ1 โ ํด โ1 ํ ( ํถ โ1 + ํ ํด โ1 ํ ) โ1 ํ ํด โ1 , (14) which, in our case, becomes: ( ํ + ํ ํ ํผํ ) โ1 = ํ โ1 โ ํ โ1 ํ ํ ( ํผ โ1 + ํํ โ1 ํ ํ ) โ1 ํํ โ1 , (15) with ํผ the identity matrix of size ํ and ํ the design matrix: ํ = [ ํฃ ํ ( ํฑ (1) ) ํฃ ํ ( ํฑ (2) ) โฏ ํฃ ํ ( ํฑ (ํ) ) ] ํ โโ ํรํ . Note that the normalization costs of steps (i) and (i) are the same regardless of the method used. The most efficient method will therefore be the one with the lowest cost in step (i). Moreover, the DI method updates both ํ ํ (ํ ํ ) and ํ ํ (ํ ํ ) โ1 . However, to detect outliers, sinceฮ ํ ํ ํ only uses ํ ํ (ํ ํ ) โ1 , we do not need to compute ํ ํ (ํ ํ+ํ ), so we can use the ISM or WMI methods. Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 3 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection 4. Computational costs For the sake of simplifying calculations, the size of the moment matrix ํ ํ (ํ) will be referred to as ํ in this section. Appendix A reports the costs of the intermediate steps used in this section in terms of floating-point operations (flops). 4.1. Computational cost of the DI method for a rank-ํ update In this subsection, we calculate the computational cost of the DI method and provide an algorithm for it. 4.1.1. Computational cost of a matrix inversion The computational cost of inverting the moment matrix ํ ํ (ํ ํ ) of size ํ ร ํ is dependent on the algorithm em- ployed. Using an LU factorization, the computational cost of ํ ํ (ํ ํ ) โ1 is expressed as [9, Theorem 2.31]: ํ ( 2ํ 3 ) ํํํํํ .(16) Since our moment matrix ํ ํ (ํ ํ ) is symetric positive definite (spd), employing Cholesky decomposition reduces the cost to [5, Section I.A]: ํ ( 5 6 ํ 3 ) ํํํํํ .(17) 4.1.2. Rank-ํ update computational cost To perform a rank-ํ update and apply equation (11), we first need to compute ํ โ ํ=1 ํฃ ํ ( ํฑ (ํ) ) ํฃ ํ ( ํฑ (ํ) ) ํ . This in- volves performing ํ column-vector by row-vector products as described in equation (29), with a computational cost of ํรํ 2 ํํํํํ . Then, we need to sum the ํ resulting matrices of size ํ ร ํ , requiring ํ โ 1 term-by-term matrix additions with a total cost of (ํ โ 1)ํ 2 ํํํํํ . Next, we perform a term-by-term matrix addition with ํ, which incurs a cost of ํ 2 ํํํํํ . Thus, the computational cost of updating ํ to form ํ ํขํํํํกํํ = ํ + โ ํ ํ=1 ํฃ ํ ( ํฑ (ํ) ) ํฃ ํ ( ํฑ (ํ) ) ํ amounts to: ํํ 2 + (ํ โ 1)ํ 2 + ํ 2 = 2ํํ 2 ํํํํํ . Note that if we compute ํ ํขํํํํกํํ using the design ma- trix ํ = [ ํฃ ํ ( ํฑ (1) ) ํฃ ํ ( ํฑ (2) ) โฏ ํฃ ํ ( ํฑ (ํ) ) ] ํ โโ ํรํ and the relation ํ ํขํํํํกํํ = ํ + ํ ํ ํ, the computational cost remains equivalent: 1 matrix-by-matrix product accord- ing to equation (32) and 1 term-by-term addition yields a cost of: 2ํํ 2 โํ 2 +ํ 2 = 2ํํ 2 ํํํํํ . However, considering the superior optimization in Python, we will employ this computational approach for efficiency during our tests. Finally, we need to compute the inverse of ํ ํขํํํํกํํ , which is a spd matrix with the cost described in equation (17). Thus, the computational cost for the DI method (Algo- rithm 1) is ํ ( 5 6 ํ 3 ) + 2ํํ 2 ํํํํํ .(18) Algorithm 1 Direct Inversion Algorithm with ํ ( 5 6 ํ 3 ) + 2ํํ 2 ํํํํํ (18) Require: Matrix ํ with size ํ ร ํ , and vectors ํฃ ํ ( ํฅ (ํ) ) ,ํ = 1...ํ to add to ํ; 1: Construct the design matrix ํ = [ ํฃ ํ ( ํฑ (1) ) ํฃ ํ ( ํฑ (2) ) โฏ ํฃ ํ ( ํฑ (ํ) ) ] ํ โโ ํ รํ ; 2: ํ ํขํํํํกํํ โ ํ + ํ ํ ํ;โณ 2ํํ 2 ํํํํํ 3: Compute ํ โ1 ํขํํํํกํํ using Cholesky decomposition; โณ ํ ( 5 6 ํ 3 ) ํํํํํ 4: Output ํ โ1 ํขํํํํกํํ ; 4.2. Computational cost of the ISM method for a rank-N update In this subsection, we calculate the computational cost of the ISM method and provide an algorithm for it. 4.2.1. Computational cost of the numerator There are three ways to compute the numerator of the SMโs formula given by equation (13) and recalled below: ( ํ + ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ) โ1 = ํ โ1 โ ํ โ1 ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ํ โ1 1 + ํฃ ํ (ํฑ) ํ ํ โ1 ํฃ ํ (ํฑ) . - Compute the outer product first, i.e., ํ โ1 ( ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ ) ํ โ1 . This involves 1 column- vector by row-vector product (equation (29)) and 2 matrix products (equation (32)), totaling: ํ 2 +2(2ํ 3 โ ํ 2 ) = 4ํ 3 โ ํ 2 ํํํํํ . - Compute left to right. This involves 1 product matrix by column-vector (equation (31)), 1 column-vector by row-vector product (equation (29)), and 1 matrix product (equation (32)), totaling: 2ํ 2 โํ +ํ 2 +2ํ 3 โ ํ 2 = 2ํ 3 + 2ํ 2 โ ํ ํํํํํ . - Compute the matrices-vectors products first, i.e., ( ํ โ1 ํฃ ํ (ํฑ) ) ( ํฃ ํ (ํฑ) ํ ํ โ1 ) . This involves 2 matrix- vector products (equations (31) and (30)), and 1 column-vector by row-vector product (equation (29)), totaling: 2(2ํ 2 โํ )+ํ 2 = 5ํ 2 โ2ํ ํํํํํ . Moreover, since ํ is symmetric, ํ โ1 ํฃ ํ (ํฑ) = ( ํฃ ํ (ํฑ) ํ ํ โ1 ) ํ , so we only have 1 matrix-vector product to compute. Thus, the computational cost is 3ํ 2 โ ํ ํํํํํ . The most effective way to compute the numerator of equation (13) is to first perform the matrix-vector products, resulting in a computational cost of 5ํ 2 โ 2ํ ํํํํํ , and in our spd case, a cost of 3ํ 2 โ ํ ํํํํํ .(19) 4.2.2. Computational cost of the denominator Since we have already computed ํ โ1 ํฃ ํ (ํฑ) during the numerator computation, we only need to perform 1 row- vector by column-vector product (equation (28)), and 1 addition which amounts to 2ํ โ 1 + 1 flops, i.e., 2ํ ํํํํํ .(20) Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 4 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection 4.2.3. Sherman-Morrison computational cost The numerator of equation (13) gives us an ํ รํ matrix, and the denominator is real. Thus, there is a term-by-term division that costs ํ 2 ํํํํํ . However, if we perform the division after computing ํ โ1 ํฃ ํ (ํฑ) and before the column- vector by row-vector product of the numerator, we only divide a vector of size ํ , reducing the division cost from ํ 2 to ํ . Finally, we have a term-by-term matrix subtraction that costs ํ 2 ํํํํํ . Thus, the computational cost of the SMโs formula for a spd matrix is (19) + (20) + ํ + ํ 2 = 3ํ 2 โ ํ + 2ํ + ํ + ํ 2 flops, i.e., 4ํ 2 + 2ํ ํํํํํ .(21) 4.2.4. Iterative Sherman-Morrison computational cost The SM formula provides a method for computing the inverse of a matrix that has been modified by a rank-1 update. Specifically, if the matrix is updated with ํ outer products of the form ํฃ ํ (ํฑ)ํฃ ํ (ํฑ) ํ , the inverse can be obtained by applying the SM formula iteratively ํ times. As a result, the computational cost with a rank-ํ update using the ISM method for a spd matrix (Algorithm 2) is 4ํํ 2 + 2ํํ ํํํํํ .(22) Algorithm 2 Iterative Sherman-Morrison Algorithm with 4ํํ 2 + 2ํํ ํํํํํ (22) Require: The inverse of the matrix ํ: ํ โ1 with size ํ รํ , and vectors ( ํฃ ํ ( ํฅ (ํ) )) ํ=1...ํ โโ ํ to add to ํ; 1: Initialize ํ โ1 ํขํํํํกํํ โ ํ โ1 ; 2: for ํ = 1 to ํ do 3: ํโ ํ โ1 ํขํํํํกํํ ร ํฃ ํ ( ํฅ (ํ) ) ;โณ 2ํ 2 โ ํ ํํํํํ 4: ํโ 1 + ํฃ ํ ( ํฅ (ํ) ) ํ ร ํ;โณ 2ํ ํํํํํ 5: ํ ํ โ ํโํ;โณ ํ ํํํํํ 6: ํ โ1 ํขํํํํกํํ โ ํ โ1 ํขํํํํกํํ โ ํ ํ ร ํ ํ ;โณ 2ํ 2 ํํํํํ 7: end for 8: Output ํ โ1 ํขํํํํกํํ ; 4.3. Computational cost of the WMI method for a rank-N update In this subsection, we calculate the computational cost of the WMI method and provide an algorithm for it. The WMI is given by equation (15) is recalled below: ( ํ + ํ ํ ํผํ ) โ1 = ํ โ1 โ ํ โ1 ํ ํ ( ํผ โ1 + ํํ โ1 ํ ํ ) โ1 ํํ โ1 . The first step is to compute the matrix productsํํ โ1 ํ ํ , with a computational cost of equation (32): ํถ ํรํ ํ รํ + ํถ ํรํ ํ รํ = 2ํํ 2 โ ํํ โโโโโ ํ =ํรํ โ1 +2ํ 2 ํ โ ํ 2 โโโโโ ํ รํ ํ = 2ํํ 2 + (2ํ 2 โ ํ)ํ โ ํ 2 ํํํํํ . Next, we perform a term-by-term addition with ํผ, which has a computational cost of ํ 2 ํํํํํ . Then, we need to invert the resulting summed matrix, which by construction is a spd ํ ร ํ matrix. The computa- tional cost of this inversion step is ํ ( 5 6 ํ 3 ) ํํํํํ (17). Since we have already computed ํ = ํ ร ํ โ1 , and given that ํ is spd, then ํ โ1 รํ ํ = (ํ รํ โ1 ) ํ = ํ ํ . Thus, we only have 2 matrix products to compute, that is 2ํ 2 ํ โ ํํ โโโโโ ํ=ํ ํ รํํํฃ +2ํํ 2 โ ํ 2 โโโโโ ํรํ = (2ํ โ 1)ํ 2 + (2ํ 2 โ ํ)ํ ํํํํํ . Finally, a term-by-term subtraction is performed, which costs ํ 2 ํํํํํ . Combining these steps, the overall computational cost for applying the WMI method (Algorithm 3) is: 2ํํ 2 + (2ํ 2 โ ํ)ํ โ ํ 2 + ํ 2 + ํ ( 5 6 ํ 3 ) + (2ํ โ 1)ํ 2 + (2ํ 2 โ ํ)ํ + ํ 2 = 4ํํ 2 + (4ํ 2 โ 2ํ)ํ + ํ ( 5 6 ํ 3 ) ํํํํํ (23) Algorithm 3 Woodbury Matrix Identity Algorithm with 4ํํ 2 + (4ํ 2 โ 2ํ)ํ + ํ ( 5 6 ํ 3 ) ํํํํํ (23) Require: The inverse of the matrix ํ: ํ โ1 with size ํ รํ , and vectors ( ํฃ ํ ( ํฅ (ํ) )) ํ=1...ํ โโ ํ to add to ํ; 1: Initialize ํผ โโ ํรํ and construct the design matrix ํ = [ ํฃ ํ ( ํฑ (1) ) ํฃ ํ ( ํฑ (2) ) โฏ ํฃ ํ ( ํฑ (ํ) ) ] ํ โโ ํ รํ ; 2: ํ โ ํํ โ1 ;โณ 2ํํ 2 โ ํํ ํํํํํ 3: ํโ ํผ + ํ ํ ํ ;โณ 2ํ 2 ํ ํํํํํ 4: Compute ํ โ1 using Cholesky decomposition; โณ ํ ( 5 6 ํ 3 ) ํํํํํ 5: ํโ ํ ํ ํ โ1 ;โณ 2ํ 2 ํ โ ํํ ํํํํํ 6: ํ โ1 ํขํํํํกํํ โ ํ โ1 โ ํํ ;โณ 2ํํ 2 ํํํํํ 7: Output ํ โ1 ํขํํํํกํํ ; 5. Comparison of the computational costs of the DI, ISM, and WMI methods We recall that the size of the moment matrix ํ ํ (ํ ํ ) is ํ รํ , and that we want to addํ new data points intoํ ํ (ํ ํ ). The computational cost of each method is summarized in Table 1. The goal of this section is to find the number of new data ํ to add such that the DI cost given by equation (18) is lower than that of the ISM method given by equation (22) or that of the WMI method given by equation (23). 5.1. Direct inversion method vs. Iterative Sherman-Morrison method The condition on ํ under which the DI method is prefer- able to the ISM method is obtained by comparing their computational cost given by equations (18) and (22), re- spectively. The condition is obtained by isolating ํ in the following inequality: ํ ( 5 6 ํ 3 ) + 2ํํ 2 < 4ํํ 2 + 2ํํ โบ 5 6 ํ 2 < (2ํ + 2)ํ โบ ํ > 5ํ 2 12 ( ํ + 1 ) (24) Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 5 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection Table 1 Computational cost of the three compared methods (in flops) DI (Algorithm 1)ISM (Algorithm 2)WMI (Algorithm 3) Computational cost ํ ( 5 6 ํ 3 ) + 2ํํ 2 (18)4ํํ 2 + 2ํํ (22)4ํํ 2 + (4ํ 2 โ 2ํ)ํ + ํ ( 5 6 ํ 3 ) (23) (a) DI (subsection 4.1) vs. ISM (subsection 4.2)(b) DI (subsection 4.1) vs. WMI (subsection 4.3) Figure 1: Theoretical representation of thresholds from which choosing DI over ISM or WMI 5.2. Direct inversion method vs. Woodbury matrix identity method We next compare the DI method with the WMI method: ํ ( 5 6 ํ 3 ) + 2ํํ 2 < 4ํํ 2 + (4ํ 2 โ 2ํ)ํ + ํ ( 5 6 ํ 3 ) โบ 5 6 ํ 3 < 2ํํ 2 + (4ํ 2 โ 2ํ)ํ + 5 6 ํ 3 โบ 5 6 ํ 3 + 4ํ ํ 2 + 2(ํ 2 โ ํ )ํ โ 5 6 ํ 3 > 0(25) Here we have a cubic equation that we need to solve with a fixed s. An empirical approximation is given in subsection 6.1. 5.3. Woodbury matrix identity method vs. Iterative Sherman-Morrison method For ํ โฅ 1 and ํ โฅ 1, the computation cost of the ISM method (equation (22)) is generally lower than that of the WMI method (equation (23)). However, these figures do not account for memory efficiency (read/write operations), nor for the optimization of matrix computations compared to vector calculations and successive iterations, which intro- duce significant overheads. The measurements carried out in subsection 6.2 show that the WMI method quickly becomes more efficient than the ISM method. 6. Optimal method selection based on the number ํ of new data 6.1. Theoretical guide for method selection based on new data volume ํ The results shown in Figure 1 are obtained by computing equations (24) and (25) for all ํ โJ1,4200K and ํ โ J1,10000K. Figure 1a shows that one must choose the DI method over the ISM method when ํ > ํ โ2.4002 โ 5ํ โ12. Figure 1b shows that one must choose the DI method over the WMI method when ํ > ํ โ3.7506. We name this relation the empirical DI over WMI threshold: ํ > ํ 3.7506 (26) Table 2 summarizes the different theoretical thresholds identified for choosing the DI method over the ISM or WMI methods. 6.2. Experimental results The following results were obtained on a laptop equipped with an Intel Core Ultra 5 135U CPU and 32 GB of RAM, running Windows 11. All results can be reproduced using the Python code and data available on GitHub 3 . We generate a random dataset of ํ = 2000 samples in โ 1287 . This simulates 2000 vectors ํฃ ํ (ํฑ) with ํฑ โโ ํ=8 and polynomials of degree at most ํ = 5: ํ = ํ ํ (ํ) = ( ํ + ํ ํ ) = ( 8 + 5 5 ) = 1287. For every ํ in 1,2,3,4,5,10,20,30,40,50,100,200, 300,400,500,750,1000, we compute the moment matrix ํ ํ (ํ ํ ) and its inverse ํ ํ (ํ ํ ) โ1 with the first ํ = ํ โํ samples. Then we update the inverse of the moment matrix with the last ํ samples using the three different methods and report the average execution times on the ํํ = 200 simulations in Table 3. Note that for the ISM method, we reduced ํํ to 50 when ํ > 50. We also computed, for each method, the error of the inverse moment matrix:ํ ํ = โํผโํ ํขํํํํกํํ รํ โ1 ํขํํํํกํํ,ํ โ ํน , where ํ โ1 ํขํํํํกํํ,ํ corresponds to the updated inverse mo- ment matrix ํ โ1 ํขํํํํกํํ found by the method ํ โ DI, ISM, 3 Link to the GitHub repository: https://github.com/fgrivet/cost-trade- offs-matrix-inversion-update Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 6 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection Table 2 Theoretical thresholds for choosing the DI method over the ISM or WMI methods DI over ISMDI over WMI Empirical DI over WMI ํ5ํ 2 โ(12(ํ + 1)) (24) 5 6 ํ 3 + 4ํ ํ 2 + 2(ํ 2 โ ํ )ํ โ 5 6 ํ 3 > 0 (25, to solve with ํ fixed)ํ โ3.7506 (26) Table 3 Average execution time (in seconds) for a rank-ํ update of the moment matrix ํ ํ (ํ ํ ) of size 1287 ร 1287 Smallest time โ Second best time ํDIISMWMI 10.085946380.0077431450.0077865 20.084216940.011417140.007551792 30.080868380.014240550.005989846 40.079419040.021481580.005903258 50.074378320.018835380.004278867 100.083148580.062156030.008458068 200.080200280.099598040.008605508 300.089235320.19004650.01317246 400.080280410.19666030.01141213 500.084503160.29697780.01532174 1000.082192190.48513110.02341882 2000.085731090.96126090.04289328 3000.093682111.4530230.06708244 4000.090565861.7950420.09021928 5000.097773042.2074090.1229474 7500.10828573.8672140.2134477 10000.1189393.8297420.3350951 Table 4 Errors for a rank-ํ update of the moment matrix ํ ํ (ํ ํ ) of size 1287 ร 1287 with conditioning ํํํํ(ํ ํ (ํ ํ )), ํ = ํ โ ํ Smallest error โ Second best error ํDIISMWMIํํํํ(ํ ํ (ํ ํ )) 11.580182e-131.582999e-131.583327e-1383.9304 21.593879e-132.9023511.605027e-1384.0449 31.594071e-134.3589411.595547e-1384.2276 41.593055e-135.3335181.609318e-1385.1681 51.59972e-136.1500921.604706e-1385.412 101.598138e-139.5442791.610103e-1386.4186 201.606954e-1314.105321.628296e-1388.7121 301.607787e-1317.766591.641526e-1391.1043 401.607618e-1320.821281.662739e-1393.5837 501.61661e-1323.757831.677276e-1396.0594 1001.618503e-1338.348131.773381e-13104.694 2001.625134e-1371.47272.144143e-13136.435 3001.656789e-13120.57272.788295e-13219.31 4001.654148e-13210.00064.830752e-13325.742 5001.655513e-13417.02481.316689e-12604.421 7501.681712e-135.943877e+1021941.933.44460e+18 10001.694056e-132.559507e+11503624.64.4279e+18 WMI and โ.โ ํน is the Frobenius norm [3] given by: โํดโ ํน = โ โ ํ,ํ |ํ ํ,ํ | 2 (27) The errors for each method and each value of ํ are reported in Table 4. We observe that the errors associated with the WMI and ISM methods increase rapidly, when ํ > 500, indicating numerical instability. This behavior is largely explained by the poor conditioning ( โ 10 18 ) of the empirical moment matrix constructed from limited data (1250 and 1000 samples, respectively). According to Vershynin [12], and Wainwright [14], in- creasing the number of samples improves conditioning and numerical stability. So we repeated the experiment with a larger sample size: ํ = 15000. The errors for all methods stabilize around 10 โ13 across all tested values of ํ (up to ํ = 1500) and matrix sizes up to ํ = 2000, confirming that the large errors observed in the small-sample regime are primarily due to ill conditioning of the moment matrix. For ISM method, the error still increases with ํ for both ํ = 2000 and ํ = 15000, reflecting the accumulation of floating-point rounding errors inherent to successive rank- one updates [4]. However, with sufficiently many samples, the error growth remains bounded and no longer exhibits the rapid instability observed for small sample sizes. The numerical stability of the Sherman-Morrison and Woodbury updates has been analyzed previously in [16]. 6.3. Comparison with theoretical results The goal of this subsection is to validate the theoretical thresholds identified in subsection 6.1, Table 2 and the values measured in subsection 6.2, Table 3. In Table 5, we present the theoretical and experimental thresholds for selecting the DI method over the ISM or WMI methods, with the matrix size ํ fixed at 1287. First, we observe that the empirical DI over WMI approximation, given by equation (26) as 343.145, is very close to the theoretical threshold of equation (25), which is 343.250. Second, as predicted in subsection 5.3, the WMI method proves to be more efficient than the ISM method even for small values of ํ ; specifically, starting from ํ = 2. When comparing the thresholds for choosing the DI method over the ISM or WMI methods, we observe that the threshold (24) for the ISM method is significantly over- estimated. The theoretical prediction is 535.834, whereas experimentally, it falls between 10 and 20. Conversely, the threshold (25) for the WMI method is substantially under- estimated. The theoretical prediction is 343.250, but ex- perimental results place it between 400 and 500, which is approximately ํ 3 . These discrepancies arise because our analysis considered only the number of operations (flops) and did not account for the efficiency of memory access Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 7 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection Table 5 Thresholds for choosing the DI method over the ISM or WMI methods with ํ = 1287: ํ = 8 and ํ = 5 DI over ISMDI over WMI Empirical DI over WMI ํ ํกโํํํํํกํํํํ (24) = 5ร1287 2 12(1287+1) โ 535.834 (25) โ 343.250(26) = 1287 3.7506 โ 343.145 ํ ํํฅํํํํํํํํกํํ Between 10 and 20Between 400 and 500 Figure 2: Fastest method for each (s, k) pair based on execution time (read/write operations) or the intrinsic Python optimizations for matrix calculations compared to vector calculations or iterative loops. 6.4. Experimental guide for method selection based on new data volume ํ We repeated the experiment conducted in subsection 6.2. This time, we used ํ = 2000 samples inโ ํ with ํ ranging in 10,20,50,100,250,500,750,1000 and ํํ = 200 sim- ulations. However, for the ISM method, when ํ > 500 and ํ > 50, we reduced ํํ to 50. Figure 2 illustrates the fastest method for performing the rank-ํ update for each pair (ํ ,ํ). A green point signifies that ISM (Algorithm 2) is the fastest method, a blue point indicates that WMI (Algorithm 3) is the fastest method, and a red point denotes that the DI (Algorithm 1) is the fastest method. Firstly, we notice that for a small matrix size (ํ = 10), the DI method outperforms the other two. Secondly, as expected, ISM is faster for a rank-1 update. Then the WMI method is faster for a rank-ํ update when ํ < ํ 3 , and finally, for a large update, the DI method is superior. Our recommendation, for a Python CPU implementa- tion, is hence as follows: - Use the ISM method (Algorithm 2) for a rank-1 up- date. - Use the WMI method (Algorithm 3) for a rank-ํ update when ํ โค ํ 3 , where ํ is the matrix size. - Use the DI method (Algorithm 1) if ํ > ํ 3 . 7. Conclusion This paper compares three methods for updating the inverse of a matrix after a rank-ํ update, assuming the orig- inal inverse is known. Although motivated by Christoffel- function-based outlier detection, the analysis applies more broadly to problems involving rank-ํ updates of spd invert- ible matrices. Theoretical thresholds for method selection are derived and summarized in Table 2. They are validated empirically with comprehensive experiments that account for implementation effects. This note includes practical guidance for Python CPU implementations. With ํ denoting the matrix dimension and ํ the update rank, it recommends the following rule: - For ํ = 1, use ISM (Algorithm 2). - For ํ โค ํ 3 , use WMI (Algorithm 3). - For ํ > ํ 3 , use DI (Algorithm 1). Although the theoretical computational costs derived in this work are expressed in terms of the matrix dimension ํ and the update rank ํ, and are therefore independent of any specific programming language or hardware platform, the practical recommendations and proposed selection rule are empirically validated only for Python implementations executed on CPU. In particular, the quantitative thresholds identified in this study depend on memory-access patterns characteristic of standard Python numerical linear algebra Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 8 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection libraries, and should not be assumed to transfer directly to other languages, libraries, or hardware architectures. Several directions for future work emerge from this study. First, it would be of interest to extend the empirical comparison to other programming languages and computa- tional platforms, including compiled implementations such as C++ and GPU-accelerated environments, where differ- ences in parallelism and memory hierarchies may signifi- cantly affect performance. Second, beyond inverse update strategies, a fundamental limitation of Christoffel-function- based anomaly detection lies in the size ํ of the moment matrix, which grows rapidly with the ambient dimension. Future work could therefore focus on reducing the effective dimension of this matrix in order to lower computational cost and enable CF-based scoring to scale to high-dimensional data streams. CRediT authorship contribution statement Florian Grivet: Writing - original draft, Writing - review and editing, Conceptualization, Formal analysis, Methodology, Software, Data curation, Investigation, Val- idation. Louise Travรฉ-Massuyรจs: Writing - review and editing, Supervision, Funding acquisition, Resources, Con- ceptualization, Validation. Data availability All code and datasets are available on GitHub: github.com/fgrivet/cost-trade-offs-matrix-inversion-update. All the tests conducted in this paper can be reproduced. The README file explains how to use the code. Declaration of generative AI and AI-assisted technologies in the manuscript preparation process During the preparation of this work the authors used Mistral in order to check the text grammar. After using this tool/service, the authors reviewed and edited the content as needed and hence take full responsibility for the content of the submitted article. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential compet- ing interests: Louise Travรฉ-Massuyรจs reports financial sup- port was provided by AI Interdisciplinary Institute ANITI. Acknowledgments This work has benefited from the AI Interdisciplinary Institute ANITI funded by the France 2030 program under the Grant agreements nยฐANR-19-P3IA-0004 and nยฐANR- 23-IACL-0002. The authors would like to thank Jean-Bernard Lasserre (LAAS-CNRS) and Didier Henrion (LAAS-CNRS) for their insightful discussions and contributions to the Christoffel function method used in this paper. A. Computational costs of generic elementary products In this appendix, we detail the computational costs of the intermediate steps employed in section 4, expressed in terms of flops. Notations โ We adopt the following generic notation: scalars are denoted by lowercase letters (e.g., ํ), vectors by boldface lowercase letters (e.g., ํ), and matrices by uppercase letters (e.g., ํด). The dimensions of vectors and matrices are represented by ํ, ํ and ํ. The computational costs are denoted by ํถ ํ ํ , where the superscript ํ indicates the dimension of the left element in the multiplication, and the subscript ํ indicates the dimension of the right element. A.1. Computational cost of multiplying a row-vector by a column-vector ํ 1 โฏํ ํโ1 ํ ํ [] ํ 1 โฎ ํ ํโ1 ํ ํ โก โข โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฅ โฆ ํค ร = Figure 3: Representation of the product of a row-vector by a column-vector Let ํ โโ 1รํ , ํ โโ ํร1 and ํค โโ such that ํค = ํร ํ (Figure 3). ํค = ํ โ ํ=1 ํ ํ รํ ํ , i.e. ํ multiplications and ํโ1 additions, which provides the following cost: ๎ฏ 1รํ ํร1 = 2ํ โ 1 ํํํํํ .(28) A.2. Computational cost of multiplying a column-vector by a row-vector ํ 1 ํ 2 โฎ ํ ํ โก โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฆ ํ 1 โฏํ ํโ1 ํ ํ [] ํ 1,1 ํ ํ,ํโ1 โก โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฆ ร = Figure 4: Representation of the product of a column-vector by a row-vector Let ํ โโ ํร1 , ํ โโ 1รํ and ํ โโ ํรํ such that ํ = ํ ร ํ (Figure 4). For every element of ํ , we have ํ ํ,ํ = ํ ํ ร ํ ํ , i.e. 1 multiplication for ํ ร ํ elements, which provides the following cost: ๎ฏ ํร1 1รํ = ํ ร ํ ํํํํํ .(29) When considering two vectors of equal length ํ, the compu- tational cost of their product is ๎ฏ ํร1 1รํ = ํ 2 ํํํํํ . Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 9 of 10 Cost Trade-offs in Matrix Inversion Updates for Streaming Outlier Detection A.3. Computational cost of multiplying a row-vector by a matrix ํ 1 ํ 2 โฏํ ํ [] ํต 1,1 ํต 1,ํ โฎ โฎ ํต ํ,1 ํต ํ,ํ โก โข โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฅ โฆ ํค 1 ํค 2 โฏํค ํ [] ร = Figure 5: Representation of the product of a row-vector by a matrix Let ํ โโ 1รํ , ํต โโ ํรํ and ํฐ โโ 1รํ such that ํฐ = ํ ร ํต (Figure 5). For every element of ํฐ, we have ํค ํ = ํ โ ํ=1 ํ ํ ร ํต ํ,ํ , i.e. ํ multiplications and ํ โ 1 additions for ํ elements, which equals to ํ ร (2ํ โ 1), which provides the following cost: ๎ฏ 1รํ ํรํ = 2ํํ โ ํ ํํํํํ .(30) When considering a square matrix of size ํ ร ํ, the computational cost of this product is ๎ฏ 1รํ ํรํ = 2ํ 2 โ ํ ํํํํํ . A.4. Computational cost of multiplying a matrix by a column-vector ํด 1,1 โฏํด 1,ํ ํด ํ,1 โฏํด ํ,ํ โก โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฆ ํ 1 ํ 2 โฎ ํ ํ โก โข โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฅ โฆ ํค 1 ํค 2 โฎ ํค ํ โก โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฆ ร = Figure 6: Representation of the product of a matrix by a column-vector Let ํด โโ ํรํ , ํ โโ ํร1 and ํฐ โโ ํร1 such that ํฐ = ํด ร ํ (Figure 6). For every element of ํฐ, we have ํค ํ = ํ โ ํ=1 ํด ํ,ํ ร ํ ํ , i.e. ํ multiplications and ํ โ 1 additions for ํ elements, which equals to ํ ร (2ํ โ 1), which provides the following cost: ๎ฏ ํรํ ํร1 = 2ํํ โ ํ ํํํํํ .(31) When considering a square matrix of size ํ ร ํ, the computational cost of this product is ๎ฏ ํรํ ํร1 = 2ํ 2 โ ํ ํํํํํ A.5. Computational cost of multiplying a matrix by a matrix ํด 1,1 โฏํด 1,ํ ํด ํ,1 โฏํด ํ,ํ โก โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฆ ํต 1,1 ํต 1,2 โฏํต 1,ํ โฎ โฎ ํต ํ,1 ํต ํ,2 โฏํต ํ,ํ โก โข โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฅ โฆ ํ 1,1 โฏํ 1,ํ ํ ํ,1 ํ ํ,2 โฏํ ํ,ํ โก โข โข โข โข โฃ โค โฅ โฅ โฅ โฅ โฆ ร = Figure 7: Representation of the product of a matrix by a matrix Let ํด โโ ํรํ , ํต โโ ํรํ and ํ โโ ํรํ such that ํ = ํด ร ํต (Figure 7). For every element of ํ , we have ํ ํ,ํ = ํ โ ํ=1 ํด ํ,ํ ร ํต ํ,ํ , i.e. ํ multiplications and ํโ1 additions for ํรํ elements, which equals to ํํร(2ํโ1), which provides the following cost: ๎ฏ ํรํ ํรํ = 2ํํํ โ ํํ ํํํํํ .(32) When considering squares matrices of size ํ ร ํ, the computational cost of this product is ๎ฏ ํรํ ํรํ = 2ํ 3 โํ 2 ํํํํํ . References [1] Bartlett, M.S., 1951. An Inverse Matrix Adjustment Arising in Discriminant Analysis. The Annals of Mathematical Statistics 22, 107 โ 111. doi:10.1214/aoms/1177729698. [2] Ducharlet, K., Travรฉ-Massuyรจs, L., Lasserre, J.B., Le Lann, M.V., Miloudi, Y., 2025. Leveraging the christoffel function for outlier detection in data streams. International Journal of Data Science and Analytics 20, 2021โ2037. doi:10.1007/s41060-024-00581-2. [3] Golub, G.H., Van Loan, C.F., 2013. Matrix Computations - 4th Edition. Johns Hopkins University Press, Philadelphia, PA. doi:10. 1137/1.9781421407944. [4] Higham, N.J., 2002. Accuracy and stability of numerical algorithms. SIAM. doi:10.1137/1.9780898718027. [5] Krishnamoorthy, A., Menon, D., 2013. Matrix inversion using cholesky decomposition, in: 2013 Signal Processing: Algorithms, Architectures, Arrangements, and Applications (SPA), p. 70โ72. URL: https://ieeexplore.ieee.org/document/6710599. [6] Lasserre, J.B., Pauwels, E., 2016. Sorting out typicality with the inverse moment matrix sos polynomial, in: Proceedings of the 30th International Conference on Neural Information Processing Systems, Curran Associates Inc., Red Hook, NY, USA. p. 190โ198. URL: https://dl.acm.org/doi/10.5555/3157096.3157118. [7] Lasserre, J.B., Pauwels, E., 2019. The empirical christoffel function with applications in data analysis. Advances in Computational Mathematics 45, 1439โ1468. doi:10.1007/s10444-019-09673-1. [8] Lasserre, J.B., Pauwels, E., Putinar, M., 2022. The Christof- felโDarboux Kernel for Data Analysis. 1st ed., Cambridge University Press. doi:10.1017/9781108937078. [9] Lu, J., 2024. Numerical matrix decomposition. arXiv:2107.02579. [10] Lu, T., Wang, L., Zhao, X., 2023. Review of anomaly detection algorithms for data streams. Applied Sciences 13. doi:10.3390/ app13106353. [11] Sherman, J., Morrison, W.J., 1950. Adjustment of an Inverse Matrix Corresponding to a Change in One Element of a Given Matrix. The Annals of Mathematical Statistics 21, 124 โ 127. doi:10.1214/aoms/ 1177729893. [12] Vershynin, R., 2018. High-dimensional probability: An introduction with applications in data science. volume 47. Cambridge university press. [13] Vu, Mai Trang, Bachoc, Franรงois, Pauwels, Edouard, 2022. Rate of convergence for geometric inference based on the empirical christoffel function. ESAIM: PS 26, 171โ207. doi:10.1051/ps/2022003. [14] Wainwright, M.J., 2019. High-dimensional statistics: A non- asymptotic viewpoint. volume 48. Cambridge university press. doi:10.1017/9781108627771. [15] Woodbury, M., of Statistics, P.U.D., 1950. Inverting Modified Ma- trices. Memorandum Report / Statistical Research Group, Princeton, Department of Statistics, Princeton University. [16] Yip, E.L., 1986. A note on the stability of solving a rank-p modifica- tion of a linear system by the shermanโmorrisonโwoodbury formula. SIAM Journal on Scientific and Statistical Computing 7, 507โ513. doi:10.1137/0907034. [17] Zhou, P., 2025. A survey of streaming data anomaly detection in network security. PeerJ Comput. Sci. 11, e3066. doi:10.7717/ peerj-cs.3066. Grivet and Travรฉ-Massuyรจs: Published in Array. DOI: 10.1016/j.array.2026.100737Page 10 of 10