Why can't this implementation get the results of LTS papers?

[tcoln]

I ran this code, but it couldn't get the results in the LTS paper.

[GillesVandewiele]

Which dataset did you test it on, and did you use the same hyper-parameter configuration as in the paper, with a batch_size of 1 and SGD as optimizer?

[tcoln]

Yes, of course, I test the first n datasets in LTS paper， and I have set the same hyper parameters as in the paper.

[GillesVandewiele]

So you also explicitly set the batch_size to 1? Which is not the default value of the ShapeletModel and is not mentioned in the paper. Perhaps be a bit more concrete and show a code snippet, else it will be difficult to help you. I managed to reconstruct the results of the LTS paper pretty well.

   clf = ShapeletModel(n_shapelets_per_size=shapelet_dict, 
                                    max_iter=max_it, verbose_level=1, batch_size=1,
                                    optimizer='sgd', weight_regularizer=reg)

[tcoln]

Thank you, I have explicitly set the batch_size=1, but it stiil doesn't work.

[GillesVandewiele]

On which dataset? And what is your result and what is the result reported in the paper? Please show me the code, if possible, you use to gather the result.

[rtavenar]

Hi @tcoln

Thanks for this feedback. I agree with @GillesVandewiele that hyper parameter setting is crucial, and that maybe with your code, we could give a hand. To the best of my knowledge, this implements what is in the paper. But there is one issue with this: on the small UCR/UEA datasets, there is a high variance in the results one can get using this kind of Stochastic Gradient Descent based algorithms. For example, if you compare results from the paper to those from the accompanying webpage or those from the timeseriesclassification.com website, you will see important differences.

Finally, I am unsure if the initialization of the weight is documented in the paper, and if it is implemented accordingly in tslearn. I will have to check that, but I'm on holidays at the moment, so I will do that later.

PS: also, in tslearn, we use hard min pooling instead of soft min pooling, but this should have very limited impact on the performance.

[tcoln]

I test on BeetleFly dataset as follows:

if __name__ == "__main__":
    params = util.read_params() # read parameters in LTS paper
    for i, param in params.items():
        dataset_name = param['dname']
        if dataset_name in ['Adiac', 'Beef']:
                continue
        k = param['K']
        L = param['L']
        R = param['R']
        lamda = param['lamda']
        epochs = param['maxIter']
        batch = 1

        X_train, y_train, X_test, y_test = UCR_UEA_datasets().load_dataset(dataset_name)
        X_train = TimeSeriesScalerMeanVariance().fit_transform(X_train)
        X_test = TimeSeriesScalerMeanVariance().fit_transform(X_test)
        n_ts = len(y_train)
        ts_sz = X_train.shape[1]
        n_classes = len(set(y_train))

        t0 = time.time()
        n_shapelets_per_size = grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_classes, k, L, R)
        clf = SerializableShapeletModel(n_shapelets_per_size=n_shapelets_per_size,
                        max_iter=epochs,
                        batch_size= batch,
                        weight_regularizer= lamda,
                        verbose_level=0)
        clf.fit(X_train, y_train)
        pred = clf.predict(X_test)
        acc = numpy.sum(y_test == pred)*1.0/len(y_test)
        print ('%s\t%d\t%d\t%d\t%d\tK=%.2f\tL=%.3f\tR=%d\tlamda=%.3f\tepochs=%d\tacc=%f\ttime=%d'%(dataset_name,n_ts,len(y_test),ts_sz,n_classes,k,L,R,lamda,epochs, acc, time.time()-t0))

BeetleFly 20 20 512 2 K=0.15 L=0.125 R=1 lamda=0.010 epochs=5000 acc=0.650000 time=181

But the accuracy is only 0.65, do you know why?

[GillesVandewiele]

Hmm, that is indeed a large difference... This could be due to 2 reasons:

1. The datasets are very small... Classifying one instance more or less correctly results in a difference of 5% in the accuracy score.

