GPU Acceleration

[hpswalters]

First, Thank you to the team for developing such an incredible tool. In the past, I've done my statistical forecasting in R, then migrated to multiprocessing calling R from Python. It is nice to have everything in one place...it makes the pipelines much easier. In addition, sktime offers so many additional capabilities that I'm just beginning to get my mind around.

Forecasting at scale (several thousand products), quickly is always a primary concern. Is there any way to leverage the GPU for the traditional forecasting algorithms (ARIMA, ETS, BATS/TBATS, etc) using sktime?

[fkiraly]

Depends - if you want to apply the algorithms to hierarchical data - as you mention, for thousands of products - or multivariate data, univariate algorithms like the ones mentionend above broadcast across variables and hierarchy instances.

This broadcasting can be parallelized via abstract backends, currently joblib and dask, and secondary backends of joblib such as ray or spark. This can be set via set_config (there is a known bug that prevents you from setting it properly on the current version 0.25.0, to be fixed in the upcoming release today/tomorrow). Further, if you want to add your favourite backend, this can be easily done in the sktime.utils.parallel module.

I wonder whether one could use GPUs for instance by using the ray backend of joblib,
https://docs.ray.io/en/latest/train/user-guides/using-gpus.html
or a GPU integration of spark: https://github.com/NVIDIA/spark-rapids

However, this will probably not work out of the box, as it seems that the models themselves need to be implemented with GPU acceleration in mind. In particular, the aforemetiond models (ARIMA, ETS, etc) are not currently implemented based on pytorch, RAPIDS, etc - so, the above may or may not be useful re GPU - but CPU parallelism via joblib is definitely supported.

GPU is also typically sth for neural networks, and sktime implements a pytorch adapter for forecasting, although the scope of model support is work in progress - in case you want to contribute, additions are welcome. Here, some of the models in the deep learning module can be trained on hierarchical data, @benHeid can perhaps give more input.

(I have myself not used GPU acceleration for hierarchical forecasting, only CPU - therefore, if you have any insights to add, I would be curious)

[fkiraly]

PS @hpswalters, I updated my response after a bit of research. Since you are looking at email, I am not sure if it gets resent with the update.

[benHeid]

Since the current torch based forecasters are quiet new they are not using GPU. However, this should be easy to add by calling .to(torch. device('cuda')) at different points if cuda is available and requested by the user.

[fkiraly]

well, sounds like a "good first issue"!

Probably needs to be done in the base pytorch adapter, and the for its descedants?

[hpswalters]

Thanks for you fast response. Unfortunately, I know of no instances where classical time series are accelerated on a GPU; however I’m very early in my search. Currently I multiprocess breaking dataset into groups that are then forecasted in process pools. From the looks of things, there seems to be a good deal of multi- core use already built-in since running with a pool of 2 frequently maxes my threadripper (32 physical cores).  If you are interested I can provide some high-level stats on how fast things are runningOn Jan 22, 2024, at 6:40 PM, Franz Király ***@***.***> wrote: Depends - if you want to apply the algorithms to hierarchical data, or multivariate data, univariate algorithms like the ones mentionend above broadcast across variables and hierarchy instances. This broadcasting can be parallelized via abstract backends, currently joblib and dask, and secondary backends of joblib such as ray or spark. Though, that's CPU, not GPU. GPU is typically sth for neural networks, and sktime implements a pytorch adapter for forecasting, although the scope of model support is work in progress - in case you want to contribute. If you know existing implementations of classical models for GPU, pointers would be appreciated - or ideas to abstract backends, or a direct contribution in sktime, too. Generally, accelerated versions of standard models are fiddlier, but sktime can easily abstract dependencies and backends, since we manage these on estimator level, not package level. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: ***@***.***>

[fkiraly]

If you are interested I can provide some high-level stats on how fast things are running

That would be nice!

The following should work for any forecaster, parallelization using CPU workers (not on 0.25.0, the bug will be fixed in 0.25.1)

myforecaster = MyFavouriteClass(params)

myforecaster.set_config(**{
    "backend:parallel": "joblib",  # set backend here
    "backend:parallel:params": {n_jobs=-1},  # pass params to backed, e.g., to joblib.Parallel
})

