タイムワープ編集距離

時系列のデータ分析において、タイム ワープ編集距離(TWED) は、離散時系列のペア間の類似度 (または非類似度) の尺度であり、弾性という物理的な概念を使用して 2 つの系列の時間単位の相対的な歪みを制御します。他の距離尺度 (DTW (動的時間ワーピング) や LCS (最長共通部分列問題) など)と比較すると、TWED はメトリックです。計算時間の複雑さは ですが、コリドーを使用して検索空間を減らすことで、特定の状況で大幅に削減できます。メモリ空間の複雑さはまで削減できます。これは、2009 年に P.-F. Marteau によって初めて提案されました。 ${\displaystyle O(n^{2})}$ ${\displaystyle O(n)}$

${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{p},B_{1}^{q})=Min{\begin{cases}\delta _{\lambda ,\nu }(A_{1}^{p-1},B_{1}^{q})+\Gamma (a_{p}^{'}\to \Lambda )&{\rm {delete\ in\ A}}\\\delta _{\lambda ,\nu }(A_{1}^{p-1},B_{1}^{q-1})+\Gamma (a_{p}^{'}\to b_{q}^{'})&{\rm {match\ or\ replacement}}\\\delta _{\lambda ,\nu }(A_{1}^{p},B_{1}^{q-1})+\Gamma (\Lambda \to b_{q}^{'})&{\rm {delete\ in\ B}}\end{cases}}}$
一方

${\displaystyle \Gamma (\alpha _{p}^{'}\to \Lambda )=d_{LP}(a_{p}^{'},a_{p-1}^{'})+\nu \cdot (t_{a_{p}}-t_{a_{p-1}})+\lambda }$
${\displaystyle \Gamma (\alpha _{p}^{'}\to b_{q}^{'})=d_{LP}(a_{p}^{'},b_{q}^{'})+d_{LP}(a_{p-1}^{'},b_{q-1}^{'})+\nu \cdot (|t_{a_{p}}-t_{b_{q}}|+|t_{a_{p-1}}-t_{b_{q-1}}|)}$
${\displaystyle \Gamma (\Lambda \to b_{q}^{'})=d_{LP}(b_{p}^{'},b_{p-1}^{'})+\nu \cdot (t_{b_{q}}-t_{b_{q-1}})+\lambda }$

再帰は 次 の ように初期化されます。${\displaystyle \delta _{\lambda ,\nu }}$
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{0},B_{1}^{0})=0,}$
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{0},B_{1}^{j})=\infty \ {\rm {{for\ }j\geq 1}}}$
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{i},B_{1}^{0})=\infty \ {\rm {{for\ }i\geq 1}}}$
${\displaystyle a'_{0}=b'_{0}=0}$

Pythonラッパーを使用したCでのTWEDアルゴリズムの実装は^[1]で入手可能である。

TWEDは、[1]で入手可能なTime Series Subsequence Search Pythonパッケージ（略してTSSEARCH）にも実装されています。

TWEDのR実装は、状態やイベントのシーケンス、より一般的には離散シーケンスデータのマイニング、記述、視覚化を行うRパッケージであるTraMineRに統合されています。^[2]

さらに、cuTWED は、G. Wright (2020) による改良アルゴリズムを使用する TWED のCUDAアクセラレーション実装です。この方法はメモリ内で線形であり、大規模に並列化されています。cuTWED は CUDA C/ C++で記述されており、Python バインディングが付属しており、Marteau のリファレンス C 実装用の Python バインディングも含まれています。

 numpyを npとして インポートする


def  dlp ( A ,  B ,  p = 2 ): 
    cost  =  np . sum ( np . power ( np . abs ( A  -  B ),  p )) 
    return  np . power ( cost ,  1  /  p )