2. Sometimes, the train and test sets are partitioned rather carefully since samples can be correlated. In many papers, they will just randomly shuffle the data and re-partition the data without taking that property into account (which causes some label leakage). One good example is the Wine dataset, which causes rather low accuracies using the original splits, while the accuracy increases significantly by re-partitioning. Probably, multiple spectrographs were taken per grape/wine bottle, and therefore it is important to have the same group of measurements in only either the train or only the test set.

3. I managed to get ~0.9 on the BeetleFly dataset using the original partioning. Do you use the following config?

nr_shap, l, r, reg, max_it = 0.15, 0.125, 1, 0.01, 5000

EDIT: 4) After inspecting the file again, I notice that I do not use the classifier of the ShapeletModel from tslearn. I just use the discovered shapelets, calculate the distances and feed those distance matrices to a sci-kit learn classifier :). Perhaps that's where the difference comes from...

A crappy version of the script that managed to produce these results can be found at https://github.com/IBCNServices/GENDIS/blob/master/gendis/experiments/lts_vs_gendis.py

Seems like I also adapted the grabocka_params_to_shapelet_size_dict method:

def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
    base_size = int(l * ts_sz)
    d = {}
    for sz_idx in range(r):
        shp_sz = base_size * (sz_idx + 1)
        d[shp_sz] = n_shapelets
    return d

Which is then called as:

shapelet_dict = grabocka_params_to_shapelet_size_dict(
    X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
)

[tcoln]

Hi @tcoln

Thanks for this feedback. I agree with @GillesVandewiele that hyper parameter setting is crucial, and that maybe with your code, we could give a hand. To the best of my knowledge, this implements what is in the paper. But there is one issue with this: on the small UCR/UEA datasets, there is a high variance in the results one can get using this kind of Stochastic Gradient Descent based algorithms. For example, if you compare results from the paper to those from the accompanying webpage or those from the timeseriesclassification.com website, you will see important differences.

Finally, I am unsure if the initialization of the weight is documented in the paper, and if it is implemented accordingly in tslearn. I will have to check that, but I'm on holidays at the moment, so I will do that later.

PS: also, in tslearn, we use hard min pooling instead of soft min pooling, but this should have very limited impact on the performance.

Hi, @rtavenar
The initializations of shapelets and logistic regression are implemented as the same as described in the LTS paper. And the page of the accompanying webpage has been closed by the author?

[tcoln]

Hmm, that is indeed a large difference... This could be due to 2 reasons:

1. The datasets are very small... Classifying one instance more or less correctly results in a difference of 5% in the accuracy score.
2. Sometimes, the train and test sets are partitioned rather carefully since samples can be correlated. In many papers, they will just randomly shuffle the data and re-partition the data without taking that property into account (which causes some label leakage). One good example is the Wine dataset, which causes rather low accuracies using the original splits, while the accuracy increases significantly by re-partitioning. Probably, multiple spectrographs were taken per grape/wine bottle, and therefore it is important to have the same group of measurements in only either the train or only the test set.
3. I managed to get ~0.9 on the BeetleFly dataset using the original partioning. Do you use the following config?

nr_shap, l, r, reg, max_it = 0.15, 0.125, 1, 0.01, 5000

A crappy version of the script that managed to produce these results can be found at https://github.com/IBCNServices/GENDIS/blob/master/gendis/experiments/lts_vs_gendis.py

Seems like I also adapted the grabocka_params_to_shapelet_size_dict method:

def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
    base_size = int(l * ts_sz)
    d = {}
    for sz_idx in range(r):
        shp_sz = base_size * (sz_idx + 1)
        d[shp_sz] = n_shapelets
    return d

Which is then called as:

shapelet_dict = grabocka_params_to_shapelet_size_dict(
    X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
)

Yes, as you found, I use the true number of shapelets (K * number of ts) as in the paper, but it still can't achieve a reasonable result.

[GillesVandewiele]

I just ran the code listed below the get the following results (test accuracy of 0.95):