# then continue using the forecaster on hierarchical data
myforecaster.fit(y, fh=[1,2,3])
myforecaster.predict()

[fkiraly]

Unfortunately, I know of no instances where classical time series are accelerated on a GPU; however I’m very early in my search.

I suppose on could do the same as with the various accelerated versions of sklearn, though it would probaly be a resource intensive project...

I think though that, architecturally, unlike for the sklearn variants, there would not need to be a separate package for sktime, because - as said (in an edited version of my response abbove), sktime manages python dependencies and backends on estimator level. This architecture was necessary due to the high fragmentation of the time series space, and it turns out to be beneficial in issues such as parallelization or acceleration.

[hpswalters]

Looks like everything is coming through. Getting a little late here so I will provide some code snips tomorrow of how I run things and then you can help me understand how to integrate some of the performance metrics (the code you provided). Right now I’m just looking at processing timeOn Jan 22, 2024, at 7:21 PM, Franz Király ***@***.***> wrote: PS @hpswalters, I updated my response after a bit of research. Since you are looking at email, I am not sure if it gets resent with the update. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were mentioned.Message ID: ***@***.***>

[fkiraly]

then you can help me understand how to integrate some of the performance metrics (the code you provided)

Sure! I just notice that performance metrics cannot be parallelized right now since the configs are not linked to the parallelization point, but that should be only a small change to make.

[fkiraly]

added parallelization for performance metrics here: #5813

[hpswalters]

Breaking some of this down so I hope it is understandable. I'm a code slinger, not a developer--I code to get the job done faster.
Not everything is included here to keep it short--just the main stuff.
If you have any questions, don't hesitate to ask

#forecasters listed as lambda functions so I can use list comprehension and persist data for further analysis
model_mapping = {
    1: (lambda: skES(trend="additive", seasonal="additive", sp=12),'sktime exponential smoothing additive trend additive season'),
    2: (lambda: skES(trend="mul", seasonal="additive", sp=12),'sktime exponential smoothing multiplicative trend additive season'),
    3: (lambda: skES(trend="add", seasonal="mul", sp=12),'sktime exponential smoothing additive trend multiplicative season'),
    4: (lambda: skES(trend="mul", seasonal="mul", sp=12),'sktime exponential smoothing multiplicative trend multiplicative season'), 
    5: (lambda: skES(trend="mul", seasonal="mul", sp=12, use_boxcox=True),'sktime exponential smoothing multiplicative trend multiplicative season boxcox transform'),     
    6: (lambda: skAutoETS(n_jobs=1),'sktime auto exponential smoothing'),
    7: (lambda: skAutoArima(start_p=0, max_p=3, sp=12, seasonal=True, n_jobs=1,suppress_warnings=True),'sktime auto arima with seasonal component'),
    8: (lambda: skAutoArima(start_p=0, max_p=3, sp=1, seasonal=False, n_jobs=1,suppress_warnings=True),'sktime auto arima with no seasonal component'),
    9: (lambda: skAutoArima(start_p=0, max_p=3, n_jobs=1,suppress_warnings=True), 'sktime auto arima no seasonal spec'),
    10: (lambda: skProphet(),'sktime FB prohet'),
    11:(lambda: skSTL(sp=12), 'sktime STL'),
    12:(lambda: skTheta(sp=12), 'sktime Theta'),
    13:(lambda: sk_nixAutoETS(), 'sktime interface to Nixtla autoETS'),
    14:(lambda: sk_nixAutoCES(),'sktime interface to Nixtla autoCES'),
    15:(lambda: sk_nixAutoTheta(),'sktime interface to Nixtla autoTheta'),
    16:(lambda: sk_nixAutoArima(),'sktime interface to Nixtla autoArima'),
    17:(lambda: skPoly(degree=2), 'sktime polynomial trend'),
    18:(lambda: skMultiSeason(season_length=[3,12]), 'sktime multiple season'),
    19:(lambda: skAutoEnsForecaster(forecasters=[('sktime auto exponential smoothing',skAutoETS(n_jobs=1)),('sktime auto arima with no seasonal component',skAutoArima(start_p=0, max_p=3, sp=1, seasonal=False, n_jobs=1,suppress_warnings=True))]),'sktime auto ensemble with auto ets and auto arima'),
    20:(lambda: skAutoEnsForecaster(forecasters=[('sktime auto exponential smoothing',skAutoETS(n_jobs=1)),('sktime auto arima with no seasonal component',skAutoArima(start_p=0, max_p=3, sp=1, seasonal=False, n_jobs=1,suppress_warnings=True)),('sktime Theta', skTheta(sp=12))]),'sktime auto ensemble with auto ets, auto arima, auto Theta'),
    21:(lambda: skAutoEnsForecaster(forecasters=[('sktime auto exponential smoothing',skAutoETS(n_jobs=1)),('sktime auto arima with no seasonal component',skAutoArima(start_p=0, max_p=3, sp=1, seasonal=False,                                                                                                                                                                       n_jobs=1,suppress_warnings=True)),('sktime fb prophet', skProphet())]),'sktime auto ensemble with auto ets, auto arima, fb prophet'),
    22:(lambda:skEns())
    }
