{'fraction_non_alphanumeric': 0.05054573570927661, 'fraction_numerical': 0.022099149479438262, 'mean_word_length': 4.236723163841808, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistence B-Spline Grids: Stable Vector Representation of Persistence Diagrams Based on Data Fitting

Many attempts have been made in recent decades to integrate machine learning (ML) and topological data analysis. A prominent problem in applying persistent homology to ML tasks is finding a vector representation of a persistence diagram (PD), which is a summary diagram for representing topological features. From the perspective of data fitting, a stable vector representation, namely, persistence B-spline grid (PBSG), is proposed based on the efficient technique of progressive-iterative approximation for least-squares B-spline function fitting. We theoretically prove that the PBSG method is stable with respect to the metric of 1-Wasserstein distance defined on the PD space. The proposed method was tested on a synthetic data set, data sets of randomly generated PDs, data of a dynamical system, and 3D CAD models, showing its effectiveness and efficiency


Introduction
Inferring the topological structures of data at different scales is one major task in the advent of the big-data era. Persistent homology (Edelsbrunner et al., 2002;Zomorodian and Carlsson, 2005) provides a robust tool to capture topological features in data. These homological features in a dimension are summarized in a persistence diagram (PD) or a barcode (Ghrist, 2008). Persistent homology has also been applied successfully to deal with various practical problems in different scientific and engineering fields, such as biochemistry (Xia and Wei, 2015;Xia et al., 2018) and information science (De Silva and Ghrist, 2007;Rieck et al., 2018).
The machine learning (ML) community has shown a growing interest in combining traditional ML methods and the current techniques of computational topology (Chazal et al., 2015;Chen and Quadrianto, 2016;Khrulkov and Oseledets, 2018). In ML applications, classification tasks require comparing the similarity between two PDs. Therefore, the distance metrics of PDs, i.e., the Wasserstein distance and the bottleneck distance (Cohen-Steiner et al., 2007), play an important role. Although some efficient algorithms (Kerber et al., 2017) were proposed to compute the distance metrics, it is limited to applying ML methods and statistics on the analysis of PDs. One of the alternatives is to extend the representation of topological features into the form of vectors that have the property of stability with respect to the distance of PDs. A vector representation has the advantages of easy distance computation and being available for more ML tools than distance matrix form. For instance, a vectorizing representation of PD is available to support vector machine (SVM), a commonly used ML tool that requires vectors as input. This context, therefore, provides a clear motivation to represent PDs by stable vectors.

Related Work
Many approaches (Pachauri et al., 2011;Carriere et al., 2015;Kališnik, 2018;Bubenik, 2015;Reininghaus et al., 2015;Adams et al., 2017;Chevyrev et al., 2020) have been proposed to transform PDs into vectors. Some attempts have been developed when studying practical applications. For example, while studying the cortical thickness measurements of the human cortex for the study of Alzheimer's disease, Pachauri et al. (2011) rasterized PDs on a regular grid to vectorize the concentration map for SVM training. Some studies (Di Fabio and Ferri, 2015;Adcock et al., 2016;Kališnik, 2018) represented PDs by using an algebraic polynomial or function. For example, Kališnik (2018) identified tropical coordinates on the barcode space and proved the stability of this representation with respect to the barcode metrics. Some related studies (Mileyko et al., 2011) developed statistical measurements of the PDs, such as persistence landscapes (Bubenik, 2015), the functional representation satisfying the proposed ∞-landscape stability and p-landscape stability with the same condition as the stability of the p-Wasserstein distance (Cohen-Steiner et al., 2010). In this paper, we will compare the performance among our methods, persistence landscapes (PLs) (Bubenik, 2015), persistence images (PIs) (Adams et al., 2017), persistence scale space kernel (PSSK) (Reininghaus et al., 2015), persistence weighted Gaussian kernel (PWGK) (Kusano et al., 2018), sliced Wasserstein kernel (SWK) (Carriere et al., 2017), and persistence bag-of-words (PBoW) (Zieliński et al., 2019) in some ML classification tasks. Therefore, the related works on these competitive methods are reviewed.
Persistence images: PIs proposed by Adams et al. (2017) combined the idea of the kernel framework in Reininghaus et al. (2015) and the idea of "pixel" counting in Rouse et al. (2015). Specifically, the product of a Gaussian distribution and a weighted function corresponding to each PD point in the birth-persistence coordinates are added up to form a persistence surface. By integrating the persistence surface on each pixel divided on the domain, a grayscale image is produced, which is often reshaped to be a vector. In particular, the Gaussian function tends to the Dirac delta if the variance σ approaches zero, and PIs degenerate into the weighted summation in Rouse et al. (2015). Among their advantages, PIs allow users to assign a weight on each point in a PD, and provide an efficient and easily understandable approach to vectorize PDs for ML tasks. Similar to the case of the kernel in Reininghaus et al. (2015), PIs have 1-Wasserstein stability, and the persistence surface of a PI is additive, introducing difficulties in proving the stability results with respect to p-Wasserstein distance when 1 < p ≤ ∞.
Kernel methods: Kernel methods have been widely developed for finding an appropriate measurement of PD distance. The main idea of the framework of kernel methods is to transform PDs into probability density functions so that the distance between the PDs can be evaluated. Inspired by a heat diffusion problem with a Dirichlet boundary condition, Reininghaus et al. (2015) proposed the PSSK. PSSK is a framework of an additive multi-scale kernel, and has 1-Wasserstein stability. Kusano et al. (2018) summarized previous studies based on kernel functions and developed the PWGK that vectorizes PDs by using kernel embedding of measures into reproducing kernel Hilbert spaces. The stability of PWGK was proved with respect to the Hausdorff distance between data. The SWK for PDs was proposed by Carriere et al. (2017). The points on PDs are projected on a series of lines passing through the origin, and the slopes of these lines are adjusted. By SWK, the distance of PDs is measured as the distances between the projected points on these lines. SWK is a negative semi-definite kernel, and it has 1-Wasserstein stability. Furthermore, other kernel methods (Le and Yamada, 2018;Padellini and Brutti, 2017;Carrière and Bauer, 2019) are being developed to be incorporated into kernel-based ML models.
Persistence bag-of-words: Considering the considerable time consumption of kernel methods in practice caused by explicit computation of kernel matrix, Zieliński et al. (2019) proposed PBoW to transform a PD into a finite-dimension vector. This method extends the conventional bag-of-words to deal with PDs. The codebook, a finite-size representation, is generated by conducting clustering and obtaining the centers of clusters on the diagram into which all PDs of training samples are consolidated. The vector representation is produced by counting the number of points closest to the codeword in the codebook in a given order. In a mild condition, the authors gave a stable form of PBoW that has 1-Wasserstein stability. However, the obtained vectorizing representation via PBoW depends on the training samples and the clustering methods.

Contributions of the Present Study
From the perspective of data fitting and approximation, the weighted kernel framework can be thought to interpolate with basis functions the "importance" values. To transform a PD into an explicit vector derived from a B-spline control grid, we proposed a data fitting framework based on the persistence B-spline function that approximates the "importance" values, named eminence functions. The information of PDs is encoded into the finitely dimensional vectors produced by the control grid of the function, so that one can reconstruct precisely the persistence surface from the vectors. Notably, while in PI (Adams et al., 2017), the persistence surface is discretized into pixels to generate the vectors. Then, some information is "lost" during the discretization. The vectors produced by our method preserve the whole information of the persistence B-spline function, thereby improving effectiveness and efficiency.
To efficiently compute the persistence B-spline function, the technique of progressiveiterative approximation for least squares B-spline function fitting (LSPIA) (Deng and Lin, 2014) is used to robustly generate the vectors. We also proved that the vectorizing representations have stability with respect to the 1-Wasserstein distance under some conditions. Generally, the proposed framework for PD representation is (1) stable with respect to the Wasserstein distance in different choices of eminence values, (2) efficient to compute the control grid to form the vector, and (3) flexible to select the "importance" values (local eminence values and global eminence values) for different practical tasks.

Background
In this section, we provide brief introductions on B-spline function and persistent homology, related to the core concepts of the proposed method.

B-Spline Function
First, B-spline is briefly introduced, and readers can refer to Piegl and Tiller (2012) for details. A B-spline of order n is a piecewise polynomial function of degree n − 1 with a variable t. A series of non-decreasing sequence {· · · ξ i−1 , ξ i , ξ i+1 , · · · } denotes knots at which the pieces of polynomials join. Given a sequence of knots, the B-splines of order 1 are defined by The B-splines of kth-order are defined in the way of recursion by and thereby, nth-order B-spline depends on (k −1)-th B-spline. The B-spline has the weight property, i.e., i B i,n (t) = 1, for all n ≥ 1. Moreover, the continuity of B i,n (t) depends on the order n and the multiplication. Specifically, for a B-spline basis function of order n, if r knots are coincident, then the first n − r − 1 derivatives of B i,n (t) are continuous across that knot. If the knots are distinct, then the first n − 2 derivatives are continuous across each knot. And if the distance between each knot and its predecessor, i.e., ξ i+1 − ξ i , is the same, the knot vector and the corresponding B-splines are called uniform. A (p, q)-order uniform B-spline function, of which the degree is (p − 1, q − 1), is defined by where {P ij }'s are the control points in the space with the size h × h, and the pth-order B-spline basis B i,p (s) and qth-order B-spline basis B j,q (t) are defined, respectively, on the knot sequences on [0, 1] which makes the property of endpoint interpolation hold. Note that, in our implementation, the bi-cubic uniform B-spline function, i.e., p = q = 4, is used for data fitting.

Persistent Homology
Persistent homology Edelsbrunner et al., 2002) aims to infer the topological structures of the underlying shape from a point cloud sampled from it. Persistent homology is based on formalism in algebraic topology, namely, homology (Hatcher, 2002;Kaczynski et al., 2006), which describes the "holes" with dimensions of a topological space in the language of algebra. Intuitively, a homology class in the k-th homology group of a topological space represents a k-dimensional "hole". To introduce the concept of PD, we intuitively explain the filtration based on the Vietoris-Rips (V-R) complex, the meaning of PDs, and the stability theorems of PD. The textbook  provides further information.
Let X = {p 1 , · · · , p n } be a finite point set sampled from a shape in a metric space with the metric d. A ball with its center p i and radius r is denoted as B(p i , r) = {x : d(x, p i ) ≤ r}, i = 1, · · · , n. For a fixed r ≥ 0, a k-simplex is defined as a subset τ = {p i 0 , · · · , p i k } of X which satisfies that the balls of any two points intersect, i.e., B(p i j 1 , r) ∩ B(p i j 2 , r) = ∅ for any p i j 1 , p i j 2 ∈ τ . The set of these simplexes with different dimensions forms a simplicial complex, called the V-R complex of X with radius parameter r, denoted as V R(X, r). For radius parameters r j 1 ≤ r j 2 , V R(X, r j 1 ) is a subcomplex of V R(X, r j 2 ), and a nested sequence of simplicial complexes F(X) = {V R(X, r j ) | r j ≥ 0} is obtained, called a V-R complex filtration.
As parameter r increases, some homological invariants emerge, and some disappear because of the appearance of simplexes. Mathematically, for r j 1 ≤ r j 2 , the inclusion between V R(X, r j 1 ) and V R(X, r j 2 ) induces a homomorphism κ r j 2 r j 1 from the k-th homology group H k (V R(X, r j 1 )) to H k (V R(X, r j 2 )). The k-th persistent homology of the filtration F(X) is defined by the collection of homology groups and the induced homomorphisms, specifically demonstrated in Zomorodian and Carlsson (2005). Intuitively, a k-th homological class is born at r j 1 if it emerges in H k (V R(X, r j 1 )), and it dies at r j 2 if it disappears in H k (V R(X, r j 2 )). Its lifetime r j 2 − r j 1 is called the persistence of the homological class. A k-th PD is a multi-set of two-dimensional points P D k = {(x i , y i ) | i ∈ I} that summarizes the birth-death pairs (x i , y i ) = (r j 1 , r j 2 ) of the k-th homological class. On the tame condition in Zomorodian and Carlsson (2005), we can assume that all PDs have finite cardinality |I| < ∞.
To measure the similarity of two k-th PDs denoted as P D (1) and P D (2) , distance metrics are introduced. Let ∆ = {(x, x) | x ≥ 0} be the diagonal point set with infinite multiplicity, and b : P D (1) ∪ ∆ → P D (2) ∪ ∆ be a bijection. The bottleneck distance between two PDs is defined by and for 1 ≤ p ≤ ∞, the p-Wasserstein distance is defined by where b ranges over all bijections from P D (1) ∪ ∆ to P D (2) ∪ ∆. One of the most important properties of PDs is that the transformation from data to a PD is continuous with respect to the distance metrics. As proved in Cohen-Steiner et al. (2007); Chazal et al. (2014), the bottleneck distance between two k-th PDs from point sets X and Y is controlled by the

