diff --git a/pes_to_spec/model.py b/pes_to_spec/model.py
index ff640e142cfdef3ebded288f58fb15d62da7c8f7..0df98b8ff2b8c62936fb6a0adc23c28bfab593f3 100644
--- a/pes_to_spec/model.py
+++ b/pes_to_spec/model.py
@@ -8,7 +8,7 @@ from autograd import grad
 from scipy.signal import fftconvolve
 from scipy.signal import find_peaks_cwt
 from scipy.optimize import fmin_l_bfgs_b
-from sklearn.decomposition import PCA
+from sklearn.decomposition import PCA, IncrementalPCA
 from sklearn.pipeline import Pipeline, FeatureUnion
 from sklearn.base import TransformerMixin, BaseEstimator
 from sklearn.base import RegressorMixin
@@ -436,7 +436,7 @@ def _fit_estimator(estimator, X: np.ndarray, y: np.ndarray):
 
 class MultiOutputWithStd(MetaEstimatorMixin, BaseEstimator):
 
-    def __init__(self, estimator, *, n_jobs=8):
+    def __init__(self, estimator, *, n_jobs=-1):
         self.estimator = estimator
         self.n_jobs = n_jobs
 
@@ -525,11 +525,11 @@ class Model(TransformerMixin, BaseEstimator):
         x_model_steps = list()
         x_model_steps += [('select', SelectRelevantLowResolution(channels, tof_start, delta_tof))]
         if n_nonlinear_kernel > 0:
-            x_model_steps += [('fex', Pipeline([('prepca', PCA(n_pca_lr, whiten=True)),
+            x_model_steps += [('fex', Pipeline([('prepca', IncrementalPCA(n_pca_lr, whiten=True, batch_size=n_pca_lr*2)),
                                                 ('nystroem', Nystroem(n_components=n_nonlinear_kernel, kernel='rbf', gamma=None, n_jobs=8)),
                                                 ]))]
         x_model_steps += [
-                          ('pca', PCA(n_pca_lr, whiten=True)),
+                          ('pca', IncrementalPCA(n_pca_lr, whiten=True, batch_size=n_pca_lr*2)),
                           ('unc', UncertaintyHolder()),
                           ]
         self.x_model = Pipeline(x_model_steps)
@@ -593,9 +593,12 @@ class Model(TransformerMixin, BaseEstimator):
 
         Returns: Smoothened high resolution spectrum.
         """
+        print("Fitting x ...")
         x_t = self.x_model.fit_transform(low_res_data)
+        print("Fitting y ...")
         y_t = self.y_model.fit_transform(high_res_data, smoothen__energy=high_res_photon_energy)
         #self.fit_model.set_params(fex__gamma=1.0/float(x_t.shape[0]))
+        print("Fitting linear model ...")
         self.fit_model.fit(x_t, y_t)
 
         # calculate the effect of the PCA