#quick and dirty function to reference the model id and return the forecaster and description
def get_model_by_id(model_id):
    model_entry = model_mapping.get(model_id)
    if model_entry:
        model_func, model_desc = model_entry
        return model_func(), model_desc
    else:
        raise ValueError(f"Model ID {model_id} is not defined.")

#function to actually perform the fit, forecast, and error calculation
#I'm calculating the average error out to lag3 based on the 3-period holdout from the test set
def fit_forecast_error(model_def,train,test,fh,desc):
    model_def.fit(train)
    sk_forecast = model_def.predict(fh=fh)
    rmse = np.sqrt(mean_squared_error(test, sk_forecast))
    return sk_forecast, rmse

#passing in the training and test set as well as forecast horizon and the pool id (i)
def forecast_and_evaluate(i,train,test,fh):
    forecast_rows = []
    rmse_rows = []
    average=test['PeriodDemand'].mean() 
    model_ids_to_use = [1,2,3,4,6,7,8,10,11,12,17,18,19,20,21]
    for model_id in model_ids_to_use:
        model_def, desc = get_model_by_id(model_id)
        horizon = 0
        try:
            fit_forecast, rmse = fit_forecast_error(model_def, train, test, fh, desc)
            
            # Accumulate forecast data
            for forecast_value in fit_forecast['PeriodDemand']:
                horizon += 1
                forecast_rows.append({'index': i, 'algoIdx': model_id, 'horizon': horizon, 'desc': desc, 'forecast': round(forecast_value,1)})
            
            # Accumulate RMSE data
            rmse_rows.append({'index': i, 'algoIdx': model_id,'desc': desc, 'rmse': round(rmse,1), 'CofV': round(rmse/average,4), 'calc_error': None})

        except Exception as e:
            print(f"Error in forecasting for index {i}, model {desc}: {e}")
            # Handle error by adding appropriate number of rows
            for _ in fh:
                forecast_rows.append({'index': i,'algoIdx': model_id, 'horizon': fh,'desc': desc, 'forecast': None})
            rmse_rows.append({'index': i,'algoIdx': model_id, 'desc': desc, 'rmse': None, 'CofV': None, 'calc_error': e})
        
    df_forecasts = pd.DataFrame(forecast_rows)
    df_rmse_errors = pd.DataFrame(rmse_rows)
    return df_forecasts, df_rmse_errors,i
#main program
if __name__ == "__main__":  

#pull the data from my database
#break each time series into training and test sets 

#Begin forecasting using multi-processing for forecasts and evaluation
    cores_to_use = 2
    all_forecasts = []
    all_rmses = []
    task_start_times = {}
    time_start_forecast = time.perf_counter()   

    def collect_result(result):
        if not hasattr(collect_result, "counter"):
            collect_result.counter = 0  # Initialize the counter attribute
        collect_result.counter += 1   # Increment the counter
        forecasts_df, rmses_df, index = result
        all_forecasts.append(forecasts_df)
        all_rmses.append(rmses_df)
        if (collect_result.counter % 10)==0:
            end_time = time.perf_counter()   
            elapsed_time = end_time - task_start_times[index]
            print(f"{collect_result.counter} products have been run through the forecasting engine. Time taken: {elapsed_time:.2f} seconds")
   #The multiprocessing part
    with Pool(cores_to_use) as pool:
        fh = list(range(1, len(dict_test_data[1]) + 1))
        results = []
        for i in range(1, total_records + 1):
            task_start_times[i] = time.perf_counter()  
           #using async 
            result = pool.apply_async(forecast_and_evaluate, (i, dict_train_data[i], dict_test_data[i], fh), callback=collect_result)
            results.append(result)

        # Wait for all tasks to complete
        for r in results:
            r.wait()

    # Concatenate results as before
    consolidated_forecasts_df = pd.concat(all_forecasts, ignore_index=True)
    consolidated_rmses_df = pd.concat(all_rmses, ignore_index=True)
            
    timeForecastAndError = time.perf_counter()
    TotalTime = timeForecastAndError - timeStart
    forecastAndErrorTime = timeForecastAndError - timeTrainTestSets

