All places where inv XTX was used has been substituted by the

linalg.solve approach
Doctrings updated accordingly with a note about over-subsription

nrt/fit_methods.py

@@ -82,6 +82,14 @@ def rirls(X, y, M=bisquare, tune=4.685,
     according to weight function and tuning parameter.
     Basically a clone from `statsmodels` that should be much faster.
 
+    Note:
+        For best performances of the multithreaded implementation, it is
+        recommended to limit the number of threads used by MKL or OpenBLAS to 1.
+        This avoids over-subscription, and improves performances.
+        By default the function will use all cores available; the number of cores
+        used can be controled using the ``numba.set_num_threads`` function or
+        by modifying the ``NUMBA_NUM_THREADS`` environment variable
+
     Args:
         X (np.ndarray): 2D (n_obs x n_features) design matrix
         y (np.ndarray): 1D independent variable
@@ -176,6 +184,14 @@ def ccdc_stable_fit(X, y, dates, threshold=3):
     2. first observation / RMSE < threshold
     3.  last observation / RMSE < threshold
 
+    Note:
+        For best performances of the multithreaded implementation, it is
+        recommended to limit the number of threads used by MKL or OpenBLAS to 1.
+        This avoids over-subscription, and improves performances.
+        By default the function will use all cores available; the number of cores
+        used can be controled using the ``numba.set_num_threads`` function or
+        by modifying the ``NUMBA_NUM_THREADS`` environment variable
+
     Args:
         X ((M, N) np.ndarray): Matrix of independant variables
         y ((M, K) np.ndarray): Matrix of dependant variables
@@ -209,9 +225,7 @@ def ccdc_stable_fit(X, y, dates, threshold=3):
             # each iteration
             y_ = y_sub[-jdx:]
             X_ = X_sub[-jdx:]
-            XTX = np.linalg.inv(np.dot(X_.T, X_))
-            XTY = np.dot(X_.T, y_)
-            beta_sub = np.dot(XTX, XTY)
+            beta_sub = np.linalg.solve(np.dot(X_.T, X_), np.dot(X_.T, y_))
             resid_sub = np.dot(X_, beta_sub) - y_
 
             # Check for stability
@@ -252,6 +266,14 @@ def roc_stable_fit(X, y, dates, alpha=0.05, crit=0.9478982340418134):
     The implementation roughly corresponds to the fit of bfastmonitor
     with the history option set to 'ROC'.
 
+    Note:
+        For best performances of the multithreaded implementation, it is
+        recommended to limit the number of threads used by MKL or OpenBLAS to 1.
+        This avoids over-subscription, and improves performances.
+        By default the function will use all cores available; the number of cores
+        used can be controled using the ``numba.set_num_threads`` function or
+        by modifying the ``NUMBA_NUM_THREADS`` environment variable
+
     Args:
         X ((M, N) np.ndarray): Matrix of independant variables
         y ((M, K) np.ndarray): Matrix of dependant variables
@@ -262,6 +284,7 @@ def roc_stable_fit(X, y, dates, alpha=0.05, crit=0.9478982340418134):
         crit (float): Critical value corresponding to the chosen alpha. Can be
             calculated with ``_cusum_rec_test_crit``.
             Default is the value for alpha=0.05
+
     Returns:
         beta (numpy.ndarray): The array of regression estimators
         residuals (numpy.ndarray): The array of residuals
@@ -300,9 +323,8 @@ def roc_stable_fit(X, y, dates, alpha=0.05, crit=0.9478982340418134):
         # Subset and fit
         X_stable = _X[stable_idx:]
         y_stable = _y[stable_idx:]
-        XTX = np.linalg.inv(np.dot(X_stable.T, X_stable))
-        XTY = np.dot(X_stable.T, y_stable)
-        beta[:, idx] = np.dot(XTX, XTY)
+        beta[:, idx] = np.linalg.solve(np.dot(X_stable.T, X_stable),
+                                       np.dot(X_stable.T, y_stable))
         fit_start[idx] = _dates[stable_idx]
 
     residuals = np.dot(X, beta) - y