[LR] TRAIN Accuracy = 1.0
[LR] TRAIN Logloss = 0.6611554461961124
[LR] TEST Accuracy = 0.95
[LR] TEST Logloss = 0.681617564392368

It is now up to you to find the difference between the provided script and yours :). Also, as mentioned by @rtavenar the variance on these datasets can be rather large, so you better take the average over multiple runs...

import time
from collections import Counter
import warnings; warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.metrics import accuracy_score, log_loss
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

from tslearn.shapelets import ShapeletModel
from tslearn.datasets import UCR_UEA_datasets


def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
    base_size = int(l * ts_sz)
    d = {}
    for sz_idx in range(r):
        shp_sz = base_size * (sz_idx + 1)
        d[shp_sz] = n_shapelets
    return d

def fit_lr(X_distances_train, y_train, X_distances_test, y_test, out_path):
    lr = GridSearchCV(
            LogisticRegression(random_state=1337), 
            {
              'penalty': ['l1', 'l2'], 
              'C': [10**i for i in range(-2, 6)] + [5**i for i in range(-2, 6)],
              'class_weight': [None, 'balanced']
            }
        )
    lr.fit(X_distances_train, y_train)
    
    hard_preds = lr.predict(X_distances_test)
    proba_preds = lr.predict_proba(X_distances_test)
    
    hard_preds_train = lr.predict(X_distances_train)
    proba_preds_train = lr.predict_proba(X_distances_train)

    print("[LR] TRAIN Accuracy = {}".format(accuracy_score(y_train, hard_preds_train)))
    print("[LR] TRAIN Logloss = {}".format(log_loss(y_train, proba_preds_train)))
    print("[LR] TEST Accuracy = {}".format(accuracy_score(y_test, hard_preds)))
    print("[LR] TEST Logloss = {}".format(log_loss(y_test, proba_preds)))

    hard_preds = pd.DataFrame(hard_preds, columns=['prediction'])
    proba_preds = pd.DataFrame(proba_preds, columns=['proba_{}'.format(x) for x in set(list(y_train) + list(y_test))])

    hard_preds.to_csv(out_path.split('.')[0]+'_lr_hard.csv')
    proba_preds.to_csv(out_path.split('.')[0]+'_lr_proba.csv')

def lts_discovery(X_train, y_train, X_test, y_test,  nr_shap, l, r, reg, max_it, shap_out_path, pred_out_path, timing_out_path):
    # Fit LTS model, print metrics on test-set, write away predictions and shapelets
    shapelet_dict = grabocka_params_to_shapelet_size_dict(
            X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
    )
    
    clf = ShapeletModel(n_shapelets_per_size=shapelet_dict, 
                        max_iter=max_it, verbose_level=1, batch_size=1,
                        optimizer='sgd', weight_regularizer=reg)

    start = time.time()
    clf.fit(
        np.reshape(
            X_train, 
            (X_train.shape[0], X_train.shape[1], 1)
        ), 
        y_train
    )
    learning_time = time.time() - start

    print([len(x) for x in clf.shapelets_])
    print(clf.get_weights())

    print('Learning shapelets took {}s'.format(learning_time))

    X_distances_train = clf.transform(X_train)
    X_distances_test = clf.transform(X_test)

    fit_lr(X_distances_train, y_train, X_distances_test, y_test, pred_out_path)

hyper_parameters_lts = {
	'BeetleFly': 				[0.15, 0.125, 1, 0.01, 5000],
}

result_vectors = []
data_loader = UCR_UEA_datasets()