Method and Algorithm
In this section, we propose a method for transforming a PD into a finite-dimensional vector. As illustrated in Fig. 1, PD points are transformed from the birth-death coordinate to the birth-persistence coordinate and normalized. Then, an eminence value is assigned to each transformed 2D point. Subsequently, the eminence values are fitted by a uniform bi-cubic B-spline function. Finally, a finite-dimensional vector is obtained by concatenating the control grid of the B-spline function. A PD with a certain dimension k is a multi-set above the diagonal, i.e., where I = {1, 2, · · · , M } is an index set. To present our method, some assumptions should be made as follows: Assumption 1 (1) All PDs we considered contain a finite number of points over the diagonal, i.e., |I| < ∞; (2) All points in all PDs are bounded, 1 and denote m as the maximum of the difference between upper and lower bounds of the birth coordinates and the persistence, i.e., letting x min = min l∈I (x l ) and x max = max l∈I (x l ), These assumptions are based on the theory of a tame persistent homology introduced in Section 2.2 and practical observation. This section is organized as follows. In Section 3.1, the coordinate transformation and eminence function are introduced. In Section 3.2, the persistence B-spline grid, the key concept of our framework, is defined. In Section 3.3, the algorithm based on the approximation technique LSPIA is summarized for the vectorization of PDs.

