PERF: Improve performance of multivarite OLS

Add short circuits when applicable to multivariate OLS
General clean up of small issues

closes #158

Author: bashtage

linearmodels/asset_pricing/model.py

@@ -2,7 +2,6 @@
 Linear factor models for applications in asset pricing
 """
 import numpy as np
-from numpy.linalg import pinv
 from patsy.highlevel import dmatrix
 from patsy.missing import NAAction
 from scipy.optimize import minimize
@@ -16,7 +15,7 @@
 from linearmodels.compat.numpy import lstsq
 from linearmodels.iv.data import IVData
 from linearmodels.utility import (AttrDict, WaldTestStatistic, has_constant,
-                                  matrix_rank, missing_warning)
+                                  missing_warning)
 
 
 def callback_factory(obj, args, disp=1):
@@ -103,9 +102,9 @@ def _validate_data(self):
             raise ValueError('portfolios must not contains a constant or equivalent.')
         if has_constant(f)[0]:
             raise ValueError('factors must not contain a constant or equivalent.')
-        if matrix_rank(f) < f.shape[1]:
+        if np.linalg.matrix_rank(f) < f.shape[1]:
             raise ValueError('Model cannot be estimated. factors do not have full column rank.')
-        if matrix_rank(p) < p.shape[1]:
+        if np.linalg.matrix_rank(p) < p.shape[1]:
             raise ValueError('Model cannot be estimated. portfolios do not have full column rank.')
 
     @property
@@ -212,13 +211,13 @@ def fit(self, cov_type='robust', debiased=True, **cov_config):
         fc = np.c_[np.ones((nobs, 1)), f]
         rp = f.mean(0)[:, None]
         fe = f - f.mean(0)
-        b = pinv(fc) @ p
+        b = np.linalg.pinv(fc) @ p
         eps = p - fc @ b
         alphas = b[:1].T
 
         nloading = (nfactor + 1) * nportfolio
         xpxi = np.eye(nloading + nfactor)
-        xpxi[:nloading, :nloading] = np.kron(np.eye(nportfolio), pinv(fc.T @ fc / nobs))
+        xpxi[:nloading, :nloading] = np.kron(np.eye(nportfolio), np.linalg.pinv(fc.T @ fc / nobs))
         f_rep = np.tile(fc, (1, nportfolio))
         eps_rep = np.tile(eps, (nfactor + 1, 1))  # 1 2 3 ... 25 1 2 3 ...
         eps_rep = eps_rep.ravel(order='F')
@@ -251,7 +250,7 @@ def fit(self, cov_type='robust', debiased=True, **cov_config):
 
         # Return values
         alpha_vcv = vcv[:nportfolio, :nportfolio]
-        stat = float(alphas.T @ pinv(alpha_vcv) @ alphas)
+        stat = float(alphas.T @ np.linalg.pinv(alpha_vcv) @ alphas)
         jstat = WaldTestStatistic(stat, 'All alphas are 0', nportfolio, name='J-statistic')
         params = b.T
         betas = b[1:].T

linearmodels/panel/data.py

@@ -1,8 +1,7 @@
 from itertools import product
 
 import numpy as np
-import pandas as pd
-from pandas import DataFrame, Panel, Series
+from pandas import DataFrame, Panel, Series, MultiIndex, get_dummies, Categorical
 
 from linearmodels.compat.numpy import lstsq
 from linearmodels.compat.pandas import (is_categorical,
@@ -32,10 +31,10 @@ class _Panel(object):
     def __init__(self, df):
         self._items = df.columns
         index = df.index
-        self._major_axis = pd.Series(index.levels[1][index.labels[1]]).unique()
-        self._minor_axis = pd.Series(index.levels[0][index.labels[0]]).unique()
-        self._full_index = pd.MultiIndex.from_product([self._minor_axis,
-                                                       self._major_axis])
+        self._major_axis = Series(index.levels[1][index.labels[1]]).unique()
+        self._minor_axis = Series(index.levels[0][index.labels[0]]).unique()
+        self._full_index = MultiIndex.from_product([self._minor_axis,
+                                                    self._major_axis])
         new_df = df.reindex(self._full_index)
         self._frame = new_df
         i, j, k = len(self._items), len(self._major_axis), len(self.minor_axis)
@@ -45,12 +44,12 @@ def __init__(self, df):
     @classmethod
     def from_array(cls, values, items, major_axis, minor_axis):
         index = list(product(minor_axis, major_axis))
-        index = pd.MultiIndex.from_tuples(index)
+        index = MultiIndex.from_tuples(index)
         i, j, k = len(items), len(major_axis), len(minor_axis)
         values = np.swapaxes(values.copy(), 0, 2).ravel()
         values = np.reshape(values, ((j * k), i))
 
-        df = pd.DataFrame(values, index=index, columns=items)
+        df = DataFrame(values, index=index, columns=items)
         return cls(df)
 
     @property
@@ -82,7 +81,7 @@ def convert_columns(s, drop_first):
         s = s.astype('category')
 
     if is_categorical(s):
-        out = pd.get_dummies(s, drop_first=drop_first)
+        out = get_dummies(s, drop_first=drop_first)
         out.columns = [str(s.name) + '.' + str(c) for c in out]
         return out
     return s
@@ -169,18 +168,18 @@ def __init__(self, x, var_name='x', convert_dummies=True, drop_first=True, copy=
             except ImportError:
                 pass
 
-        if isinstance(x, Series) and isinstance(x.index, pd.MultiIndex):
+        if isinstance(x, Series) and isinstance(x.index, MultiIndex):
             x = DataFrame(x)
         elif isinstance(x, Series):
             raise ValueError('Series can only be used with a 2-level MultiIndex')
 
         if isinstance(x, (Panel, DataFrame)):
             if isinstance(x, DataFrame):
-                if isinstance(x.index, pd.MultiIndex):
+                if isinstance(x.index, MultiIndex):
                     if len(x.index.levels) != 2:
                         raise ValueError('DataFrame input must have a '
                                          'MultiIndex with 2 levels')
-                    if isinstance(self._original, (pd.DataFrame, PanelData, pd.Series)):
+                    if isinstance(self._original, (DataFrame, PanelData, Series)):
                         for i in range(2):
                             index_names[i] = x.index.levels[i].name or index_names[i]
                     self._frame = x
@@ -388,9 +387,9 @@ def general_demean(self, groups, weights=None):
         if not isinstance(groups, PanelData):
             groups = PanelData(groups)
         if weights is None:
-            weights = PanelData(pd.DataFrame(np.ones((self._frame.shape[0], 1)),
-                                             index=self.index,
-                                             columns=['weights']))
+            weights = PanelData(DataFrame(np.ones((self._frame.shape[0], 1)),
+                                          index=self.index,
+                                          columns=['weights']))
         weights = weights.values2d
         groups = groups.values2d.astype(np.int64, copy=False)
 
@@ -417,11 +416,11 @@ def demean_pass(frame, weights, root_w):
             return frame
 
         # Swap out the index for better performance
-        init_index = pd.DataFrame(groups)
+        init_index = DataFrame(groups)
         init_index.set_index(list(init_index.columns), inplace=True)
 
         root_w = np.sqrt(weights)
-        weights = pd.DataFrame(weights, index=init_index.index)
+        weights = DataFrame(weights, index=init_index.index)
         wframe = root_w * self._frame
         wframe.index = init_index.index
 
@@ -619,7 +618,7 @@ def dummies(self, group='entity', drop_first=False):
         axis = 0 if group == 'entity' else 1
         labels = self._frame.index.labels
         levels = self._frame.index.levels
-        cat = pd.Categorical(levels[axis][labels[axis]])
-        dummies = pd.get_dummies(cat, drop_first=drop_first)
+        cat = Categorical(levels[axis][labels[axis]])
+        dummies = get_dummies(cat, drop_first=drop_first)
         cols = self.entities if group == 'entity' else self.time
         return dummies[[c for c in cols if c in dummies]].astype(np.float64, copy=False)

linearmodels/system/_utility.py

@@ -1,9 +1,6 @@
 import numpy as np
 import pandas as pd
-from numpy import cumsum, diag, eye, zeros
-from numpy.linalg import inv
-
-from linearmodels.utility import matrix_rank
+from numpy.linalg import inv, matrix_rank
 
 
 def blocked_column_product(x, s):
@@ -85,9 +82,19 @@ def blocked_inner_prod(x, s):
     """
     k = len(x)
     widths = list(map(lambda m: m.shape[1], x))