def  twed ( A 、 timeSA 、 B 、 timeSB 、 nu 、 _lambda ): 
"""指定された時系列 A と B のタイム ワープ編集距離 (TWED) を計算します。""" # [distance, DP] = TWED(A、 timeSA、 B、 timeSB、 lambda、 nu) # # A := 時系列 A (例: [ 10 2 30 4]) # timeSA := 時系列 A のタイムスタンプ (例: 1:4) # B := 時系列 B # timeSB := 時系列 B のタイムスタンプ# lambda := 削除操作のペナルティ# nu := 弾性パラメーター - 距離測定には nu >=0 が必要# 参照: # Marteau、P.、F. (2009)。「タイム ワープ編集距離と剛性調整による時系列マッチング」。# IEEE Transactions on Pattern Analysis and Machine Intelligence。 31 (2): 306–318. arXiv:cs/0703033 # http://people.irisa.fr/Pierre-Francois.Marteau/    
    
    
    
    
    
    
    
    
    
    
    
    

    # 入力引数をチェック
    if  len ( A )  !=  len ( timeSA ): 
        print ( "A の長さは timeSA の長さと等しくありません" ) 
        return  None ,  None

    if  len ( B )  !=  len ( timeSB ): 
        print ( "B の長さは timeSB の長さと等しくありません" ) 
        return  None ,  None

     nu  <  0の場合: 
        print ( "nu は負です" ) 
        return  None ,  None

    # パディングを追加
    A  =  np . array ([ 0 ]  +  list ( A )) 
    timeSA  =  np . array ([ 0 ]  +  list ( timeSA )) 
    B  =  np . array ([ 0 ]  +  list ( B )) 
    timeSB  =  np . array ([ 0 ]  +  list ( timeSB ))

    n  =  len ( A ) 
    m  =  len ( B ) 
    # 動的計画法
    DP  =  np . zeros (( n ,  m ))

    # DP行列を初期化し、最初の行と列を無限大に設定します。
    DP [ 0 ,  :]  =  np . inf 
    DP [:,  0 ]  =  np . inf 
    DP [ 0 ,  0 ]  =  0

    # 最小コストを計算します
    for  i  in  range ( 1 ,  n ): 
        for  j  in  range ( 1 ,  m ): 
            # さまざまな操作のコストを計算して保存します
            C  =  np . ones (( 3 ,  1 ))  *  np . inf 
            # A での削除
            C [ 0 ]  =  ( 
                DP [ i  -  1 ,  j ] 
                +  dlp ( A [ i  -  1 ],  A [ i ]) 
                +  nu  *  ( timeSA [ i ]  -  timeSA [ i  -  1 ]) 
                +  _lambda 
            ) 
            # B での削除
            C [ 1 ]  =  ( 
                DP [ i ,  j  -  1 ] 
                +  dlp ( B [ j  -  1 ],  B [ j ]) 
                +  nu  *  ( timeSB [ j ]  -  timeSB [ j  -  1 ]) 
                +  _lambda 
            ) 
            # 両方の時系列でデータ ポイントを保持
            C [ 2 ]  =  ( 
                DP [ i  -  1 ,  j  -  1 ] 
                +  dlp ( A [ i ],  B [ j ]) 
                +  dlp ( A [ i  -  1 ],  B [ j  -  1 ]) 
                +  nu  *  ( abs ( timeSA [ i ]  -  timeSB [ j ])  +  abs (時間SA [ i  -  1 ]  - 時間SB [ j -  1 ])) 
            ) 
            # 最小コストの操作を選択し、DP行列を更新します。DP 
            [ i , j  ] =  np  . min ( C ) distance 
    =  DP  [ n -  1  , m  -  1  ]を
    返します。 距離、 DP

最もコスト効率の高いパスを見つけるためにバックトラッキングします。

def  backtracking ( DP ): 
"""最もコスト効率の良いパスを計算します。""" # [ best_path ] = BACKTRACKING (DP) # DP := TWED 関数の DP 行列    
    
    

    x  =  np .形状( DP ) 
    i  =  x [ 0 ]  -  1 
    j  =  x [ 1 ]  -  1

    # パスのインデックスは反対方向に保存されます
    # path = np.ones((i + j, 2 )) * np.inf; 
    best_path  =  []

    steps  =  0だが
    i != 0またはj != 0の場合 : best_path . append (( i - 1 , j - 1 ))      
             

        C  =  np . 1 (( 3 ,  1 ))  *  np . inf

        # 両方の時系列のデータポイントを保持
        C [ 0 ]  =  DP [ i  -  1 ,  j  -  1 ] 
        # A で削除
        C [ 1 ]  =  DP [ i  -  1 ,  j ] 
        # B で削除
        C [ 2 ]  =  DP [ i ,  j  -  1 ]

        # コストが最も低いインデックスを見つける
        idx  =  np . argmin ( C )

        if  idx  ==  0 : 
            # 両方の時系列のデータポイントを保持
            i  =  i  -  1 
            j  =  j  -  1 
        elif  idx  ==  1 : 
            # A で削除
            i  =  i  -  1 
            j  =  j 
        else : 
            # B で削除
            i  =  i 
            j  =  j  -  1 
        steps  =  steps  +  1

    best_path . append (( i  -  1 ,  j  -  1 ))

    best_path.reverse () はbest_path [ 1 : ]を返します
     