Coordinate Transformation and Eminence Functions
As the transformation adopted in Adams et al. (2017), PDs are transformed into the birth-persistence coordinates and normalized into a unit region. According to Assumption 1, m is the maximum of the difference between upper and lower bounds of the birth coordinates and the persistence in all PDs we considered, and the lower bound of the birth coordinates in all PDs is denoted as c. 2 The coordinate transformation φ : The normalized birth-persistence coordinate is denoted by (s l , t l ) = φ(x l , y l ) for (x l , y l ) ∈ P D (6).
To measure the significance of a birth-death pair in a PD, an eminence value is assigned to each pair by defining different eminence functions. The principles of designing an eminence function depend on the understanding of different persistence in a practical task. Generally, the birth-death pair close to the diagonal is thought to be noise and marked with a small eminence value. By contrast, the pair with large persistence is marked with a large eminence value. Therefore, to design the eminence function, we propose two alternatives, namely, local eminence function and global eminence function.
Local Eminence Function: Given a birth-death pair (x, y), the local eminence function assigns to each PD point a local eminence value that depends only on its persistence w = y − x. To make the stability hold in Section 4, the univariate local eminence function should satisfy the following conditions: 1. Reduced 0-th PD is used to overcome the emergence of an infinitely persistent feature. 2. In particular, in classification tasks, m and c depend on the PDs of training data. For the test data, the PDs are truncated by m and c.
Global Eminence Function: To consider persistence of a PD point and the distribution of all points on a PD, the global eminence function assigns a global eminence value to a PD point. Given a PD, the global eminence function has the form where K(x, y) is a basis function and (x l , y l ) ∈ P D. For stability in Section 4, the global eminence function is supposed to follow the principles: (1) the basis function K(x, y) is bounded and differentiable, i.e., K(x, y) ∞ < ∞; (2) the first-order partial derivatives of K(x, y) with respect to x and y are bounded, i.e., sup (x,y)∈R 2 ∇K(x, y) 2 < ∞.

Persistence B-Spline Grid
For each transformed coordinate (s l , t l ) in a PD, an eminence value z l is assigned through a local eminence function, i.e., or through a global eminence function, i.e., z l = E G,P D (x l , y l ).
Then, the eminence values are fitted by a uniform bi-cubic B-spline function with h × h size of the control points as follows: where h is the grid size satisfying h ≥ 4, and the uniform cubic B-spline basis B i (s) and B j (t) are both defined on the uniform knot sequence To obtain the uniform bi-cubic B-spline function (12) that fits the eminence values, the following least-squares problem should be solved: It is equivalent to solving the normal equation, where Z = ( z 11 , z 12 , · · · , z 1h , z 21 , z 22 , · · · , z 2h , · · · , z h1 , z h2 , · · · , z hh ) T , and the matrix B is the collocation matrix of the B-spline basis at {(s l , t l ) | l ∈ I}. By arranging the B-spline basis in a vector and denoting B i (s)B j (t) as B (i−1)×h+j (s, t), i.e., the collocation matrix B at {(s l , t l ), l ∈ I} can be written as We define some terminologies as follows.
Definition 1 (Persistence B-spline function, PBSG, and vector of PBSG) (1) The uniform bi-cubic B-spline function (12) which fits the eminence values is called the persistence B-spline function.
(2) The control grid of the persistence B-spline function (12) is called the persistence B-spline grid (PBSG).
(3) The vector Z generated by arranging PBSG with the lexicographic order, shown in (16), is called the vector of PBSG.
In the persistence B-spline function (12), the size of its control grid (PBSG) is h × h. We call h or h × h the size of the PBSG.

Definition 2 (Local eminence vector and global eminence vector)
(1) The vector Z (17) is called a local eminence vector, if its elements are taken from a local eminence function by Eq. (10).
(2) The vector Z (17) is called a global eminence vector, if its elements are taken from a global eminence function by Eq. (11).
The local and global eminence vectors are both called eminence vectors.
If a PBSG is generated by fitting a local eminence vector, we call it a local PBSG. Similarly, if a PBSG is produced by fitting a global eminence vector, we call it a global PBSG.

