How to make this matrix efficiently?

Question

Originally posted on stats.stackexchange, I'll pair the post down to something a bit more general.

Suppose I have vectors $\{\mathbf{\delta}, \mathbf{x}_1, \ldots, \mathbf{x}_J\}$, where $\delta \in \mathbb{R}^{J}$ and $\mathbf{x}_i \in \mathbb{R}^{I}$. Furthermore, let $\mathbf{A} \in \mathbb{R}^{I\times J}$ be such that $$ \mathbf{A} = \begin{pmatrix} \delta_1 \mathbf{x}_1 \\ \delta_2 \mathbf{x}_2 \\ \vdots \\ \delta_J \mathbf{x}_J \end{pmatrix} $$

Is there any sequence of BLAS-optimized operations that can make $\mathbf{A}$ from my set of vectors?

To be concrete, I understand that if I have the matrix $$ \mathbf{B} = \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_J \end{pmatrix} $$ then $$ \text{diag}(\delta) \, \mathbf{B} := \begin{pmatrix} \delta_1 & 0 & \cdots & 0 \\ 0 & \delta_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \delta_j \end{pmatrix} \mathbf{B} = \mathbf{A} $$ But it's not clear to me if making $\text{diag}(\delta)$ is an efficient operation, or even if the multiplication between $\text{diag}(\delta)$ and $\mathbf{B}$ is efficient.

Edit: I noticed that a similar, fortran, specific question was asked and answered. However that question addresses efficient ways to multiply diagonal matrices (in fortran, no less). Whereas, my question references the use of diagonal matrices as one way in which the solution can be approached. My thoughts are that this question is a bit more general than the aforementioned one.

Answer 1

Suppose I have vectors $\{\mathbf{\delta}, \mathbf{x}_1, \ldots, \mathbf{x}_J\}$, where $\delta \in \mathbb{R}^{J}$ and $\mathbf{x}_i \in \mathbb{R}^{I}$. Furthermore, let $\mathbf{A} \in \mathbb{R}^{I\times J}$ be such that $$ \mathbf{A} = \begin{pmatrix} \delta_1 \mathbf{x}_1 \\ \delta_2 \mathbf{x}_2 \\ \vdots \\ \delta_J \mathbf{x}_J \end{pmatrix} $$

There's some confusion about notation here. As written, $\mathbf{A} \in \mathbb{R}^{IJ}$; that is, is has one column and $IJ$ rows, instead of $I$ rows and $J$ columns.

Is there any sequence of BLAS-optimized operations that can make $\mathbf{A}$ from my set of vectors?

What I think you're proposing essentially reduces to

\begin{align} \begin{pmatrix} \mathbf{x}_{1} & \ldots & \mathbf{x}_{J} \end{pmatrix} \begin{pmatrix} \delta_1 & 0 & \cdots & 0 \\ 0 & \delta_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \delta_j \end{pmatrix}, \end{align}

followed by reshaping the result. I would not do that.

Instead, dscal each of the $\mathbf{x}_{i}$ by their respective $\delta_{i}$ in-place. As DougLipinski points out, this is basically just a for loop (or do loop in Fortran). A tuned BLAS implementation might do some loop unrolling.

To be concrete, I understand that if I have the matrix $$ \mathbf{B} = \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \\ \vdots \\ \mathbf{x}_J \end{pmatrix} $$ then $$ \text{diag}(\delta) \, \mathbf{B} := \begin{pmatrix} \delta_1 & 0 & \cdots & 0 \\ 0 & \delta_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \delta_j \end{pmatrix} \mathbf{B} = \mathbf{A} $$ But it's not clear to me if making $\text{diag}(\delta)$ is an efficient operation, or even if the multiplication between $\text{diag}(\delta)$ and $\mathbf{B}$ is efficient.

The memory movement isn't worth it; also, the multiplication you've proposed isn't conformal, so it would not work. Likely, you intended to do what I show above.