Topological Feature Vectors for Chatter Detection in Turning Processes

TL;DR

This paper introduces a topological data analysis approach using persistence diagrams to detect chatter in turning processes, achieving high accuracy and fast computation suitable for real-time applications.

Contribution

It proposes a novel topological feature-based method for chatter detection that outperforms traditional wavelet and EEMD methods in accuracy and efficiency.

Findings

Achieved up to 97% accuracy in chatter detection.

Topological features outperform traditional methods in certain configurations.

Runtime reduced to under a second enabling online detection.

Abstract

Machining processes are most accurately described using complex dynamical systems that include nonlinearities, time delays, and stochastic effects. Due to the nature of these models as well as the practical challenges which include time-varying parameters, the transition from numerical/analytical modeling of machining to the analysis of real cutting signals remains challenging. Some studies have focused on studying the time series of cutting processes using machine learning algorithms with the goal of identifying and predicting undesirable vibrations during machining referred to as chatter. These tools typically decompose the signal using Wavelet Packet Transforms (WPT) or Ensemble Empirical Mode Decomposition (EEMD). However, these methods require a significant overhead in identifying the feature vectors before a classifier can be trained. In this study, we present an alternative approach based on featurizing the time series of the cutting process using its topological features. We first embed the time series as a point cloud using Takens embedding. We then utilize Support Vector Machine, Logistic Regression, Random Forest and Gradient Boosting classifier combined with feature vectors derived from persistence diagrams, a tool from persistent homology, to encode chatter's distinguishing characteristics. We present the results for several choices of the topological feature vectors, and we compare our results to the WPT and EEMD methods using experimental turning data. Our results show that in two out of four cutting configurations the TDA-based features yield accuracies as high as 97%. We also show that combining Bezier curve approximation method and parallel computing can reduce runtime for persistence diagram computation of a single time series to less than a second thus making our approach suitable for online chatter detection.

Figures

Tables

Table 1. Table 1 : The number of tagged datasets for each overhang case.

Overhang length (cm (inch))	Stable	Mild chatter	Chatter	Total
5.08 (2)	19	8	12	39
6.35 (2.5)	9	4	3	16
8.89 (3.5)	9	3	2	14
11.43 (4.5)	13	4	5	22

Table 2. Table 2 : Feature matrix for persistence landscapes λ 1 subscript 𝜆 1 \lambda_{1} and λ 3 subscript 𝜆 3 \lambda_{3} corresponding to persistence diagrams X 1 subscript 𝑋 1 X_{1} through X n subscript 𝑋 𝑛 X_{n} . The entries in the cells are the values of each of the features y i , j b i superscript subscript 𝑦 𝑖 𝑗 subscript 𝑏 𝑖 y_{i,j}^{b_{i}} where i 𝑖 i is the persistence landscape number, j 𝑗 j is the persistence diagram number, while b i subscript 𝑏 𝑖 b_{i} is the feature number corresponding to the i 𝑖 i th landscape.

Persistence Diagrams	Label	${\lambda }_1$				${\lambda }_3$
$X_1$	1	$y_{1,1}^1$	$y_{1,1}^2$	$\dots$	$y_{1,1}^{\left|{𝐛}_1\right|}$	$y_{3,1}^1$	$y_{3,1}^2$	$\dots$	$y_{3,1}^{\left|{𝐛}_3\right|}$
$X_2$	0	$y_{1,2}^1$	$y_{1,2}^2$	$\dots$	$y_{1,2}^{\left|{𝐛}_1\right|}$	$y_{3,2}^1$	$y_{3,2}^2$	$\dots$	$y_{3,2}^{\left|{𝐛}_3\right|}$
$⋮$	$⋮$	$⋮$	$⋮$		$⋮$	$⋮$	$⋮$		$⋮$
$X_n$	1	$y_{1,n}^1$	$y_{1,n}^2$	$\dots$	$y_{1,n}^{\left|{𝐛}_1\right|}$	$y_{3,n}^1$	$y_{3,n}^2$	$\dots$	$y_{3,n}^{\left|{𝐛}_3\right|}$

Table 3. Table 3 : Feature matrix for persistence images. The notation I i , j k superscript subscript 𝐼 𝑖 𝑗 𝑘 I_{i,j}^{k} for each feature uses k 𝑘 k to indicate the persistence diagram index, while i , j 𝑖 𝑗 i,j are row and column numbers, respectively, in the persistence image.

Persistence Diagrams	Label	Persistence Image
$X_1$	1	$I_{1,1}^1$	$\dots$	$I_{1,100}^1$	$I_{2,1}^1$	$\dots$	$I_{2,100}^1$	$\dots$	$I_{100,1}^1$	$\dots$	$I_{100,100}^1$
$X_2$	0	$I_{1,1}^2$	$\dots$	$I_{1,100}^2$	$I_{2,1}^2$	$\dots$	$I_{2,100}^2$	$\dots$	$I_{100,1}^2$	$\dots$	$I_{100,100}^2$
$⋮$	$⋮$	$\ddots$
$X_n$	1	$I_{1,1}^p$	$\dots$	$I_{1,100}^p$	$I_{2,1}^p$	$\dots$	$I_{2,100}^p$	$\dots$	$I_{100,1}^p$	$\dots$	$I_{100,100}^p$

Table 4. Table 4 : Feature matrix for path signatures for n 𝑛 n persistence diagrams and using the first λ 1 subscript 𝜆 1 \lambda_{1} and second λ 2 subscript 𝜆 2 \lambda_{2} persistence landscapes.

Diagrams

Label

${\lambda }_1$

${\lambda }_2$

$X_1$

1

$S_{1,1}^1$

$S_{1,1}^2$

$S_{1,1}^{1,1}$

$S_{1,1}^{1,2}$

$S_{1,1}^{2,1}$

$S_{1,1}^{2,2}$

$S_{2,1}^1$

$S_{2,1}^2$

$S_{2,1}^{1,1}$

$S_{2,1}^{1,2}$

$S_{2,1}^{2,1}$

$S_{2,1}^{2,2}$

$X_2$

0

$S_{1,2}^1$

$S_{1,2}^2$

$S_{1,2}^{1,1}$

$S_{1,2}^{1,2}$

$S_{1,2}^{2,1}$

$S_{1,2}^{2,2}$

$S_{2,2}^1$

$S_{2,2}^2$

$S_{2,2}^{1,1}$

$S_{2,2}^{1,2}$

$S_{2,2}^{2,1}$

$S_{2,2}^{2,2}$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$⋮$

$X_n$

1

$S_{1,n}^1$

$S_{1,n}^2$

$S_{1,n}^{1,1}$

$S_{1,n}^{1,2}$

$S_{1,n}^{2,1}$

$S_{1,n}^{2,2}$

$S_{2,n}^1$

$S_{2,n}^2$

$S_{2,n}^{1,1}$

$S_{2,n}^{1,2}$

$S_{2,n}^{2,1}$

$S_{2,n}^{2,2}$

Table 5. Table 5 : Comparison of runtimes (seconds) for embedding parameters and persistence diagram computation of all overhang cases with parallel and serial computing.

	5.08 cm (2 inch)		6.35 cm (2.5 inch )		8.89 cm (3.5 inch)		11.43 cm (4.5 inch)
	Parallel	Serial	Parallel	Serial	Parallel	Serial	Parallel	Serial
Persistence Diagram	9420	84346	3448	23570	2073	11319	4819	37617

Table 6. Table 6 : Runtime (seconds) for embedding parameter computation and persistence diagram computation of a single time series with different methods.

	Embedding Parameters	Persistence Diagram
		Number of points$\approx 1000$	Number of points$\approx 100$			Number of points$\approx 300$
Overhang Distances	Delay and Dimension	Subsampled Point Cloud	Greedy Permutation $n_{perm}=100$	Bézier $r=1$ $spg=100$ (Serial)	Bézier $r=1$ $spg=100$ (Parallel)	Greedy Permutation $n_{perm}=300$	Bézier $r=3$ $spg=100$ (Serial)	Bézier $r=3$ $spg=100$ (Parallel)
5.08 cm (2 inch)	242.63	266.95	0.2	106.00	0.93	3.62	569.78	6.21
6.35 cm (2.5 inch)	191.10	208.95	0.29	106.15	0.85	3.79	538.98	6.19
8.89 cm (3.5 inch)	166.79	296.21	0.18	96.00	0.72	3.87	541.62	6.48
11.43 cm (4.5 inch)	200.51	276.38	0.17	113.86	0.77	3.53	600.29	6.99

Table 7. Table 7 : Runtime (seconds) for performing classification with a single time series for TDA based methods and signal decomposition based ones.

	Topological Data Analysis				Signal Decomposition
Overhang Distances	Persistence Landscapes	Template Functions	Carlsson Coordinates	Persistence Images	WPT	EEMD
5.08 cm (2 inch)	1.01	0.97	0.97	0.97	0.03	0.52
6.35 cm (2.5 inch)	0.92	0.90	0.90	0.87	0.08	0.65
8.89 cm (3.5 inch)	0.81	0.81	0.76	0.76	0.09	0.70
11.43 cm (4.5 inch)	0.87	0.81	0.81	0.81	0.06	0.52

Table 8. Table 8 : Comparison of results for each method where WPT is the Wavelet Packet Transform, and EEMD stand for Ensemble Empirical Mode Decomposition.

Overhang Length cm (inch)	Persistence Landscapes	Persistence Images	Template Functions	Carlsson Coordinates	Kernel Method	Persistence Paths	WPT	EEMD
5.08 (2)	$96.8\%$	$96.4\%$	$91.5\%$	$93.6\%$	$74.5\%$*	$83.0\%$	$93.9\%$	$84.2\%$
6.35 (2.5)	$88.6\%$	$85.8\%$	$89.3\%$	$86.3\%$	$58.9\%$	$84.2\%$	$100.0\%$	$78.6\%$
8.89 (3.5)	$92.2\%$	$93.0\%$	$83.9\%$	$95.7\%$	$87.0\%$	$85.9\%$	$84.0\%$	$90.7\%$
11.43 (4.5)	$68.6\%$	$72.5\%$	$65.1\%$	$72.2\%$	$59.3\%$	$70.0\%$	$87.5\%$	$79.1\%$

Table 9. Table 9 : Closed-form expressions for the first and second level signature paths for landscape functions.

Signature	Equation
$S_1$	$\begin{array}{c}{\int }_{x_1}^{x_n}𝑑P_t^1={\int }_{x_1}^{x_n}𝑑t=x_n-x_1\end{array}$
$S_2$	$\begin{array}{c}{\int }_{x_1}^{x_n}𝑑P_t^2={\int }_{x_1}^{x_2}\frac{y_2-y_1}{x_2-x_1}𝑑t+\dots +{\int }_{x_i}^{x_{i+1}}\frac{y_{i+1}-y_i}{x_{i+1}-x_i}𝑑t+\dots +{\int }_{x_{n-1}}^{x_n}\frac{y_n-y_{n-1}}{x_n-x_{n-1}}\end{array}$
$S_{1,1}$	$\begin{array}{c}\int \int 𝑑P_t^1𝑑P_t^1={\int }_{x_1}^{x_n}\left({\int }_{x_1}^{t_2}𝑑t_1\right)𝑑t_2={\int }_{x_1}^{x_n}\left(t_2-x_1\right)𝑑t_2={\left(\frac{t_2^2}{2}-x_1{t}_2\right)|}_{t_2=x_1}^{t_2=x_n}\end{array}$
$S_{1,2}$
$S_{2,1}$
$S_{2,2}$