Algorithm for PBSG
As stated to obtain the solution of the least-squares fitting problem (14), the normal equation (15) should be solved. However, the coefficient matrix of the normal equation (15) can be either non-singular or singular. When it is singular or near singular, special treatments should be used to obtain a robust solution. In this article, we use LSPIA (Deng and Lin, 2014;Lin et al., 2018) to effectively and efficiently solve the linear system (15). Regardless of whether the coefficient matrix is singular, LSPIA can converge to the solution of the linear system (15) in a unified format, without any special treatment. When the coefficient matrix is singular, an infinite number of solutions are possible, and LSPIA converges to the solution with minimum Euclidean norm (Lin et al., 2018).
LSPIA begins with an initial B-spline function where the initial control points are all zero, i.e., z be the difference between the eminence values (10)(11) and the corresponding values S (0) (s l , t l ) on the function S (0) , and take the initial differences for the control points as where the setting of the parameter C is discussed later. Then, adding the difference for the control points (20) to the control points z ij , we obtain the new control points which then generate a new function S (1) (s, t). Suppose that the k-th function S (k) (s l , t l ) is obtained. We calculate δ and the (k + 1)st B-spline function is generated as In general, the iteration will stop when the given iteration time N is reached, or the difference error is less than a threshold ε, i.e., For the convenience to show the stability of PBSGs in Section 4, we set a maximum iteration time N to stop the LSPIA iteration.
The parameter C in (20) and (21) is a positive constant. If the LSPIA algorithm converges (Lin and Zhang, 2013;Deng and Lin, 2014;Lin et al., 2018). In our implementation, we set C = max 1≤i,j≤h M l=1 B i (s l )B j (t l ) to speed up the convergence rate of LSPIA. Based on the result of convergence of a general LSPIA format, the convergence of this algorithm is proved in Appendix 8. Finally, we summarize the algorithm to compute the PBSG by LSPIA in Algorithm 1.

Stability of Vectorizing Representation
In this section, we prove the stability of PBSG with respect to 1-Wasserstein distance between PDs by the following steps: (1) We first show the stability of local and global eminence values and the transformed coordinates.
(2) A core lemma on the LSPIA iteration method is shown to bridge the vectors of eminence values, transformed coordinates, and PBSG.
(3) The stability of local and global PBSG is induced.
be a given PD and P D (2) be the perturbed PD, respectively. In the 1-Wasserstein distance between these two PDs, the mapping b is the matching bijection that achieves the infimum, satisfying, l )), l = 1, 2, · · · , M. Therefore, we have, In addition, the vectors of transformed coordinates from P D (i) , i = 1, 2, given by (7), are denoted as, respectively. The eminence vectors of P D (i) , given by (17), are denoted as Z (i) , i = 1, 2, respectively. And the vectors of PBSG (16) generated from P D (i) are represented by Z (i) , i = 1, 2, respectively. To prove the stability of PBSG, we use the assumptions mentioned in Section 3 (Assumption 1).

Stability of Eminence Values and Transformed Coordinates
Stability of local eminence values: In the case that the eminence values are assigned via a local eminence function (8), the eminence value is determined only by the persistence of the PD point. Thus, we have z l = E L,P D (x l , y l ) = E L (y l −x l ). The following lemma shows that the difference between the local eminence vectors, i.e., Z (1) and Z (2) , is controlled by the multiplication of the 1-Wasserstein distance between P D (i) , i = 1, 2 and a constant.
Lemma 3 Supposed that the local eminence function E L (8) is differentiable and has a bounded first-order derivative, i.e., sup x∈R |E L (x)| < ∞, we have where sup |E L | denotes sup x∈R |E L (x)|.
Proof It follows by using Lemma 9 in Appendix 7 that for 1 ≤ l ≤ M , where sup x∈R |E L (x)| := sup |E L |. Because for any a, b ∈ R, |a−b| ≤ |a|+|b| ≤ 2 max{|a|, |b|} holds, we have where P Consider the difference of Z (1) and Z (2) measured by 1-norm, and it finally holds that Stability of global eminence values: For a global eminence function determined by all points in a PD, as shown in (9), the following lemma shows that the difference of the global eminence vectors under 1-norm, i.e., Z (1) − Z (2) 1 , is controlled by the multiplication of the 1-Wasserstein distance between P D (i) , i = 1, 2 and a constant.
Lemma 4 Supposed that basis function K(x, y) of the global eminence function in (9) is bounded and differentiable with bounded first-order partial derivatives, we have , and m is given in Assumption 1.
The proof of Lemma 4 is attached in Appendix 7. Stability of coordinate transformation: Recall that birth-death coordinates are transformed to the birth-persistence coordinates through (7). In the following lemma, we show that the vector of the transformed coordinates is stable with respect to the 1-Wasserstein distance between PDs. Specifically, the difference between the vectors of the transformed coordinates under 1-norm is controlled by the multiplication of the 1-Wasserstein distance between P D (i) , i = 1, 2 and a constant.
Lemma 5 For the vectors of transformed coordinates u (i) from P D (i) , i = 1, 2, through the transformation in (7), we have where m is given in Assumption 1.

Stability on LSPIA
In this section, we provide a lemma based on the matrix form of LSPIA when the iteration time N is finite. The LSPIA iterative format shown in (21) can be rewritten in a matrix form, introduced in Lin et al. (2018), as follows: where Z (k) represents the vector of PBSG after k iterations of LSPIA, is a diagonal matrix, C is an iterative parameter (a constant), Z is the eminence vector given in (16), B is the collocation matrix shown in (18), and O is a zero vector. The iterative format in (33) is equivalent to where E is an identity matrix. In the following, we will estimate the difference between the vectors of PBSG Z from P D (i) , i = 1, 2, under 2-norm after k + 1 iterations of LSPIA.
Lemma 6 Let B (i) be the collocation matrices (18) at the transformed coordinate u (i) , i = 1, 2, respectively. For the vectors of PBSG Z where α = 2(h − 3) is the bound of |B i (s)| given in Lemma 12, h is the size of PBSG, C is the iterative parameter given in (23), and E represents the local eminence function E L,P D (8) or global eminence function E G,P D (9).
The proof of Lemma 6 is provided in Appendix 7.

Stability of PBSG
In this section, by using the lemmas proved above, we provide stability results of the vector of local and global PBSG based on a finite iterations N of LSPIA.
The following theorem shows that the vector of local PBSG is stable with respect to the 1-Wasserstein distance between PDs.
Theorem 7 (Stability for Local Eminence Function) Assume that the iterative parameter C satisfies (23) and C ≥ M . For the local eminence vectors Z (i) , produced from a local eminence function E L , and the vectors of PBSG Z where the stability coefficient in which α = 2(h − 3) is the bound of |B i (s)| in Lemma 12, h is the size of PBSG, m is given in Assumption 1, and sup |E L | = sup x∈R |E L (x)|.
Analogously, the following theorem shows that the vector of global PBSG is stable with respect to the 1-Wasserstein distance between PDs.
Theorem 8 (Stability for Global Eminence Function) Assume that the iterative parameter C satisfies (23) and C ≥ M 2 . For the global eminence vectors Z (i) , produced from a global eminence function E G,P D (x, y), and the vectors of PBSG Z where the stability coefficient in which α = 2(h − 3) is the bound of |B i (s)| in Lemma 12, h is the size of PBSG, m is given in Assumption 1, K ∞ = max (x,y)∈R 2 |K(x, y)|, and sup ∇K 2 = sup (x,y)∈R 2 ∇K 2 .
Eqs. (37) and (39) explicitly give upper bounds of the stability coefficients C loc and C glob . However, it is difficult to induce the theoretical supremum for the stability coefficients. In Section 5.2, we estimated the stability coefficients by computing