-    cum_width = cumsum([0] + widths)
+    s_is_diag = np.all((s - np.diag(np.diag(s))) == 0.0)
+
+    w0 = widths[0]
+    homogeneous = all([w == w0 for w in widths])
+    if homogeneous and not s_is_diag:
+        # Fast path when all x have same number of columns
+        # Slower than diag case when k is large since many 0s
+        x = np.hstack(x)
+        return x.T @ x * np.kron(s, np.ones((w0, w0)))
+
+    cum_width = np.cumsum([0] + widths)
     total = sum(widths)
-    out = zeros((total, total))
+    out = np.zeros((total, total))
 
     for i in range(k):
         xi = x[i]
@@ -97,7 +104,7 @@ def blocked_inner_prod(x, s):
         out[sel_i, sel_i] = prod
 
     # Short circuit if identity
-    if np.all((s - diag(diag(s))) == 0.0):
+    if s_is_diag:
         return out
 
     for i in range(k):
@@ -183,7 +190,7 @@ def blocked_full_inner_product(x, s):
 def inv_matrix_sqrt(s):
     vecs, vals = np.linalg.eigh(s)
     vecs = 1.0 / np.sqrt(vecs)
-    out = vals @ diag(vecs) @ vals.T
+    out = vals @ np.diag(vecs) @ vals.T
     return (out + out.T) / 2
 
 