Table 10. Table 10 : Persistence Landscape Results with Support Vector Machine (SVM) classifier. Landscape number column represents which landscapes was used to extract features.

	2 inch		2.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$96.8\pm 1.5\%$	$96.7\pm 0.5\%$	$87.0\pm 4.8\%$	$86.9\pm 2.4\%$
2	$86.6\pm 1.7\%$	$91.5\pm 0.8\%$	$87.0\pm 3.8\%$	$86.0\pm 2.0\%$
3	$89.1\pm 2.7\%$	$93.4\pm 0.5\%$	$84.7\pm 4.7\%$	$86.6\pm 1.9\%$
4	$90.5\pm 1.9\%$	$92.9\pm 0.4\%$	$85.8\pm 2.6\%$	$86.1\pm 1.4\%$
5	$90.5\pm 2.3\%$	$92.2\pm 0.7\%$	$88.4\pm 1.5\%$	$85.2\pm 0.6\%$
	3.5 inch		4.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$92.2\pm 5.4\%$	$93.2\pm 1.8\%$	$66.8\pm 5.6\%$	$78.5\pm 3.0\%$
2	$78.3\pm 4.3\%$	$85.2\pm 2.1\%$	$65.3\pm 4.3\%$	$69.3\pm 2.3\%$
3	$82.6\pm 4.3\%$	$83.6\pm 2.0\%$	$68.6\pm 4.8\%$	$72.6\pm 1.5\%$
4	$84.3\pm 6.8\%$	$84.3\pm 2.8\%$	$67.6\pm 5.2\%$	$71.8\pm 2.6\%$
5	$85.7\pm 4.8\%$	$84.3\pm 3.1\%$	$66.8\pm 5.5\%$	$73.0\pm 3.4\%$

Table 11. Table 11 : Persistence Landscape Results with Logistic Regression (LR) classifier. Landscape number column represents which landscapes was used to extract features.

	2 inch		2.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$95.7\pm 1.6\%$	$98.3\pm 0.4\%$	$82.6\pm 4.2\%$	$91.3\pm 2.5\%$
2	$85.5\pm 1.9\%$	$93.7\pm 1.0\%$	$86.3\pm 4.7\%$	$91.5\pm 2.1\%$
3	$88.5\pm 1.1\%$	$93.8\pm 0.6\%$	$84.7\pm 4.3\%$	$92.0\pm 2.6\%$
4	$89.4\pm 1.6\%$	$94.8\pm 0.8\%$	$86.5\pm 5.0\%$	$90.8\pm 2.4\%$
5	$88.3\pm 1.7\%$	$94.7\pm 0.6\%$	$86.0\pm 3.4\%$	$91.1\pm 2.7\%$
	3.5 inch		4.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$91.3\pm 3.4\%$	$96.4\pm 1.8\%$	$63.1\pm 3.5\%$	$82.0\pm 3.4\%$
2	$82.2\pm 6.6\%$	$90.7\pm 2.6\%$	$63.1\pm 3.5\%$	$82.0\pm 3.4\%$
3	$87.8\pm 6.4\%$	$91.4\pm 2.8\%$	$64.2\pm 4.5\%$	$84.1\pm 1.1\%$
4	$90.4\pm 4.3\%$	$91.8\pm 1.5\%$	$63.1\pm 5.5\%$	$82.9\pm 2.4\%$
5	$85.7\pm 4.8\%$	$93.4\pm 1.9\%$	$65.9\pm 4.6\%$	$83.0\pm 2.2\%$

Table 12. Table 12 : Persistence Landscape Results with Random Forest (RF) classifier. Landscape number column represents which landscapes was used to extract features.

	2 inch		2.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$96.1\pm 1.7\%$	$100.0\pm 0.0\%$	$86.7\pm 5.6\%$	$100.0\pm 0.0\%$
2	$87.7\pm 2.2\%$	$100.0\pm 0.0\%$	$88.6\pm 4.0\%$	$100.0\pm 0.0\%$
3	$88.0\pm 1.9\%$	$100.0\pm 0.0\%$	$87.4\pm 2.4\%$	$100.0\pm 0.0\%$
4	$88.8\pm 1.2\%$	$100.0\pm 0.0\%$	$86.7\pm 5.8\%$	$100.0\pm 0.0\%$
5	$89.6\pm 2.2\%$	$100.0\pm 0.0\%$	$86.7\pm 5.8\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$91.3\pm 3.9\%$	$100.0\pm 0.0\%$	$65.1\pm 7.7\%$	$100.0\pm 0.0\%$
2	$80.9\pm 7.1\%$	$100.0\pm 0.0\%$	$66.9\pm 2.9\%$	$100.0\pm 0.0\%$
3	$83.0\pm 6.3\%$	$100.0\pm 0.0\%$	$66.8\pm 4.9\%$	$100.0\pm 0.0\%$
4	$80.4\pm 5.6\%$	$100.0\pm 0.0\%$	$63.6\pm 4.2\%$	$100.0\pm 0.0\%$
5	$85.2\pm 5.9\%$	$100.0\pm 0.0\%$	$64.7\pm 7.1\%$	$100.0\pm 0.0\%$

Table 13. Table 13 : Persistence Images Results with pixel size = 0.1.

	2 inch		2.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$82.5\pm 1.4\%$	$82.7\pm 0.6\%$	$79.3\pm 4.2\%$	$81.4\pm 2.9\%$
LR	$80.2\pm 2.0\%$	$82.4\pm 0.9\%$	$77.0\pm 7.0\%$	$82.4\pm 4.5\%$
RF	$96.4\pm 1.1\%$	$100.0\pm 0.0\%$	$85.8\pm 5.3\%$	$100.0\pm 0.0\%$
GB	$95.9\pm 1.5\%$	$100.0\pm 0.0\%$	$84.4\pm 4.4\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$82.6\pm 5.5\%$	$83.2\pm 2.5\%$	$66.9\pm 5.8\%$	$63.7\pm 2.9\%$
LR	$82.6\pm 5.1\%$	$85.7\pm 2.0\%$	$62.7\pm 5.9\%$	$65.5\pm 3.3\%$
RF	$93.0\pm 2.9\%$	$100.0\pm 0.0\%$	$72.5\pm 6.4\%$	$100.0\pm 0.0\%$
GB	$91.3\pm 4.3\%$	$100.0\pm 0.0\%$	$68.5\pm 2.8\%$	$100.0\pm 0.0\%$

Table 14. Table 14 : Persistence Images Results with pixel size = 0.05.

	2 inch		2.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$82.0\pm 2.5\%$	$83.4\pm 1.2\%$	$81.2\pm 5.1\%$	$80.2\pm 2.8\%$
LR	$80.3\pm 2.3\%$	$79.5\pm 1.0\%$	$68.6\pm 7.3\%$	$75.2\pm 3.4\%$
RF	$96.0\pm 1.4\%$	$100.0\pm 0.0\%$	$85.8\pm 4.1\%$	$100.0\pm 0.0\%$
GB	$95.5\pm 1.6\%$	$100.0\pm 0.0\%$	$85.6\pm 3.4\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$85.2\pm 4.0\%$	$81.8\pm 2.5\%$	$63.4\pm 5.2\%$	$65.6\pm 2.4\%$
LR	$80.9\pm 5.9\%$	$85.9\pm 3.0\%$	$63.6\pm 4.5\%$	$64.8\pm 2.1\%$
RF	$90.9\pm 3.0\%$	$100.0\pm 0.0\%$	$70.7\pm 3.1\%$	$100.0\pm 0.0\%$
GB	$90.4\pm 4.7\%$	$100.0\pm 0.0\%$	$70.0\pm 5.4\%$	$100.0\pm 0.0\%$

Table 15. Table 15 : Template Function Results for H 0 subscript 𝐻 0 H_{0} diagrams.

	2 inch		2.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$82.4\pm 1.8\%$	$80.6\pm 1.2\%$	$82.1\pm 6.5\%$	$84.8\pm 2.2\%$
LR	$81.7\pm 2.1\%$	$83.7\pm 1.2\%$	$78.1\pm 7.4\%$	$85.9\pm 2.3\%$
RF	$87.6\pm 4.0\%$	$100.0\pm 0.0\%$	$86.5\pm 2.5\%$	$100.0\pm 0.0\%$
GB	$86.7\pm 2.5\%$	$99.5\pm 0.5\%$	$84.9\pm 5.7\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$83.5\pm 6.7\%$	$81.4\pm 3.5\%$	$62.4\pm 5.5\%$	$66.0\pm 2.8\%$
LR	$71.7\pm 5.2\%$	$90.0\pm 1.8\%$	$69.3\pm 4.6\%$	$70.9\pm 2.2\%$
RF	$78.7\pm 5.3\%$	$100.0\pm 0.0\%$	$72.4\pm 3.9\%$	$100.0\pm 0.0\%$
GB	$77.8\pm 10.0\%$	$100.0\pm 0.0\%$	$72.2\pm 5.5\%$	$100.0\pm 0.0\%$

Table 16. Table 16 : Template Function Results for H 1 subscript 𝐻 1 H_{1} diagrams.

	2 inch		2.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$84.0\pm 1.9\%$	$85.9\pm 1.5\%$	$84.0\pm 7.3\%$	$86.4\pm 2.9\%$
LR	$87.5\pm 1.7\%$	$90.1\pm 0.7\%$	$75.8\pm 6.9\%$	$92.9\pm 2.3\%$
RF	$90.5\pm 1.6\%$	$100.0\pm 0.0\%$	$89.3\pm 4.6\%$	$100.0\pm 0.0\%$
GB	$91.5\pm 1.2\%$	$100.0\pm 0.0\%$	$85.3\pm 4.7\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$83.0\pm 4.9\%$	$82.0\pm 2.4\%$	$64.7\pm 3.6\%$	$64.8\pm 1.8\%$
LR	$78.7\pm 10.0\%$	$100.0\pm 0.0\%$	$64.1\pm 4.4\%$	$81.8\pm 2.8\%$
RF	$81.3\pm 8.5\%$	$100.0\pm 0.0\%$	$61.4\pm 5.9\%$	$100.0\pm 0.0\%$
GB	$83.9\pm 5.5\%$	$100.0\pm 0.0\%$	$65.1\pm 6.1\%$	$100.0\pm 0.0\%$

Table 17. Table 17 : Carlsson Coordinates Results.

	2 inch		2.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$87.8\pm 2.0\%$	$87.1\pm 1.2\%$	$72.1\pm 7.1\%$	$79.7\pm 3.9\%$
LR	$93.1\pm 1.8\%$	$92.9\pm 0.9\%$	$86.3\pm 6.2\%$	$100.0\pm 0.0\%$
RF	$93.0\pm 2.0\%$	$100.0\pm 0.0\%$	$69.5\pm 4.8\%$	$70.9\pm 2.6\%$
GB	$93.6\pm 1.7\%$	$100.0\pm 0.0\%$	$84.7\pm 5.0\%$	$100.0\pm 0.0\%$
	3.5 inch		4.5 inch