Experiments
This section is organized as follows. In Section 5.1, we assessed the proposed vectorizing representation with different choices of parameters. In Section 5.2, we estimated the upper bounds of the stability coefficients on data sets of randomly generated PDs. In Section 5.3, we compare the performance of our representation and other methods, namely PIs, PLs, SWK, PWGK, PSSK, PBoW, in classification tasks to distinguish the dynamical systems and classify CAD models. In Section 5.4, we explained the differences between PI and our method, and conducted a comparison experiment to show the better performance of our method.
In all experiments, the 0th or 1st PDs were computed by constructing the filtration of V-R complexes of a point cloud in the space equipped with Euclidean metric by using the Python package Ripser (Bauer, 2017). In classification tasks, five commonly used ML classifiers were used, namely, k-nearest neighbor classifier (kNN), random forest (RF) in Breiman (2001), gradient boosted decision trees (GBDT) in Friedman (2001), logistic regression (LR), and linear support vector machine (LSVM) in Fan et al. (2008). All experiments were run on a PC with an Intel Core(TM) i7-7700 CPU@3.60 GHz×8. The data sets, MATLAB, and Python source code can be found on https://github.com/ ZC119/PB.

Effect of Parameter Choice in PBSG
The vector of PBSG relies on multiple parameters, including the choice of eminence function and its parameter, the iteration time N of LSPIA, and the size of PBSG h. In this section, we evaluated the effect of parameter choice of PBSG on a synthetic data set containing five categories, and measured the classification accuracy as a quantitative index that indicates the influence of different choices of parameters.
As illustrated in Figure 2, the data set contains a circle with radius 0.4; two concentric circles with radii 0.2 and 0.4, respectively; two disjoint circles each with radius 0.2; a cluster of points sampled at random in the unit square; and two clusters of points sampled at random separately in two squares with edge length 0.5. We generated 50 images for each type, and sampled 1, 000 points for each image. The sampled data were perturbed by adding Gaussian noises (η = 0.025). For each classifier, data were split into 70% training set and 30% test set. After performing 100 trials, the average classification accuracy was computed. On each classifier, cross-validation was used to select the optimal parameters of classifiers. To measure the fitting precision of the persistence B-spline function (12), we used the root mean square (RMS) error as follows: where n is the number of points in the given PD, and δ (k) l is presented in (21). Given a collection of PDs, the average RMS error and the maximum RMS error were used to measure the fitting precision on the collection of PDs.
Choice of local eminence function: To choose a desirable local eminence function, we tested several commonly used candidates, including the linear function E L,P D (x, y) = E L (y − x) = b(y − x), the arctan function E L (y − x) = arctan[b(y − x) 2 ], the exponential function E L (y − x) = e b(y−x) − 1, where b is a parameter, and the sigmoidal function where the parameter b means a translation on the standard sigmoidal function. With the sigmoidal function, the eminence values of the PD points with relatively small and medium persistence varied considerably lot, whereas the change of those with markedly large persistence was insignificant. Moreover, by using the sigmoidal function, the average classification accuracy of five classifiers was higher than that of other candidates. Therefore, the sigmoidal function was selected as a local eminence function.
Choice of global eminence function: As shown in (9), the global eminence function relies on a basis function and the persistence of PD points. In our framework, we assume that a two-dimensional basis function is the multiplication of two one-dimensional basis functions, i.e., K(x, y; a 1 , a 2 ) = ϕ(x; a 1 ) · ϕ(y; a 2 ), where a 1 and a 2 are tunable parameters. In the pretest, we set the parameters a 1 = a 2 . We tested commonly employed basis functions, including the Gaussian basis function ϕ(x; a) = e −ax 2 , multiquadric function ϕ(x; a) = 1 + ax 2 , inverse multiquadric function ϕ(x; a) = 1/ √ 1 + ax 2 , and sinc-Gaussian function in Zhang and Ma (2011) as follows: where a is a tunable parameter. Through the pretest, the vectors of PBSG with the sinc-Gaussian function as a global function had superior performance over all other candidates on five classifiers. In addition, Zhang and Ma (2011) showed the C ∞ continuity and the properties such as interpolatory property, partition of unity, and local support under a small tolerance. According to these properties, the basis function and its derivative are easily verified to be bounded. Therefore, we selected the sinc-Gaussian function (43) as the basis function in the global eminence function.
Effect of the parameters in eminence functions: In the pretests, we selected the sigmoidal function in (42) as the local eminence function and sinc-Gaussian function in (43) as the basis function of the global eminence function. We also evaluated the effect of different tunable parameters b in local eminence function and a in the basis function of global eminence function on the performance of classification. Classification was conducted by setting the size of PBSG to be 20 × 20 and the LSPIA iteration time N = 100. The parameter b in the local eminence function was set to be −10, −1, 1, and 10. The parameter a in the global eminence function was 0.1, 1, 5, and 10. The average classification accuracy was measured on five classifiers. The average RMS error was less than 10 −2 and the maximum RMS error was less than 10 −1 . The experimental results are shown in Table 1. The classification accuracy on all classifiers was close to 100%, demonstrating that the classifiers could distinguish the categories of synthetic data. The fluctuation of accuracy on every classifier under different parameters of local and global eminence function was small (less than 1.5%). The classification accuracy was not sensitive to the parameter in the eminence function in the given range, i.e., b ∈ [−10, 10], and a ∈ [0.1, 10].  the iteration time was increased from 100 to 1,000 at increments of 100 with the size of PBSG 20 × 20. The average RMS error was less than 10 −2 and the maximum RMS error was less than 10 −1 . For the effect of iteration time on accuracy, as shown in Fig. 3, the accuracy of the vectors of local PBSG was unchanged after 500 iterations of LSPIA on the classifier GBDT, and the accuracy on the other four classifiers was stable, as shown in Fig. 3(a). The accuracy of global PBSGs was also insensitive to the increase of iteration, as shown in Fig. 3(c).
Effect of the size of PBSG: Analogously, we set b = 10 and a = 1. The PBSG size increased from 10 × 10 to 100 × 100 at increments of 10 on each axis. The average and maximum RMS errors were less than 10 −2 and 10 −1 , respectively. The PBSG size v.s. accuracy curves demonstrate the effect of PBSG size (Fig. 3). The accuracy of local PBSGs fluctuated on the GBDT classifier as the PBSG size increased, and the accuracy on the kNN classifier dropped from 80 × 80 PBSG, while the accuracy was stable on RF, LR, and linear SVM, as shown in Fig. 3(b). Moreover, as shown in Fig. 3(d), the accuracy of global PBSGs on the kNN classifier began to decrease on the 90 × 90 PBSG, whereas the accuracy remained stable on the other classifiers.
In general, according to the results on the synthetic dataset, as LSPIA iteration and PBSG size increased, the classification accuracy of local and global PBSGs on the RF, LR, and linear SVM classifiers was more stable than that on the kNN and GBDT classifiers. The accuracy did not fluctuate much on five classifiers as the iteration and PBSG size ascended from 100 iterations and a 10 × 10 PBSG. A reasonable explanation for this result is that the increase of iteration leads to the reduction of RMS errors and thus robust vectors of PBSG. The increase in the size of PBSG causes a growth in vector dimension and a decrease in RMS error. Consequently, we set the LSPIA iteration to be 100 and the PBSG size to be 20 × 20 in the subsequent experiments.

