Re: Feature proposal: weighted gram matrix / sandwich product

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]