Classifier	Test Set	Training Set	Test Set	Training Set
SVM	$95.2\pm 4.5\%$	$95.2\pm 1.9\%$	$68.1\pm 6.5\%$	$68.6\pm 2.4\%$
LR	$90.9\pm 9.2\%$	$93.6\pm 1.4\%$	$70.8\pm 4.0\%$	$72.6\pm 3.0\%$
RF	$95.7\pm 4.3\%$	$100.0\pm 0.0\%$	$72.2\pm 4.4\%$	$100.0\pm 0.0\%$
GB	$92.2\pm 4.3\%$	$100.0\pm 0.0\%$	$71.9\pm 3.3\%$	$100.0\pm 0.0\%$

Table 18. Table 18 : Kernel Method Results with LibSVM package

	2 inch	2.5 inch	3.5 inch	4.5 inch
Kernel Scale ($\sigma$)	Test Set Score and Deviation
0.2	*	*	$30.4\pm 6.2\%$	*
0.25	$74.5\%$**	$58.9\pm 28.5\%$	$87.0\pm 3.6\%$	$59.3\pm 9.6\%$

Table 19. Table 19 : Path signature method results obtained with SVM classifier. Landscape number represent which landscapes was used to compute signatures.

	2 inch		2.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$83.0\pm 2.9\%$	$84.8\pm 0.9\%$	$82.7\pm 5.3\%$	$82.9\pm 2.0\%$
2	$79.2\pm 3.1\%$	$80.2\pm 1.0\%$	$84.2\pm 3.5\%$	$80.4\pm 1.4\%$
3	$78.6\pm 1.8\%$	$79.4\pm 0.7\%$	$79.1\pm 2.9\%$	$80.8\pm 1.1\%$
4	$79.3\pm 2.5\%$	$79.5\pm 0.8\%$	$80.6\pm 5.6\%$	$78.6\pm 2.3\%$
5	$0.0\pm 0.0\%$	$0.0\pm 0.0\%$	$80.9\pm 5.1\%$	$78.8\pm 1.2\%$
	3.5 inch		4.5 inch
Landscape Number	Test Set	Training Set	Test Set	Training Set
1	$81.2\pm 6.3\%$	$82.6\pm 2.2\%$	$70.0\pm 6.4\%$	$71.4\pm 2.2\%$
2	$82.9\pm 7.2\%$	$81.8\pm 2.4\%$	$64.1\pm 5.4\%$	$67.0\pm 2.4\%$
3	$81.2\pm 10.5\%$	$82.2\pm 3.8\%$	$67.7\pm 6.6\%$	$65.8\pm 2.3\%$
4	$85.9\pm 4.7\%$	$80.8\pm 1.6\%$	$64.5\pm 4.2\%$	$65.4\pm 1.2\%$
5	$82.4\pm 9.1\%$	$82.0\pm 3.1\%$	$64.8\pm 6.9\%$	$65.4\pm 2.1\%$

Equations

$$R_{ϵ}(K,d)=\{\sigma \subset K|\underset{(x,y)\in \sigma }{\max }d(x,y)<ϵ\},$$

$$R_1\subseteq R_2\subseteq \dots \subseteq R_m$$

$${\Phi }_{\varphi ,\beta }$$