Bound Estimation of Stability Coefficients
In Section 4, we proved the stability of the vector of PBSG with respect to 1-Wasserstein distance. The stability coefficients (37) and (39) mainly rely on the iterative parameter C, the number of iterations N , and the size of the PBSG h. In this section, the upper bounds of stability coefficients were estimated using numerical experiments performed on the data sets of randomly generated PDs. Specifically, two data sets were generated, each of which includes 50 randomly generated PDs and the perturbed 50 PDs. In the first collection of PDs, 100 points above the diagonal were randomly generated in each PD. In the second collection of PDs, 1, 000 points in each PD were randomly generated.
Given P D (1) and the perturbed P D (2) , let Z (i) (N ) , i = 1, 2, be the vectors of PBSG generated after N iteration of LSPIA, respectively. In the experiments, we estimated the stability coefficients by computing /W 1 (P D (1) , P D (2) ), using the same eminence functions as those in Section 5.1, and setting the parameter b = 1 and a = 100 in the local and global eminence functions, respectively. On the first collection of PDs that contain 100 points, we set C = M in local PBSG and C = M 2 in global PBSG. When N ≥ 1000, the maximum RMS error of local PBSGs was not more than 0.3. When N ≥ 10 5 , the maximum RMS error of global PBSGs was not more than 1.3. On the second collection of PDs that contain 1, 000 points, we set C = max 1≤i,j≤h M l=1 B i (s l )B j (t l ) in both local and global PBSGs. In this setting, when N ≥ 50, the maximum RMS error of local PBSGs was less than 0.9, and when N ≥ 100, the maximum RMS error of global PBSGs was less than 3.0.
As shown in Fig. 4, the averages of stability coefficient C loc and C glob were computed with error bars. In all of the experiments, the stability coefficients were all less than 10.0. After reaching a small RMS error, under a fixed PBSG size h, the average stability coefficients and their standard deviations were nearly unchanged with an increment of N . Meanwhile, under a fixed iteration N , the average stability coefficients and their standard deviations increased with the increment of PBSG size h from 20 to 50 with an increment of 10. Specifically, in Fig. 4 (a-c), the average stability coefficient when h = 40 was approximately two times more than that when h = 20. In Fig. 4 (d), the average stability coefficient when h = 40 was nearly five times more than that when h = 20.
In conclusion, the experiment demonstrates that after the persistence B-spline function reached a small RMS error, the practical stability coefficients did not sharply change as the  iteration and the PBSG size increased. Under a fixed PBSG size, the stability coefficient did not change much with the increment of iteration. Under a fixed iteration time, the stability coefficient was proportional to the PBSG size.

Applications on Classification
To compare with other methods including computing distance matrix of p-Wasserstein distance, PIs, PLs, kernel methods (PSSK, PWGK, SWK), and PBoW, some experiments were conducted on a data set generated from a discrete dynamical system and a 3D CAD model data set. In these experiments, software (Kerber et al., 2017) was used to compute the Wasserstein distances and the bottleneck distance. For PLs, the Persistence Landscapes Toolbox (Bubenik, 2015) was used. For PIs, the open-source MATLAB codes (Adams et al., 2017) were used. The PIs were computed with resolutions of 20 × 20, and the parameter σ in PIs was adjusted to guarantee the best result.

