I've been investigating problems where our MKL-based sparse matrices are massively underperforming scipy.sparse. For example:

  operation           storage memory        time
0    matvec  scipy.sparse csc      0  0.00211215
1    matvec  quantcore.matrix      0   0.0266283

This is a matrix with 3e6 rows and 3 columns

It seems like having a small number of columns makes MKL perform quite poorly. I'm not sure why that's the case. But, it may be worth having a check and just falling back to scipy.sparse in narrow cases like this. This kind of narrow case may actually be the dominant use case for sparse matrices because they will be a small component of a SplitMatrix.

This might fix https://github.com/Quantco/quantcore.glm/issues/323. I think we were passing in the number of rows into matvec when we mean to pass in the number of columns. But maybe I'm misunderstanding what's going on.

The function signature for matvec is

https://github.com/Quantco/quantcore.matrix/blob/9ef54c6cb21e8d8063c0968fe47c300b79d3af4b/src/quantcore/matrix/ext/categorical.pyx#L61-L62

but before we were passing in the number of rows as last argument.

I see the build script for linux uses jemalloc with disable-tls: "./autogen.sh --disable-cxx --with-jemalloc-prefix=local --with-install-suffix=local --disable-tls --disable-initial-exec-tls",

However, the source distribution in PyPI doesn't run that script when installing it through pip, relying instead on whatever jemalloc it finds when it tries to compile. If, for example, one tries to install tabmat from source through pip, it will later on fail to import, complaining about an error with jemalloc:

cannot allocate memory in static TLS block

The installation wia conda froge was getting stuck in the "Solving environment" part so I tried to install with pip, given that the package is available on Pypi. pip install glum runs in seconds, but then I am unable to import stuff from it, with the following error:

In [1]: from glum import GeneralizedLinearRegressor
---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-1-0284693fe484> in <module>
----> 1 from glum import GeneralizedLinearRegressor

~/anaconda3/lib/python3.9/site-packages/glum/__init__.py in <module>
      1 import pkg_resources
      2 
----> 3 from ._distribution import TweedieDistribution
      4 from ._glm import GeneralizedLinearRegressor
      5 from ._glm_cv import GeneralizedLinearRegressorCV

~/anaconda3/lib/python3.9/site-packages/glum/_distribution.py in <module>
      6 import numpy as np
      7 from scipy import sparse, special