for dataset in hyper_parameters_lts:
    print(dataset)

    # Load the training and testing dataset (features + label vector)
    X_train, y_train, X_test, y_test = data_loader.load_dataset(dataset)

    X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1]))
    X_test = np.reshape(X_test, (X_test.shape[0], X_test.shape[1]))

    # Map labels to [0, .., C-1]
    map_dict = {}
    for j, c in enumerate(np.unique(y_train)):
        map_dict[c] = j
    y_train = pd.Series(y_train).map(map_dict).values
    y_test = pd.Series(y_test).map(map_dict).values

    # Get the best hyper-parameters for LTS (from original paper)
    nr_shap, l, r, reg, max_it = hyper_parameters_lts[dataset]

    # LTS Discovery & write away results to disk
    lts_discovery(X_train, y_train, X_test, y_test, nr_shap, l, r, reg, max_it,
            'results/lts_vs_genetic/{}_learned_shapelets_{}.txt'.format(dataset, int(time.time())), 
            'results/lts_vs_genetic/{}_learned_shapelets_predictions_{}.csv'.format(dataset, int(time.time())), 
            'results/lts_vs_genetic/{}_learned_runtime_{}.csv'.format(dataset, int(time.time()))
    )

[tcoln]

Yes， only on some datasets the code can achieve the results of LTS paper.

[tcoln]

I just ran the code listed below the get the following results (test accuracy of 0.95):

[LR] TRAIN Accuracy = 1.0
[LR] TRAIN Logloss = 0.6611554461961124
[LR] TEST Accuracy = 0.95
[LR] TEST Logloss = 0.681617564392368

It is now up to you to find the difference between the provided script and yours :). Also, as mentioned by @rtavenar the variance on these datasets can be rather large, so you better take the average over multiple runs...

import time
from collections import Counter
import warnings; warnings.filterwarnings('ignore')

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.metrics import accuracy_score, log_loss
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV

from tslearn.shapelets import ShapeletModel
from tslearn.datasets import UCR_UEA_datasets

def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
base_size = int(l * ts_sz)
d = {}
for sz_idx in range(r):
shp_sz = base_size * (sz_idx + 1)
d[shp_sz] = n_shapelets
return d

def fit_lr(X_distances_train, y_train, X_distances_test, y_test, out_path):
lr = GridSearchCV(
LogisticRegression(random_state=1337),
{
'penalty': ['l1', 'l2'],
'C': [10i for i in range(-2, 6)] + [5i for i in range(-2, 6)],
'class_weight': [None, 'balanced']
}
)
lr.fit(X_distances_train, y_train)

hard_preds = lr.predict(X_distances_test)
proba_preds = lr.predict_proba(X_distances_test)

hard_preds_train = lr.predict(X_distances_train)
proba_preds_train = lr.predict_proba(X_distances_train)

print("[LR] TRAIN Accuracy = {}".format(accuracy_score(y_train, hard_preds_train)))
print("[LR] TRAIN Logloss = {}".format(log_loss(y_train, proba_preds_train)))
print("[LR] TEST Accuracy = {}".format(accuracy_score(y_test, hard_preds)))
print("[LR] TEST Logloss = {}".format(log_loss(y_test, proba_preds)))

hard_preds = pd.DataFrame(hard_preds, columns=['prediction'])
proba_preds = pd.DataFrame(proba_preds, columns=['proba_{}'.format(x) for x in set(list(y_train) + list(y_test))])

hard_preds.to_csv(out_path.split('.')[0]+'_lr_hard.csv')
proba_preds.to_csv(out_path.split('.')[0]+'_lr_proba.csv')

def lts_discovery(X_train, y_train, X_test, y_test, nr_shap, l, r, reg, max_it, shap_out_path, pred_out_path, timing_out_path):
# Fit LTS model, print metrics on test-set, write away predictions and shapelets
shapelet_dict = grabocka_params_to_shapelet_size_dict(
X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
)

clf = ShapeletModel(n_shapelets_per_size=shapelet_dict, 
                    max_iter=max_it, verbose_level=1, batch_size=1,
                    optimizer='sgd', weight_regularizer=reg)

start = time.time()
clf.fit(
    np.reshape(
        X_train, 
        (X_train.shape[0], X_train.shape[1], 1)
    ), 
    y_train
)
learning_time = time.time() - start

print([len(x) for x in clf.shapelets_])
print(clf.get_weights())