$$=(\beta ({\mathbf{x}}_t),\beta ({\mathbf{x}}_{t+\tau },\beta ({\mathbf{x}}_{t+2\tau },\dots ,\beta ({\mathbf{x}}_{t+(d-1)\tau }))),$$

$$p(t)=\sum _{i=0}^3(1-t{)}^{3-i}{t}^i{p}_{i}0<t<1,$$

$$\text{centering}L(p_0,p_1,p_2,p_3)=\sum _{i=1}^l||p(t_i)-x_i|{|}^2,\text{@add@centering}$$

$$\displaystyle\sum_{k=1}^{l}2\Bigg{(}\sum_{i=0}^{3}\binom{3}{i}(1-t_{k})^{3-i}t_{k}^{i}p_{i}-x_{k}\Bigg{)}\Bigg{(}\frac{\partial p(t_{k})}{\partial p_{i}}\Bigg{)}=0$$

$$\displaystyle\leavevmode\resizebox{208.13574pt}{}{$\frac{\partial p(t_{k})}{\partial p_{i}}=(1-t_{k})^{3}\hat{e}_{1}+3(1-t_{k})^{2}t_{k}\hat{e}_{2}+3(1-t_{k})t_{k}^{2}\hat{e}_{3}+t_{k}^{3}\hat{e}_{4}=0$},$$

$$\Bigg{(}\sum_{k=1}^{l}A_{k}\Bigg{)}p=\sum_{k=1}^{l}b_{k}\quad\text{or}\quad Ap=b.$$

$$A_k=\left[\begin{array}{cccc}(1-t_k{)}^6& 3(1-t_k{)}^5{t}_k& 3(1-t_k{)}^4{t}_k^2& (1-t_k{)}^3{t}_k^3\\ 3(1-t_k{)}^5{t}_k& 9(1-t_k{)}^4{t}_k^2& 9(1-t_k{)}^3{t}_k^3& 3(1-t_k{)}^2{t}_k^4\\ 3(1-t_k{)}^4{t}_k^2& 9(1-t_k{)}^3{t}_k^3& 9(1-t_k{)}^2{t}_k^4& 3(1-t_k)t_k^5\\ (1-t_k{)}^3{t}_k^3& 3(1-t_k{)}^2{t}_k^4& 3(1-t_k)t_k^5& t_k^6\end{array}\right]$$

$$p=\left[\begin{array}{cccc}{p}_0^1& p_0^2& \dots & p_0^d\\ p_1^1& p_1^2& \dots & p_1^d\\ p_2^1& p_2^2& \dots & p_2^d\\ p_3^1& p_3^2& \dots & p_3^d\end{array}\right]b_k=\left[\begin{array}{c}(1-t_k{)}^3\\ 3(1-t_k{)}^2{t}_k\\ 3(1-t_k)t_k^2\\ t_k^3\end{array}\right]{\left[\begin{array}{c}{x}_k^1\\ x_k^2\\ ⋮\\ x_k^d\end{array}\right]}^T,$$

$$d(\overrightarrow{l_0{l}_1},\overrightarrow{m_0{m}_1})=\underset{l\in \overrightarrow{l_0{l}_1},m\in \overrightarrow{m_0{m}_1}}{\min }d(l,m),$$

$$f(l,m)=d(l(s),m(t){)}^2=||l(s)-m(t)|{|}^2,$$

$$W_{\infty }(X_1,X_2)=\underset{\eta :X_1\to X_2}{\inf }\underset{x_1\in X_1}{\sup }||x_1-\eta (x_1)|{|}_{\infty },$$

$$g_{(b,d)}(x)=\{\begin{array}{ll}0& \text{if x∈(b,d)}\\ x-b& \text{if x∈(b,2b+d]}\\ -x+d& \text{if x∈(2b+d,d)}\end{array}$$

$$T(b_i,d_i)=(b_i,d_i-b_i)=(b_i,p_i),$$

$$D_k(x,y)=\frac{1}{2\pi {\sigma }^2}{e}^{-[(x-b_k{)}^2+(y-p_k{)}^2]/2{\sigma }^2}.$$

$$W(k)=W(b_k,p_k)=\{\begin{array}{ll}0& \text{ if }{p}_k\le 0;\\ \frac{p_k}{b}& \text{ if }0<p_k<b;\\ 1& \text{ if }{p}_k\ge b.\end{array}$$

$$S(x,y)=\sum _{k\in T(X)}W(k)D_k(x,y).$$

$$I_{i,j}(S)=\iint Sdxdy,$$

$$\left[\begin{array}{cccc}{I}_{1,1}& I_{1,2}& \dots & I_{1,100}\\ ⋮& \ddots & & \\ I_{100,1}& & & I_{100,100}\end{array}\right],$$

$$\begin{array}[]{rl}f_{1}(X)&=\sum{b_{i}(d_{i}-b_{i})},\\ f_{2}(X)&=\sum{(d_{\rm max}-d_{i})(d_{i}-b_{i})},\\ f_{3}(X)&=\sum{b_{i}^{2}(d_{i}-b_{i})^{4}},\\ f_{4}(X)&=\sum{(d_{\rm max}-d_{i})^{2}(d_{i}-b_{i})^{4}},\\ f_{5}(X)&=\max\{(d_{i}-b_{i})\}\end{array}$$

$${\kappa }_{\sigma }(X,Y)=\frac{1}{8\pi \sigma }\sum _{z_1ϵX,z_2ϵY}\exp \left(-\frac{||z_1-z_2|{|}^2}{8\sigma }\right)-\exp \left(-\frac{||z_1-{\hat{z}}_2|{|}^2}{8\sigma }\right),$$

$${\kappa }_{\sigma }=\left[\begin{array}{cccc}{\kappa }_{\sigma }(X_1,X_1)& {\kappa }_{\sigma }(X_1,X_2)& \dots & {\kappa }_{\sigma }(X_1,X_N)\\ ⋮& \ddots & & \\ {\kappa }_{\sigma }(X_N,X_1)& & & {\kappa }_{\sigma }(X_N,X_N)\end{array}\right].$$

$$\gamma (t)={\gamma }_t=[{\gamma }_t^1,{\gamma }_t^2,\dots ,{\gamma }_t^d],$$

$$S(\gamma {)}_{a,t}^i={\int }_a^{t}d{\gamma }_s^i={\gamma }_t^i-{\gamma }_a^i,(a<s<t);$$

$$S(\gamma {)}_{a,t}^{i,j}={\int }_a^{t}S(\gamma {)}_{a,s}^{i}d{\gamma }_s^k={\int }_a^t{\int }_a^{s}d{\gamma }_r^{i}d{\gamma }_s^j,(a<r<s<t);$$

$$S(\gamma {)}_{a,t}^{i,j,\dots ,k}={\int }_a^{t}S(\gamma {)}_{a,s}^{i,j,\dots ,k-1}d{\gamma }_s^k={\int }_a^t\dots {\int }_a^{t_2}d{\gamma }_{t_1}^i\dots d{\gamma }_{t_k}^k,(a<t_1<t_2<\dots <t).$$

$${\gamma }_t({\lambda }_k(t))=[t,{\lambda }_k(t)].$$

$${\lambda }_k(t)=\{\begin{array}{ll}\frac{y_{i+1}-y_i}{x_{i+1}-x_i}t+\frac{y_{i+1}(x_{i+1}-x_i)+x_{i+1}(y_i-y_{i+1})}{x_{i+1}-x_i},& \text{ for }i\in \{1,2,\dots ,n\}\text{ and }t\in [x_1,x_n],\\ 0,& \text{ otherwise }.\end{array}$$

$$P=[P_t^1,P_t^2]=\left[t,\frac{y_{i+1}-y_i}{x_{i+1}-x_i}t+\frac{y_{i+1}(x_{i+1}-x_i)+x_{i+1}(y_i-y_{i+1})}{x_{i+1}-x_i}\right],$$

$$dP=[d{P}_t^1,d{P}_t^2]=[dt,\frac{y_{i+1}-y_i}{x_{i+1}-x_i}dt].$$

$$S_1={\int }_{x_1}^{x_n}d{P}_t^1={\int }_{x_1}^{x_n}dt=x_n-x_1$$

Full text

Keywords: Chatter detection, machine learning, topological data analysis, persistence homology, featurization, Bézier curves

1 Introduction

Chatter is characterized by excessive vibrations that can negatively affect the surface finish or shorten the lifetime of cutting tools during machining operations such as turning and milling. Consequently, identification and mitigation of chatter has become a prominent research topic in recent decades. Some of the challenges associated with chatter identification is that it is a complex phenomenon that depends on several factors including the dynamic properties of the tool and the workpiece. Therefore, as these properties vary during the cutting process, the results of predictive models become invalid, thus necessitating a data-based approach for more reliable chatter detection. Motivated by this goal, many studies in the literature have focused on extracting chatter features from signals obtained using sensors mounted on the cutting center. Most of these studies are based on analyzing the spectrum of force or acceleration signals often in combination with machine learning techniques.

For example, Thaler et al. utilized short time Fourier transform with quadratic discriminant analysis to detect chatter in a band sawing process [1]. In addition, Yesilli et al. extracted features from experimental turning signals using Fast Fourier Transform (FFT), Power Spectral Density (PSD) and the Auto-correlation Function (ACF) [2]. However, the two most common methods for analyzing cutting signals are the Wavelet Packet Transform (WPT) and the Empirical Mode Decomposition (EMD). Chen and Zheng used WPT to isolate the part of the signal carrying chatter markers, and they utilized it to compute both time and frequency domain features [3] that can be used as feature vector for machine learning. Similarly, Ji et al. used EMD to decompose the acceleration signal from the tool vibration into Intrinsic Mode Functions (IMFs) [4]. They then identified the IMFs that carry chatter information called the informative IMFs and they used them to define feature vectors for chatter identification. Li et al. utilizes coarse-grained entropy rate to detect chatter in turning [5]. In addition, Liu et al. applied EMD to servo motor current signal to compute the energy and the kurtosis of the informative IMFs as features for chatter [6]. They then used Support Vector Machine (SVM) algorithm to train a classifier for chatter detection. SVM is the most widely used machine learning algorithm for chatter classification from sensory signals. Some less commonly used methods include logistic regression [7], deep belief network [8], back propagation neural network [9], convolutional neural networks [10], and Hidden Markov Model [11, 12]. Zuo et al. proposed an approach which is based on spiking neural networks for fault detection in bearings [13]. Although WPT and EEMD are widely used for detecting chatter, these methods have some limitations that preclude them from being adopted as general chatter detection tools. To elaborate, Yesilli et al. recently showed that identifying appropriate feature vectors using these two methods is signal-dependent and it requires skilled operators [14]. Further they also showed that the accuracy of the trained classifier is highly sensitive to changes in the dynamic properties of the tool-workpiece system [14]. In addition, to these the features obtained from signal decomposition methods does not provide extensive information about the data set. Therefore, we propose a novel approach based on Topological Data Analysis (TDA) which removes the requirement of preprocessing. This increases the level of automation of feature extraction. In addition, we can have more information about the structure of the data using the features obtained with TDA.

Topological Data Analysis (TDA) [15, 16, 17, 18] is a promising tool for generating feature vectors for chatter detection comes from a new field with many mature computational tools. TDA, and more specifically persistent homology, provides a quantifiable way for describing the topological features in a signal [19]. Specifically, by embedding the sensory signal into a point cloud, it is then possible to use persistent homology to produce a multiscale summary of the topological features of the signal thus enabling the analysis of the underlying dynamical system. The homology classes that correspond to the embedded signal are often reported using a planar diagram that shows how long each topological feature persisted. The application of TDA tools to machining dynamics has only been recently explored [20, 21, 22]. Specifically, Ref. [21] and [23] show that maximum persistence—a single number from the persistence diagram—can be used to ascertain the stability of simulated data from a stochastic turning model. Khasawneh et al. incorporated more information from the persistence diagram by extracting $5$ features including Carlsson coordinates ([24]) and the maximum persistence [22], see Section 1.1 for more details on featurizing persistence diagrams. The resulting feature vector in combination with SVM was used to train a chatter classifier, and it was applied to simulated deterministic and stochastic turning data with success rates as high as $97\%$ in the deterministic case. In addition, Yesilli et al. utilized Carlsoon Coordinates and Template Functions [25] to diagnose chatter in milling simulations and show that these two featurization methods are noise robust [26].

However, despite the active work in the literature on featurizing persistence diagrams, all prior studies on chatter detection with TDA have utilized only a small fraction of the persistence diagram for constructing a feature vector. Further, these publications only studied simulated signals and no sensory signals from actual cutting tests have been tested. Therefore, this work aims to collect and summarize state-of-the-art featurization tools for persistence diagrams, and apply them for the first time for chatter classification using actual experimental signals obtained from an accelerometer mounted on the cutting tool during a turning process. The methods that we investigate for featurizing the resulting persistence diagrams and classifying chatter time series include persistence landscapes [27], Carlsson coordinates [24], persistence images [28], an example kernel method [29], and path signatures of persistence landscapes [30]. Moreover, we provide the run time for each featurization method including the runtime for persistence diagram computation which constitutes the majority of the total computation time. To reduce the runtime for persistence diagram computation, we utilize Bézier curve approximation method [31], greedy permutation [32] and parallel computing.

The paper is organized as follows. In Section 1.1 we provide a brief literature survey of featurization methods for persistence diagrams. In Section 2 we describe the experimental apparatus as well as how the data was processed and tagged. Section 3 provides the needed background material on the tools we use from TDA. The method we use for featurizing the persistence diagrams and how we compute persistence diagrams are explained in Section 4. Section 5 presents the classification results and runtime comparisons, while the concluding remarks can be found in Section 6.

1.1 An overview for featurization of persistence diagrams

Persistence diagrams can be compared to each other using metrics such as the Wasserstein and Bottleneck distances. This allows creating similarity or distance matrices between different persistence diagrams that can be used for machine learning. For example, in previous studies on machine learning applications with TDA, these two distances were utilized as similarity measures between persistence diagrams for classification algorithms [33]. However, working directly with persistence diagrams is difficult in part due to the non-uniqueness of geodesics which leads to non-uniqueness of the Fréchet mean for a collection of diagrams [34, 35, 36]. This necessitates using methods for extracting feature vectors that live in a Euclidean space from the persistence diagrams, which enables using traditional machine learning tools with topological features. The following paragraph provides a list of some of these methods.

Adcock et al. describe a subring of polynomial functions on persistence diagrams and a convenient system of free generators [24]. However, their construction is not continuous with respect to the bottleneck distance, and it is not applicable to the case of infinite persistence diagrams which are relevant when working with fractal structure. Other tools for extracting feature vectors from persistence diagrams are based on functional summaries [37], which turn persistence diagrams into functions. Perhaps the most common functional summary of persistence diagrams is persistence landscapes [27] which represent persistence diagrams in the form of piecewise continuous functions. Although averages of persistence diagrams are not well-defined, persistence landscapes have well-established statistical properties [38] and their averages and pairwise distances are defined which enables combining them with traditional machine learning tools. Another functional summary of persistence diagrams is persistence images  [39, 28], which are closely related to Chen et al. persistence intensity functions [40]. Persistence images are also related to prior works on the size function which predates persistence [41, 42]. These images are obtained by placing Gaussians at each point in the diagram, and then utilizing a histogram of the image for featurization and machine learning.

Another featurization method is based on path signatures of persistence landscapes [30]). Specifically, signatures for $n$-dimensional paths of bounded variation can be computed and used as features for machine learning [43]. Therefore, by composing persistence landscapes, which can be viewed as one type of persistence path embedding, with path signatures, it is possible to extract mesh-free feature vectors from the underlying persistence diagram. However, path signatures do not come with stability guarantees for all input paths.

Persistence diagrams can also be featurized using a general class of functions called template functions [25]. These functions are only required to be continuous and compactly supported. The specific realization of these functions that we work with in this paper are Chebyshev interpolating polynomials [44].

Kernel-based methods have also been used for featurizing persistence diagrams [29, 45, 46, 47, 48, 49, 50]). However, instead of testing all of the available options, we choose the persistence scale space kernel given by [29] as a representative of this group of methods. This kernel is defined by treating the persistence diagram as a sum of Dirac deltas at each point in the persistence diagram, and using this as the initial condition for a heat diffusion problem. This allows obtaining a closed form solution for this kernel.

2 Experimental Setup and Signal Processing

This section describes the details of the cutting tests 2.1, the data collection, signal tagging 2.3, and time series conditioning 2.2.

2.1 Cutting tests

Figure 1 shows the experimental setup which features a Clasuing-Gamet $33$ cm ($13$ inch) engine lathe instrumented with three accelerometers. Two PCB 352B10 uniaxial accelerometers were used to measure the $x$ and $y$ accelerations of the boring rod, see Fig. 1. These uniaxial accelerometers have frequency range between 2 and $10000$ Hz and sensivity of 1.02 mV/(m/$s^{2}$). They were mounted $3.81$ cm ($1.5$ inch) away from the cutting tool to keep them safe from the cutting debris. Further, a PCB 356B11 triaxial accelerometer was attached to the tool holder with superglue. Triaxial accelerometer has frequency range between $2$ to $10000$ Hz for $y$ and $z$ axes and it has range of $2$ to $7000$ Hz for $x$ axis. It has also sensitivity of 1.02 mV/(m/$s^{2}$). The acceleration signals were collected with a sampling rate of $160$ kHz using Matlab’s data acquisition toolbox and a NI USB-6356 data acquisition box, which has 8 analog inputs with 1.25MS/s/ch sampling rate and 16-bit resolution, see Section 2.2 for more information on the applied signal processing.

The used cutting tool is a $0.04$ cm ($0.015$ inch) radius Titanium nitride coated insert which was attached to an S10R-SCLCR3S boring bar, where the latter is part of the Grizzly T10439 carbide insert boring bar set. This boring bar set uses the $80^{\circ}$ diamond cutting inserts and has negative $7^{\circ}$ end and side cutting angle. 6061 aluminum cylindrical workpiece with Brinell hardness of 95 is used in the cutting tests.

The cutting tests were performed by varying the depth of cut and the cutting speed while holding the stiffness of the boring rod constant. Some of the cutting depths are 0.0127 cm, 0.0508 cm and 0.0635 cm (0.005 inch, 0.02 inch and 0.025 inch). Each vibration signal is labeled with its cutting speed and depth of cut and available in the Mendeley repository [51]. The feed rate for all cutting tests was kept constant and it is 0.00508 cm/rev (0.002 inch/rev). Data was collected from four sets of cutting configurations where each configuration corresponds to a different boring rod stiffness/eigenfrequency. Experimental data is available online in raw and processed format [51]. The stiffness of the rod was varied by changing the overhang length of the tool, which is the distance between the heel of the boring rod and the back surface of the tool holder. Four different overhang lengths were considered: $5.08$ cm ($2$ inch), $6.35$ cm ($2.5$ inch), $8.89$ cm ($3.5$ inch) and $11.43$ cm ($4.5$ inch). Increasing the overhang length leads to a stiffer rod and to higher chatter frequencies since the latter are close to the eigenfrequency of the rod. Similarly, smaller overhang length leads to a more compliant system, and to lower chatter frequencies.

2.2 Signal Processing

The experimental data was sampled at the rate of $160$ kHz, which is an oversampling by a factor of $16$ of the maximum accelerometer measuring frequency of $10$kHz. Since no in-line analog filter was used, this oversampling is necessary to avoid aliasing effects by utilizing a digital filter before downsampling the signal to $10$kHz. The filter we used is Butterworth low pass filter with order $100$ from Matlab’s signal processing toolbox. The resulting filtered and downsampled signals were the ones used in tagging the data as described in Section 2.3.

2.3 Data Tagging

Upon examining the acceleration signals, we found that the $x$-axis signals from the tri-axial accelerometer (see Fig. 1) contained the best Signal-to-Noise Ratio (SNR). The other acceleration signals contained redundant information albeit with lower SNR; therefore, only the $x$-axis signal from the tri-axial accelerometer was used for tagging the data. Both time and frequency domain information were simultaneously considered when tagging the data. Specifically, in the time domain, we separated the different parts of the raw time series according to amplitude, i.e., portions with large amplitude, low amplitude, and with sudden increase in amplitude. In the frequency domain, frequencies under $5$ kHz in the downsampled signal were taken into account.

We established specific criteria to label the raw data which we describe below. We also provide the number of tagged time series for each overhang case in Table 1. Figure 2a shows that during the same cutting test it was possible for different parts of the time series to be tagged differently. A signal was tagged chatter-free (stable) if it has a low amplitude in both time and frequency domains. These signals have the highest peaks at the spindle rotation frequencies in the Fourier spectrum [52], as shown in Fig. 2b. On the other hand, we characterized intermediate chatter signals as having low amplitude in the time domain but high amplitude in the frequency domain, see Fig. 2c. In contrast, chatter signals have large amplitude in both time and frequency domains, as shown in Fig. 2a and the spectrum superimposed on Figs. 2b–2d. For intermediate chatter and chatter signals, the highest peaks are around the chatter frequencies in the frequency domain. Any case that did not fit with any of the above categories was labeled as unknown. Figure 2d shows that these cases are generally composed of time series that have sudden increase and decrease in amplitude in the time domain combined with low amplitude in the frequency domain.

The data labeling was verified by spot-checking the tagged data against pictures taken during the cutting tests. Figure 3 shows one such comparison for a sample tagged signal. The figure shows that the surface finish characteristics agree with the assigned tagging, thus confirming the validity of the labeling.

3 Topological Data Analysis

Topological Data Analysis (TDA) extracts information by investigating the shape of the data. We study persistent homology which is a powerful tool of Topological Data Analysis (TDA) to extract features from the persistence diagrams and use them in supervised machine learning algorithms. Experimental data is embedded using Takens embedding theorem [53] and 1-D persistent homology is investigated for feature matrix generation. In this section, we briefly explain persistent homology, and one can refer to [15, 16, 17, 18, 54, 55] for detailed information about TDA and persistent homology.

Persistent homology provides a compact tool for studying the topology of data embedded in a Euclidean space which is often called a point cloud. The resulting shape information is represented by a two dimensional plot called the persistence diagram. We can obtain a persistence diagram for different shape characteristics of interest. For instance, if we are interested in the connectivity of the points in the point cloud, then we can represent that information in a [math]-dimensional ([math]-D) persistence diagram. Alternatively, if we are interested in loops in the point cloud, then we compute the $1$-D persistence diagram. For voids, we go to the $2$-D persistence diagram and so on. In this study we only extracted features from the [math]-D and $1$-D persistence, so in the following we will introduce the basic idea for obtaining the $1$-D persistence, i.e., for representing loops that emerge and disappear as the point cloud is thickened. The process for obtaining the [math]-D persistence is similar, only instead of considering loops, we would need to track the connectivity of the points as the point cloud is uniformly thickened.

Consider the point cloud shown in Fig. 5a. Then, we start to thicken the point cloud, i.e., expand disks with radius $\epsilon$ around each data point. As $\epsilon$ is increased, disks can start to intersect. The intersection of two disks forms an edge as shown by the two edges in Fig. 5b. Increasing $\epsilon$ further can lead to three disks intersecting thus forming a triangle which we fill in as shown by the nine triangles in Fig. 5c. At some values of $\epsilon$, some disk intersections will lead to cycles. The time (here the $\epsilon$ value) at which a cycle appears is called the birth time of the cycle. Figure 5d shows three example cycles numbered 1,2 and 3 with birth time $b_{1}=b_{2}=b_{3}$.

As we continue thickening the disks in the point cloud, more disks will intersect leading to more triangles filling in, and at some point, some cycles may fill in. The time at which a cycle disappears is called its death time. For example, Figs. 5e–g show the death times of cycles $1$–$3$, respectively. The information about the birth and death of cycles is succinctly summarized in a persistence diagram. In this diagram, each point corresponds to the paired birth and death times of a cycle. For example, Fig. 5h shows the tuples $(b_{1},d_{1})$, $(b_{2},d_{2})$, and $(b_{3},d_{3})$ corresponding to the birth and death times of the cycles $1$–$3$ shown in Fig. 5d. The cycle that persists the longest is characterized by the highest point above the diagonal, and is called the maximum persistence. In this example, cycle $3$ leads to the maximum persistence in the persistence diagram.

3.1 Simplicial complexes

Let $\{u_{0},\ldots,u_{k}\}\in\mathbb{R}^{d}$ be a set of data points, and the vectors defined between these data points ($u_{1}$-$u_{0}$,$u_{2}$-$u_{0}$,…,$u_{k}$-$u_{0}$) are linearly independent. A geometric *k-*simplex, $\sigma$ is a set of all points in $\mathbb{R}^{d}$ such that $\sum_{j=0}^{k}\lambda_{j}u_{j}$ where $\sum_{j=0}^{n}\lambda_{j}=1$ and $\lambda_{j}\geq 0$ for all j. Figure 14 provides illustrations for $0,1,$ and $2-$simplex. Each data point on a point cloud is represented as 0-dimensional simplex and they are called vertices. When two vertices are connected, an edge is formed and it is 1-dimensional simplex. Connection of three vertices will form $2-$simplex which is a triangle. Simplicies spanned by any subset of $u_{0},\ldots,u_{k}$ are the faces of $\sigma$. In general, $n-$simplex contains $n+1$ vertices, and the set of these simplices, are called geometric simplicial complexes, $K$, if the following two conditions are satisfied [54]:

1. If $\sigma\in K$, then faces of $\sigma$ are also in $K$,
2. If two $n-$simplex, $\sigma_{1}$ and $\sigma_{2}$ are in $K$, then the intersection of them is either common face or empty. The dimension of the simplicial complex, $K$ is equal to the largest dimension of its simplices.

3.2 Persistent Homology

The simplicial complex $K$ is used to compute homology $H_{n}(K)$ in different dimension to identify the shape of the data. 0 dimensional homology, $H_{0}(K)$ represents connected components and one dimensional homology, $H_{1}(K)$ represents loops, while two dimensional homology $H_{2}(K)$ represents voids. In persistent homology, the simplicial complex $K$ is not fixed and it is varying over time.

We start expanding disks centered at data points of point cloud and let $\epsilon$ be the radius of these disks. As $\epsilon$ increases, n dimensional simplicies are formed as shown in Fig. 5. Intersection of two disks forms an edge (1-simplex), while a triangle (2-simplex) is formed when three disk intersect with each other (see Fig. 14). Each $\epsilon$ will result in different simplicial complexes, and they can be approximated using filtration functions. In this study, we used a Python package that employs Rips complex whose definition is given as

[TABLE]

where $d$ is the distance between vertices of simplicial complex $\sigma$. Let {$\epsilon_{1}<\epsilon_{2}<\ldots<\epsilon_{m}$ } be the set of varying radius of the disk. These radii form Rips complexes such that

[TABLE]

where $R_{j}=R_{\epsilon_{j}}(K,d)$. Then, we can choose a specific dimension $n$ to identify the shape of the data along the simplicial complexes $R_{j}$. For instance, if a loop is seen first in $R_{i}$, this is called birth time ($b=\epsilon_{i}$). When it disappears in $R_{j}$, this is denoted as death time ($d=\epsilon_{j}$, where $d>b$). That allows us to generate persistence diagrams, $D$ for a given point cloud and selected persistent homology dimension. Figure 5h provides an example of a persistence diagram. The horizontal axis represents the birth time, while the vertical axis is for the death time. All points on the diagram is on the upper diagonal and the points with larger lifetime $(d-b)$ is far away from the diagonal point, and these points includes major information about the shape and the structure of the data. In this study, we mostly focused on the one dimensional persistent homology to extract features.

4 Method

The method that we develop for chatter detection using topological features can be summarized using Fig. 6. The parts of the pipeline related to data collection, processing, and labeling were described in Sections 2.1–2.3, respectively. In this section, we describe the rest of the steps shown in Fig. 6.

Recall that the cutting tests are composed of four different overhang configurations, and that each configuration includes a different number of time series that correspond to different labels, rotational speeds, and depths of cut. Therefore, the time series are grouped with respect to these three parameters and they are normalized to have zero mean and unit variance. This normalization reduces the effect of large feature values on smaller ones ([56]).

The next step is to split long time series into smaller pieces to reduce the computation time needed for finding the delay reconstruction parameters (see Section 4.1.1), and for obtaining the persistence diagrams (see Section 3.2). Upon finding the appropriate embedding parameters, the data is embedded using delay reconstruction, also known as Takens’ delay embedding, see Section 4.1.1. The resulting point cloud is then used to compute the corresponding persistence diagrams using two different approaches: 1)traditional way where persistence diagrams are computed with Ripser package for Python, 2) approximating to point cloud with Bezier curves and computing persistence diagrams using line segments generated with Bézier curves ([31]). Second method is a different approach introduced in [31], and it uses Bézier curves to approximate to reduce time to compute persistence diagrams. Section 4.1 explains both approaches for persistence diagram computation and sections 4.2.1–4.2.5 describe different methods in the literature for featurizing the resulting diagrams to obtain a feature vector in Euclidean space that can be used with existing machine learning tools such as SVM. We mainly featurize and use the resulting $1$-dimensional persistent homology $H_{1}$ except when we use the template function approach where we also utilize the [math]-dimensional persistence $H_{0}$ ([25]).