Classification of Dynamical System Data with Different Parameters
Dynamical systems can simulate some natural phenomena. Different parameters in a dynamical system lead to different types of behavior. In this experiment, PLs, PIs, PBSGs, and PDs were tested to observe their performances when classifying a complicated dynamical system with different parameters.
The dynamical system proposed by Lindstrom (2002) describes a discrete food chain model defined by where The variables X, Y and Z are related to the trophic levels of the food chain system. Osipenko (2006) studied the attractor of the system with fixed parameters M 1 = 1.0, M 2 = M 3 = 4.0 and variable parameter M 0 ∈ [3.00, 3.65]. Referring to Chapter 17 in Osipenko (2006), we set M 0 = 3.0, 3.3, 3.48, 3.54, 3.57, 3.532, 3.571, 3.3701, 3.4001, respectively, Thus, nine classes of results were obtained with various parameters M 0 . For each class, the initial values (X 0 , Y 0 , Z 0 ) were generated randomly in the region (1, 2) × (0, 1) × (0, 1), and the system was iterated 2,000 times to obtain the points in 3D Euclidean space. Each combination of parameters was repeated 50 times, and eventually 450 results were generated from the system. In the experiment, the data were randomly split into a 70%-30% training-test data set, and the classification was conducted 10 times so that the average accuracy and its error bars were computed. In Fig. 5, we visualize an instance of a 3D point cloud of the dynamical system, the corresponding PD, the persistence B-spline function, and PBSG. Comparison of kNN with different metrics: By computing k nearest neighbors of a target point under a metric of the vector space, the kNN classifier determines the label of the target through voting according to those of k nearest neighbors. The metric of vector space plays an important role. To clarify the effect on the accuracy of the kNN classifier with different metrics, the performance of p-Wasserstein (W p ) distance matrix, PLs, PIs, local and global PBSGs was obtained under the metrics L 1 , L 2 , and L ∞ , respectively. For the method of computing the p-Wasserstein distance matrix, the chosen metrics L 1 , L 2 , and L ∞ represent 1-Wasserstein, 2-Wasserstein, and the bottleneck distance, respectively. For other methods, the metrics L 1 , L 2 , and L ∞ mean the Manhattan distance, Euclidean distance, and Chebyshev distance, respectively. As shown in Table 2, the classification accuracy on the kNN classifier and the time cost are given regarding the W p distance matrix, local PBSG (b = 1), global PBSG (a = 10), PIs (σ = 10 −3 ), and PLs. The vectors of global  PBSG performed the best on the kNN classifier with L 1 , L 2 , and L ∞ metrics. PIs, PLs, and PBSGs were computed efficiently, whereas the direct computation of the W p distance matrix required considerable time. For PIs and PBSGs, they performed better on the kNN classifier by setting the metrics to be Euclidean distance. Local PBSGs had relatively low accuracy on the kNN classifier with the L ∞ metric. For the global PBSGs, the difference of accuracy was within approximately 5.0% among kNN classifiers with L 1 , L 2 , and L ∞ . Comparison among PBSGs and other competitors: We first measured the accuracy of local and global PBSGs under different parameters on five classifiers, shown in Table 3. For the performance on different classifiers, global PBSGs (a = 1) had the best performance on four classifiers. Local PBSGs also performed well on RF, GBDT, LR, and LSVM classifiers. The accuracy did not fluctuate much by tuning the parameter of the eminence functions.
To compare local and global PBSGs with other competitive methods, the classification accuracy on the SVM classifier was measured, shown in Table 4. In experiments, the hyper-parameters of SVM were determined by cross-validation, and for PWGK, the parameters K and ρ were selected to be 100 and 0.01, respectively. For the performance of different methods on SVM, global PBSGs (a = 1) performed the best. The performance of SWK was close to that of global PBSGs. Local PBSGs performed better than did other methods except for SWK.

Classification of 3D CAD Models
PDs can represent topological features with different geometric scales. In this section, the vectorization of PDs is applied to classify some 3D CAD models. The benchmark data set ModelNet40 (Wu et al., 2015) was used, which contains 12,311 CAD models classified into 40 object categories. Initially, we computed the global PBSGs from the point clouds of all the models by uniform sampling, and trained the classifiers GBDT and SVM by the training data. The accuracy of PBSG on the full benchmark was approximately 35.2% on the GBDT classifier and approximately 28.0% on the SVM classifier to classify shapes from 40 object categories. Given that similar shapes may have different semantics, the geometrically and topologically relying methods could not recognize those shapes. For example, the four categories, namely, bottle, bowl, cup, and flowerpot, have similar geometric and topological structures. We, therefore, selected seven categories, namely, piano, person, airplane, radio, guitar, wardrobe, and sofa, with 1,970 models in total. The 1st PDs of each category were different so that the topologically relying methods could distinguish them.
As shown in Table 5, global PBSGs achieved the best performance on five classifiers, and local PBSGs performed better than PIs on kNN, RF, GBDT, and LSVM. In addition, as shown in Table 6, the accuracy of global PBSGs reached the best performance in all compared methods, 1.0% more than the second-best (SWK). As for the time cost, PBoW took the least time (approximately 200 seconds). Local and global PBSGs were also efficient. In the process of PBSG computation, approximately 200 and 400 seconds were needed to compute the local and global eminence values, respectively, and approximately 200 seconds were needed in generating the PBSGs via LSPIA. However, kernel methods (SWK, PWGK, and PSSK) needed over 2 hours in training and testing stages because the construction of kernel matrices was time-consuming. In conclusion, the global PBSGs were efficient, reaching the best accuracy in the classification of 3D shapes.

Comparison between PBSGs and PIs
From the perspective of data fitting, PBSG and PI both fit the weighting values on the birth-death pairs in a PD with a function. Specifically, while PI produces a quasiinterpolation surface (Cheney, 1995) by fitting the weights based on Gaussian kernels, PBSG uses a uniform bi-cubic B-spline function to fit the weights. However, the vectorization of the two methods is different. On the one hand, PI discretizes the domain and integrates on each pixel to produce a vector, thereby losing some information when discretization. On the other hand, PBSG uses the control coefficients of the B-spline function to form a vector. Given that a uniform bi-cubic B-spline function is entirely determined by its control coefficients, the vectorization of PBSG does not lose any information. Loosely speaking, the relationship between PBSG and PI is similar to that between vectorized and rasterized images.
To compare the PIs and PBSGs in practical cases, we designed some experiments on the data sets of the 3D dynamical system and 3D CAD models. In these experiments, we set the basis function K(x, y) in the global eminence function of PBSG as the Gauss function: where σ is a parameter and σ 2 is the variance. Thus, the weights generated at points on a PD are the same in PBSG and PI. For convenience, we denote by PBSG Gauss a PBSG with Gaussian function as the basis function. On experiments of the 3D dynamical system, the resolution of PIs and the size of PBSG Gauss were both fixed as 20 × 20, and the parameters σ in PI and PBSG Gauss were set as 10 −1 , 10 −2 , 10 −3 , and 10 −4 , respectively. On the experiments of the 3D CAD data set, while the parameters σ were fixed to be 10 −2 in PI and PBSG Gauss , the resolution of PIs and the size of PBSG Gauss increased from 20 × 20 to 40 × 40 at increments of 10. In PBSG Gauss , the iteration of LSPIA was set to be 100, and the maximum RMS errors were less than 0.05.  On the one hand, as shown in Table 7, with the dimensions of PIs and PBSG Gauss fixed, on the data set of the 3D dynamical system, the classification accuracy of PBSG Gauss was higher than that of PI on all classifiers as the parameter σ changed. As shown in Table 8, on the 3D CAD data set, PBSG Gauss performed better than did PI in most cases as the parameter σ varied. On the other hand, as listed in Table 9, with the parameter σ unchanged, on the data set of the 3D dynamical system, PBSG Gauss had better performance than did PI on all classifiers as the PBSG size (the resolution for PIs) increased from 20×20  to 40 × 40. In Table 10, the accuracy of PBSG Gauss on the 3D CAD model data set was higher than that of PI on kNN, RF, GBDT, and LSVM as the resolution increased. In conclusion, when the weights of PBSG were the same as that of PI, the performance of PBSG Gauss was better than PI in most cases. It should mention that the PBSG with Eq. (43) as the basis function performs better than PBSG Gauss .

