Re: Feature proposal: weighted gram matrix / sandwich product
Christian Lorentzen via NumPy-Discussion <[email protected]>
| Newsgroups | gmane.comp.python.numeric.general |
|---|---|
| Message-ID | <[email protected]> |
Hi Kevin As I wrote in my original mail, BLAS has no routine for this operation, you need at least 2 BLAS (or numpy) operatios and I'd like to avoid the drawback of those solutions because this is a very performance critical and memory bound operation (X can be very large). The undertaking of implementing a GEMM-like sandwich product exploiting the symmetry, simd and thread-parallelism is difficult, but doable as demonstrated, e.g., here: https://github.com/Quantco/tabmat/blob/main/src/tabmat/ext/dense_helpers-tmpl.cpp <https://github.com/Quantco/tabmat/blob/main/src/tabmat/ext/dense_helpers-tmpl.cpp>. The question remains: Is this something that numpy could provide? Best Christian On 04.09.2025 14:19, Kevin Sheppard wrote: > A lighter alternative would be for a downstream project to use > Cython + SciPy to create a small binary module that mimics the numba > code but uses LAPACK directly. It isn't particularly difficult to > write the equivalent code, and we use a lot of this in statsmodels. > You then would only need to distribute a relatively small binary, and > the rest would be already included in the install. > > Kevin > > > On Tue, Sep 2, 2025 at 9:02 PM Christian Lorentzen > <[email protected]> wrote: > > As nice as numba is, it’s a very heavy dependency that many high > profile libraries want to avoid. > > Best > Christian > >> Am 18.08.2025 um 14:25 schrieb Kevin Sheppard >> <[email protected]>: >> >> >> This is a good place for numba if performance and memory use are >> a concern. >> >> Good NumPy >> 125 ms ± 2.42 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) >> >> Naive numba (memory efficient) >> 386 ms ± 1.83 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) >> >> Parallel numba (memory efficient, cpu bounded) >> 46.8 ms ± 1.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops >> each) >> >> Here is a gist that produced the timings: >> https://gist.github.com/bashtage/39c50b223580b9bc3f763a3b1be3e478 >> >> All were on a 12 core Ryzen. >> >> Kevin >> >> On Mon, Aug 18, 2025 at 12:35 PM Christian Lorentzen via >> NumPy-Discussion <[email protected]> wrote: >> >> If W is just a vector, i.e. my main use case, than what you >> write is what I wanted to imply and actually how I >> implemented it in several cases. >> The point is that X.T * weights[None, :] still doubles memory >> (same size as X) and this solution doesn't exploit the symmetry. >> >> Best >> Christian >> On 17.08.2025 22:13, Kevin Jacobs wrote: >>> Why not: >>> >>> XTWX = (X.T * weights[None, :]) @ X >>> >>> if `weights` is a vector? Allocating `diag(weights)` is the >>> big memory and CPU hog at least for my generally long and >>> skinny X matrices. >>> >>> -Kevin >>> >>> >>> On Fri, Aug 15, 2025 at 10:08 AM Christian Lorentzen via >>> NumPy-Discussion <[email protected]> wrote: >>> >>> Dear numpy community >>> >>> I would like to propose a mew feature in np.linalg to >>> compute weighted gram matrix (sandwich product), see >>> https://github.com/numpy/numpy/issues/29559. >>> >>> >>> Proposed new feature or change: >>> >>> The solvers for many important statistical models like >>> Generalized Linear Models often need to compute |X.T @ W >>> @ X| where |X| is a 2-dimensional array of features >>> (features in columns, observations in rows) and |W| is >>> very often a diagonal array of (current) weights. This >>> computation is usually the main computational bottleneck. >>> >>> It would be great if numpy could provide an efficient >>> implementation of it, e.g. >>> |np.linalg.sandwicht_product(X, weight=w)|. >>> >>> Why numpy? It seems even more unrealistic to me, to get >>> it into BLAS implementations. Numpy has support for SIMD >>> via Highways. >>> >>> >>> Computational alternatives >>> >>> Drawbacks of |(X.T * diag(W)) @ X|: >>> >>> * Additional memory allocation to compute |X.T @ >>> diag(W)| , same size as (usually large) |X|. >>> * The result is symmetric, but this fact is not used. >>> So at least a factor of 2 is possible. >>> Note that without weights, |X.T @ X| uses BLAS syrk. >>> But the weights are crucial. >>> >>> Drawback of |Z = np.sqrt(diag(W))[:, None] * X| and then >>> |Z @ Z|: >>> >>> * Additional memory allocation for |Z|, same size as >>> (usually large) |X|. >>> * Taking square roots for possibly large |diag(W)| >>> >>> >>> Additional information >>> >>> https://github.com/Quantco/tabmat has an implementation >>> of it with XSIMD. >>> >>> >>> Best >>> Christian Lorentzen >>> >>> _______________________________________________ >>> NumPy-Discussion mailing list -- [email protected] >>> To unsubscribe send an email to >>> [email protected] >>> https://mail.python.org/mailman3//lists/numpy-discussion.python.org >>> Member address: [email protected] >>> >>> >>> _______________________________________________ >>> NumPy-Discussion mailing list [email protected] >>> To unsubscribe send an email [email protected] >>> https://mail.python.org/mailman3//lists/numpy-discussion.python.org >>> Member address:[email protected] >> _______________________________________________ >> NumPy-Discussion mailing list -- [email protected] >> To unsubscribe send an email to [email protected] >> https://mail.python.org/mailman3//lists/numpy-discussion.python.org >> Member address: [email protected] >> _______________________________________________ NumPy-Discussion mailing list -- [email protected] To unsubscribe send an email to [email protected] https://mail.python.org/mailman3//lists/numpy-discussion.python.org Member address: [email protected]