@@ -261,7 +268,7 @@ def _verify_constraints(self):
     def _compute_transform(self):
         r = self._ra
         c, k = r.shape
-        m = eye(k) - r.T @ inv(r @ r.T) @ r
+        m = np.eye(k) - r.T @ inv(r @ r.T) @ r
         vals, vecs = np.linalg.eigh(m)
         vals = np.real(vals)
         vecs = np.real(vecs)

linearmodels/system/covariance.py

@@ -181,15 +181,27 @@ def __init__(self, x, eps, sigma, full_sigma, gls=False, debiased=False, constra
         self._name = 'Heteroskedastic (Robust) Covariance'
 
         k = len(x)
-        weights = inv(sigma) if gls else eye(k)
-        bigx = blocked_diag_product(x, weights)
         nobs = eps.shape[0]
-        e = eps.T.ravel()[:, None]
-        bigxe = bigx * e
-        m = bigx.shape[1]
-        xe = zeros((nobs, m))
-        for i in range(nobs):
-            xe[i, :] = bigxe[i::nobs].sum(0)[None, :]
+
+        if gls:
+            weights = inv(sigma)
+            bigx = blocked_diag_product(x, weights)
+            e = eps.T.ravel()[:, None]
+            bigxe = bigx * e
+            m = bigx.shape[1]
+            xe = zeros((nobs, m))
+            for i in range(nobs):
+                xe[i, :] = bigxe[i::nobs].sum(0)[None, :]
+        else:
+            # Do not require blocking when not using GLS
+            k_tot = sum(map(lambda a: a.shape[1], x))
+            xe = empty((nobs, k_tot))
+            loc = 0
+            for i in range(k):
+                offset = x[i].shape[1]
+                xe[:, loc:loc+offset] = x[i] * eps[:, i:i+1]
+                loc += offset
+
         self._moments = xe
 
     def _xeex(self):