4.1 Persistence Diagram Computation

In this section, we briefly explain how to embed a time series to higher dimension and to compute persistence diagrams using Bézier curve approximation method. The steps for persistence diagram computation are summarized as shown in Fig. 7.

4.1.1 Delay Reconstruction

Takens’ theorem lays a theoretical framework for studying deterministic dynamical systems ([53]). It states that in general we can obtain an embedding of the attractor of a deterministic dynamical system from a one-dimensional recording of a corresponding trajectory. This embedding is a smooth map $\Psi:M\to N$ between the manifolds $M$ and $N$ that diffeomorphically maps $M$ to $N$.

Specifically, assume we have an observation function $\beta(\mathbf{x}):M\to\mathbb{R}$, where for any time $t\in\mathbb{R}$ the point $\mathbf{x}$ lies on an $m$-dimensional manifold $M\subseteq\mathbb{R}^{d}$. While in practice we do not have the the flow of the system for an time $t\in\mathbb{R}$ given by $\phi^{t}(\mathbf{x}):M\times\mathbb{R}\to M$, the observation function implicitly captures the time evolution information according to $\beta(\phi^{t}(x))$, typically in the form of the one-dimensional, discrete and equi-spaced time series $\{\beta_{n}\}_{n\in\mathbb{N}}$.

Takens’ theorem states that by choosing an embedding dimension $d\geq 2m+1$, where $m$ is the dimension of a compact manifold $M$, and a time lag $\tau>0$, then the map $\Phi_{\phi,\beta}:M\to\mathbb{R}^{d}$ given by

[TABLE]

is an embedding of $M$, where $\phi^{d-1}$ is the composition of $\phi$ $d-1$ times and $\mathbf{x}_{t}$ is the value of $\mathbf{x}$ at time $t$.