[fkiraly]

I'm a code slinger, not a developer

Not sure what that means, is that a D&D thing?

Some comments come to mind:

• you can pass iterables for fh and it will produce forecasts at all the horizons in the iterable
• are you aware of the evaluate utility, from sktime.forecasting.model_evaluation? It does what your code does, more or less afaik, and you can specify parallel backends
• would you be intersted to contribute to the benchmarking module? Some parts of it are current work in progress for forecasting.

[hpswalters]

Hmm. The code lost the indenting. I can attach as a file if that helps

[fkiraly]

The way you can do it is using this:

(I'm copying a screenshot from stackoverflow, because typing it out would display as a code block instead of the markdown source)

I've changed it in the code above, you can view the source by clicking on "edit"

[hpswalters]

I'm a code slinger, not a developer Not sure what that means, is that a D&D thing? Not a D&D thing. It's a term I've heard used to describe somebody like me who codes to get faster. I'm not in Python every day although I try to touch it enough so that I don't forget too much. At present, I'm attempting to automate and consolidate workflows I have developed--the forecasting is one part. I will have to look into the evaluate utility Yes, I would be interested in contributing, although, as I have said, my skills leave something to be desired

…

On Tue, Jan 23, 2024 at 4:14 PM Franz Király ***@***.***> wrote: I'm a code slinger, not a developer Not sure what that means, is that a D&D thing? Some comments come to mind: - you can pass iterables for fh and it will produce forecasts at all the horizons in the iterable - are you aware of the evaluate utility, from sktime.forecasting.model_evaluation? It does what your code does, more or less afaik, and you can specify parallel backends - would you be intersted to contribute to the benchmarking module? Some parts of it are current work in progress for forecasting. — Reply to this email directly, view it on GitHub <#5809 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHDHTABIM32XPOSIKSKYAX3YQAR5BAVCNFSM6AAAAABCFSZRECVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DEMRVGY3DS> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[fkiraly]

It's a term I've heard used to describe somebody like me
who codes to get faster.

I know, I was just joking! 😄

Yes, I would be interested in contributing, although, as I have said, my
skills leave something to be desired

Don't know about that, looks decent enough

[hpswalters]

Thanks for the tip!

…

On Thu, Jan 25, 2024 at 3:46 PM Benedikt Heidrich ***@***.***> wrote: Since the current torch based forecasters are quiet new they are not using GPU. However, this should be easy to add by calling .to(torch. device('cuda')) at different points if cuda is available and requested by the user. — Reply to this email directly, view it on GitHub <#5809 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AHDHTAAMBJQKE4FMAIBRUTTYQLAB5AVCNFSM6AAAAABCFSZRECVHI2DSMVQWIX3LMV43SRDJONRXK43TNFXW4Q3PNVWWK3TUHM4DENBZHAYTM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[yarnabrina]

This is after 2 months, but I think statsforecast supports GPU. I don't have GPU to test myself, but quoting from their documentation. We have them exposed through adapter, may be you can try at least those with GPU for better performance. It's not going to work with all models though.

This gives me an idea: can we consider adding a tag for estimators which supports GPU, and may be one which needs GPU? The needs one probably don't exist yet, but will do after Hugging Face adapter is done and concrete estimators are added with any big model.

[fkiraly]

hm, sure, we can add that tag. Would it make sense to have some more granular tagging for the type of GPU backend?

[yarnabrina]

May be, but my GPU knowledge is next to zero so can not comment.