arXiv:2604.21698v1 [cs.CL] 23 Apr 2026
Fixation Sequences as Time Series: A Topological Approach to Dyslexia Detection Marius Huber∗ David R. Reich† Lena A. Jäger Department of Computational Linguistics, University of Zürich, Switzerland {marius.huber,davidrobert.reich,lenaann.jaeger}@uzh.ch
Abstract
features are robust to noise and coordinate-invariant, and useful in data analysis, visualization and machine learning tasks [7, 10]. In recent years, persistent homology, a key tool within TDA, has been increasingly applied to a wide variety of data sets (for accessible introductions, see [8, 19]). Informally, computation of persistent homology is obtained from a so-called filtration on a data set, which can be thought of as a way of “sweeping” the data set with respect to some relevant feature. Persistent homology tracks topological features (the “shape”) of a data set across various stages of a filtration, producing a provably stable representation of the topology and geometry of the underlying data set [12].
Persistent homology, a method from topological data analysis, extracts robust, multi-scale features from data. It produces stable representations of time series by applying varying thresholds to their values (a process known as a filtration). We develop novel filtrations for time series and introduce topological methods for the analysis of eyetracking data, by interpreting fixation sequences as time series, and constructing “hybrid models” that combine topological features with traditional statistical features. We empirically evaluate our method by applying it to the task of dyslexia detection from eye-tracking-while-reading data using the Copenhagen Corpus, which contains scanpaths from dyslexic and non-dyslexic L1 and L2 readers.1 Our hybrid models outperform existing approaches that rely solely on traditional features, showing that persistent homology captures complementary information encoded in fixation sequences. The strength of these topological features is further underscored by their achieving performance comparable to established baseline methods. Importantly, our proposed filtrations outperform existing ones.
If the data set at hand is a time series, that is, a sequence of values indexed by time, persistent homology may be computed from a filtration that is defined on the time series itself, and it enables the characterization of geometric and dynamical structures intrinsic to temporal signals [21, 35, 42]. In this context, a filtration may be defined, for example, by thresholding the values the time series attains, and gradually letting that threshold sweep the entirety of the time series. For an illustration, see Figure 1, which shows a time series (in red) being swept by a horizontal line (in blue); see Example 2.1 for a full discussion of this. The resulting persistent homology can 1 Introduction thus capture properties of the time series such as its amTopological data analysis (TDA) has emerged as a power- plitude and periodic behavior. ful framework for extracting features from complex data The main contribution of this work is twofold: first, we sets that capture not only more traditional, statistics- extend existing topological methods for time series classibased properties of a data set, but also its shape. These fication, and second, we introduce the above framework to eye-tracking data and empirically evaluate whether TDAbased features add additional information that is not captured by traditional statistical or hand-crafted features. In existing work, persistent homology of a time series is usually computed with respect to a so-called horizontal filtration that sweeps the time series by horizontal lines
∗ Supported by the Swiss National Science Foundation (grant no. 209413). † Supported by the Swiss National Science Foundation (grant IZCOZ0 220330/1 (EyeNLG)) 1 All code to run experiments is available at https://github. com/DiLi-Lab/topological-dyslexia-detection and archived at https://doi.org/10.5281/zenodo.19238279.
1
2
of increasing heights [14, 33]. In particular, this horizontal filtration is agnostic to the temporal coordinate of the time series. In order to capture the temporal information contained in a time series, we introduce non-horizontal filtrations which explicitly depend on the temporal coordinate. In the case where the time series is a fixation sequence, these non-horizontal filtrations enable sensitivity to, for instance, the total number of fixations in a sequence, the time elapsed between consecutive fixations, or the presence and number of regressions. Second, and as a case study introducing persistent homology to eyetracking data, we consider the task of dyslexia detection from reading data. Concretely, we extend existing models for dyslexia detection to contain the TDA-based features, in order to investigate whether persistent homology can provide features which are complementary to more traditional, hand-crafted ones.
Background and Related Work
We now review how persistent homology may be used for time series classification, and then provide relevant background on dyslexia detection.
2.1
Persistent Homology for Time Series Classification
We provide a conceptual review of persistent homology in time series classification, as described, for instance, in [33] or [14]. For a more technical discussion, we refer to Appendix A. In the following, we use the term time series to refer to a finite sequence ((t1 , x1 ), . . . , (tn , xn )) of n pairs of real numbers, where the t1 < t2 < · · · < tn represent points in time, and x1 , x2 , . . . , xn denote the values that the time series attains at the respective points in time. In the context of eye-tracking data, for example, a fixation sequence naturally gives rise to a time series where the t-values represent the onset of a fixation and the x-values represent either the horizontal or vertical coordinate of the corresponding fixation. The former is expected to be more sensitive to eye movements within a single line than to those across lines, and vice versa for the latter.2 In our experiments, we therefore work with both time series extracted from a single fixation sequence (see Section 4.2 for details). Persistent homology [18, 20] is a topological tool that, broadly speaking, “sweeps” a time series and produces a so-called persistence diagram that captures characterizing features of the time series. One approach to this is to represent a time series (t1 , x1 ), . . . , (tn , xn ) by the graph in the (t, x)-plane that is obtained by drawing a straight line connecting the point (t1 , x1 ) with (t2 , x2 ), one connecting (t2 , x2 ) with (t3 , x3 ), and so on. Persistent homology is obtained by sweeping this graph by a horizontal line moving from bottom to top, and recording the evolution of the number of line segments that constitute the portion of the graph lying below the horizontal line.
To extract topological features from a fixation sequence of an individual reading a text, we first translate the information contained in it into a pair of time series: one containing the horizontal coordinates of the fixations indexed by the time of their onset, and the other containing the vertical coordinates. This separation of horizontal and vertical movements is done because our topological approach is suitable only for time series where the values are numbers (as opposed to vectors). In the context of eye movements, however, this separation also enables persistent homology to capture reading behavior constrained to a single line (from the horizontal movements) as well as that happening across lines (from vertical movements); this is particularly useful when dealing with fixation sequences stemming from reading of multiline stimuli, as is the case in the present article. From these two time series, we compute persistent homology, and combine the resulting topological features with the traditional ones obtained from existing models for dyslexia detection. We find that when adding these topological features, the resulting models outperform their counterparts that use only traditional features. Moreover, we find that even when removing all traditional features, the topological features on their own are sufficient to achieve near-stateof-the-art classification performance in dyslexia detection. This suggests that persistent homology captures aspects of fixation sequences that are complementary to those captured by traditional methods. For the above, we train and evaluate various pipelines on the Copenhagen Corpus of Eye Tracking Recordings from Natural Reading of Danish Texts (CopCo) [4, 24, 40], a large eye-tracking data set specifically designed for dyslexia research.
Example 2.1 (Sweeping the time series corresponding to a fixation sequence by horizontal lines). Consider the 2 One could generalize the notion of time series to one that captures both the horizontal and the vertical coordinate of a fixation. This, however, would yield a so-called multivariate time series, a format that is unsuitable for the methods of time series classification outlined here. Note also that with time series stemming from other domains such as finance or audio signals, we would not necessarily expect one coordinate to be more informative than others.
2
0 2 4 6 8 10 12 14 onset [100 × ms]
birth 0 2 4 6 8 10 12 14 onset [100 × ms]
12 10 8 6 4 2 0
birth 0 2 4 6 8 10 12 14 onset [100 × ms]
x-coord. [100 × px]
12 10 8 6 4 2 0
x-coord. [100 × px]
x-coord. [100 × px]
x-coord. [100 × px]
12 10 8 6 4 2 0
12 10 8 6 4 2 0
death
0 2 4 6 8 10 12 14 onset [100 × ms]
Figure 1: Graph corresponding to the time series from Example 2.1 (leftmost panel) and three stages of sweeping of the graph by a horizontal line for increasing heights (second through fourth panels). The number of line segments constituting the portion of the graph lying under the horizontal line changes from one to two and back to one as the height of the horizontal line changes from x = 200.0 to x = 400.0 and to x = 600.0, respectively. Table 1: An example of a fixation sequence.
below it now consists of a line segment near the point (366, 100.0) in the (t, x)-plane; see the second panel from Onset [ms] x-coordinate [px] y-coordinate [px] the left of Figure 1, where the horizontal line at height x = 200.0 is indicated in blue, and where the portion of 0 963.1 533.0 the graph lying below it is indicated by the thickened seg244 210.7 95.4 ment. As the horizontal line moves further upwards, we 366 100.0 130.8 see that a new line segment is “born” as the line passes 726 457.2 124.3 the value x = 281.6: at height, for example, x = 400.0 984 569.6 120.7 the portion of the graph lying under the horizontal line 1415 281.6 136.6 now consists of two separate line segments; see the second panel from the right of Figure 1. This newly born horizontal movements of the fixation sequence given in segment “dies” as the height of the horizontal line passes Table 1 (a procedure analogous to the following one ap- x = 569.6, since it gets absorbed into the3 line segment plies to the vertical movements). The corresponding time that was born earlier at height x = 100.0. Correspondingly, the number of segments is reduced back to one at series is specified by the sequence height, for example, x = 600.0; see the rightmost panel ((t1 , x1 ), . . . , (t6 , x6 )) = ((0, 963.1), (244, 210.7), of Figure 1. Note that this behavior is a consequence of the fact that the time series returns to a value near (366, 100.0), (726, 457.2), x = 200 toward its end (indicating a within-line regres(984, 569.6), (1415, 281.6)). sion of the reader). Indeed, if no within-line regressions Its graph is obtained by connecting consecutive (t, x)- were present, the graph in Figure 1 would lack the local pairs by straight lines; it is pictured in the leftmost panel maximum responsible for the death event. of Figure 1, where fixations are indicated by blue dots, Sweeping the graph corresponding to a time series thus and dotted lines indicate the t- and x-entries of the pairs yields a sequence of the number of segments present at that constitute the time series. Sweeping this graph starts each height of the horizontal line. This information is out with a horizontal line that is situated at a height recorded in a persistence diagram [17]; we illustrate this that is theoretically infinitely far below the graph. For a concept in the following example, and refer to Appendix A horizontal line placed at any height that is smaller than for details. the global minimum of the time series, the portion of the graph lying under that line is empty. As the height 3 One defines that “older” connected line segments absorb of the horizontal line passes the global minimum value “younger” ones (and not vice versa) for technical reasons; withx = 100.0, this changes: for a horizontal line at height, out this definition, computation of persistent homology would not for example, x = 200.0, the portion of the graph lying lead to a well-defined result [18, Chapter VII.1]. 3
1200
diagram from a time series, and is known as ordinary persistent homology. In our experiments, we also use extended persistent homology [13]. In a nutshell, this is obtained by first sweeping the graph corresponding to a time series by horizontal lines as above, and then sweeping the same graph in the opposite direction—from top to bottom—to capture additional information. We point out that extended persistent homology produces three persistence diagrams from a time series, and refer to Appendix A for details.
1000
death
800 600 400 200 0
In order to make persistence diagrams usable in machine learning pipelines, one must transform them into vectors of a fixed length. Methods for doing so abound [6, 9, 11, 16], and we choose to use the method known as persistence image [1]. In a nutshell, a persistence image is obtained from a persistence diagram by first applying a shearing that maps the diagonal of the persistence diagram to a horizontal line. The resulting diagram is then converted into a “heat map” by placing a bivariate Gaussian distribution with a given standard deviation (also referred to as bandwidth in the present context) at each point of the diagram. This yields a single-channel image commonly referred to as persistence surface; see Figure 3 for an illustration of the process of obtaining a persistence surface from a persistence diagram. Finally, the persistence image is obtained from the persistence surface by discretizing it into a pixel grid of a resolution prescribed by the user. We refer to Appendix B for details on persistence images, and conclude this section by pointing out that a persistence image created as above depends on the bandwidth and resolution used; in applications, these are treated as hyperparameters.
0 200 400 600 800 10001200 birth Figure 2: Persistence diagram corresponding to the time series from Example 2.1, obtained from sweeping by horizontal lines. The point located at (281.6, 569.6) records the line segment that is born at x = 281.6 and dies at x = 569.6. Example 2.2 (Persistence diagram). Consider the time series given in Example 2.1. By definition, a persistence diagram is obtained by placing a point at the coordinates (b, d) in R2 , where b and d record the height of the horizontal line at which a line segment is born and dies, respectively. In Example 2.1, the first line segment is born at height x = 100.0, whereas the second one is born at x = 281.6 and dies at x = 569.6. While the former never dies and thus does not have a death time, the birth and death times of the latter are recorded by including the point (281.6, 569.6) in the persistence diagram; see Figure 2. Additionally, a persistence diagram is endowed with the diagonal line which indicates the region of the persistence diagram where birth time equals death time; the greater the vertical distance of a point corresponding to a line segment to this diagonal line, the greater the difference between its death and birth time, and hence the longer this segment survives for. Line segments surviving for a long time are deemed to represent more “significant” features of the time series than short-lived ones. A persistence diagram thus provides a summary of the evolution of the line segments as one sweeps a time series by horizontal lines.
2.2
Predicting Dyslexia from Tracking-While-Reading Data
Eye-
Dyslexia is one of the most common learning disabilities that affects reading, impacting an estimated 9–12% of the population [27, 44]. For children with dyslexia, identifying the condition early is crucial for helping them succeed and keep pace in school [22, 45, 46]. Eye-tracking technology has emerged as a promising tool in this context. It offers a good balance between cost, invasiveness and data We point out that persistence diagrams computed from quality, and as the technology continues to improve and a time series like above typically contain many more become more affordable, it is increasingly realistic to use points than that pictured in Figure 2. Indeed, the “spar- it on a large scale. A single eye tracker in a school could sity” of that persistence diagram is a consequence of the serve multiple screening purposes, for example, evaluatfact that the underlying time series is very short. ing reading proficiency in a second language [3]. These Remark 2.3. The process above produces one persistence developments raise the question of whether dyslexia can 4
death–birth
death
D 7→ π(D)
birth
π(D) 7→ ρπ(D)
birth
Figure 3: Schematic showing the transformation of a persistence diagram D in (birth, death)-coordinates (left) to (birth, death-birth)-coordinates (middle) by means of a shearing π, and the resulting persistence surface ρπ(D) (right). conditional probability P (cj | Tj,k ), whereas in reader level aggregation, we aim to model P cj {Tj,k }k . The performance of our classification models will be evaluated using the Area Under the Receiver Operating Characteristic curve (ROC AUC). The ROC AUC provides a robust measure of discriminatory power, as it considers the trade-off between true positive and false positive rates across various classification thresholds. In particular, ROC AUC is less sensitive to class imbalance than accuracy, which is relevant in our case given the unequal distribution of dyslexic and non-dyslexic participants in the data set we used (see Section 5.1).
be detected from eye movements. Researchers have been exploring links between dyslexia and eye-movement behavior since the 1980s [32, 34, 37]. While early studies focused on deficits in eye-movement control, later research shifted attention toward difficulties in phonological decoding as the core cause of dyslexia. Importantly, eye movements are closely tied to phonological decoding [25] throughout life [5, 29, 30, 38, 39], which provides a strong theoretical basis for using eye-movement data to infer dyslexia. Indeed, a growing body of experimental research [2, 4, 23, 26, 28, 31, 36, 41, 43] shows that machine learning approaches can effectively screen for dyslexia using features derived from eye movements. Here, we use the work of [4] and [36] as reference methods.
3
Problem Setting
4
We empirically study the usefulness of casting fixation sequences as time series, and combining topological features computed from them with traditional eye-movement features. This is to investigate whether the former capture information complementary to that captured by the latter. We do this for the problem of dyslexia detection from eye-tracking-while-reading data for two levels of aggregation: trial-level and reader-level. For each individual j, their reading behavior during the reading of stimulus text k is recorded as a sequence Tj,k = (t1 , x1 , y1 ), . . . , (tMj,k , xMj,k , yMj,k )
Method
We inject our topological features extracted from fixation sequences into the two models introduced in [4] and [36] (the “baseline models”). We compare the resulting “hybrid models” trained with these features against the two baseline models, and against models that take into account only the topological features (the “TSH-models”, where TSH stands for “time series homology”). We perform this for both trial-level and reader-level aggregation, and including/excluding data stemming from non-native readers, which results in four experiment settings (here, a trial refers to a participant’s reading of one stimulus text; see Section 5.1 for details).
of Mj,k fixations, where tm denotes the onset of the m-th fixation and (xm , ym ) denotes its location. In our training data set D = {(Tj,k , cj )}j,k , each fixation sequence Tj,k is We now introduce an extension of the sweeping method associated with a binary label cj ∈ {0, 1}, indicating the discussed in Section 2.1, and then present our dyslexia deabsence (cj = 0) or diagnosis (cj = 1) of dyslexia in the tection pipeline for each of the three model classes menreader. In trial-level aggregation, we aim to model the tioned above. 5
x-coord. [100 × px]
x-coord. [100 × px]
12 10 8 6 4 2 0 0 2 4 6 8 10 12 14 onset [100 × ms]
the line. As in the case of the horizontal filtration, the birth and death times of the line segments being born and dying along the way are recorded in a persistence diagram. For an illustration of this, see Figure 4. In that figure, the graph corresponding to the time series is swept by a line of slope c ≈ −0.3 moving from left to right (indicated in blue), and the portion of the graph lying to the left of that line is indicated by the thickened segments.
12 10 8 6 4 2 0 0 2 4 6 8 10 12 14 onset [100 × ms]
Crucially, the sloped filtration is no longer agnostic to the time coordinate.6 For eye movements, this makes the filtration sensitive to, for instance, the total number of fixations. For the sigmoid and arctan filtrations, one replaces the sloped line by a curve described by the logistic and the arctan functions, respectively, and lets that curve move from the left to the right while keeping track of the number of line segments lying to the left of the curve, as before. These filtrations can “see” backward and forward in the time coordinate thanks to their curved shape, thus capturing, for instance, the time elapsed between two fixations.
Figure 4: Two stages of sweeping of the graph corresponding to a time series by a line of slope c ≈ −0.3 moving from left to right (indicated in blue). The number of line segments constituting the portion of the graph lying to the left of the line changes from two to one as the line moves from the left panel to the right one.
4.1
Computation of Topological Features and Non-Horizontal Filtrations
To compute the topological features of a time series, we use the horizontal filtration as described in Section 2.1.4 A shortcoming of that method is that this process is agnostic to the time coordinate of the time series. Indeed, the birth and death times of line segments are dependent only on the values of the local extrema of the time series, but not on the point in time at which these local extrema occur.5 We therefore introduce three new filtrations, which we call sloped, sigmoid and arctan filtrations. We illustrate the sloped filtration with an example, and refer the reader to Appendix C for technicalities, and details on the sigmoid and arctan filtrations.
4.2
Pipeline for Dyslexia Detection M
j,k Given a fixation sequence Tj,k = {(ti , xi , yi )}i=1 stemming from the k-th trial of individual j and consisting of Mj,k fixations, we wish to predict the binary label cj ∈ {0, 1}, indicating the absence or diagnosis of dyslexia. For the hybrid models combining traditional features with topological ones, we first compute the feature vecBL tor vj,k from Tj,k via the relevant baseline model. For the topological feature vector, we split up the fixation Mj,k x := {(ti , xi )}i=1 sequence Tj,k into the sequences Tj,k
Example 4.1 (Sloped filtration of a time series). Consider the time series from Example 2.1, whose corresponding graph is pictured in the left of Figure 1. To start, one places a line with slope c ̸= 0 far enough to the left in the (t, x)-plane so that the graph corresponding to the time series lies entirely to the right of that line (in experiments, c is a tunable hyperparameter). One then lets the sloped line move horizontally to the right while keeping track of the evolution of the number of line segments that constitute the portion of the graph that lies to the left side of
M
y j,k := {(ti , yi )}i=1 and Tj,k containing just the x- and ycoordinates of the fixation sequence, respectively. We iny x terpret each of Tj,k and Tj,k as a time series, min-maxscale their values, and compute persistent homology with respect to either the horizontal, sloped, sigmoid or arctan filtration; for the latter three filtrations, we treat the slope c as a tunable hyperparameter. This yields one x or three persistence diagrams (PDs) for each of Tj,k and y Tj,k , depending on whether ordinary or extended persistence is computed (cf. Remark 2.3). This results in two or six diagrams, respectively, for each initial fixation se-
4 The term “filtration” is simply a more technical term for a sweeping method. 5 As an example of this, observe that inverting the time coordinate of a time series and computing persistent homology using the horizontal filtration yields the same result as computing persistent homology of the original time series.
6 Indeed, inverting the time coordinate of the time series and using the sloped filtration will result in a persistence diagram that is different from the one computed from the non-inverted time series; cf. Footnote 5.
6
Time series x Tj,k PD(s) Fixation sequence Tj,k
PI(s)
flatten concat
Time series y Tj,k
PCA
PI(s)
PD(s)
flatten Feature vector TSH vj,k
TSH Figure 5: Computation of the topological feature vector vj,k from a fixation sequence. A fixation sequence Tj,k is y x split up into two time series Tj,k and Tj,k , after which the persistence diagrams (PDs) and the corresponding persistence images (PIs) are computed for each. The persistence images are each flattened, concatenated and TSH reduced via PCA to yield the feature vector vj,k .
and a SVM for the TSH-models.
quence Tj,k . Each persistence diagram is vectorized into a persistence image (PI) of a fixed resolution. The persistence images are flattened, min-max-scaled and concatenated into a single vector. Using a resolution of (n, m) for the persistence images, this resulting vector has length 2·1·n·m = 2nm if ordinary persistence is used, and length 2·3·n·m = 6nm if extended persistence is used; in our experiments, we set n = m = 50 throughout. Pixel values of persistence images are highly correlated, so we reduce this vector using PCA to obtain the topological feature vector TSH vj,k corresponding to the fixation sequence Tj,k ; in our experiments, we use PCA with 250 and 3·250 = 750 components for ordinary and extended persistence, respectively, as this number of components is empirically sufficient to capture the information stored in the pixel values. TSH , see For a schematic depiction of the computation of vj,k BL TSH into a Figure 5. Finally, we concatenate vj,k and vj,k single feature vector vj,k , containing both traditional and topological features of Tj,k . For trial-level aggregation, we pass vj,k into the classifier specified by the baseline model for prediction of cj,k ∈ {0, 1}; for the baseline introduced in [4], this is a support vector machine (SVM), whereas for the one introduced in [36], both a SVM and a random forest (RF) are used. For reader-level aggregation, we compute the mean of all vectors vj,1 , vj,2 , . . . stemming from individual j before passing the resulting vector into the respective classifier for prediction of cj ∈ {0, 1}.
5
Experiments
In the following, we describe the data set and evaluation protocol used in our experiments; all code to run them is available at https://github.com/DiLi-Lab/ topological-dyslexia-detection and archived at https://doi.org/10.5281/zenodo.19238279.
5.1
Data Set
The Copenhagen Corpus of Eye-Tracking Recordings from Natural Reading (CopCo; [24]) is a corpus of eyetracking data collected during natural reading of Danish texts. These consisted of 46 transcriptions of Danish speeches and 12 articles from Danish Wikipedia, which were presented to participants in multiline portions containing at most 10 lines. We refer to a participant’s reading of one such portion as a trial. Participants read a varying number of texts according to their own speed, and reading data was collected tracking the right eye, using an EyeLink 1000 Plus eye tracker with a sampling rate of 1000 Hz. The data set comprises recordings from 58 participants, including 22 L1 (native) speakers of Danish not diagnosed with dyslexia, 19 L1 speakers diagnosed with dyslexia, and 17 L2 speakers of Danish without dyslexia.7
For the baseline and TSH-models, we pass just the fea7 The non-L1 speakers stem from different L1 backgrounds. BL TSH ture vector vj,k and vj,k , respectively, into the classifier; Moreover, as in [40], we refer to the non-L1 speakers as L2 speakthe latter is either a SVM or a RF for the baseline models, ers for simplicity, even though Danish is L3 or L4 for some of the 7
Table 2: Search spaces used for randomized hyperparameter tuning.
We discarded the data from four out of the 22 L1 readers without dyslexia due to poor calibration or a diagnosed attention disorder. Moreover, the data from one of the L1 readers with dyslexia was discarded due to the absence of a dyslexia screening result, leaving data from a total of 53 readers. The fixation sequences used in our models were taken directly from CopCo. These were previously extracted using the DataViewer software by SR Research, which applies a velocity- and acceleration-based saccade detection method [24]. We discarded those trials consisting of fewer than five fixations. Each participant was assigned a binary label based on the original CopCo annotations: participants diagnosed with dyslexia were labeled as “dyslexic”, and all others as “non-dyslexic”. This results in a data set consisting of 4653 fixation sequences (out of which 1134 or 24.37% are labeled as “dyslexic”) on the trial-level, and consisting of 53 readers (out of which 18 or 33.96% are labeled as “dyslexic”) on the readerlevel.
5.2
Model class
Domain/values
∗
hybrid, TSH
Slope c Bandwidth σ †
[−4, −0.5] ∪ [0.5, 4] [10−3 , 10−1 ]
TSH
Kernel‡ C (linear kernel)† C (rbf kernel)† γ (rbf kernel)†
{linear, rbf} [10−2 , 101 ] [10−1 , 102 ] [10−4 , 10−2 ]
∗ Sampled uniformly with respect to angle † Sampled log-uniformly ‡ Sampled uniformly
Table 3: Mean ROC AUC scores by aggregation level and model class. Subscripts indicate standard deviation. Mean ROC AUC score Aggregation
Model class
Trial-level
Reader-level
Evaluation Protocol
We train the pipeline introduced in Section 4.2 for each of the hybrid, baseline and TSH-models, and in each of the settings specified by aggregation level and inclusion/exclusion of non-L1 readers. We evaluate the pipelines via k-fold nested cross-validation, where k = 5 for trial-level aggregation and k = 10 for reader-level aggregation, using Area Under the Receiver Operating Characteristic curve (ROC AUC) as the scoring function. In both the trial-level and reader-level aggregation, the folds are stratified according to dyslexia label. In the trial-level aggregation, the folds are additionally such that data from a given reader appears in precisely one fold, ensuring that the model always makes predictions about unseen readers. 5.2.1
Hyperparameter
Including L2
Excluding L2
Hybrid Baseline TSH
0.880.07 0.850.12 0.810.17
0.890.12 0.870.17 0.870.16
Hybrid Baseline TSH
0.920.10 0.890.08 0.880.17
0.990.02 0.970.04 0.950.07
on the parameters of their respective classifiers; we tune these according to their original implementations. For the TSH-models, the pipeline depends on the slope parameter c ̸= 0, the bandwidth σ > 0, and on the hyperparameters specifying the SVM classifier; we tune these via randomized search according to the distributions given in Table 2. Note that we excluded values between −0.5 and 0.5 for the slope c for computational stability. In all instances of randomized search, 200 iterations are performed.
Hyperparameter Tuning
For the hybrid models, the pipeline described in Section 4.2 depends on the slope parameter c ̸= 0 specifying the filtration, the bandwidth σ > 0 of the persistence image, and the parameters of the classifier stemming from the respective baseline model; we tune these via randomized search according to the distributions for the slope and bandwidth given in Table 2, and the distributions for the parameters of the classifier as given in the respective original implementation. The baseline models depend only
5.3
Results
In Table 3, we report the classification performance of the model classes described in Section 4. For each model class, we select the best-performing model and report its mean and standard deviation of the ROC AUC across all folds, and we highlight in bold the best performance in each of the four settings resulting from including/excluding L2-readers and aggregating on trial/reader-level (for a detailed breakdown of perfor-
participants.
8
mances, see Table D.1 in Appendix D). The results show that the hybrid models consistently outperform the baseline and the TSH-models, and, moreover, that the TSHmodels perform on par with the baseline models, or come close to doing so. Finally, the results demonstrate that our novel, non-horizontal filtrations mostly outperform the existing horizontal one.
6
tures derived from fixation sequences are a viable and effective addition to dyslexia detection pipelines, and likely to be useful in analysis of eye-tracking data in general. Future work should focus on validating these findings on larger and more diverse data sets as well as on tasks other than dyslexia detection, and on interpreting the topological features to better connect them to cognitive theories of reading and clinical practice.
Discussion
References
This study assessed the value of introducing novel filtrations for time series and of using these filtrations to extract topological features from fixation sequences with persistent homology. We provide strong empirical evidence that the proposed topological features encode information about reading behavior and are complementary to traditional features: across all experimental settings (trial-/reader-level and including/excluding nonnative readers), the hybrid models that combine topological and traditional features consistently achieved the highest ROC AUC scores. Moreover, seeing that each of the highest-scoring models makes use of one of the newly introduced filtrations, we empirically showed that these provide improvements over existing filtrations. The competitive performance of the TSH-models further indicates that topological representations of fixation time series are themselves highly informative. Despite abstracting away from explicit linguistic alignment and relying only on geometric and temporal structure, these models perform on par with established baselines. This suggests that dyslexia-related differences manifest not only in more traditional, aggregate reading measures, but also in higherorder temporal and structural patterns of eye movements that are captured by persistent homology. The introduction of non-horizontal filtrations appears particularly relevant, as they explicitly incorporate temporal ordering and spacing of fixations, properties that are central to reading dynamics. Performance gains are more pronounced at the reader level, which is expected given the aggregation of multiple trials into a more stable behavioral profile. At the same time, the very high reader-level ROC AUCs, especially when excluding L2 readers, should be interpreted cautiously due to the limited number of participants. The observed drop in performance when including L2 readers highlights population heterogeneity as an important factor. Overall, the results demonstrate that topological fea-
[1] Henry Adams, Tegan Emerson, Michael Kirby, Rachel Neville, Chris Peterson, Patrick Shipman, Sofya Chepushtanova, Eric Hanson, Francis Motta, and Lori Ziegelmeier. Persistence images: A stable vector representation of persistent homology. Journal of Machine Learning Research, 18(8):1–35, 2017. URL http://jmlr.org/papers/v18/16-337.html. [2] Thomais Asvestopoulou, Victoria Manousaki, Antonis Psistakis, Ioannis Smyrnakis, Vassilios Andreadakis, Ioannis M. Aslanides, and Maria Papadopouli. Dyslexml: Screening tool for dyslexia using machine learning, 2019. URL https://arxiv. org/abs/1903.06274. [3] Yevgeni Berzak, Boris Katz, and Roger Levy. Assessing language proficiency from eye movements in reading. In Marilyn Walker, Heng Ji, and Amanda Stent, editors, Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long Papers), pages 1986–1996, New Orleans, Louisiana, June 2018. Association for Computational Linguistics. doi: 10.18653/v1/N18-1180. URL https://aclanthology.org/N18-1180/. [4] Marina Björnsdóttir, Nora Hollenstein, and Maria Barrett. Dyslexia prediction from natural reading of Danish texts. In Tanel Alumäe and Mark Fishel, editors, Proceedings of the 24th Nordic Conference on Computational Linguistics (NoDaLiDa), pages 60– 70, Tórshavn, Faroe Islands, May 2023. University of Tartu Library. URL https://aclanthology.org/ 2023.nodalida-1.7/. [5] Hazel I. Blythe. Developmental changes in eye movements and visual information encoding associated with learning to read. Current Directions in Psychological Science, 23(3):201–207, 2014. doi: 10. 9
1177/0963721414530145. URL https://doi.org/ 10.1177/0963721414530145.
10.1145/1064092.1064133. URL https://doi.org/ 10.1145/1064092.1064133.
[6] Peter Bubenik. Statistical topological data analysis [13] David Cohen-Steiner, Herbert Edelsbrunner, and using persistence landscapes. Journal of Machine John Harer. Extending persistence using poincaré Learning Research, 16(3):77–102, 2015. URL http: and lefschetz duality. Foundations of Computa//jmlr.org/papers/v16/bubenik15a.html. tional Mathematics, 9(1):79–103, Feb 2009. ISSN 1615-3383. doi: 10.1007/s10208-008-9027-z. URL [7] Gunnar Carlsson. Topology and data. https://doi.org/10.1007/s10208-008-9027-z. Bull. Amer. Math. Soc. (N.S.), 46(2):255– 308, 2009. ISSN 0273-0979,1088-9485. doi: [14] Sebastiano Cultrera di Montesano, Herbert Edels10.1090/S0273-0979-09-01249-X. URL https: brunner, Monika Henzinger, and Lara Ost. Dynami//doi.org/10.1090/S0273-0979-09-01249-X. cally maintaining the persistent homology of time series. In Proceedings of the 2024 Annual ACM-SIAM [8] Gunnar Carlsson and Mikael Vejdemo-Johansson. Symposium on Discrete Algorithms (SODA), pages Topological data analysis with applications. Cam243–295, Philadelphia, PA, 2024. SIAM. ISBN 978-1bridge University Press, Cambridge, Cambridge, 61197-791-2. doi: 10.1137/1.9781611977912.11. URL UK, 2022. ISBN 978-1-108-83865-8. doi: 10.1017/ https://doi.org/10.1137/1.9781611977912.11. 9781108975704. URL https://doi.org/10.1017/ 9781108975704. [15] Department of Computer Science, ETH Zürich. [9] Mathieu Carrière, Steve Y. Oudot, and Maks OvsLecture notes: Topological data analysis (tda25), janikov. Stable topological signatures for points n.d. URL https://ti.inf.ethz.ch/ew/courses/ on 3d shapes. In Proceedings of the Eurographics TDA25/Script.pdf. Accessed: 2026-03-23. Symposium on Geometry Processing, SGP ’15, page 1–12, Goslar, DEU, 2015. Eurographics Association. [16] Barbara Di Fabio and Massimo Ferri. Compardoi: 10.1111/cgf.12692. URL https://doi.org/10. ing persistence diagrams through complex vectors. 1111/cgf.12692. In Vittorio Murino and Enrico Puppo, editors, Image Analysis and Processing — ICIAP 2015, pages [10] Frédéric Chazal and Bertrand Michel. An introduc294–305, Cham, 2015. Springer International Pubtion to topological data analysis: Fundamental and lishing. ISBN 978-3-319-23231-7. doi: 10.1007/ practical aspects for data scientists. Frontiers Ar978-3-319-23231-7 27. URL https://doi.org/10. tif. Intell., 4:667963, 2021. doi: 10.3389/FRAI.2021. 1007/978-3-319-23231-7_27. 667963. URL https://doi.org/10.3389/frai. 2021.667963. [17] Edelsbrunner, Letscher, and Zomorodian. Topological persistence and simplification. Discrete & Com[11] Frédéric Chazal, Brittany Terese Fasy, Fabrizio putational Geometry, 28(4):511–533, Nov 2002. ISSN Lecci, Alessandro Rinaldo, and Larry Wasserman. 1432-0444. doi: 10.1007/s00454-002-2885-2. URL Stochastic convergence of persistence landscapes https://doi.org/10.1007/s00454-002-2885-2. and silhouettes. In Proceedings of the Thirtieth Annual Symposium on Computational Geometry, SOCG’14, page 474–483, New York, NY, USA, [18] Herbert Edelsbrunner and John L. Harer. Computational topology. American Mathematical Soci2014. Association for Computing Machinery. ISBN ety, Providence, RI, 2010. ISBN 978-0-8218-4925-5. 9781450325943. doi: 10.1145/2582112.2582128. URL doi: 10.1090/mbk/069. URL https://doi.org/10. https://doi.org/10.1145/2582112.2582128. 1090/mbk/069. An introduction. [12] David Cohen-Steiner, Herbert Edelsbrunner, and John Harer. Stability of persistence diagrams. [19] Robert Ghrist. Barcodes: the persistent topology of data. Bull. Amer. Math. Soc. (N.S.), 45(1): In Proceedings of the Twenty-First Annual Sympo61–75, 2008. ISSN 0273-0979,1088-9485. doi: 10. sium on Computational Geometry, SCG ’05, page 1090/S0273-0979-07-01191-3. URL https://doi. 263–271, New York, NY, USA, 2005. Association org/10.1090/S0273-0979-07-01191-3. for Computing Machinery. ISBN 1581139918. doi: 10
[20] Robert Ghrist. Elementary applied topology. Cre- [26] Appadurai Jothi Prabha and Renta Chintala Bharatespace Independent Publishing Platform, North gavi. Prediction of dyslexia from eye movements Charleston, SC, September 2014. ISBN 978using machine learning. IETE Journal of Re1502880857. search, 68(2):814–823, 2022. doi: 10.1080/03772063. 2019.1622461. URL https://doi.org/10.1080/ [21] Marian Gidea and Yuri Katz. Topological data 03772063.2019.1622461. analysis of financial time series: Landscapes of crashes. Physica A: Statistical Mechanics and [27] Slavica K. Katusic, Robert C. Colligan, William J. Barbaresi, Daniel J. Schaid, and Steven J. Jacobits Applications, 491:820–834, 2018. ISSN 0378sen. Incidence of reading disability in a population4371. doi: https://doi.org/10.1016/j.physa.2017. based birth cohort, 1976–1982, rochester, 09.028. URL https://www.sciencedirect.com/ minn. Mayo Clinic Proceedings, 76(11):1081–1092, science/article/pii/S0378437117309202. Nov 2001. ISSN 0025-6196. doi: 10.4065/76.11.1081. URL https://doi.org/10.4065/76.11.1081. [22] Jonathan Glazzard. The impact of dyslexia on pupils’ self-esteem. Support for Learning, [28] Anna Laurinavichyute, Anastasiya Lopukhina, and 25(2):63–69, 2010. doi: https://doi.org/10. David Robert Reich. Automatic detection of dyslexia 1111/j.1467-9604.2010.01442.x. URL https: based on eye movements during reading in Rus//nasenjournals.onlinelibrary.wiley.com/ sian. In Wanxiang Che, Joyce Nabende, Ekatedoi/abs/10.1111/j.1467-9604.2010.01442.x. rina Shutova, and Mohammad Taher Pilehvar, editors, Proceedings of the 63rd Annual Meeting of [23] Patrick Haller, Andreas Säuberli, Sarah Kiener, the Association for Computational Linguistics (VolJinger Pan, Ming Yan, and Lena Jäger. Eye-tracking ume 2: Short Papers), pages 59–66, Vienna, Ausbased classification of Mandarin Chinese readers tria, July 2025. Association for Computational Linwith and without dyslexia using neural sequence guistics. ISBN 979-8-89176-252-7. doi: 10.18653/ models. In Sanja Štajner, Horacio Saggion, Daniel v1/2025.acl-short.5. URL https://aclanthology. Ferrés, Matthew Shardlow, Kim Cheng Sheang, Kai org/2025.acl-short.5/. North, Marcos Zampieri, and Wei Xu, editors, Proceedings of the Workshop on Text Simplification, Ac- [29] Mallorie Leinenger. Survival analyses reveal how cessibility, and Readability (TSAR-2022), pages 111– early phonological processing affects eye movements 118, Abu Dhabi, United Arab Emirates (Virtual), during reading. Journal of Experimental Psychology: December 2022. Association for Computational LinLearning, Memory, and Cognition, 45(7):1316–1344, guistics. doi: 10.18653/v1/2022.tsar-1.10. URL July 2019. doi: 10.1037/xlm0000648. https://aclanthology.org/2022.tsar-1.10/. [30] Sara V. Milledge and Hazel I. Blythe. The changing role of phonology in reading development. Vi[24] Nora Hollenstein, Maria Barrett, and Marina sion, 3(2):23, 2019. ISSN 2411-5150. doi: 10. Björnsdóttir. The copenhagen corpus of eye track3390/vision3020023. URL https://www.mdpi.com/ ing recordings from natural reading of Danish texts. 2411-5150/3/2/23. In Nicoletta Calzolari, Frédéric Béchet, Philippe Blache, Khalid Choukri, Christopher Cieri, Thierry [31] Mattias Nilsson Benfatto, Gustaf Öqvist Seimyr, Jan Declerck, Sara Goggi, Hitoshi Isahara, Bente MaeYgge, Tony Pansell, Agneta Rydberg, and Christer gaard, Joseph Mariani, Hélène Mazo, Jan Odijk, and Jacobson. Screening for dyslexia using eye tracking Stelios Piperidis, editors, Proceedings of the Thirduring reading. PLOS ONE, 11(12):1–16, 12 2016. teenth Language Resources and Evaluation Conferdoi: 10.1371/journal.pone.0165508. URL https:// ence, pages 1712–1720, Marseille, France, June 2022. doi.org/10.1371/journal.pone.0165508. European Language Resources Association. URL [32] Richard K Olson, Reinhold Kliegl, and Brian J https://aclanthology.org/2022.lrec-1.182/. Davidson. Dyslexic and normal readers’ eye move[25] International Dyslexia Association. Dyslexia baments. Journal of Experimental Psychology: Human sics, n.d. URL https://dyslexiaida.org/ Perception and Performance, 9(5):816–825, October dyslexia-basics/. Accessed: 2026-03-23. 1983. doi: 10.1037/0096-1523.9.5.816. 11
[33] Lara Ost, Sebastiano Cultrera di Montesano, and reading. Journal of Experimental Psychology: LearnHerbert Edelsbrunner. Banana Trees for the Pering, Memory, and Cognition, 24(2):476–497, March sistence in Time Series Experimentally. In Os1998. doi: 10.1037/0278-7393.24.2.476. win Aichholzer and Haitao Wang, editors, 41st International Symposium on Computational Geometry [40] David R. Reich, Shuwen Deng, Marina Björnsdóttir, Lena Jäger, and Nora Hollenstein. Reading does not (SoCG 2025), volume 332 of Leibniz International equal reading: Comparing, simulating and exploitProceedings in Informatics (LIPIcs), pages 71:1– ing reading behavior across populations. In Nicoletta 71:13, Dagstuhl, Germany, 2025. Schloss Dagstuhl Calzolari, Min-Yen Kan, Veronique Hoste, Alessan– Leibniz-Zentrum für Informatik. ISBN 978dro Lenci, Sakriani Sakti, and Nianwen Xue, ed3-95977-370-6. doi: 10.4230/LIPIcs.SoCG.2025. itors, Proceedings of the 2024 Joint International 71. URL https://drops.dagstuhl.de/entities/ Conference on Computational Linguistics, Language document/10.4230/LIPIcs.SoCG.2025.71. Resources and Evaluation (LREC-COLING 2024), [34] George Th. Pavlidis. Do eye movements pages 13586–13594, Torino, Italia, May 2024. ELRA hold the key to dyslexia? Neuropsycholoand ICCL. URL https://aclanthology.org/ gia, 19(1):57–64, 1981. ISSN 0028-3932. doi: 2024.lrec-main.1187/. https://doi.org/10.1016/0028-3932(81)90044-0. URL https://www.sciencedirect.com/science/ [41] Luz Rello and Miguel Ballesteros. Detecting readarticle/pii/0028393281900440. ers with dyslexia using machine learning with eye tracking measures. In Proceedings of the 12th In[35] Jose A. Perea and John Harer. Sliding windows ternational Web for All Conference, W4A ’15, New and persistence: An application of topological methYork, NY, USA, 2015. Association for Computing ods to signal analysis. Foundations of ComputaMachinery. ISBN 9781450333429. doi: 10.1145/ tional Mathematics, 15(3):799–838, Jun 2015. ISSN 2745555.2746644. URL https://doi.org/10.1145/ 1615-3383. doi: 10.1007/s10208-014-9206-z. URL 2745555.2746644. https://doi.org/10.1007/s10208-014-9206-z. [42] Lee M. Seversky, Shelby Davis, and Matthew Berger. [36] Peter Raatikainen, Jarkko Hautala, Otto Loberg, On time-series topological data analysis: New data Tommi Kärkkäinen, Paavo Leppänen, and Paavo and opportunities. In 2016 IEEE Conference on Nieminen. Detection of developmental dyslexia Computer Vision and Pattern Recognition Workwith machine learning using eye movement data. shops (CVPRW), pages 1014–1022, Las Vegas, NV, Array, 12:100087, 2021. ISSN 2590-0056. doi: USA, 2016. IEEE. doi: 10.1109/CVPRW.2016.131. https://doi.org/10.1016/j.array.2021.100087. URL https://www.sciencedirect.com/science/ [43] Soroosh Shalileh, Dmitry Ignatov, Anastasiya article/pii/S2590005621000345. Lopukhina, and Olga Dragoy. Identifying dyslexia [37] Keith Rayner. Do faulty eye movements cause dyslexia? Developmental Neuropsychology, 1(1):3– 15, 1985. doi: 10.1080/87565648509540294. URL https://doi.org/10.1080/87565648509540294.
in school pupils from eye movement and demographic data using artificial intelligence. PLOS ONE, 18(11):1–26, 11 2023. doi: 10.1371/journal. pone.0292047. URL https://doi.org/10.1371/ journal.pone.0292047.
[38] Keith Rayner, Sara C. Sereno, Mary F. Lesch, and Alexander Pollatsek. Phonological codes are auto- [44] Sally E. Shaywitz, Bennett A. Shaywitz, Kenneth R. Pugh, Robert K. Fulbright, R. Todd Constable, matically activated during reading: Evidence from W. Einar Mencl, Donald P. Shankweiler, Alvin M. an eye movement priming paradigm. Psychological Liberman, Pawel Skudlarski, Jack M. Fletcher, Science, 6(1):26–32, 1995. doi: 10.1111/j.1467-9280. Leonard Katz, Karen E. Marchione, Cheryl Lacadie, 1995.tb00300.x. URL https://doi.org/10.1111/ Christopher Gatenby, and John C. Gore. Funcj.1467-9280.1995.tb00300.x. tional disruption in the organization of the brain [39] Keith Rayner, Alexander Pollatsek, and Katherine S. for reading in dyslexia. Proceedings of the National Binder. Phonological codes and eye movements in Academy of Sciences, 95(5):2636–2641, 1998. doi: 12
10.1073/pnas.95.5.2636. URL https://www.pnas. at least) h. We point out that, in practice, the funcorg/doi/abs/10.1073/pnas.95.5.2636. tion f is often given as f = g|Γ(T ) for some function g : R2 → R. For ease of exposition, we will assume now [45] Joseph K. Torgesen. Individual differences that f is continuous and generic in the sense that no two in response to early interventions in reading: local extrema of f map to the same value under f . ObThe lingering problem of treatment resisters. serve that we have T f,h = ∅ for small enough values of Learning Disabilities Research & Practice, 15 h. Indeed, this is true for any value of h that is strictly (1):55–64, 2000. doi: 10.1207/SLDRP1501\ smaller than the global minimum of f , which we denote 6. URL https://www.tandfonline.com/doi/abs/ by fˇ. Similarly, we have T = Γ(T ) for any h ≥ fˆ, where f,h 10.1207/SLDRP1501_6. fˆ is the largest value that f attains on Γ(T ). Moreover, [46] Frank R. Vellutino, Jack M. Fletcher, Margaret J. we have that Tf,h1 ⊆ Tf,h2 whenever h1 ≤ h2 . Put toSnowling, and Donna M. Scanlon. Specific reading gether, {Tf,h }+∞ h=−∞ forms a sequence of nested subsets disability (dyslexia): what have we learned in of the graph of T that starts out as the empty set and the past four decades? Journal of Child Psy- eventually equals the entire graph Γ(T ). This means that chology and Psychiatry, 45(1):2–40, 2004. doi: {Tf,h }+∞ h=−∞ forms a filtration of Γ(T ). If h is set to some https://doi.org/10.1046/j.0021-9630.2003.00305.x. value smaller than fˇ and starts increasing, Tf,h ceases to URL https://acamh.onlinelibrary.wiley.com/ be empty as soon as h passes the value fˇ from below. doi/abs/10.1046/j.0021-9630.2003.00305.x. At that value of h, the number of connected components of Tf,h changes from zero to one.8 This logic generalizes from the global minimum value fˇ of f to any local miniA Persistent Homology of Time mum of f : whenever the parameter h passes the value of a local minimum of f from below, a new connected comSeries ponent of Tf,h is born. Conversely, whenever h passes the value of a local maximum f from below, the number of In this section, we review persistent homology in the conconnected components of Tf,h decreases by one, since in text of time series. In doing so, we will limit ourselves this case an existing connected component dies as it is abto those aspects of persistent homology that are relesorbed by one that was born earlier. By a similar line of vant to the setting of this paper, and we refer the reader −∞ argument, we see that the sequence T f,h h=+∞ forms a to [15, 18, 20] for details. In the present setting, we mean by a time series a piece- filtration of Γ(T ) (as h decreases from +∞ to −∞). The f,h increases by one wise linear function T : I → R, where I ⊆ R is a closed in- number of connected components of T whenever h passes a local maximum of f from above, and terval. Furthermore, we require that T have only finitely it decreases by one if it passes a local minimum of f from many points of non-linearity in its domain. Note that any above. such time series is univariate by definition. In practice, Persistent homology keeps track of the evolution of the a time series is specified by the finite sequence of tuples n 2 connected components of the f -sublevel and f -superlevel {(ti , xi )}i=1 ⊆ R , where xi := T (ti ), I = [t1 , tn ], and sets of a time series that we described in the previous t2 < t3 < · · · < tn−1 are the points of non-linearity of T paragraph by means of certain diagrams. Namely, given in the interior of I. The graph of a time series is defined 2 a time series T : I → R, and some function f : Γ(T ) → R, as Γ(T ) := {(t, T (t)) | t ∈ I} ⊆ R . the thresholding parameter h is first set to −∞ and then Given a time series T : I → R, a function f : Γ(T ) → R increased up to +∞. In this first phase, we record the and a number h ∈ R, the f -sublevel set of T at h is defined births and deaths of connected components of Tf,h by as means of a tuple (b, d) for each connected component, Tf,h := {(t, T (t)) | f (t, T (t)) ≤ h, t ∈ I} ⊆ Γ(T ). where b and d denote the birth and death of a component, respectively. Once h reaches +∞, the process is Symmetrically, the f -superlevel set of T at h is defined as reversed. That is, we now let h decrease back to −∞ in T f,h := {(t, T (t)) | f (t, T (t)) ≥ h, t ∈ I} ⊆ Γ(T ). 8 Here and in the following, a connected component is defined as one connected “piece”; in Example 2.1, connected components are
In words, Tf,h (resp. T f,h ) consists of that portion of the the line segments constituting the portion of the graph lying below graph of T on which f attains values of at most (resp. the horizontal line. 13
this second phase and record (birth, death)-tuples along the way as before. Note that the connected component corresponding to fˇ (the one that has the smallest birth time of all) never dies in the first phase. While its death time is thus technically d = +∞, we treat the first birth in the second phase as its death, and hence this component has death time d = fˆ.9 Consequently, we discard the event where h passes fˆ from above from the set of birth events of the second phase. The resulting collection of all (birth, death)-tuples is recorded in three persistence diagrams, each of whose horizontal and vertical axes measure birth and death time, respectively. The first and second of these diagrams are populated with the (birth, death)tuples stemming from the first and second phase, respectively, while the third diagram is populated with the single (birth, death)-tuple equaling (fˇ, fˆ). Moreover, each of the diagrams is endowed with the diagonal indicating the region where b = d, the idea being that the further a point is from this diagonal, the longer the lifetime (defined as the difference between death time and birth time) of its corresponding connected component is. These diagrams are called the ordinary, relative and essential diagram of f , respectively. Note that in the first phase, the death time of any connected component is larger than its birth time, whereas in the second phase the opposite is true. Hence all points in the ordinary and relative diagram are above and below the diagonal, respectively. The collection of the ordinary, relative and essential diagrams is referred to as extended persistent homology of a time series. In contrast, ordinary persistent homology of a time series is the result that is obtained from stopping the above process after the first phase, which produces just the ordinary diagram. In that diagram, the point corresponding to the component with birth time b = fˇ is either omitted, or recorded on a dashed line atop the rest of the points that serves as a placeholder for infinite death time. We conclude this section by providing an example of an extended persistence diagram (see Example 2.2 for an example of an ordinary persistence diagram).
part of Figure A.1, Tf,h (resp. T f,h ) is given by the portion of Γ(T ) that lies below (resp. above) the dashed line defined by x = h, and hence, for the choice of h depicted, consists of two (resp. three) connected components. As h passes the value fˇ = a2 from below, the connected component corresponding to the local extremum of T marked by the square is born. By definition, the death time of this component is d = fˆ, and we record this putting a point at (b, d) = (fˇ, fˆ) in the essential diagram of f . Similarly, the local minimum of T at t = e1 gives rise to a connected component of Tf,h with (b, d) = (e1 , d2 ), and we record this tuple in the ordinary diagram of f . Proceeding like this yields the persistence diagram given in the right part of Figure A.1, where we combined the essential, ordinary and relative diagram into one. For future reference, we call the filtration from Example A.1 the horizontal filtration of a time series [14, 33].
B
Persistence Diagrams and Machine Learning
In Appendix A, we discussed how (extended or ordinary) persistent homology produces one or three diagrams associated with a time series. In order to make these diagrams usable for machine learning, one must turn these diagrams into a format suitable for processing by machine learning algorithms, a procedure known as vectorization. Vectorization methods for persistence diagrams abound, and we limit ourselves to reviewing that known as persistence image, as introduced in [1]. n Given a persistence diagram D = {(bi , di )}i=1 ⊆ R2 consisting of (birth, death)-pairs, one starts out by discarding any points from D that have at least one infinite component. One then applies the linear transformation π : R2 → R2 given by π(x, y) := (x, y − x); this is a shearing mapping the diagonal of the persistence diagram D to the horizontal axis of R2 . Define the persistence surface of π(D) as
Example A.1 (Extended persistence diagram). Let ρπ(D) : R2 → R T : I → R be a time series and consider the function X (x, y) 7→ w(d) · φd,σ2 (x, y), f : Γ(T ) → R given by projection onto the vertical axis, d∈π(D) that is, f (t, T (t)) = T (t) for t ∈ I. The filtrations {Tf,h }h∈R and T f,h h∈R may be thought of as sweeping where w : R2 → R≥0 is a weighting function and Γ(T ) by a horizontal line at vertical height h. In the left φ 2 denotes the probability density function of a twod,σ dimensional normalized Gaussian with mean d ∈ R2 and 9 In technical terms, this choice stems from the fact that the 2 connected component with birth time fˇ dies in relative homology variance σ ∈ R>0 . In words, ρπ(D) is obtained by placing at time fˆ during the second phase of the filtration. a Gaussian at each point of the diagram π(D) (weighted 14
x
d
fˆ = d2
Γ(T )
d2
b2
b2
c2
c2 x=h
e2 fˇ = a2 a1 b1 c1
d1
e1
e2 a2
t
a 2 e2
c2
b2
d1
b
Figure A.1: Graph of a time series T with the horizontal filtration indicated by the line x = h (left), and the combination of the essential, ordinary and relative persistence diagrams (right). according to the position of the point) and summing these up. See Figure 3 for an illustration of this process. Finally, the persistence image is obtained by discretizing a relevant subdomain of ρπ(D) into a grid of pixels and assigning each pixel the integral of ρπ(D) over that pixel, resulting in an n × m-matrix, where n and m denote the height and width of the pixel grid, respectively. The persistence image thus created depends on the parameters σ > 0 and (n, m) ∈ Z2≥1 ; these are referred to as the bandwidth and resolution, respectively. In practice, the persistence image is usually flattened into a vector of length nm before feeding it to machine learning algorithms. We conclude this section by pointing out that, in practice, the weighting function w is often chosen to be w(x, y) = 1 or w(x, y) = y for (x, y) ∈ R2 . In the former case, all points of a persistence diagram are weighted equally, whereas in the latter case, points of D are effectively weighted according to their lifetimes.
This motivates the definition and use of non-horizontal filtrations, that is, filtrations of T that explicitly depend on the time coordinate. In the following, we define a general family of nonhorizontal filtrations for time series, from which we will pick three explicit filtrations for our experiments. To that end, let T : I → R be a time series, and let F : J → [0, 1] be continuous and bijective, where J ⊆ R is some interval (possibly of infinite length). Letting T̂ and Ť denote the global maximum and minimum values that T attains, respectively, we rescale and shift the function F to obtain G := (T̂ − Ť )F + Ť : J → [Ť , T̂ ].
(1)
By construction, the maximum and minimum values of G coincide with those of T , and, moreover, G is continuous and bijective, too. In particular, G is strictly monotonic. Defining Gh (t) := G(t − h), this implies that the graphs of Gh sweep out the strip R × [Ť , T̂ ] as h ranges from −∞ to +∞. Given a point (t, x) ∈ R × [Ť , T̂ ], this in turn guarantees that there exists a unique value h ∈ R C Non-Horizontal Filtrations for such that Gh (t) = x. In particular, this applies to any point (t0 , T (t0 )) ∈ Γ(T ), and we define fF : Γ(T ) → R Time Series as the function that maps a point (t0 , T (t0 )) ∈ Γ(T ) to the unique value of h ∈ R such that Gh (t0 ) = T (t0 ). As discussed in Appendix A, one can construct a filtration The function fF thus defines a filtration of Γ(T ), with of the graph Γ(T ) ⊆ R2 of a time series T : I → R from respect to which we can compute (ordinary or extended) a function f : Γ(T ) → R. Letting f be projection onto persistent homology of T . We call fF : Γ(T ) → R the the vertical axis of R2 yields the well-known horizontal filtration associated to F : J → [0, 1]. filtration of Γ(T ). Note that the horizontal filtration of a time series is sensitive only to the values of T at its local Example C.1 (Associated filtration). Let T : I → R be extrema, but not to the point in time at which they occur. a time series with global minimum and maximum values 15
t 7→ t + 1/2. Following the above, we have that Gh (t) = (T̂ − Ť )(t + 1/2 − h) + Ť . Given any point (t0 , T (t0 )) ∈ Γ(T ), one can check that Gh (t0 ) = T (t0 ) ⇔ h = t0 −
T (t0 ) − Ť T̂ − Ť
x-coord. [100 × px]
F : [−1/2, 1/2] → [0, 1]
x-coord. [100 × px]
Ť and T̂ , respectively, and consider the function defined by
12 10 8 6 4 2 0 0 2 4 6 8 10 12 14 onset [100 × ms]
12 10 8 6 4 2 0 0 2 4 6 8 10 12 14 onset [100 × ms]
Figure C.2: Sweeping curves defined by the logistic (left) and arctan (right) functions. While these curves look very similar to the naked eye, their asymptotic behaviors are quite different.
1 + . 2
This value of h is, by definition, the filtration value fF (t0 , T (t0 )). This filtration is obtained by sweeping out the strip R × [Ť , T̂ ] with a line of slope T̂ − Ť . Note this time series, one replaces the sloped line by a curve that, unlike the horizontal filtration of T , this filtration described by the logistic and the arctan functions, respecexplicitly depends on the time-coordinate of a point in tively, and lets that curve move from the left to the right Γ(T ). as before. The sloped filtration can be thought of as a tilted verWith the above, we are now in a position to introduce the non-horizontal filtrations that we will use in our ex- sion of the usual horizontal filtration of a time series. The idea behind the choice of the functions σc and τc in Defperiments. inition C.2 is that these functions—unlike λc —can “look Definition C.2. Let T : I → R be a time series. Given ahead” when sweeping out the strip R × [Ť , T̂ ] ⊆ R2 in a parameter c ∈ R \ {0}, the sloped filtration of T is the which Γ(T ) is supported; see Figure C.2 for illustrations filtration associated to the function of these filtrations. While these filtrations look very similar to the naked λc : [−1/2c, 1/2c] → [0, 1] eye, we choose to use them both in our experiments bet 7→ ct + 1/2. cause the logistic and the arctan functions converge to their asymptotic extrema at quite different rates. The sigmoid filtration of T is the filtration associated to The computation from Example C.1 generalizes to an the function explicit formula for the filtration level of a point of Γ(T ) in terms of the function F . σc : R → [0, 1] Proposition C.3. Let T : I → R be a time series with 1 t 7→ . global minimum and maximum values Ť and T̂ , respec1 + exp(−4ct) tively, and let F : J → [0, 1] be continuous and bijective, The arctan filtration of T is the filtration associated to where J ⊆ R is an interval. Then the filtration level of a the function point (t0 , T (t0 )) ∈ Γ(T ) under the filtration associated to F is given by τc : R → [0, 1] T (t0 ) − Ť 1 1 fF (t0 , T (t0 )) = t0 − F −1 . t 7→ arctan(cπt) + . π 2 T̂ − Ť Crucially, note that the functions λc , σc and τc in the above definition are continuous and bijective, and hence are valid functions to construct an associated filtration of a time series from. For an example illustrating the sloped filtration, see Example 4.1. For the sigmoid and arctan filtrations of 16
Proof. By definition, the filtration level of (t0 , T (t0 )) ∈ Γ(T ) is given as the solution to Gh (t0 ) = T (t0 ), which, in turn, is equivalent to (T̂ − Ť )F (t0 − h) + Ť = T (t0 ). Re0 )−Ť arranging, this is equivalent to F (t0 − h) = T (t . Note T̂ −Ť 0 )−Ť that T (t ∈ [0, 1], and hence the preceding equation is T̂ −Ť
equivalent to h = t0 − F −1
T (t0 )−Ť T̂ −Ť
1. We have that λ′c (0) = σc′ (0) = τc′ (0) = c, so that c can be regarded as a parameter that controls the slope of these functions in a consistent manner. In practice, the parameter c ∈ R \ {0} is a hyperparameter that can be tuned according to the problem setting.
, as desired.
This proposition allows us to provide explicit formulas for the sloped, sigmoid and arctan filtration of a time series. Corollary C.4. Let T : I → R be a time series with global minimum and maximum values Ť and T̂ , respectively, and let (t0 , T (t0 )) ∈ Γ(T ). Given a parameter c ∈ R \ {0}, the sloped, sigmoid and arctan filtrations of T are given by 1 1 T (t0 ) − Ť − , fλc (t0 , T (t0 )) = t0 − c 2 T̂ − Ť ! 1 T̂ − T (t0 ) , fσc (t0 , T (t0 )) = t0 + log 4c T (t0 ) − Ť and
2. As a consequence of Corollary C.4, we have that fσc (t, T (t)) = fτc (t, T (t)) = ±∞ whenever t ∈ I is such that T (t) ∈ {Ť , T̂ }. To avoid infinite filtration values in practice, we “pad” the range of the time series T by a percentage of its actual range. This amounts to replacing the values of Ť and T̂ in Equation (1) by Ť −ε(T̂ −Ť ) and T̂ +ε(T̂ −Ť ), respectively, for some small value of ε > 0.
D
Experiment Results
In Table D.1, we report for each model the mean and standard deviation of the ROC AUC across all k folds (where k = 5 for trial-level aggregation and k = 10 for reader-level aggregation). We denote by “TSH” and “BL” the models using respectively. time series homology and the baseline models, respecProof. The sloped filtration is, by definition, the one asso- tively. Moreover, we use subscripts to denote the filtraciated to the function λc (t) = ct + 1/2, t ∈ [−1/2c, 1/2c], tion used in the topological models, and to denote the 1 whose inverse is given by λ−1 c (x) = c (x − 1/2), x ∈ [0, 1]. name of the baseline models. For example, the model Hence, by Proposition C.3, the filtration is given by name “BLBjö +TSHhorizontal ” denotes the hybrid model combining the features obtained via the baseline from [4] 1 1 T (t0 ) − Ť with those obtained via the TSH-model using the horizon− fλc (t0 , T (t0 )) = t0 − , c 2 T̂ − Ť tal filtration. For the hybrid and TSH-models, we subdivide the result according to use of ordinary and extended as claimed. For the sigmoid filtration, note that the inverse of σc is persistence. We highlight in bold the best performing 1 given by σc−1 (x) = − 4c log x1 − 1 , x ∈ [0, 1]. From this, model in each of the four settings resulting from including/excluding L2-readers and aggregating on trial/readerwe obtain that level. ! 1 T̂ − Ť fσc (t0 , T (t0 )) = t0 + log −1 . 4c T (t0 ) − Ť 1 T (t0 ) − Ť 1 − fτc (t0 , T (t0 )) = t0 − tan π , cπ 2 T̂ − Ť
Rearranging the argument of the logarithm yields the claimed formula. Finally, for the arctan filtration, note that the inverse 1 tan (π(x − 1/2)), x ∈ [0, 1]. of τc is given by τc−1 (x) = cπ From this, we obtain that 1 T (t0 ) − Ť 1 fτc (t0 , T (t0 )) = t0 − tan π − , cπ 2 T̂ − Ť as claimed. We conclude this section with a couple of remarks. 17
Table D.1: Mean ROC AUC scores by model and aggregation level. Subscripts indicate standard deviation. Mean ROC AUC score
READER-LEVEL
TRIAL-LEVEL
Model name
Including L2
Excluding L2
Ord. persistence
Ext. persistence
Ord. persistence
Ext. persistence
BLBjö +TSHhorizontal BLBjö +TSHsloped BLBjö +TSHsigmoid BLBjö +TSHarctan BLRaa-RF +TSHhorizontal BLRaa-RF +TSHsloped BLRaa-RF +TSHsigmoid BLRaa-RF +TSHarctan BLRaa-SVC +TSHhorizontal BLRaa-SVC +TSHsloped BLRaa-SVC +TSHsigmoid BLRaa-SVC +TSHarctan BLBjö BLRaa-RF BLRaa-SVC TSHhorizontal TSHsloped TSHsigmoid TSHarctan
0.830.11 0.820.12 0.820.12 0.810.12 0.820.12 0.850.09 0.840.10 0.850.10 0.870.07 0.880.08 0.870.08 0.870.08
0.810.12 0.810.12 0.800.12 0.810.12 0.750.13 0.770.12 0.790.10 0.830.10 0.870.07 0.880.07 0.880.07 0.850.11
0.880.13 0.890.12 0.880.13 0.880.13 0.860.16 0.860.17 0.860.18 0.860.18 0.860.19 0.870.19 0.860.21 0.870.20
0.880.12 0.880.12 0.880.12 0.880.13 0.830.15 0.830.16 0.820.17 0.840.18 0.880.18 0.880.18 0.860.21 0.870.20
0.800.19 0.800.17 0.790.17 0.760.17
0.800.19 0.810.17 0.790.16 0.770.17
0.850.18 0.870.16 0.850.18 0.830.15
0.850.19 0.870.16 0.850.19 0.840.16
BLBjö +TSHhorizontal BLBjö +TSHsloped BLBjö +TSHsigmoid BLBjö +TSHarctan BLRaa-RF +TSHhorizontal BLRaa-RF +TSHsloped BLRaa-RF +TSHsigmoid BLRaa-RF +TSHarctan BLRaa-SVC +TSHhorizontal BLRaa-SVC +TSHsloped BLRaa-SVC +TSHsigmoid BLRaa-SVC +TSHarctan BLBjö BLRaa-RF BLRaa-SVC TSHhorizontal TSHsloped TSHsigmoid TSHarctan
0.890.07 0.770.14 0.840.15 0.800.15 0.840.16 0.850.18 0.890.12 0.870.11 0.890.11 0.840.13 0.920.10 0.870.11
0.750.08 0.800.15 0.840.13 0.760.11 0.640.10 0.790.16 0.810.12 0.860.09 0.900.13 0.860.13 0.880.12 0.900.09
0.960.05 0.920.10 0.920.08 0.900.16 0.860.09 0.820.11 0.880.14 0.800.20 0.980.03 0.950.07 0.990.02 0.990.02
0.870.13 0.880.19 0.940.05 0.910.09 0.850.12 0.890.06 0.930.10 0.940.05 0.970.07 0.950.07 0.990.02 0.990.02
0.810.16 0.820.14 0.850.12
0.840.21 0.870.17 0.870.18
0.890.08 0.880.12 0.880.10 0.840.16 0.810.18 0.830.21 0.880.17
0.930.08 0.940.08 0.970.04 0.840.16 0.840.12 0.820.11 0.830.17
18
0.920.09 0.950.07 0.930.10 0.800.19
0.870.11 0.900.12 0.920.05 0.870.11