For noise-free data of infinite precision any time lag $\tau$ can be used; however, in practice, the choice of $\tau$ can influence the resulting embedding. In this paper $\tau$ was found by using the method of Least Median of Squares (LMS) ([57]) combined with the magnitude of the Fast Fourier Transform (FFT) of the signal. Specifically, we obtain the FFT spectrum and we identify the maximum significant frequency in the signal using LMS. We then use Nyquist’s sampling criterion to choose the delay value according to the inequality described in [58]. This approach yielded reasonable delay values in comparison to the standard mutual information function approach ([59]) where the mutual information function is plotted for several values of $\tau$ and the first dip in the plot indicates the $\tau$ value to use. This is because (1) the mutual information function is not guaranteed to have a minimum, thus leading to a failed selection of $\tau$ and, (2) the identification of the first true dip, if it exists, is not easy to automate especially in non-smooth plots.

The embedding dimension $d\in\mathbb{N}$ is computed by using the False Nearest Neighbor (FNN) approach ([60, 61]). In the FNN approach, the delay reconstruction is applied to the time series using increasing dimensions. The distances between neighboring points in one dimension are re-computed when the points are embedded into the next higher dimension. By keeping track of the percent of points that appear to be neighbors in a low dimension but are farther apart in a higher dimension (termed false neighbors), it is possible to identify a threshold that indicates that the attractor has been sufficiently unfolded. Applying FNN to all of the time series yielded the values in the range $d\in\{1,2,\ldots,10\}$, depending on the time series being reconstructed. Upon identifying $\tau$ and $d$ for each time series, delay reconstruction was used to embed the signal into a point cloud $P\subseteq\mathbb{R}^{d}$. The shape of the resulting point cloud was then quantified using persistence, as described in Section 3. We study five different methods to extract features from the resulting persistence diagrams as we show in Sections 4.2.1–4.2.5.

4.1.2 Persistence Diagram Computation with Bézier Curve Approximation

Bézier curves are introduced by Pierre Bézier and it has been used in computer aided design softwares and path planning applications for robots ([62, 63, 64, 65]). Recently, [31] utilized Bézier curves to speed up the persistence diagram computation. In this approach, the first step is to divide point cloud into groups. User should define number of samples per group ($spg$). Then, the Bézier curves is fit to each group individually. Second step is to generate line segments using these curves, and the number of line segments per curve ($r$) is also selected by user. Pairwise distance matrix between the line segments is computed and it is given to $Ripser$ as input, and it will provide the persistence diagrams as output based on selected maximum number of homology dimension. This section explains two of three main steps, which are fitting a bezier curve and computing the distance matrix between the line segments, and the effect of $spg$ and $r$ on approximation of persistence diagrams.

Fitting a Bezier Curve: In this implementation, we use cubic Bézier curve expression defined as

[TABLE]

where $t$ is the parametrization variable and $p_{i}$ is called the control points of the curve. For cubic Bézier curves, we totally have four different control points. The first one and the last one should match with the first and last point of the group of samples which we fit Bézier curve, respectively. Solution for these control points is obtained by using least squares error method and the expression is given as

[TABLE]

where $x_{i}$ is the data points of point cloud. Solution of $\nabla L=0$ provides the control points. The expression for $\nabla L=0$ can be rewritten such that

[TABLE]

where $t_{k}$ represents varying parametrization variable along a Bézier curve. The above expression can also be written in matrix form as

[TABLE]

The equations for $A_{k}$, $p$ and $b_{k}$ are defined as

[TABLE]

where $d$ is the dimension of the data set. To illustrate the Bézier curve fitting, we embed a sinusoidal signal to dimension two. Then, we divide the samples into groups and fit bezier curves. Figure 8 represents the control points and fitted lines on the groups.

Computing pairwise distance matrix: Bézier curves are split into intervals, and the number of intervals ($r$) are selected by the user. The end points of the intervals are connected to each other and line segments are generated. The next step is to compute distance matrix between the line segments. Lets assume that $\overrightarrow{l_{0}l_{1}}$ and $\overrightarrow{m_{0}m_{1}}$ represent two lines. The distance between these two lines is defined as follows

[TABLE]

where $l_{0}$, $l_{1}$, $m_{0}$ and $m_{1}$ are the end points of the two lines. The distance is computed by minimizing the function given as

[TABLE]

where $s$ and $t$ are parametrization variables. [31] solved this problem using gradient descent algorithm, however we employ the simplicial homology global optimisation (SHGO) ([66]). After finding the parameters making function $f$ minumum, we use them values to compute the distance between two lines. This is repeated for all combinations between line segments. We only compute the upper diagonal of the distance matrix, since it is symmetric. Then, it is given to $Ripser$ package as input to obtain persistence diagrams.

Effect of parameter selection: The parameters sample per group ($spg$) and the number of line segments ($r$) define how well we approximate to persistence diagrams of a time series. We choose a time series belong to 8.89 cm (3.5 inch) overhang distance, and we computed persistence diagrams of that time series with $r=1,3,5,7$ and $spg=100$. In Fig. 9, we plot both diagrams computed with $Ripser$ in blue and approximated diagrams in green color. All diagrams represents the first dimensional persistence ($H_{1}$) and they belong $50^{th}$ time series in the 8.89 cm (3.5 inch) case. Figure 9 shows that approximated diagrams converge to diagram obtained from $Ripser$ as we increase the number of line segments ($r$). We also validate this statement by computing Bottleneck distances whose definition is given as ([67]),

[TABLE]

where $X_{1}$ and $X_{2}$ represent two different diagrams and $\eta$ is the bijections between the points of the diagrams. Figure 10 shows the Bottleneck distances between the diagrams and it is seen that increasing $r$ results in smaller distances for the 8.89 cm (3.5 inch) case which means that larger values of $r$ approximates to blue one better. In addition, this type of behavior can be observed when we decrease the $spg$, since it will lead to larger number of groups and line segments.

4.2 Feature Extraction Using Persistence Diagrams

In this section, we explain how to extract features from persistence diagrams using persistence landscapes, persistence images, Carlsson Coordinates, kernel for persistence diagrams and signature path of persistence path’s signatures.

4.2.1 Persistence Landscapes

Persistence landscapes are functional summaries of persistence diagrams [68]. They are obtained by rotating persistence diagrams by $45\degree$ clockwise, and drawing isosceles right triangles for each point in the rotated diagram [37], see Fig. 11 where the landscape functions are represented by $\lambda_{k}$. Given a persistence diagram, we define the piecewise linear functions [68]

[TABLE]

where $b$ and $d$ correspond to birth and death times, respectively. Figure 11 shows that there are several landscape functions $\lambda_{k}(x)$ indexed by the subscript $k\in\mathbb{N}$. For example, the first landscape function $\lambda_{1}(x)$ is obtained by connecting the topmost values of all the functions $g_{(b_{i},d_{i})}(x)$ [68]. If we connect the second topmost components of $g_{(b,d)}(x)$, we obtain the second landscape function $\lambda_{2}$. The other landscape functions are obtained similarly. Note that the landscape functions are also piecewise linear functions.

Featurization of persistence landscapes: The persistence landscapes—$\lambda_{k}(x)$ where $x$ corresponds to the birth time—were computed using the persistence diagrams obtained from each of the embedded acceleration signals. Although these persistence landscapes can be utilized to featurize the persistence diagram, there is no one way to define these features. In this paper, we choose to extract a feature vector from the persistence landscapes by defining (a) a set of length $|\mathbb{K}|$ of the landscapes $\{\lambda_{k}\}_{k\in\mathbb{K}}$ where $\mathbb{K}\subset\mathbb{N}$ to work with, and (b) for the $k$th landscape, a mesh of non-empty, distinct birth times $\mathbf{b}_{k}=\{x_{i}\in\mathbb{R}\}$ where the corresponding values of the death times $\mathbf{d}_{k}=\{\lambda_{k}(x_{i})\mid x_{i}\in\mathbf{b}_{k}\}$ constitute the entries of the feature vector for the $k$th landscape. We can then combine features from all $|\mathbb{K}|$ landscapes to obtain the full feature vector $\mathbf{d}=\{\mathbf{d}_{k}\}_{k\in\mathbb{K}}$ that can be used with the machine learning algorithms.

Although the choice of $\mathbb{K}$, the set of landscapes to use, can be optimized using cross validation for example, in this study we set $\mathbb{K}=\{1,2,\ldots,5\}$ since it gives good results for our data. Similarly, the mesh may also be optimized in a similar way; however, this a more difficult task due to the infinite domain of $\mathbf{b}$, so we define our mesh as follows and as shown in Fig. 12.

Let $\lambda_{i,j}$ be the $i$th landscape corresponding to the $j$th persistence diagram from a training set in a supervised learning setting. Fix $i$ and overlay the chosen landscape functions corresponding to all of the persistence diagrams in the training set. Fig. 12 provides an example to this process and it selects second landscapes to extract features. Now project all the points that define the linear pieces of each of the landscape functions onto the birth axis. The red dots in Fig. 12 represents these projected points. Sort the projected points in ascending order and remove duplicates. The resulting set of points is our mesh $\mathbf{b}_{i}$ with length $|\mathbf{b}_{i}|$ for the $i$th landscape. We can repeat the same process to get the feature vector for all the $|\mathbb{K}|$ landscapes and construct the overall feature vector $\mathbf{b}$. We emphasize that a separate mesh is computed for each selected landscape number, and that the number of features will generally vary for each landscape function. Now to pull the features out of a given landscape function, we need to evaluate that function at the mesh points. Computationally, we efficiently accomplish that using piecewise linear interpolation functions.

Upon extracting the features from the persistence landscapes, we construct a feature matrix which represents all the tagged feature vectors. For instance, Table 2 shows an example feature matrix obtained from the first and second landscapes corresponding to each of the $n$ persistence diagrams in the training set. This table shows data with two labels: [math] for no chatter, and $1$ for chatter. It also denotes each feature with $y_{i,j}^{b_{i}}$ where $i\in\{1,3\}$ is the landscape number, $j\in\{1,2,\ldots,n\}$ is the corresponding persistence diagram number, while the superscript $b_{i}\in\{1,2,\ldots,|\mathbf{b}_{i}|\}$ is the feature number corresponding to the $i$th landscape. These feature matrices can then be used with supervised machine learning algorithms for example to train a classifier.

For our data, we computed the persistence diagrams and the corresponding first five persistence landscapes for each overhang case. We then randomly split the resulting landscapes into a training set ($67\%$) and test set ($33\%$), created feature matrices for each of the first five landscapes separately, and used SVM with ’rbf’ kernel, Logistic Regression and Random Forest algorithms to obtain and test a classifier. We repeated the split-train-test process $10$ times, and every time new meshes were computed from the training sets and these same meshes were used with the corresponding test sets. The mean accuracy and the standard deviation of the classification computed from $10$ iterations individually using each of the first $5$ landscapes can be found in Tables 10–12 in the appendix. In the results section we utilize the results with the highest accuracy for each of the overhang cases from Tables 10–12 when comparing the different TDA-based featurization methods.

4.2.2 Persistence Images

Persistence images are another functional summary of persistence diagrams [28, 37]. The first step in converting a persistence diagram $X=\{(b_{i},d_{i})\mid i\in\{1,2,\ldots,|X|\}\}$ to persistence images is to define the linear transformation

[TABLE]

which transforms the persistence diagram from the birth-death coordinates to birth-lifetime coordinates. Let $D_{k}(x,y):\mathbb{R}^{2}\to\mathbb{R}$ be the normalized symmetric Gaussian centered at $(b_{k},p_{k})$ with standard deviation $\sigma$ according to