----> 8 from tabmat import MatrixBase, StandardizedMatrix
      9 
     10 from ._functions import (

~/anaconda3/lib/python3.9/site-packages/tabmat/__init__.py in <module>
----> 1 from .categorical_matrix import CategoricalMatrix
      2 from .constructor import from_csc, from_pandas
      3 from .dense_matrix import DenseMatrix
      4 from .matrix_base import MatrixBase
      5 from .sparse_matrix import SparseMatrix

~/anaconda3/lib/python3.9/site-packages/tabmat/categorical_matrix.py in <module>
    171 from .ext.split import sandwich_cat_cat, sandwich_cat_dense
    172 from .matrix_base import MatrixBase
--> 173 from .sparse_matrix import SparseMatrix
    174 from .util import (
    175     check_matvec_out_shape,

~/anaconda3/lib/python3.9/site-packages/tabmat/sparse_matrix.py in <module>
      4 from scipy import sparse as sps
      5 
----> 6 from .ext.sparse import (
      7     csc_rmatvec,
      8     csc_rmatvec_unrestricted,

ImportError: /home/mathurin/anaconda3/lib/python3.9/site-packages/tabmat/ext/../../tabmat.libs/libjemalloclocal-691a3dac.so.2: cannot allocate memory in static TLS block

googling did not help, is there a way to make the pip-installed version work ?

• Allow SplitMatrix to be constructed from another SplitMatrix.
• Allow inputs of SplitMatrix to be 1-d
• Implement __getitem__ for column subset
• Also had to implement column subsetting for CategoricalMatrix
• __repr__ uses the __repr__ method of components instead of str()

ToDo:

• [ ] FIX BUG WITH _split_col_subsets (first confirm that it's a bug)
• [ ] Add testing for new features

Checklist

• [ ] Added a CHANGELOG.rst entry

Currently, CategoricalMatrix does not provide an easy way to drop a column. We are required to include a category for every row in the dataset, but in an unregularized setting, it is nice to sometimes drop one column.

Something sort of like this is already implemented by the cols parameter to the matrix vector and sandwich functions.

Requirements: pip install libsvmdata

The following script gives me a segfault:

from libsvmdata import fetch_libsvm
from glum import GeneralizedLinearRegressor

X, y = fetch_libsvm("rcv1.binary")
clf = GeneralizedLinearRegressor(alpha=0.01, fit_intercept=False,
                                 family="gaussian")
clf.fit(X, y)

Output:

In [1]: %run glum_segfault.py
Dataset: rcv1.binary
[1]    271745 segmentation fault (core dumped)  ipython

I'm using glum 2.0.3

@qb3

Bumps pypa/cibuildwheel from 2.2.2 to 2.3.0.

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting @dependabot rebase.

We are currently observing issues when using quantcore.matrix in conjunction with onnx and onnxruntime on MacOS. The call to python -c 'import onnx; import quantcore.matrix.ext.dense; import onnxruntime' fails with a bus error or segfault whereas the call DYLD_INSERT_LIBRARIES=$CONDA_PREFIX/lib/libjemalloc.dylib python -c 'import onnx; import quantcore.matrix.ext.dense; import onnxruntime' passes just fine. This indicates that using an unnamespaced jemalloc may be problematic here as the following traceback indicates:

collecting ... Process 6259 stopped
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=2, address=0x4efffffff7)
    frame #0: 0x000000013c0f3704 libjemalloc.2.dylib`je_free_default + 240
libjemalloc.2.dylib`je_free_default:
->  0x13c0f3704 <+240>: str    x20, [x8, w9, sxtw #3]
    0x13c0f3708 <+244>: ldr    w8, [x19, #0x200]
    0x13c0f370c <+248>: sub    w9, w8, #0x1              ; =0x1
    0x13c0f3710 <+252>: str    w9, [x19, #0x200]
Target 0: (python) stopped.
(lldb) bt
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=2, address=0x4efffffff7)
  * frame #0: 0x000000013c0f3704 libjemalloc.2.dylib`je_free_default + 240
    frame #1: 0x0000000142745010 onnxruntime_pybind11_state.so`std::__1::__hash_table<std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, std::__1::__unordered_map_hasher<std::__1::type_index, std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, pybind11::detail::type_hash, pybind11::detail::type_equal_to, true>, std::__1::__unordered_map_equal<std::__1::type_index, std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, pybind11::detail::type_equal_to, pybind11::detail::type_hash, true>, std::__1::allocator<std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > > > >::__rehash(unsigned long) + 76
    frame #2: 0x0000000142744dd0 onnxruntime_pybind11_state.so`std::__1::pair<std::__1::__hash_iterator<std::__1::__hash_node<std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, void*>*>, bool> std::__1::__hash_table<std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, std::__1::__unordered_map_hasher<std::__1::type_index, std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, pybind11::detail::type_hash, pybind11::detail::type_equal_to, true>, std::__1::__unordered_map_equal<std::__1::type_index, std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > >, pybind11::detail::type_equal_to, pybind11::detail::type_hash, true>, std::__1::allocator<std::__1::__hash_value_type<std::__1::type_index, std::__1::vector<bool (*)(_object*, void*&), std::__1::allocator<bool (*)(_object*, void*&)> > > > >::__emplace_unique_key_args<std::__1::type_index, std::__1::piecewise_construct_t const&, std::__1::tuple<std::__1::type_index const&>, std::__1::tuple<> >(std::__1::type_index const&, std::__1::piecewise_construct_t const&, std::__1::tuple<std::__1::type_index const&>&&, std::__1::tuple<>&&) + 480
    frame #3: 0x00000001427427dc onnxruntime_pybind11_state.so`pybind11::detail::generic_type::initialize(pybind11::detail::type_record const&) + 396
    frame #4: 0x0000000142751688 onnxruntime_pybind11_state.so`pybind11::class_<onnxruntime::ExecutionOrder>::class_<>(pybind11::handle, char const*) + 140
    frame #5: 0x00000001427513f8 onnxruntime_pybind11_state.so`pybind11::enum_<onnxruntime::ExecutionOrder>::enum_<>(pybind11::handle const&, char const*) + 52
    frame #6: 0x000000014272b5c8 onnxruntime_pybind11_state.so`onnxruntime::python::addObjectMethods(pybind11::module_&, onnxruntime::Environment&) + 296
    frame #7: 0x0000000142734e68 onnxruntime_pybind11_state.so`PyInit_onnxruntime_pybind11_state + 340
    frame #8: 0x000000010019f994 python`_imp_create_dynamic + 2412
    frame #9: 0x00000001000b40f8 python`cfunction_vectorcall_FASTCALL + 208
    frame #10: 0x000000010016bfd8 python`_PyEval_EvalFrameDefault + 30088

My suggestion would be to add an output to the jemalloc-feedstock as described in https://github.com/conda-forge/jemalloc-feedstock/issues/23 that comes with a prefixed version of the library.

Bumps google-github-actions/setup-gcloud from 0.2.0 to 0.2.1.

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting @dependabot rebase.

Dependabot will merge this PR once CI passes on it, as requested by @xhochy.

Updating the flake8 config to match the new flake8 config from glm_benchmarks.

Changes:

• changed linter according to the issue
• added simple docstrings to public functions (most functions were in the main matrix classes)
• preceded function names with underscores if the functions were only being used internally
• added “no docstrings in magic function” flake8 error to list of ignores (didn’t seem helpful)
• added # noqa in places where flake8 errors were just creating issues in unhelpful places

Closes #45

Checklist

• [ ] Added a CHANGELOG.rst entry

Bumps pypa/cibuildwheel from 2.11.3 to 2.11.4.

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting @dependabot rebase.

This throws an error:

import numpy as np
import tabmat
from scipy.sparse import csc_matrix

rng = np.random.default_rng(seed=123)
X = rng.standard_normal(size=(100,20))
Xd = tabmat.DenseMatrix(X[:,:10])
Xs = tabmat.SparseMatrix(csc_matrix(X[:,10:]))
Xm = tabmat.SplitMatrix([Xd, Xs])
Xm.sandwich(np.ones(X.shape[0]))

---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-2-91ba52e4f568> in <module>
      8 Xs = tabmat.SparseMatrix(csc_matrix(X[:,10:]))
      9 Xm = tabmat.SplitMatrix([Xd, Xs])
---> 10 Xm.sandwich(np.ones(X.shape[0]))

~/anaconda3/envs/py39/lib/python3.9/site-packages/tabmat/split_matrix.py in sandwich(self, d, rows, cols)
    287             idx_i = subset_cols_indices[i]
    288             mat_i = self.matrices[i]
--> 289             res = mat_i.sandwich(d, rows, subset_cols[i])
    290             if isinstance(res, sps.dia_matrix):
    291                 out[(idx_i, idx_i)] += np.squeeze(res.data)

~/anaconda3/envs/py39/lib/python3.9/site-packages/tabmat/dense_matrix.py in sandwich(self, d, rows, cols)
     62         d = np.asarray(d)
     63         rows, cols = setup_restrictions(self.shape, rows, cols)
---> 64         return dense_sandwich(self, d, rows, cols)
     65 
     66     def _cross_sandwich(

src/tabmat/ext/dense.pyx in tabmat.ext.dense.dense_sandwich()

Exception: 

Compare against this:

Xd = tabmat.DenseMatrix(X[:,:10].copy())
Xs = tabmat.SparseMatrix(csc_matrix(X[:,10:]))
Xm = tabmat.SplitMatrix([Xd, Xs])
Xm.sandwich(np.ones(X.shape[0]))

(No error)

I find it convenient to be able to check directly inside a python shell.

In [1]: import tabmat

In [2]: tabmat.__version__
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-2-aae82a909ca3> in <module>
----> 1 tabmat.__version__

AttributeError: module 'tabmat' has no attribute '__version__'

this is tabmat 3.0.7 installed via PyPI

I think we've discussed this, but I don't remember the conclusion and can't find an issue now.

We recommend from_pandas as the way "most users" should construct tabmat objects. from_pandas then guesses which columns should be treated as categorical. I think it would be really nice to have Patsy-like formulas as an alternative, since

1. R users (including many economists) like using formulas, and
2. It's easy to infer from a Patsy formula which columns are categorical, which are sparse (generally interactions with categoricals), and which are dense (everything else), so this could remove some of the guesswork from tabmat and improve performance.

I'm not sure how feasible this would be, since Patsy is a sizable library that allows for fairly sophisticated formulas and it would be quite an endeavor to replicate all of the functionality. A few ways of doing this would be

1. Don't change any code, but document how Patsy can already be used to construct a dataframe that can then be passed to tabmat / glum. Warn that this involves creating a large dense matrix as an intermediate. See Twitter discussion: https://twitter.com/esantorella22/status/1447980727820296198
2. Have tabmat call patsy.dmatrix with "return_type = 'dataframe'", then call tabmat.from_pandas on the resulting pd.DataFrame. That would not be any more efficient than (1), but would just save the user a little typing and the need to install patsy. On the down side, it adds a dependency and may force creation of a very large dense matrix.
3. Support very simple patsy-like formulas without having patsy as a dependency or reproducing its full functionality. That would allow the user to designate which columns should be treated as categorical in a more natural way. See Twitter discussion: https://twitter.com/esantorella22/status/1447981081358184461
4. Make it so that any Patsy formula can be used to create a tabmat object -- I'm not sure how. Might be hard.