Discussion
In this paper, we proposed a stable and flexible framework of PD vectorization based on a persistence B-spline function using the technique of LSPIA, i.e., PBSG, and proved the 1-Wasserstein stability of the vector of PBSG. The key idea of vectoring PDs is to fit the eminence values of points in a PD by a B-spline function, and reshape its control grid as a finitely dimensional vector. The LSPIA technique guarantees the robustness of PBSG computation.
Experiments were conducted to evaluate PBSGs and compare their performance with other competitive methods. As the experiments showed, local and global PBSGs were insensitive to the setting of LSPIA iteration and the size of PBSG. In practical cases, the stability coefficients did not change sharply as the iteration of LSPIA and the PBSG size increased. In the applications of classification, global PBSG achieved the best performance compared with other competitive approaches, and local PBSG also had good performance.
As for the differences between PBSG and PI, the PI surface can be considered a quasiinterpolation surface based on Gaussian basis. The vectorization of PI relies on the discretization of the domain, potentially causing information loss of the PI surface. How-ever, PBSG preserves the complete information of a uniform bi-cubic B-spline function. The experiments showed that PBSG performed better than did PI in most classification tasks. Moreover, the framework of PBSGs provides flexibility for users to select a variety of eminence functions obeying the satisfied principles, especially when dealing with diverse real-world applications.

Appendix: Lemmas and Proofs in Section 4
The following lemmas are employed in Section 4. Lemma 9 is an extended version of the mean value theorem in calculus.
Proof It follows by the mean value theorem, the Cauchy-Schwarz inequality, and the extreme value theorem. Refer to Rudin (1976).
Lemma 10 gives the relationship of norms in a finitely dimensional space.
Then there exist finite positive constants C m and C M such that C m ||x|| p 1 ≤ ||x|| p 2 ≤ C M ||x|| p 1 , for any x ∈ R n . Especially, the following inequities hold: Proof It follows immediately from Corollary 5.4.5 in Horn and Johnson (2013), and the special cases are given in Eq. (5.4.21) in Horn and Johnson (2013).
Lemma 11 provides some properties of the matrix norm · 2 used in the proofs of stability theorems in Section 4.
Lemma 11 Let · 2 be the matrix norm on R m×n induced by the l 2 -norm on R n , i.e., for A = {a ij } m,n i=1,j=1 ∈ R m×n , A 2 = max x 2 =1,x∈R n Ax 2 . For A, B ∈ R m×n and x ∈ R n , the matrix norm satisfies the following properties: (1) cA 2 = |c| A 2 , for any c ∈ R; (2) A + B 2 ≤ A 2 + B 2 ; (3) Ax 2 ≤ A 2 x 2 , for any A ∈ R m×n and x ∈ R n ; (4) A 2 = ρ(A), the largest singular value of A; Proof It follows from the definition of matrix norm in Section 5.6 in Horn and Johnson (2013), Theorem 5.6.2(b), and Example 5.6.6 in Horn and Johnson (2013).
The following lemma is provided to show that the first-order derivative of B-spline basis functions defined on the uniform knot sequence (13) is bounded.
Proof Due to the non-negativity and the weight property of B-spline basis functions in Piegl and Tiller (2012), we have B i,4 (s) ≥ 0 and h i=1 B i,4 (s) = 1, for s ∈ [0, 1]. It implies that 0 ≤ B i,4 (s) ≤ 1. As given in Section 2.1, because the knots are distinct and uniform on (0, 1), the first 2-derivatives are continuous across each knot, and therefore, B i,4 (s) is C 2 -continuous on [0, 1]. According to the formula of the first-order derivative of a k-degree basis function B i,k+1 (s) that in our implementation where k = 3, the bound of the derivative is estimated on (0, 1) as where h is given by the uniform knot sequence (13), and in our case, h is the size of the control grid. Hence, we have |B i,4 (s)| ≤ α = 2(h − 3), which gives an upper bound. (y Meanwhile, it follows that For the second half of the inequity in (49), it follows by using the fundamental theorem of calculus given in Lemma 9 in Appendix 7 that where l ) and x (2) = (x (2) i − x (2) l , y i − y (2) l ). For convenience, let sup ∇K 2 = sup (x,y)∈R 2 ∇K 2 . Since x 2 ≤ x 1 for any x ∈ R 2 , it induces from (50) that In addition, y l − x (2) l is bounded by m (Assumption 1). It is induced from (49) that Notice that i − x i − y l ) l ) .

Appendix: Convergence and Computational Complexity of LSPIA
In our theoretical framework, the constant C is larger than M 2 , where M is the number of points in a PD. In practical implementation, to accelerate the convergence of LSPIA, C is set to be B T ∞ . To clarify the convergence of LSPIA, the following theorem is summarized from Lin and Zhang (2013), Deng and Lin (2014) and Lin et al. (2018).
Theorem 13 Let B be the B-spline collocation matrix given by (18), and B T B Z = B T Z be the normal equation (15) of the least-square fitting problem (14). If the spectral radius ρ(ΛB T B) ≤ 1, then whatever B T B is singular or not, the LSPIA is convergent: (1) If B T B is non-singular, the iterative solution given in (34) converges to the solution of the linear system, i.e., (B T B) −1 B T Z; Then, by multiplying B at both sides of (74), it holds that BΛB T (Bx) = λ(Bx).
Finally, the computational complexity and memory consumption per iteration step of the LSPIA method are provided. Each iteration step of the LSPIA method requires (2M + 1)h 2 multiplications and 2M h 2 additions. Meanwhile, it is required to store two h 2 × 1 vectors, Z (k) and Z (k+1) , an M × 1 vector Z, and an M × h 2 matrix B. Therefore, each iteration requires (M + 2)h 2 + M unit memory.