[TABLE]

It was shown in [69] that the persistence images method is not very sensitive to $\sigma$, which we set to $0.1$ in this study. We also define a weighting function for the points in the persistence diagram $W(k)=W(b_{k},p_{k}):(b_{k},p_{k})\in T(X)\to\mathbb{R}$ according to

[TABLE]

Note that this is not the only possible weighting function, but it satisfies the requirements needed to guarantee the stability of persistence images [28]: it vanishes along the horizontal axis, is continuous, and is piecewise differentiable. Now define the integrable persistence surface

[TABLE]

The surface $S$ can be reduced to a finite dimensional vector by defining a grid over its domain and then assigning to each box (or pixel) in this grid the integral of the surface over that pixel. For example, the value over the $i,j$ pixel in the grid is given by

[TABLE]

where the integral is performed over that entire pixel. The persistence image corresponding to the underlying persistence diagram $X$ is the collection of all of the resulting pixels.

Featurization of persistence images: Persistence images can be used for support vector machine classification [39]. The corresponding feature vector is obtained from persistence images by concatenating the pixel values typically either row-wise or column-wise. The dimension of the resulting vector depends on the choice of the pixel size, i.e., the resolution of the persistence image. For example, let $I_{i,j}$ be a pixel in the persistence image, then a persistence image of size $100\times 100$ pixels is represented by the matrix

[TABLE]

and a typical feature vector is obtained by concatenating the entries of this matrix row-wise as shown by the rows in Table 3. The table shows a feature matrix where each persistence diagram is labeled either [math] or $1$, and the corresponding feature vector is shown using entries of the form $I_{i,j}^{k}$, where $k\in\{1,2,\ldots,n\}$ is the persistence diagram index while $i,j$ are row and column numbers, respectively, in the image.

We use Python’s PersistenceImages package to featurize the cutting signals, and then randomly split the resulting images into $67\%$-$33\%$ train-test sets. The persistence images have boundaries depending on lifetime and birth time ranges. Therefore, we find the maximum lifetime and maximum birth time by checking all diagrams of a data set. These maximum values can correspond to a point with significant lifetime and significant features can be lost if we set to boundaries to those values exactly. Accordingly, we sum each value with 1 to be able to capture all important features nearby them. We trained a classifier using SVM and the ‘rbf’ kernel, Logistic Regression, Random Forest and Gradient Boosting classifiers for two different pixel sizes: $0.05$ and $0.1$. The training and testing results for each different overhang case are available in Tables 13–14 of the appendix. When we compare the classification accuracy for persistence images to the other featurization methods, we choose the best results from this table for each cutting configuration.

4.2.3 Carlsson Coordinates

Another method for featurizing persistence diagrams is Carlsson’s four Coordinates [24] with the addition of the maximum persistence [22], i.e, the highest off-diagonal point in the persistence diagram. The basic idea of Carlsson’s coordinates is to utilize polynomials that (1) respect the inherent structure of the persistence diagram, and (2) that are defined on the persistence diagrams’ off-diagonal points. Specifically, these polynomials must be able to accommodate persistence diagrams with different numbers of off-diagonal points since the persistence diagrams can vary in size even if the original datasets are of equal size. Further, the output of the coordinates must not depend on the order in which the off-diagonal points of a persistence diagram were stored. The resulting features can be computed directly from a persistence diagram $X$ according to

[TABLE]

where $d_{\rm max}$ is the maximum death time, $b_{i}$ and $d_{i}$ are, respectively, the $i$th birth and death times, and the summations and maximum are each taken over all the points in $X$.

In order to utilize Carlsson coordinates, we compute the persistence diagrams from the embedded accelerometer signals and randomly split the data into training ($67\%$) and testing ($33\%$) sets. We then computed all five coordinates for each diagram, and utilized SVM, Logistic Regression, Random Forest and Gradient Boosting to train a classifier. The feature vectors that we tested were all $\sum\limits_{i=1}^{5}{\tbinom{5}{i}}$ combinations of these coordinates, where the term inside the summation is $5$ choose $i$. This revealed which combination of features yielded the highest accuracy in each iteration. The classification results for all of the different feature vectors are reported in Table 17 in the appendix. However, in the results section we utilize the feature vectors that yielded the highest accuracy when we compare the classification results of Carlsson coordinates to the other featurization methods.

4.2.4 Kernels for Persistence Diagrams

In addition to featurization methods, many kernel methods have also been developed for machine learning on persistence diagrams [29, 45, 46, 47, 48, 49, 50]. As an example, we choose to work with the kernel introduced by [29] which is defined for two persistence diagrams $X$ and $Y$ according to

[TABLE]

where if $z=(x,y)$, then $\hat{z}=(y,x)$, and $\sigma$ is a scale parameter for the kernel that can be used to tune the approach. for this study we investigated two values for this parameter: $\sigma=0.2$ and $\sigma=0.25$.

Given either a training or a testing set $\{X_{i}\}_{i=1}^{N}$ of labeled persistence diagrams, and using Eq. (6), we can define the kernel matrix

[TABLE]

Note that given two persistence diagrams $X$ and $Y$ whose number of points is $|X|$ and $|Y|$, respectively, then the corresponding kernel $\kappa_{\sigma}(X,Y)$ can be computed in $\mathcal{O}(|X|\cdot|Y|)$ time [29]. Therefore, the computation time for kernel methods is generally high and this can complicate optimizing the tuning parameter $\sigma$. To emphasize the effect of the computational complexity, in this paper, the long runtime for the $5.08$ cm ($2$ inch) overhang case caused by its large number of samples has forced us to report the corresponding classification results for a smaller number of iterations than the other overhang cases and the other featurization approaches.

For our data, we performed a $67\%$/$33\%$ train/test split of the labeled persistence diagrams. For each of the training and testing sets, we precomputed the corresponding kernel matrices and used Python’s LibSVM [70] for classification. For almost all but the $5.08$ cm (2 inch) overhang case where only $1$ iteration was used, we repeated the split-train-test process $10$ times and recorded the average and the standard deviation of the resulting accuracies. The resulting classification accuracies are reported in Table 8. Note that [29] describes another approach for training a classifier based on a measuring the distances between two kernels in combination with a $k$-Nearest Neighbor ($k$-NN) algorithm. However we do not explore this alternative method in this work, and only perform the computations using the kernel matrix and the LibSVM library.

4.2.5 Persistence Paths’ Signatures

Persistence paths’ signatures are a recent addition to featurization tools for persistence diagrams [30]. Let $\gamma:[a,b]\to\mathbb{R}^{d}$ be the piecewise differentiable path given by

[TABLE]

where each $\gamma_{t}^{i}=\gamma^{i}(t)$ is a continuous function with $t\in[a,b]$. The first, second, and third signatures, respectively, can be defined according to the iterated integrals [43]

[TABLE]

Other signatures are defined similarly, although, the computational cost significantly increases beyond the third level of signatures. The resulting path signatures can be used in classification algorithms as features. Looking back at persistence landscape functions in Section 4.2.1, we see that the $k$th landscape function $\lambda_{k}(t)$ can be written as a two-dimensional path

[TABLE]

Therefore, we can obtain signatures from persistence landscapes and use them as features in machine learning algorithms [30]. In this paper we use up to the second level path signatures. Specifically, let $\lambda_{r,i}$ be the $r$th persistence landscape corresponding to the $i$th persistence diagram. Then the signatures that we use from the $r$th landscape function are given by $\mathbf{S}(\gamma_{t}(\lambda_{r}(t)))=[S_{r,i}^{1},S_{r,i}^{2},S_{r,i}^{1,1},S_{r,i}^{1,2},S_{r,i}^{2,1},S_{r}^{2,2}]$.

By incorporating higher order signatures, or signatures from more landscape functions we can construct a longer feature vector for classification. For example, Table 4 shows the second level feature vectors computed using the first and second landscape functions for $n$ persistence diagrams.

In our experiment, we train a classifier using $75\%$ of the data, and we test using the remaining $25\%$. We construct a feature vector that includes up to the second path signatures for each of the first five landscape functions. Table 19 shows the classification accuracies for each configuration and for each landscape function. The best results in this table were used to compare the path signatures method to the other featurization procedures in Table 8.

5 Results and Discussion

5.1 Runtime Comparison

Runtime is a criterion for comparing the different feature extraction methods. For the TDA-based methods, the total runtime required for classification includes ruuntime of three main steps: (1) obtaining the persistence diagrams, (2) obtain features or computing kernels, and (3) training and testing the corresponding classifier. However, obtaining results with serial computing takes significantly larger runtime. Therefore, we implement paralell comuting to improve the runtime. We used High Performance Computing Center (HPCC) of Michigan State University for parallel computing. It includes several supercomputers which are composed of hundreds of nodes. Each node represents a computer and it has certain number of processors and RAM capacity. Users are allowed submit multiple jobs at the same time to HPCC, and they can define number of CPU per job and memory per CPU. To compute the embedding parameters and persistence diagrams in parallel, we request 10 CPU per job and 2GB of memory per CPU. The number of jobs submitted to HPCC changes depending on the number of time series of the overhang distance.

Embedded time series are subsampled such that every $10^{th}$ point is taken into account to compute persistence diagrams. We recorded the times to complete persistence diagram computation of all overhang cases and report them in Tab. 5. It also includes times for serial computing where one persistence diagram is computed at a time. It is seen that parallel computing reduces the computation time significantly, although most part of the runtimes for parallel computing is the queue time. Parallel computing can also be performed with some workstations available in market without having need of expensive supercomputers that HPCC have. Entry level workstations with CPU having 64 cores and 512 GB of RAM can be afforded by small workshops.

Despite long computation time for persistence-based methods, we note that after obtaining the persistence diagrams, they can be saved and used in multiple TDA-based classification methods. In addition, it was observed that delay and embedding dimension parameters do not change significantly for changing time series of a same cutting configuration. Parameters for embedding can be computed in the training phase of a classifier and they will be used in the test phase. Therefore, once these diagrams are computed, the time required for featurization and classification would be a fraction of the the ones reported in Tab. 5. Another point we wish to emphasize is that the most computationally expensive step is that of training a classifier. Once a classifier is trained, which can be done offline, the effort in classifying incoming streams of data is much smaller because a much smaller set of persistence diagrams and features are needed. Therefore, we compare the run times is needed for a single time series with different methods. Table 6 provide runtimes for embedding parameters computations and persistence diagram computation with different methods.