print('Learning shapelets took {}s'.format(learning_time))

X_distances_train = clf.transform(X_train)
X_distances_test = clf.transform(X_test)

fit_lr(X_distances_train, y_train, X_distances_test, y_test, pred_out_path)

hyper_parameters_lts = {
'BeetleFly': [0.15, 0.125, 1, 0.01, 5000],
}

result_vectors = []
data_loader = UCR_UEA_datasets()

for dataset in hyper_parameters_lts:
print(dataset)

# Load the training and testing dataset (features + label vector)
X_train, y_train, X_test, y_test = data_loader.load_dataset(dataset)

X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0], X_test.shape[1]))

# Map labels to [0, .., C-1]
map_dict = {}
for j, c in enumerate(np.unique(y_train)):
    map_dict[c] = j
y_train = pd.Series(y_train).map(map_dict).values
y_test = pd.Series(y_test).map(map_dict).values

# Get the best hyper-parameters for LTS (from original paper)
nr_shap, l, r, reg, max_it = hyper_parameters_lts[dataset]

# LTS Discovery & write away results to disk
lts_discovery(X_train, y_train, X_test, y_test, nr_shap, l, r, reg, max_it,
        'results/lts_vs_genetic/{}_learned_shapelets_{}.txt'.format(dataset, int(time.time())), 
        'results/lts_vs_genetic/{}_learned_shapelets_predictions_{}.csv'.format(dataset, int(time.time())), 
        'results/lts_vs_genetic/{}_learned_runtime_{}.csv'.format(dataset, int(time.time()))
)

Hi, @GillesVandewiele ,
I have a question, the LTS algorithm has jointly learned shaplets and linear regression(LR), why you first learn shapelets and then learn LR separately?

[GillesVandewiele]

I just ran the code listed below the get the following results (test accuracy of 0.95):
[LR] TRAIN Accuracy = 1.0
[LR] TRAIN Logloss = 0.6611554461961124
[LR] TEST Accuracy = 0.95
[LR] TEST Logloss = 0.681617564392368
It is now up to you to find the difference between the provided script and yours :). Also, as mentioned by @rtavenar the variance on these datasets can be rather large, so you better take the average over multiple runs...
import time
from collections import Counter
import warnings; warnings.filterwarnings('ignore')
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import accuracy_score, log_loss
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV
from tslearn.shapelets import ShapeletModel
from tslearn.datasets import UCR_UEA_datasets
def grabocka_params_to_shapelet_size_dict(n_ts, ts_sz, n_shapelets, l, r):
base_size = int(l * ts_sz)
d = {}
for sz_idx in range(r):
shp_sz = base_size * (sz_idx + 1)
d[shp_sz] = n_shapelets
return d
def fit_lr(X_distances_train, y_train, X_distances_test, y_test, out_path):
lr = GridSearchCV(
LogisticRegression(random_state=1337),
{
'penalty': ['l1', 'l2'],
'C': [10i for i in range(-2, 6)] + [5i for i in range(-2, 6)],
'class_weight': [None, 'balanced']
}
)
lr.fit(X_distances_train, y_train)

hard_preds = lr.predict(X_distances_test)
proba_preds = lr.predict_proba(X_distances_test)

hard_preds_train = lr.predict(X_distances_train)
proba_preds_train = lr.predict_proba(X_distances_train)

print("[LR] TRAIN Accuracy = {}".format(accuracy_score(y_train, hard_preds_train)))
print("[LR] TRAIN Logloss = {}".format(log_loss(y_train, proba_preds_train)))
print("[LR] TEST Accuracy = {}".format(accuracy_score(y_test, hard_preds)))
print("[LR] TEST Logloss = {}".format(log_loss(y_test, proba_preds)))

hard_preds = pd.DataFrame(hard_preds, columns=['prediction'])
proba_preds = pd.DataFrame(proba_preds, columns=['proba_{}'.format(x) for x in set(list(y_train) + list(y_test))])

