Re: Feature proposal: weighted gram matrix / sandwich product

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 to [email protected]
https://mail.python.org/mailman3//lists/numpy-discussion.python.org
Member address: [email protected]