First column in Tab.6 represent the runtimes of computing embedding dimension and delay parameter. When these parameters computed, they can be saved and used in embedding time series In the second column of Tab. 6, the runtime of persistence diagram computation of subsampled point cloud is given. In this way, we use nearly 1000 points from the embedded time series to compute persistence diagrams. However, computation time for a persistence diagrams of a point cloud with that size is still high as seen from Tab 6. Therefore, we also employ greedy permutation subsampling and Bézier curve approximation technique. Greedy permutation option of $Ripser$ package is utilized and it is a method that subsamples the point cloud and computes the persistence diagrams with less number of points. $n_{perm}$ is a parameter that defines the number of points selected by greedy permutation algorithm. We chose $100$ and $300$ points for this option and reported runtimes for resulting persistence diagrams. In Tab. 6, we grouped the runtimes with respect to the number of points used in the corresponding method. It is seen that Bézier curve approximation with $r=1$ and $r=3$ uses approximately $100$ and $300$ points as in the case of greedy permutation. We provide runtime for both serial and parallel computing in Tab. 6. Parallel computing can only be applied to Bézier curve approximation among the methods of persistence diagram computation given in Tab.6. The reason is that persistence diagrams are obtained directly from the $Ripser$ package for other methods. However, the steps of Bézier curve approximation method, computation of coefficients for the line segments and the distance matrix between these line segments can be performed in parallel. In these two steps, a job can only compute coefficients of the lines in a single group or a distance between two lines. The number of the jobs will be equal to group number and number of combinations between lines for coefficient computation and distance matrix computation, respectively. Ideally, all jobs for a step can be computed simultaneously, if there is no queue time. Therefore, we recorded runtimes for computation of coefficients of the line segments in a single group, computation of a distance between two line segments and persistence diagram from a distance matrix individually. Then, we sum them up and reported in Tab. 6. Combining parallel computing with Bézier curve approximation reduces the runtime significantly. Moreover, it is seen that the fastest method is the greedy permutation with $n_{perm}=100$ and Bézier curve approximation computed in parallel places second. Both method is able to complete the diagram computation in less than a second, while runtime gets larger with increasing $n_{perm}$ and $r$ parameters.

Table 7 provides the times required to complete classification of a single time series. To be fair in comparison between runtime of WPT/EEMD and TDA-based methods, we assumed that the classifier is already trained and required parameters for all methods are selected. It is seen that WPT is the fastest method and EEMD places second. Runtime for TDA-based method is comparable to the ones for EEMD. Further the WPT and EEMD methods use codes that have been highly optimized, whereas the TDA-based methods are still under active research with huge future potential for improved optimization. We believe that the runtimes for TDA-based method can be further decreased with optimization.

5.2 Classification Scores

This section presents the classification accuracies for all the methods that are introduced in Section 4 and compares them to the results in Ref. [14] which uses the Wavelet Packet Transform (WPT) and the Ensemble Empirical Mode Decomposition (EEMD). The latter two methods are used for comparison since they are some of the currently most prominent methods for chatter identification using supervised learning. For persistence images, Template Functions and Carlsson Coordinates, we applied four different classifier and these are Support Vector Machine (SVM), Logistic Regression (LR), Random Forest (RF) and Gradient Boosting (GB) algorithms. All of these classifier except Gradient Boosting is used in Persistence Landscape method, while Kernel method and Persistence Paths results are obtained with LibSVM and SVM classifiers, respectively. The classification results are summarized in Table 8 where for each cutting configuration the best results of the classification algorithms for each method are included. Table 8 also includes the classification results obtained using a new TDA approach, which is not included in Section 4, based on template functions [25]. In this table, the best accuracy for each dataset is highlighted in green. Further, methods whose accuracy are within one standard deviation of the best result in the same category are highlighted in blue.

Table 8 shows that the WPT approach yields the highest classification accuracy for the $6.35$ cm ($2.5$ inch) and the $11.43$ cm ($4.5$ inch) overhang cases. However, we also see that for the $5.08$ cm ($2$ inch) and the 8.89 cm (3.5 inch) case, persistence landscape and Carlsson Coordinates yield the highest accuracy, respectively. For the $6.35$ cm ($2.5$ inch) overhang case we point out that the number of time series is small. Specifically, for this case, less than $10$ time series were divided into small pieces and used as the test set, see Table 1. Therefore, the $100\%$ classification accuracy using WPT for this case does not represent a robust result. Nevertheless, for the same case Table 8 shows that the TDA methods based on persistence landscapes, persistence images, template functions, Carlsson coordinates, and persistence paths yield better results than EEMD—a leading approach for chatter detection. For the $8.89$ cm ($3.5$ inch) case, Carlsson coordinates method yields the highest mean accuracy of $95.7\%$, placing ahead of both WPT and EEMD. Further, the other TDA-based method for this cutting configuration score classification accuracies of at least $83.9\%$. In addition, WPT and EEMD includes manual preprocessing which requires extra time and skill [14]. Therefore, the automation of these processes are not feasible, while all steps in TDA-based feature extraction can be fully automatized.

We also compare the performance of different persistence diagram computation methods whose runtime is compared in Sec.5.1. Figure 13 shows you the mean classification accuracies and error for the persistence diagrams obtained with the ways shown in Tab. 6. Subsampling point cloud at every $10^{th}$ point is the first method to compute persistence diagrams. Tab 6 shows that Bézier curve approximation has slightly larger computation time compared to greedy permutation subsampling method when it is computed in parallel. However, it is seen that for all overhang cases Bézier curve approximation method results in higher accuracy compared to greedy permutation. Also, its result are the closest results to the ones obtained from subsampled point cloud.

Increasing number of points and greedy permutation or increasing number of line segments ($r$) generated for a group does not always yield higher accuracy as we see from Fig. 13. The reason could be that increasing the number of points or the number of line segments can cause more topological noise on persistence diagrams. For example, Fig. 9 shows that the new point appears closed to significant feature, which is the point with highest lifetime, as we increase $r$ from $5$ to $7$. Therefore, this could cause small drops in accuracy for some overhang cases as seen from the Fig. 13.

Greedy permutation and Bézier curve approximation method can provide persistence diagrams in less than 1 second without optimization. For greedy permutation method, we can not apply any optimization, however coefficient computation for the line segments in Bézier curve approximation method can be optimized. This could further decrease the runtimes and opens the possibility for exploring in-situ chatter detection using TDA-based methods especially with properly optimized algorithms.

6 Conclusion

Two of the most common methods for chatter detection are based on extracting features from the time series using the WPT and the EEMD methods. However, even after obtaining a tagged time series, utilizing either of these methods requires the analysis of at least a subset of the available time series to identify which parts of the signal’s decomposition are the most informative when training a classifier. In addition to requiring highly trained individuals to perform this critical step, both WPT and EEMD are then locked into a small set of informative packets or IMFs which limits the transfer learning ability of the resulting classifiers [14].

In contrast to the WPT and EEMD methods, we use persistent homology—the most prominent analysis tool from TDA—to obtain a summary of the persistent topological features of the data. These are based on the global structure of the point cloud embedding of the acceleration signals in a turning experiment; therefore, upon obtaining a persistence diagram, there is no manual work involved in selecting the features from the persistence diagram. Since working directly with the resulting persistence diagrams is difficult, we investigated using the leading tools for feature extraction from persistence diagrams. The featurization methods that we studied are based on persistence landscapes, persistence images, Carlsson Coordiantes, a kernel method, template functions, and persistence paths’ signatures. The resulting features are then combined with several classification algorithms for training a classifier. The classification results are then computed from multiple split-train-test sets, and the resulting mean accuracies as well as the corresponding standard deviation are recorded for each method and for each cutting configuration.

Tables 8 summarizes the classification accuracy of all the TDA-based tools as well as their WPT and EEMD counterparts. In terms of the classification accuracy across the different cutting configurations we note the Carlsson coordinates and persistence landscapes approaches yield the best accuracies for two cases, and the former has the smallest runtime in comparison to the other TDA tools. On the other hand, in remaining cases, WPT yielded the best accuracies. Table 8 shows that WPT yields an accuracy of $100\%$ (with a standard deviation of zero) for the $6.35$ cm ($2.5$ inch) overhang case; however, as pointed out in Section 5 and Table 1, the size of the test set for this cutting configuration is too small which casts some doubts about the robustness of this result. Nevertheless, for this case both template functions and Carlsson coordinates still yield at least $86\%$ classification accuracy. Specifically, Table 8 shows that for the $6.35$ cm (2.5 inch) case, persistence landscapes, and Carlsson coordinates yield accuracies that are off by $13.7\%$, and $10.7\%$, respectively, from the WPT result. Similarly, for the $11.43$ cm ($4.5$ inch) case, Carlsson Coordinates and persistence images are within $13.3\%$ and $13\%$ of the EEMD result. For the 5.08 cm (2 inch) and $8.89$ cm ($3.5$ inch) overhang case, persistence landscapes and Carlsson coordinates yield the highest accuracy scoring $96.8\%$ and $95.7\%$, respectively, with tight error bounds that do not enclose the WPT accuracies.

As for runtime comparisons, we see dramatic decrease in persistence diagram computation when we computed persistence diagram in parallel as seen from Tab. 5. Parallel computing reduces the runtime for computing persistence diagram at least $82\%$. Table 7 shows that WPT is the fastest followed by EEMD; however, the reported time for TDA-based approach is comparable to the ones of EEMD. We can reduce the runtime for computation of persistence diagram of a single time series to less than 1 seconds using Greedy permutation and Bézier curve approximation method with combination of parallel computing. Fig.13 also shows that Bézier curve approximation method results in larger accuracies compared to Greedy permutation. Therefore, combination of parallel computing and Bézier curve approximation will make TDA-based approaches feasible option for online chatter detection.

To summarize, we show that persistence features are appropriate for chatter detection in cutting processes. These features have the potential to lower the barrier to entry when tagging cutting signals as chatter or chatter free because no manual pre-processing is needed before extracting and using the features in the persistence diagram. We also note that after obtaining a classifier, the time required for the classification of new incoming data will be greatly reduced, thus opening the door for future implementation of TDA methods in-situ for chatter detection and mitigation. We also believe that applying TDA based feature extraction techniques to different machining processes and machine tools can be included in future work.

Appendix A Simplicial complexes

Let $\{u_{0},\ldots,u_{k}\}\in\mathbb{R}^{d}$ be a set of data points, and the vectors defined between these data points ($u_{1}$-$u_{0}$,$u_{2}$-$u_{0}$,…,$u_{k}$-$u_{0}$) are linearly independent. A geometric *k-*simplex, $\sigma$ is a set of all points in $\mathbb{R}^{d}$ such that $\sum_{j=0}^{k}\lambda_{j}u_{j}$ where $\sum_{j=0}^{n}\lambda_{j}=1$ and $\lambda_{j}\geq 0$ for all j. Figure 14 provides illustrations for $0,1,$ and $2-$simplex. Each data point on a point cloud is represented as 0-dimensional simplex and they are called vertices. When two vertices are connected, an edge is formed and it is 1-dimensional simplex. Connection of three vertices will form $2-$simplex which is a triangle. Simplicies spanned by any subset of $u_{0},\ldots,u_{k}$ are the faces of $\sigma$. In general, $n-$simplex contains $n+1$ vertices, and the set of these simplices, are called geometric simplicial complexes, $K$, if the following two conditions are satisfied [54]:

1. If $\sigma\in K$, then faces of $\sigma$ are also in $K$,
2. If two $n-$simplex, $\sigma_{1}$ and $\sigma_{2}$ are in $K$, then the intersection of them is either common face or empty. The dimension of the simplicial complex, $K$ is equal to the largest dimension of its simplices.

Appendix B Expressions for persistence paths’ signatures

Let the $k$th landscape functions be $\lambda_{k}(x)=y$ where $x$ and $y$ represent coordinates along the birth time and the persistence axes. Since the persistence landscapes are piecewise linear functions, we can write them in closed form in terms of the nodes $\{x_{i},y_{i}\}_{i=1}^{n}$ that define the boundaries of each of their linear pieces according to

[TABLE]

The corresponding path is given by

[TABLE]

and its differential is

[TABLE]

Using the above definitions, we can derive general expressions for the first and second level signatures, respectively, according to

[TABLE]

The equations for the second level signatures is provided in Table 9.

Appendix C Classification results