Change Point Detection with Optimal Transport and Geometric Discrepancy
Nikita Pronko
stat.ME 62H30 (Primary), 62G05 (Secondary)

[
    [
    December 30, 2023
=====================


We present novel retrospective change point detection approach based on optimal transport and geometric discrepancy. 
The method does not require any parametric assumptions about distributions separated  by change points.
It can be used both for single and multiple change point detection and estimation, while the number of change points is 
either known or unknown.
This result is achieved by construction of a certain sliding window statistic  from which change points can be derived with elementary convex geometry in 
a specific Hilbert space.
The work is illustrated with computational examples, both artificially constructed and based on actual data.




  
Introduction
Change point problem was firstly mentioned by Stewhart <cit.> as a problem of industrial quality control in the year 1928. 
This problem was solved by Girshick and Rubin <cit.> in 1953. 
The optimal solutions for parametric  online formulation of the problem were provided by Shiryaev and Pollak <cit.>,<cit.>.
Asymptotically optimal posteriori change point detection method was developed by Borovkov <cit.>. 
This method requires knowledge of parametric likelihood functions, however it  can be translated to nonparametric case with empirical  likelihood <cit.>.
The fundamental results in nonparametric study of change-point problem are due Brodsky and Darkhovsky <cit.>, and Horvath et al. <cit.>. 
These methods are generally tailored  for one-dimensional data. 
Nonparametric change point inference with multidimensional data is  still an open problem. 
New methods are still being developed usually with applications of ideas from other scopes to change point data. 
For example, <cit.>, where divergence from cluster network analysis is used in order to estimate
both the number and the locations of change points or <cit.>, where multivariate version of Wilcoxon rank statistic is computed.


This work is focused on nonparametric change point detection in multivariate data.
This problem has large scope of applications including  finance <cit.>,genomics <cit.> and signal processing <cit.>.

The problem of optimal transportation were introduced by Monge <cit.> in the year 1781. The modern statement of the problem is due to Kantorovich as it was stated in the seminal paper in 1941 <cit.>.
We recommend <cit.> and <cit.> as modern references to the subject 

Vector rank and quantile functions were introduced by Serfiling <cit.> as part Depth-Outlyingness-Quantile-Rank (DOQR) paradigm. Monge-kantorovich vector ranks were introduced by Chernozhukov et al. <cit.>.

Study of discrepancy were started by H. Weyl <cit.> in 1916  with number theory as intended area of application, 
more computational aspects of discrepancy  were pioneered by Niederreiter <cit.> (quasi- and pseudo- random number generation) and Hlawka <cit.> (numeric integration). 
We say that discrepancy is geometric when we talk about continuous distributions with fixed well-identified support, this is the opposition to the case of combinatorial discrepancy which study nonuniformity of finite distributions.
The main contribution to the theory of quadratic generalized geometric discrepancy is by Hickernel <cit.>.
This is a work of a highly synthetic nature. 
This Synthetic nature is coming from the combination  of vector ranks and geometric discrepancy, 
optimal transport  and change point problems. 
All this concepts are rarely to never seen together in scientific publications, 
however in this paper they are intertwined in a singular computational methodology.

But why these results come to existence only now? 
The answer to this questions in our case is rooted into the fact that computing discrepancy used to be computationally challenging task for non-uniform probability distributions. 
However, in the recent publication <cit.> by Chernozhukov,Carlie, Galichon and Halin the Monge-Kantorovich vector rank functions were introduced. 
Vector rank function can be understood in the context of this work as homeomorphic extension of classical one-dimensional cumulative distribution functions to multidimensional spaces. 
The idea of using low-discrepancy sequences as a tool for the estimation of the vector rank functions was already present in the <cit.>. 
As the name suggests discrepancy is the natural measure of consistency of the sequence to be low-discrepancy sequence in the same way as the order statistic is the natural measure of the sequence to be ordered. 
Thus, as  classical one-dimensional ranks and order statistic is used in the construction of the classical Kolmogorov-Smirnov and Crámer-Von Mizes tests of goodness-of-fit, 
the vector rank function and discrepancy can be used in order to construct natural extension of this tests in the multidimensional context.

Thus, the main objective of the work is the exploration of vector rank discrepancy as  methodology for nonparametric high-dimensional statistic.
We select change point problem as the test subject for this inquiry.

The strict measure theoretic approach mentioned above is strengthened in this work as we speak about probability measures 
instead of probability distributions and densities whenever is possible which contradicts current trends in the statistical science. 
This can be justified by convenience of understanding  measures as points in certain multidimensional spaces bearing intrinsic geometric structure. 
On the other hand, we view all mathematical constructions in this work as theoretical prototypes of certain concrete algorithms ad data structures  
which can be  implemented in order to solve practical tasks. 
Our approach to the problem is highly inspired by <cit.> and we use convex geometric intuition whenever possible.



    acknowledgement

 Underlying research was conducted during author's work for master's degree in  "Mathematical methods of Optimization and Stochastics" program in NRU Higher School of Economics, Moscow, Russia.
 Author wants to thank head of the program Vladimir Spokoiny, his scientific advisor Andrei Sobolevski for suggesting the research topic and support, Alexey Naumov and Maxim Panov for the discussion in IITP RAS;
 Yuri Yanovich,Alexandra Suvorikova , Elena Chernousova and Alexey Kroshnin for fruitful discussion at the IUM seminar.   
 


  
Notation

This is notation for the use in the sequel.

^d is a d-dimensional euclidean space with d < ∞ .
I = [0,1] is a unit interval and I^d is a unit hypercube in ^d. 
In line with Kolmogorov's axiomatics triple (Ø,,) is a probability space, and every observation as a random variable is assumed to be  a Borel measurable map from Ø to ^d. For brevity Ω denotes the whole probability measure structure (Ø, , ).

Every random variable ξ : Ø→^d defines the probability measure P = ξ_# that is  referred to as the probability law of ξ, where ξ_# denotes a pushforward of the measure  by ξ . That is, for every Borel set B ⊂^d it holds 
that P(B) = (ξ^-1(B)). In this case the notation ξ∼ P is used. In case ξ is a random process over domain  notation ξ∼ P 
means that P defines the hierarchical  probability distribution system of ξ as a whole. 
In this case for all t,s ∈ the measure P_t is the distribution of a single observation ξ_t and P_t,s 
is the joint distribution of ξ_t and ξ_s i.e. ξ_t ∼ P_t, (ξ_t, ξ_s) ∼ P_t,s. 
ξ P means that ξ is independently identically distributed sample  (i.i.d) with probability law P.  

A set of all probability laws with finite first and second moments is denoted by ^d.
A subset of ^d  of probability measures absolutely continuous with respect to Lebesgue measure λ is denoted by ().  is the uniform distribution over I^d. It is obvious that ∈^d.In the sequel I^n is assumed to be implicitly equipped with 𝒰^n for each n ∈, which makes it into a probability space,  and  I^0 = {∅} is equipped with the counting measure.

If ξ_n is an  sequence of random elements in some metric space with distance metric ρ, when it is said that ξ converges to x in probability if for each ϵ∈_++ it
holds lim_n →∞( ρ(ξ_n, x) >  ϵ) = 0, which is denoted as  ξ_n →_ x. In case  ρ(ξ_n,x) is not a measurable function for all n the convergence in outer probability may be introduced and denoted by ξ→_∙ x. Having A^ϵ_n = {ω∈Ω : ρ(ξ_n(ω),x) > ϵ} the convergence in outer probability is equivalent to lim_n →∞inf{(B) | B ∈ : A_n^ϵ⊂ B } = 0 for all ϵ > 0.   

If ξ is a discrete time ergodic process with convergence in probability, 

    1/n∑^n_i = 1 f(ξ_i) →_∫_^d f  dP̅

for every bounded Lipschitz function f, we say that P̅ is the weak limit distribution of 
ξ.


  
Problem Formulation
In this section a brief review for common variations of change point problem is provided.

Ordered Set  is associated with time. It is possible to discuss change point problem with continuous time <cit.>. However, this work deals only with case of discrete time. It is also possible to consider a change point detection for random fields  <cit.>.

Two main types of change point problem are offline and online detection. In offline detection    is assumed to be a finite set of size || = T. 
Without loss of generality, it is assumed that = {1, …, T }. Initial observations X are treated as a time series indexed by . 
Firstly,  we discuss the case of a single change point.  

The change point of X is an unknown moment of time ∈ such that (X_t)_ t = 1 ^∼ P^1 
 and (X_t )^T_t =  + 1∼ P^2. There are two possibilities. Firstly, both subsamples may be i.i.d with  (X_t)^_t = 1P̅^1 and
  (X_t)^T_t =  + 1P̅^2.  It is also possible to consider unknowns P̅^1 and 
  P̅^2 to be ergodic weak limit distributions of  (X_t)^_t = 1 and
  (X_t)^T_t =  + 1  respectively.
  
   Then offline single change point detection problem is a hypothesis test for H_0 : P̅_1 = P̅_2 versus H_1 : P̅_1 ≠P̅_2 
   for every possible estimate  of  having 1 < < T.  Retrieving valid value of  is a related problem which will be referred as a change point estimation.
  
  In case of multiple change points  existence of up to N change moments 
  _1  < … < _n  < …_N ∈ is assumed. In simpler formulation of problem value of N is assumed to be known, while in more complicated one N is also an object of estimation. 
	 As the previous variation  the problem is a hypothesis test about distribution of samples 
  (X_t)_t = _n - 1 + 1^_n∼P̅_n where _0 = 1 and _N + 1 = T for simplicity. 
	 The problem splits into multiple hypothesis tests for H^n_0 : P̅_n = P̅_n + 1 versus H^n_1 : P̅_n ≠P̅_n + 1 for all possible estimates   of  in case of known N. 
	 For unknown N the problem is structured around testing against H_0 :P̅_1 ≠P̅_2  ≠…≠P̅_N̂ for all estimates 
	 ( N̂, ) of (N,) where possible values of N̂ are  constricted to some meaningful finite set.

  
For online change point problem  ≅ as an ordered set, and X can be thought as infinite times series unfolding in  real time.
In this case the goal is to identify change point   as soon as possible which means using minimal amount of observations X_t with t ≥.
Branding every moment of time as a change point will achieve zero delay. However, this detection procedure is unacceptable as it achieves maximal false alarm rate. 
This means that change-point detection algorithm needs to be tuned in for minimal delay with fixed false alarm rate.

  In a more theoretical framework  a single change point _0 can be considered, which implies that infinite set { X_t : t < 0 } of observations with prior distribution can be used for detection. 
  However, in more practical situation the stream of chanege ponts : → is considered, with only finite sample { X_t : _n - 1 < t ≤_n } available for detection of _n.

Generally, methods of change point detection are focused on detection of change in such properties of statistical distributions as mean, variance and spectral density.
The method presented in this work is focused on change in ergodic  weak limit distribution of data. 
Which means that change in mean and variance must be detectable, although, 
change in spectral density is not. 
Moreover, changes restricted to such properties of distribution 
as median or higher moments also must be traceable.
Another restriction of proposed solution is a demand for all measures 
P̅_n to belong to ^d.

Construction of proposed change point detection starts with  offline single detection problem.
 We restrict our attention to models with i.i.d distributed observations along one side of change point.
 In order to give meaningful characteristic to point estimates of change point this problem can be reformulated as certain probabilistic game against Nature. 
 Two  sample spaces (Ø^1, ^1, ℙ^1) and (Ø^2, ^2, ℙ^2) with corresponding families of random variables ξ^1 and ξ^2 taking values in ^d are assumed to exist. We assume that both families
  are intrinsically i.i.d which means that ξ^1  P_1 and ξ^2  P_2 and the same holds for ξ^2. This means that both sample probability spaces can be decomposed as Ω^1 = ∏^∞_t = 1Ω^1_t and 
  Ω^2 = ∏^∞_t = 1Ω^2_t  
  such that for each ω∈Ø^1  the random variable ξ_t^1(ω)  depends only on 
  ω^t.
 The common sample space is constructed by the rule

    Ω_t = Ω^1_t ⊔Ω^2_t 
    = (  {1}×Ω^1_t ) ∪( { 2 }×Ω^2_t ),
      Ω = ∏^∞_t = 1Ω_t
 
with σ-algebras generated by a set  {{1 }× A : A ∈_1}∪{{2 }× B : B ∈_2} . Note, that the only possible probability distribution on  
the set {1,2 } is  a Bernoulli   distribution Bern(r) with a parameter r ∈ [0,1]. So probability measure on Ω can be defined by ℙ(S) = rℙ_1(ι^-1_1 S) + (1 -r)ℙ_2(ι^-1_2 S), where ι_j is a natural injection defined by ω↦ (j,ω).  Nature defines value of r and distribution laws  P_1 and P_2. Random variables (b_t,X_t) are constructed by 
letting (b_t,X_t)(j, ω) = ( j, ξ^j_t(ω) ). Then, nature generates sample of T observations of (b,X), sorts it by b, then erases values of b which retains observable sample of X. Observer either claims that there were no change point which leads to a victory in case 
P_1 = P_2 or equivalently r is too close to either 0 or 1. Otherwise, the observer gives point-estimate r̂ of r which can be converted to a classical change-point estimate = r̂ T. 

This construction also can be understood as a Bayesian assumption ∼B(T,r), where B(T,r) stands for binomial distribution. 
However, the change point problem is a single element sample problem, which means that only minimal properties of the parameter can be inferred from the data. This does not change the problem except for treating  as an estimate of the mean. 

Obviously,the  observer chooses strategy which maximizes chances of victory. This can be interpreted in terms of loses L, which are equal to L(0) = min{  1 - r, r }  = L(1) in case observer claims no change point and L(r̂) =  | r̂ - r  | otherwise.
Let  L(r̂) = min1 - r̂, r̂  in case the observer gives estimate when no real change point exists.
 
 If one preserves r while increasing T the empirical distribution of X will converge to a measure μ_r = r P̅_1 + (1 - r) P̅_2, the convex combination of measures. It possible to think of μ_r as contained inside one-dimensional interval [P̅_1, P̅_2] embedded into the infinite-dimensional convex space of probability measures. In order to construct such measures from random variables we will use random mixing map I^s_Z : Y ↦ b_sY + (1 - b_s)Z in the sequel. Here b_s ∼Bern(s) and s ∈ [0,1].
  
 
In order to propose the solution strategy for  the change point problem we  introduce vector rank function  R : ^d →^d such that R_#μ_r =, a known reference distribution of simple structure. For our application  is a uniform distribution over a unit hypercube I^d. It was shown in <cit.> and <cit.> that R̂, an estimate of R , can be recovered from data without any parametric assumptions. 
 
If data is i.i.d distributed with laws P_1 and P_2 on each side of change point, when  figure 1 will  exhibit relations between probability distribution laws discussed so far. 

The diagram shows that application of vector rank function can be thought as a probabilist's "change of basis". The pragmatic value for change point problem is in elimination of unknown distribution μ_r from the model.  

It is obvious that the diagram depicted at fig. 1 commute. As the pushforward acts as a linear map of signed measures defined over two different measurable spaces. If α and β are measures, x, y ∈, B is a measurable set and f is a measurable function, then  
 
    f_#(xα + yβ)(B) = x α( f^-1(B))  + y β(f^-1(B))
      = x f_#α (B)    + y f_#β(B).
  
 Therefore, by linearity of the pushforward
 
    = R_#μ_r =   R_# ( r P_1 + (1 -r) P_2) = r R_# P_1 + (1 - r)R_# P_2 .

 Thus, the diagram is correct.

The real result of change-point estimation will depend on measurable difference between P_1 and P_2. We introduce notion of Kullback-Leibler divergence in order to quantify this difference.

 Kullback-Leibler divergence between measures P_1 and P_2 admitting densities p_1 and
and  p_2 respectively is

    KL(P_1 , P_2) = ∫_Ø_1logp_1/p_2 dP_1


Note  that while P_1,P_2  ∈^d they admit densities. Now we formulate a sufficient condition on R that ensures that our transformations preserve divergence between distributions.

If R is invertible almost everywhere P_1, when KL(R_#P_1,R_#P_2) 
 = KL( P_1,P_2)
 Let X be a random variable such that  X ∼ P_1. Then, R(X) ∼ R_# P_1 and by
  invertibility  a.e. of R  distribution R_# P_1 has density p_1 ∘ R^-1 a.e. ; similar fact is also true for P_2 .
  Then,
 
    KL( R_# P_1, R_# P_2   ) = log p_1(  R^-1  R  X) / p_2 ( R^-1  R  X )    = log p_1 (X) / p_2  ( X)   =  KL( P_1,P_2)



  
Formulation of the method
In this chapter our approach to construction of the vector rank function's estimates R̂_n is presented. The main goal of This approach is  the avoidance of approximation of the entire function R. 
Note, that even the estimation that the aim of this task is actually not a precise estimation of values R(X_i) for each observation X_i but a construction of 
an estimate that preserves the relative geometry of a sample X along its change points. 


The work <cit.>  presents continuous, semidiscrete, 
and discrete methods of estimation of R. 
The Implementation of change point detection discussed in this section is based solely on discrete approach. The reason behind this choice is computational simplicity.

Vector rank function in general can be understood as (P, mathbbU)-almost surely homeomorphisms  between supports of sample probability distribution P and reference probability law of choice 𝕌
(R is defined almost everywhere P and continuous with respect to subset topology, and admits similar continuous inverse existing almost everywhere 𝕌)
, 
such that  R_# P = 𝕌 and both depth median and symmetries are preserved. 
This symmetries  and the concept of depth median need to be additionally defined, 
which goes beyond of the scope of this paper. In our application uniform measure over unit hypercube I^d is selected as the reference U.
In case d = 1 the cumulutive distribution function (cdf) of P is also the 'vector' rank function of the probability distribution of P.


The vector rank function is not unique in general. We use Monge-Kantarovich vector rank developed in <cit.>, which is defined as optimal transport between 
P and U

    R = _R : R_#P = ∫_^d d^2(R(x),x)  dP(x).

Vector rank role of the optimal transport map R can be intuitively justified by the equivalence of the above optimization problem to the maximization of the integral

    ∫_^d⟨ R(x) - m() , x - m(P)  ⟩ dP(x),

where m() and m(P) stands for depth medians of distributions  and P respectively. 
Thus, the optimal transport preserves geometric properties of distribution P which can 
be expressed by inner products of points relative to its center of symmetry. This is what is understood as the relative geometry of the data 
,and what we try to preserve during  vector rank estimation. 

In this work we are focused on the discrete estimation of vector rank function.
Let u be an equidistributed sequence of points in I^d. That is, for every Lipschitz continuous function  f defined on I^d it holds

    lim_n →∞1/n∑^n_i = 1 f(u_i) = ∫_I^d f(x)  dx.

Then, for T observations define Y_i = R̂_T(X_i) = u_σ^*(i), where 

    σ^* = _σ∈ S^T∑^T_i = 1  X_i - u_σ(i)_2^2.

In case convergence in (<ref>) is understood as convergence a.s or even as convergence in probability the sequence u can be taken as an i.i.d  
sequence of random variables  or as an ergodic process with a uniform weak limit distribution. 
However, this  implementation is more suited for  deterministic form of u. 
Discrepancy can be understood as a natural measure of slackness of condition (<ref>) for a fixed value of T.

Classical (one-sided) geometric discrepancy of the sequence (u_i)^n_i = 1 is described by the formula

    D(u) = max_x ∈ I^d| | { u_i | 1 ≤ i ≤ n}∩ [0,x] |/T  - ∏^d_j = 1 x^j    |.

Note that (<ref>) admits representation D(u) = φ_u _∞ for a certain function φ_u. 
This idea were used in <cit.> to introduce generalized quadratic discrepancy with supremum-norm replaced by a 
certain quadratic norm in a certain Sobolev Space. The result can be expressed  as

    (D_2^κ(u) )^2 = 
    	
     =
    	∬η(x,y)  dx dy - 2/n∑^n_i = 1∫η(x,u_i)  +1/n^2∑^n_i,j = 1η(u_j,u_i),

where η is the reproducing Hilbert kernel of the Sobolev space 

    η(x,y) = ∏^d_i = 1( M + β^2( κ(x^i) + κ(y^i) + 1/2B_2( x^i - y^i  1  ) + B_1(x^i)B_1(y^i) )

with β∈ standing for  a scale parameter and κ standing for a functional parameter with square-integrable derivative  with ∫^1_0 κ = 0 and the value M defined by

    M = 1 - β^2 ∫^1_0 (κ')^2,

and B_i is the ith Bernoulli polynomial. 
Note, that in case  d = 1 statistic (<ref>) turns into Kolmogorov-Smirnov statistic and (<ref>) turns into Cr/`amer-Von Mizes statistic for uniform distribution test.

The function κ is a functional parameter defining exact form of the discrepancy. 
In this text we use star discrepancy produced by

    κ^* = 1/6 - 1/x^2

and the centred discrepancy produced by selecting

    κ^c = -1/2 B_2( x - 1/2 1)


The case of u being a low-discrepancy sequence is a particularly well suited for our application.

low-discrepancy sequence is a deterministic sequence u in I^d designed with a goal of minimizing value the D^κ_p(u_i)^n_i = 1 for each natural number n.

Definition above is rather informal. However, it is postulated firmly that for any low-discrepancy sequence u property (<ref>) holds and lim_n →∞ D^κ_p(u_i)^n_i = 1 = 0 for any choice of p and κ. 
Thus,  the rate of convergence of D^κ_p(u_i)^n_i = 1 can be understood as a measure of efficiency of a low-discrepancy sequence u_i.

In our application we use Sobol sequence with grey code implementation. Sobol sequence were introduced by Sobol <cit.> and the grey code implementation is due <cit.>. 
Detailed investigation of the nature of this sequence goes beyond the scope of this work. However, any other low-discrepancy sequence can be used.

For a Sobol sequence value D^κ_2(u_i)^n_i = 1  converges to zero as O(n^-1 + ϵ), 
where ϵ depends on d implicitly. While convergence rate of discrepancy  for a random uniform sequence  is O(1/√(n)), for d ≪ 200 value  of ϵ < 1/2, 
which makes low-discrepancy sequence rather efficient.  
However, scrambling and effective dimension techniques can increase rate of convergence <cit.>.

It can be established that the sequences u and X have no repeating values almost surely, 
As our data is assumed to come from atomless distributions. 
This means that the problem of finding Permutation σ^* in (<ref>) is the optimal assignment problem. 
The Optimal assignment problem  is the special case the linear programming, 
which can be solved in O(n^3) time by the Hungarian algorithm <cit.><cit.>. 
Amortizations based on applications of parallel programming can improve computation complexity to O(n^2 log n). 
Approximate algorithms can be used for faster computations, for example <cit.>.

Note, that even if sample X has an i.i.d distribution, then the  resulting transport Y is not  independent itself. 
However, the covariance of the elements is converging to zero as T goes to infinity. 
Let n ≠ m be two indices less or equal  to T and i,j be coordinate indices. As (X_k)^T_k = 1 is assumed to be i.i.d, it follows that Y^i_n has a discrete uniform distribution over { u^i_k }^T_k = 1. 
Then, 

    (Y^i_n,Y^j_m) =  Y^i_n Y^j_m -  Y^i_n  Y^j_m  = ∑^T_k = 1∑^T_ l ≠ ku^i_k u^j_l/T(T - 1) - ∑^T_k = 1∑^T_l = 1 u_k^i u_l^j /T^2 =
    
     =
    ∑^T_k = 1∑^T_ l ≠ ku^i_k u^j_l/T^2(T - 1) - ∑^T_k = 1 u^i_k u^j_k/T^2,

so by bounding from above and below with equidistribution property of u in the limit case

    (Y^i_n,Y^j_m) ≤∑^T_k = 1∑^T_l = 1u^i_k u^i_l/T^2(T - 1)lim_T →∞ U  U/T - 1  = 
    lim_T →∞ 1/4(T - 1)  = 0,


    (Y^i_n,Y^j_m) ≥  - ∑^T_k = 1 u^i_k u^j_k/T^2≥    - ∑^T_k = 1 u^i_k /T^2lim_T →∞ -  U /T   = 
    lim_T →∞ -  1/2T  = 0;

 where U is a uniform random variable on [0,1]. Thus, the value (Y^i_n,Y^j_m) converges to zero. As variance of Y^i_n converges to the variance of a standard uniform distribution on [0,1]  
 we will treat sample Y as uncorrelated in asymptotic context under assumption of i.i.d. distribution of X.

 It can be easily seen that that quadratic discrepancy admits a degenerate  V-statistic representation  with the kernel :
 
    ( D^κ_2(y) )^2 =  1/n^2∑^n_i,j = 1(y_i,y_j),

 assuming sample y of n elements. Properties of V-statistic produces the  asymptotic result 
 
    n ( D^κ_2(y) )^2   V = ∑^∞_i = 1λ_i Z_i^2 ∼(λ),

 where  Z (0,1) and λ are non-zero eigenvalues of the integral operator 𝒜 defined by the relation
 
    𝒜(f)(x) = ∫^1_0 f(y)(x,y) dy

 It can be shown that 𝒜 is in fact positive-semidefinite  and  trace-class, which means that all λ_i > 0 and that
 
    ∑^∞_n = 1λ_i < ∞

 Eigenvalues of 𝒜 can be approximated by a  finite collection of N numbers λ̂ with Nyström method.
 As it was shown in <cit.> the twofold approximation of cdf of V_|N = ∑^N_i = 1λ_i^N Z_i^2 is possible for fixed natural numbers N, K and parameter α∈ (0,1]:

    (V_|N < x) ≈
    ≈1/2 - ∑^K_k = 0sin(
     1/2∑^N_1 = 1arctan(2(k + 12) αλ^N_i) - (k +  1/2) α x  ) /( k + 1/2) ∏^N_i = 1√(1 +  (2(k + 1/2) αλ^N_i )^2 ).

This formula can be used for computing quantiles  and critical values of (λ) with simple one-dimensional zero-seeking algorithm.


    Case of A Single Change Point 

In this chapter an approach for detecting a single change point is presented. We impose a model assumptions that for a fixed τ change point  either satisfy τ <    < T - τ or it does not exist .
This value τ can be selected in a such way that 1 ≪τ≪ T  and be used as a sliding window bandwidth defining  a 'control chart'-like object which we refer to as diphoragram.

 empirical sliding diphoragram for change point data (X_t)^T_i = 1 is defined by

    Δ_t^T,τ = ( D^κ_2(R̂_T(X_i))^t + τ_i = t)^2 
     =  ( D^κ_2(Y_i)^t + τ_i = t)^2

and ideal  sliding diphoragram by

    Δ_t^T,τ = ( D^κ_2( R_μ_r(X_i))^t + τ_i = t)^2

 for t in range from 1 to T'_τ = T - τ.
 
 If κ is a continuous function, then  tt as T →∞.
 With this condition discrepancy is a continuous function of data. So, convergence of vector ranks proved in <cit.> implies convergence of diphoragrams.
 Note, that with kernel representation  charts admit a difference representation

    t + 1 -t =  1/τ^2(∑^t+1 + τ_i = t + 1(Y_i,Y_t + τ + 1)
     -  ∑^t+ τ_i = t (Y_i,Y_t) ).

From the computational standpoint this means that computation of the whole time-series  and  takes only quadratic time O(dT^2) in number in observation. 
Moreover, if crisp optimal transport with low-discrepancy sequence u is used in estimation of vector rank function, then A_Y = A_u, so all values can be precomputed. 

For  the application to the change point problem  consider two increasing sequences of integers  T_n and τ_n such that T_n/τ_n =  a > 1 for all n ∈ are constructed.  
Then as Y_i are assumed to be independent t and t + τ_n are also independent random variables. 

mean sliding discrepancy for set sample (Y_i)_i = 1^T_n is computed by

    = τ_n^2/T_n∑^a - 1_j = 01 + jτ_n


Note, that the mean sliding discrepancy is not the same  as the discrepancy of the whole sample. By independence, in case H_0 holds, as n →∞ 

    ∑^a_j = 1∑^∞_i = 1λ_i/a Z_i,j^2 ∼( ( λ_i/a)^a_j = 1)^∞_i = 1,

where Z (0,1). Otherwise, there will be a sliding discrepancies t sampled form the non-uniform data which goes to ∞ in probability as n →∞. 
So, the whole sum ∞. 

In case  the single change point exists the statistic  is expected to attain the minimal value at such moment of time t^* when the  empirical distribution  of (Y_t)^t^* + τ_n_t = t^* 
approaches the empirical distribution of the whole sample, as the empirical distribution of the whole sample is converging to the . Let 
       
    t^*_n = _1 ≤ t ≤ T_n - τ_n t.

It is expected that  the ratio of numbers of elements from both sides of change point  in the subsample used in the computation of t^*_n 
and in the whole sample are equal. This can be represented in the algebraic relation

    _n - t^*_n  /τ_n  = _n/T_n = r̂_n,

where θ̃_n is the produced estimate of the change point. This provides the expression for the estimate

    θ̃_n =  t^*_n/1 - a^-1.

We accept H_0 for a fixed γ-level in case  

    p_n = ( V_|N≤) <   1 - γ,
  
where the cdf is numerically estimated as in (<ref>) for some fixed parameters N and K. Otherwise reject H_0 and state that there was a change most probable at _n. 

In order to reason about properties of the estimate _n. Let ℳ(I^d) denote  space of finite signed measures other the hypercube I^d.

Space of nonuniformities 𝒟 is  defined as a quotient of the real vector spaces

    𝒟 = ℳ(I^d)/

endowed with a Hilbert space structure by inner product defined for [ν], [μ] ∈𝒟 by 

    ⟨ [ν], [μ] ⟩ =  ∫∫(x,y)   dν(x)  dμ(y)


Note, that the relation (<ref>) is well defined as for all a,b ∈

    ⟨[ν + a ], [μ + b] ⟩  =
     ∫∫(x,y)   dν(x)  dμ(y) +
     b∫∫(x,y)   dν(x)  d y + 
     + 
     a ∫∫(x,y)   d x  dμ(y) +
     ab ∫∫(x,y)   d x  d y 
     = ∫∫(x,y)   dν(x)  dμ(y),

as ∫(x,y)  d x = 0 for any value of y and  is symmetric and is indeed an inner product as  is a positive-definite kernel. 
As the line  intersects simplex of probability measures 𝒫(I^d) only in one point ( itself) the natural projection μ↦ [μ] is injective on 𝒫(I^d), 
so we denote a nonuniformity arising from each μ∈𝒫(I^d) just as μ.   
If every subsample (Y_i)_i = t^t + τ is associated with an empirical measures μ̂_t = ∑_i = t^t + τ_n δ_Y_i, then 
    t = μ̂_t  ^2_𝒟 = ⟨μ̂_t,  μ̂_t ⟩.
 
By construction of the vector rank function r R_# P_1 + (1 - r) R_# P_2 _𝒟 = 0, hence r R_# P_1 + (1 - r) R_# P_2 = 0 in 𝒟.  This produces result

    R_# P_2 = -r/1 -r R_# P_1 ⟹   R_# P_2  _𝒟^2  = r^2/(1 - r)^2 R_# P_1 ^2_𝒟⟹
    ⟹
    	r =  R_# P_2 _𝒟/ R_# P_1 _𝒟 +  R_# P_2 _𝒟

Considering that r = 1/2 it follows that [R_# P_1]  =  - [R_# P_2], so magnitude of the discrepancy will have the same distribution for sample with equal proportions of elements with distributions 
of R_# P_1 and R_# P_2. If t^* = ( 1 - a^-1 ), then it can be shown that estimate  is unbiased: 
 

    _n = ∑^T_n - τ_n_t = 1 t (t = _t' t')1 - a^-1_n
    = 
     t^* +  B^+_n
         -  B^-_n  
      /1 - a^-1_n 
       = t^*/1 - a^-1_n =   = T_n/2,

 where

    B^-_n = ∑^⌊ r(T_n - τ_n) ⌋_t = 1 t ( t^*  - t = _t't'),   
    
     B^+_n = ∑^⌊ (1 - r)(T_n - τ_n) ⌋_t = 1  t (t^*  + t = _t' t').

However, then r ≠ 1/2 the estimate r̂ is projected to  be biased towards 1/2 as the value 
B^+ will only increase and the  value B^- will only decrease  as r decreases. Otherwise, increase  of r 
increases the value B^-  too, and  decreases the value of B^+. 
In order to proof consistency of estimator r̂_n we introduce a family of sequences for each s ∈ ((1 -a^-1)^-1, 1)

    Δ̃^s_n =  t'  having   t' = _1 ≤ t ≤ T_n - τ_n| s - t/T_n(1 - a^-1)|

Then by the weak convergence of vector ranks for every s the sequence converges to a value:

    Δ̃^s_n Δ^s   = 
       1_ s < r - a^-1/1 - a^-1(s)R_# P_1 +   1_ s > r  /1 - a^-1(s)R_# P_2 +   
     
     +  
     1_r - a^-1/1 - a^-1≤ s ≤r /1 - a^-1(s) 
     (  r   - s(1 - a^-1)/a^-1 R_# P_1 +      a^-1 - r   + s(1 - a^-1)/a^-1 R_# P_2 ) _𝒟 .

Otherwise Δ̃^r_n  0.  We can assert by  structural uniformity of Δ̃_n that 
Δ̃Δ in Skorohod's topology.  Thus, r̂_n = _s Δ_n^s _s  Δ^s  = r providing convergence in probability.
This suggests that the bias of r̂_n can be bounded by some sequence β with convergence β_n  0:

    |B^+_n - B^-_n|/T_n≤β_n

Without loss of generality assume that r < 1/2. Then, we separate positive bias into mixing and non-mixing parts:

    |B^+_n - B^-_n|/T_n <  B^+_n/T_n < 
    	
     <
    	∑_t = 1^⌊(1-r)τ_n ⌋ t(t' < Δ̃^r_n )/T_n  +  ∑_t = ⌊(1-r)τ_n ⌋ + 1^T_n - τ_n - t^* t(t' < Δ̃^r_n ) /T_n.

For non-mixing part apply Markov inequality for each fixed value of Δ̃^r_n and some value ζ > 1

    (t' < Δ̃^r_n ) = ( (t')^-ζ >  (Δ̃^r_n )^- ζ) ≤( Δ̃^r_n )^ζ( (t')^-ζ).

It is possible to use inverse as Δ̃_n > 0 for each n. With this inequality  the non-mixing part can be bounded

    ∑_t = ⌊(1-r)τ_n ⌋ + 1^T_n - τ_n - t^* t(t' < Δ̃^r_n ) /T_n≤ T_n (  ( Δ̃^r_n )^ζ) ( ( Δ̃^1_n )^-ζ) 
     ≤
    ≤ C_r T_n (  ( Δ̃^r_n )^ζ).

As Δ̃^r_n → 0 with the rate O(τ_n^-1 + ε) it is  possible to specify a constant ζ in such a way that 
T_n  (  ( Δ̃^r_n )^ζ ) ↓ 0 and the expectation  ( ( Δ̃^1_n )^-ζ ) approaches some constant value by weak convergence and therefore can be bounded. 
Thus, the non-mixing bias approaches zero as n goes to infinity. Note, that all moments of discrepancy exists as it is bounded on a positive interval for each n.
This method can be extended in order to show that the estimate is asymptotically unbiased  under an assumption of change point slackness,

Sequence of diphoragrams t has slackness property iff for some constant γ∈ (0,a^-1)  and  for all n large enough there are  time points t^**_n > t^*_n such that :

    t^**_n - t^*_n/T_n≥γ,   ∑_t = 1^t^**_n t( t = _t't') - B^-_n ≤ 0.

Then, it is possible to construct a similar bound

    | B^+_n - B^-_n|/T_n≤   T_n (  ( Δ̃^r_n )^λ( ( Δ̃^γ + r_n )^-λ| Δ̃^r_n  ) ),

which converges to zero at infinity.




    alternative methods

One of the negative properties of the method described in the previous chapter is the requirement of specification of bandwidth τ, which prevents change point detection in the proximity of  the limit points 1  and T. 
However, as space of nonuniformities  𝒟  is a metric space, it is possible to determine change point by maximizing distance between two empirical measures 
dist() = μ̂^+_ - μ̂^-_^2_𝒟, where 

    μ̂^-_  = 1/∑^_t = 1δ_Y_t, 
      μ̂^+_  = 1/T - ∑^T_t = + 1δ_Y_t .

Then, by definition of the norm the distance is computed as
 
    dist() =  
      1/()^2∑^_n,m = 1𝒦(Y_n,Y_m) - 
     
       -2/(T - )∑^_n = 1∑^T  _m = + 1𝒦(Y_n,Y_m) 
       +  1/( T - )^2∑^ T _n,m = + 1𝒦(Y_n,Y_m),

leading to a change point estimate = _dist(). As empirical measures of (<ref>) will converge to some points of interval [R_# P_1, R_# P_2] as T_n goes to infinity, the estimate _n converges to true value  in probability.

Change points detection methods of this forms were explored in the work <cit.>. Thus, we will not explore it in further depth. The important properties of this statistic is that it still can be computed in O(T^2) time and that dist() is a U-statistic.

 Geometric idea of (<ref>) suggests that the value
 
    ς() =  /T   - μ̂^+__𝒟/μ̂^-__𝒟 + μ̂^+__𝒟

 approaches zero as data size grows to infinity  iff r̂ approaches the true ratio  r. Thus, if change point exists, then change point can be estimated as = _ |ς()|.
 Or alternatively compute r̂^m by applying iteration
 
    r̂^m =   μ̂^+_r^m -1_𝒟/μ̂^-_r^m-1_𝒟 + μ̂^+_r^m - 1_𝒟,

 where μ^+_r^m - 1 and μ^-_r^m - 1 are empirical measures corresponding to the estimates obtained at previous iterations. The initial value r̂^0 can be selected to be equal to 1/2.
 It is interesting, that true location of change point can be estimated by the value of the error ς(). 
 
 Assume, that >, then the value μ^-__𝒟 should not distinct dramatically from the value | μ^-__𝒟. However,  it holds that   
 | μ^+__𝒟 < | μ^+__𝒟 with high probability as μ^+_ is closer to the uniform distribution than μ^+_. 
 Thus, 0 ≈ς() > ς() with high probability.
 
 Otherwise, then < then  the statistic ς() grows similarly to the function f(x) = Cg(x) -  x/c + x, where g is monotonically decreases, as   grows to T.
 Thus 0 ≈ς() < ς() with  hight probability.
 
This suggests, that it must be possible to establish an iterative procedure for finding optimal estimate .
Note, that in case  change point does not exist the values of |ς()| will be low  for all possible values of . So, in case change point exists it is possible to establish an algorithm:





  
Multiple Change Points 
In this chapter the situation of K possible consecutive  change points = (^1,…,^K) in the data is considered. Now, r denotes a list of positive values
   
    r =  (r^1, …, r^K) = ( ^1/T, ^2 - ^1/T,…,^K - ∑^K-1_k = 1^k/T),

 which can be understood as the first K coefficients in the convex combination  of K + 1 probability distributions P_1, …, P_K + 1. That is, define

    μ_r = ∑^K_k = 1 r^k P_k  + ( 1 -  ∑^K_k = 1 r^k)P_K + 1.

Furthermore, estimates  and r̂ are treated as lists of corresponding structure. 

By definition of vector rank function R_#μ_r =.  However, in order to apply methods similar to ones developed in the previous chapter we need one more property of the model.  

K probability measures P_1, … P_K are said to be convexly independent if for all lists of K coefficient λ∈^K_+, such that ∑^K_k =1λ_k = 1, for each i equality 
P_i = ∑^K_k = 1λ_k P_k,
 holds only if λ_j = δ_i,j for each j, where δ_i,j is the Kronecker delta.

Assume that the true value of K ∈𝔎. Then, by construction K^SMA≤ K with probability α.

 It can be postulated, that 0 of 𝒟  lies in the convex hull of R_# P_i. That is  
0 ∈conv{  R_# P_k  }^K + 1_k = 1.
It suggests that there are projections of 0 to the edges of the convex polytope π_k0 ∈ [R_# P_k, R_# P_k + 1] which minimizes π_k 0   _𝒟. Hence, if τ is taken to be small enough then local minimas of Δ^T,τ will happen in the proximity of the change points with high probability. 

The problem with this method is that proportions of points from different sides of  a change point at the local minimise
are given by the relation

    π_k 0 =   R_# P_k + 1^2_𝒟 - ⟨ R_# P_k,R_# P_k + 1⟩_𝒟/d^2_𝒟(R_# P_k,R_# P_k + 1)R_# P_k + 
    	 
     + 
    	 R_# P_k^2_𝒟 
    	- ⟨ R_# P_k,R_# P_k + 1⟩_𝒟/d^2_𝒟(R_# P_k,R_# P_k + 1)R_# P_k + 1  ,

which can not be recovered from the diphoragram. For this reason we propose an iterative procedure. Let  μ̂_^k,^k + 1^n denote empirical measures of points from the interval bounded by change point estimates [^k]^n and [^k + 1]^n .


  * start with 𝒯_1 = 𝒯.

  * for each k ≤ K select t^*_k = _t ∈𝒯_kΔ̃^T, τ_t and update T_k + 1 = { t ∈𝒯_k : | t - t^*_k | > τ}.

  * make initial change point estimation with blind adjustment [^k  ]_1= t^*_k + τ/2.

  * readjust change points for a fixed number of iterations N  with nth plus one readjustment being 
 
    [^k  ]_n + 1 =  t^*_k + [λ̂^k_1 ]_n τ,

 where
 
    [λ̂^k_1 ]_n =  μ^n__k , _k + 1^2_𝒟 - ⟨μ^n__k -1 , _k ,μ^n__k , _k + 1⟩_𝒟/d^2_𝒟(μ^n__k - 1 , _k,μ^n__k , _k + 1)

 with surrogate change points being [ ^0]_n-1 = 0 and  [ ^K + 1]_n - 1 = T .

In order to approach problem of model misspecification we apply smallest accepted model (SMA) methodology. 
For a collection of K we acquire minimizing time points t^*. Then, test for change points in the intervals bounded by _i and _i + 1 with surrogate values as above by
computing a p-value approximations.

Thus, the model can be estimated while new local minimizers are being recovered and the process can be stopped as minimal accepted value of K̂ has been achieved. 


  
Computational Results
In order to conduct computational experiments the  methods discussed in the previous chapter were implemented in Python programming language with use of numpy and scipy libraries. 
We use   package in order to  generate Sobol sequences.


    Simulations with Zero Change Points

 
For simulations with zero change points we are interested in measuring statistical significance or confidence of the test statistic T which can be understood as 

    significance(T) = ( T  rejects  H_0 | H_0  is true),   confidence(T) = (T  accepts  H_0 | H_0  is true).

In simulations with zero change points H_0 obviously is true. So, for a run of n simulations we estimate confidence as 

    conf(T,n) = #{simulations with no change points detected}/n.

We run simulations without change points and vary certain fixed parameter while measuring confidence and inverse p-value for each value of parameter. The only nontrivial results were obtained for change in data dimension d. 
 
Additional experiments were conducted with growing variance and changes in covariance structure of the observations, however no dependencies were identified. This could be due to stabilizing effect of vector rank functions.



  
Simulations with One Change Point
 
  While running a simulations with one change point we are naturally interested in the estimation of statistical power of the test T which can be understood as
  
    power(T) = ( T  rejects  H_0 | H_0  is false),

and can be estimated for n simulations as 

    pow(T,n) =  #{simulations with change points detected}/n.

We investigate dependence of pow(T,n) on differences between distributions from opposite sides of the change point.
 
 
 
 The figures shows that star discrepancy outperforms symmetric discrepancy  in detecting change both in expectation and in variance. However empirical results in <cit.> indicates that in some situations symmetric discrepancy may turn out tob be a better tool.

In the situation of existing change point not only power of the test is of interest but also a precision of change point estimations. 
As it was shown in previous chapter bias of a change point estimate increases as true ratio r of  distributions  in the sample diverges from the value of 1/2. 
We provide a computational illustrations:

 
  


    Simulations with Multiple Change Points

For case of multiple change point we provide only example with diphoragrams:
  


  
Examples with financial data:
In order to provide examples, which are not artificially constructed, financial data similar to one in <cit.>. 
This data was acquired from http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.htmlData library  of  Keneth R. French. It contains mean monthly returns from five portfolios each composed of one of  five industrial sectors in the USA, which are:

	
  * Finance,
	
  * Manufacturing,
	
  * Retail, wholesale  and some services,
	
  * Utilities,
	
  * Other.

This provides data dimension of d = 5 and total number of  observations of T = 1091 as they were recorded monthly from July of the year 1925 to the may of the year 2017.   

In the original paper <cit.> the parametric Bayesian inference was used too estimate nine and seven change points, and the total number of change point was only guessed and not inferred. 
Our results differ significantly from this original results, however different time range was used.      

It can be projected that all above change points can be attributed to the important events in the economic history. For example, change point at July 1954 can be related to the end of recession of 1953, 
which itself can be explained by the change of interrelations  of the industrial sectors leading to the change in the structure of the observed distribution. 
Hence, we can describe performance of the SMA-procedure as satisfactory.



  
Conclusion and Discussion

As result of the work a collection of methods of change point detection methods was developed. All this methods are based on interaction of vector ranks and geometric discrepancy which is a novel result. Certain basic consistency result were proved for the method in its basic form designed for detecting a single change point. However, they also can be applied for detecting and estimating multiple change points. 

Computational results shows applicability of the method both for simple artificial change point problems an problems concerning actual data from applications. 
It is strictly indicated by resulting experience that the method is much more potent in situations then the distributions are not concentric. 
Empirical  evidences of our theoretical findings were also observed. 
These theoretical results include expression of relation on different sides of change point through inner product in the Hilbert space 𝒟. 
Another Theoretical result concerns dependence of estimate's bias on ratio in which change point separates sample. 

Another positive results is the discovery of 𝒟, which can be used for establishing alternative iterative  procedures defined by relations in this Hilbert space. Furthermore, this Hilbert space structure on empirical measures can be used for proving more theoretical results in the future.  

Better proofs which cover ergodic processes and provide exact rates of convergence a still need to be worked out for the methods. 
Furthermore concentration results for change point estimates might render this algorithms interesting for practical application. 
However, they are absent in the current work.


Finally, online version of method can be implemented. The only requirement for this algorithm is fast online computation of the optimal assignment problem. 
It can be projected that such algorithm can be  derived from the Sinkhorn distance regularised algorithm which were designed by Cuturi <cit.>. 