hard_preds.to_csv(out_path.split('.')[0]+'_lr_hard.csv')
proba_preds.to_csv(out_path.split('.')[0]+'_lr_proba.csv')

def lts_discovery(X_train, y_train, X_test, y_test, nr_shap, l, r, reg, max_it, shap_out_path, pred_out_path, timing_out_path):

Fit LTS model, print metrics on test-set, write away predictions and shapelets

shapelet_dict = grabocka_params_to_shapelet_size_dict(
X_train.shape[0], X_train.shape[1], int(nr_shap*X_train.shape[1]), l, r
)

clf = ShapeletModel(n_shapelets_per_size=shapelet_dict, 
                    max_iter=max_it, verbose_level=1, batch_size=1,
                    optimizer='sgd', weight_regularizer=reg)

start = time.time()
clf.fit(
    np.reshape(
        X_train, 
        (X_train.shape[0], X_train.shape[1], 1)
    ), 
    y_train
)
learning_time = time.time() - start

print([len(x) for x in clf.shapelets_])
print(clf.get_weights())

print('Learning shapelets took {}s'.format(learning_time))

X_distances_train = clf.transform(X_train)
X_distances_test = clf.transform(X_test)

fit_lr(X_distances_train, y_train, X_distances_test, y_test, pred_out_path)

hyper_parameters_lts = {
'BeetleFly': [0.15, 0.125, 1, 0.01, 5000],
}
result_vectors = []
data_loader = UCR_UEA_datasets()
for dataset in hyper_parameters_lts:
print(dataset)

# Load the training and testing dataset (features + label vector)
X_train, y_train, X_test, y_test = data_loader.load_dataset(dataset)

X_train = np.reshape(X_train, (X_train.shape[0], X_train.shape[1]))
X_test = np.reshape(X_test, (X_test.shape[0], X_test.shape[1]))

# Map labels to [0, .., C-1]
map_dict = {}
for j, c in enumerate(np.unique(y_train)):
    map_dict[c] = j
y_train = pd.Series(y_train).map(map_dict).values
y_test = pd.Series(y_test).map(map_dict).values

# Get the best hyper-parameters for LTS (from original paper)
nr_shap, l, r, reg, max_it = hyper_parameters_lts[dataset]

# LTS Discovery & write away results to disk
lts_discovery(X_train, y_train, X_test, y_test, nr_shap, l, r, reg, max_it,
        'results/lts_vs_genetic/{}_learned_shapelets_{}.txt'.format(dataset, int(time.time())), 
        'results/lts_vs_genetic/{}_learned_shapelets_predictions_{}.csv'.format(dataset, int(time.time())), 
        'results/lts_vs_genetic/{}_learned_runtime_{}.csv'.format(dataset, int(time.time()))
)

Hi, @GillesVandewiele ,
I have a question, the LTS algorithm has jointly learned shaplets and linear regression(LR), why you first learn shapelets and then learn LR separately?

You're right. I often get better performance out of the sklearn implementation, due to the hyper-parameter tuning. I can also use other, more powerful classifiers if needed as well.

[tcoln]

Hi @tcoln

Thanks for this feedback. I agree with @GillesVandewiele that hyper parameter setting is crucial, and that maybe with your code, we could give a hand. To the best of my knowledge, this implements what is in the paper. But there is one issue with this: on the small UCR/UEA datasets, there is a high variance in the results one can get using this kind of Stochastic Gradient Descent based algorithms. For example, if you compare results from the paper to those from the accompanying webpage or those from the timeseriesclassification.com website, you will see important differences.

Finally, I am unsure if the initialization of the weight is documented in the paper, and if it is implemented accordingly in tslearn. I will have to check that, but I'm on holidays at the moment, so I will do that later.

PS: also, in tslearn, we use hard min pooling instead of soft min pooling, but this should have very limited impact on the performance.

Hi, @rtavenar
I have implemented the soft min pooling. It does have limited impact on the performance. And we still can't achieve the good performance that the paper says.