function [distance, DP] = twed ( A, timeSA, B, timeSB, lambda, nu ) % [distance, DP] = TWED( A, timeSA, B, timeSB, lambda, nu ) % 指定された時系列 A と B のタイム ワープ編集距離 (TWED) を計算します% % A := 時系列 A (例: [ 10 2 30 4]) % timeSA := 時系列 A のタイムスタンプ (例: 1:4) % B := 時系列 B % timeSB := 時系列 B のタイムスタンプ% lambda := 削除操作のペナルティ% nu := 弾性パラメーター - 距離測定には nu >=0 が必要です% % コード作成者: P.-F. Marteau - http://people.irisa.fr/Pierre-Francois.Marteau/ % 入力引数をチェックif length ( A ) ~= length ( timeSA ) warning ( 'A の長さが、timeSA の長さと等しくありません' ) return end if length ( B ) ~= length ( timeSB ) warning ( 'B の長さが、timeSB の長さと等しくありません' ) return end if nu < 0 warning ( 'nu が負です' ) return end % パディングを追加A = [ 0 A ]; timeSA = [ 0 timeSA ]; B = [ 0 B ]; timeSB = [ 0 timeSB ]; % 動的計画法DP = zeros ( length ( A ), length ( B )); % DP 行列を初期化し、最初の行と列を無限大に設定DP ( 1 , :) = inf ; DP (:, 1 ) = inf ; DP ( 1 , 1 ) = 0 ; n = length ( timeSA ); m = length ( timeSB ); % i = 2 : nの場合、j = 2 : mの場合の最小コストを計算する 
    
    
    
    
    
    
    
    
    
    
    
 
    
       
        
        
    
 
       
        
        
    
 
       
        
        
    
    
       
       
       
       
 
    
       
 
    
       
       
       
 
      
      
    
       
           
             = Dlp ( A ( i ), B ( j )); % さまざまな操作のコストを計算して保存しますC = ones ( 3 , 1 ) * inf ; % A での削除C ( 1 ) = DP ( i - 1 , j ) + Dlp ( A ( i - 1 ), A ( i )) + nu * ( timeSA ( i ) - timeSA ( i - 1 )) + lambda ; % B での削除C ( 2 ) = DP ( i , j - 1 ) + Dlp ( B ( j - 1 ), B ( j )) + nu * ( timeSB ( j ) - timeSB ( j - 1 )) + lambda ; % 両方の時系列のデータポイントを保持するC ( 3 ) = DP ( i - 1 , j - 1 ) + Dlp ( A ( i ), B ( j )) + Dlp ( A ( i - 1 ), B ( j - 1 )) + ... nu * ( abs ( timeSA ( i ) - timeSB ( j )) + abs ( timeSA ( i - 1 ) - timeSB ( j - 1 ))); % コストが最小となる操作を選択し、DP 行列を更新するDP ( i , j ) =  
         
            
                 
         
            
                                
            
                                
            
                               
                        
         
            
               min ( C ); 
end end distance = DP ( n , m ); % ユーク​​リッド距離を計算する関数function [cost] = Dlp ( A, B ) cost = sqrt ( sum ( ( A - B ) .^ 2 , 2 ) ); end end        
    
 
       
 
    
       
               
    
 

最もコスト効率の高いパスを見つけるためにバックトラッキングします。

function  [path] = backtracking ( DP ) % [ path ] = BACKTRACKING ( DP ) % 最もコスト効率の高いパスを計算します% DP := TWED 関数の DP 行列x = size ( DP ); i = x ( 1 ); j = x ( 2 ); % パスのインデックスは反対方向で保存されますpath = ones ( i + j , 2 ) * Inf ; steps = 1 ; while ( i ~= 1 || j ~= 1 ) path ( steps , :) = [ i ; j ]; C = ones ( 3 , 1 ) * inf ; % 両方の時系列でデータ ポイントを保持しますC ( 1 ) = DP ( i - 1 , j - 1 ); % A で削除C ( 2 ) = DP ( i - 1 , j ); % B での削除C ( 3 ) = DP ( i , j - 1 ); % コストが最小となるインデックスを検索[ ~ , idx ] = min ( C ); switch idx case 1 % 両方の時系列でデータ ポイントを保持i = i - 1 ; j = j - 1 ; case 2 % A での削除i = i - 1 ; j = j ; case 3 % B での削除i = i ; j = j - 1 ; end steps = steps + 1 ;  
    
    
    
 
      
      
      
 
    
           
 
      
           
            
     
             
     
        
               
        
             
        
             
     
        
           
     
         
             
                
                    
                    
             
                
                    
                  
             
                
                  
                    
        
            
    終了
path ( steps , :) = [ i j ]; % パスは逆方向に計算されました。path = path ( 1 : steps , :); path = path ( end : - 1 : 1 , :);終了        
 
    
       
         
 

• Marteau, P.; F. (2009). 「時系列マッチングのための剛性調整によるタイムワープ編集距離」. IEEE Transactions on Pattern Analysis and Machine Intelligence . 31 (2): 306–318. arXiv : cs/0703033 . doi :10.1109/TPAMI.2008.76. PMID  19110495. S2CID  10049446.