# How to make this matrix efficiently?

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.

• possible duplicate of What is the best way to multiply a diagonal matrix (in fortran) – Federico Poloni Jun 18 '15 at 12:10
• If you're using a C-type language that uses row-major storage order I doubt you'll do better than for (j=0; j<J; j++) {for (i=0; i<I; I++) { A[j][i] = delta[j] * B[j][i]; }};. Storage order and cache efficiency is going to be the determining factor for speed. – Doug Lipinski Jun 18 '15 at 13:09
• This was basically my original line of thought. But, computational mathematics isn't really my field, and I figured there might be a more clever way that doing the row scaling myself... – StevieP Jun 18 '15 at 14:08
• @FedericoPoloni see edit. – StevieP Jun 18 '15 at 15:07

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.

• I suppose I wasn't explicit enough about the specification of $\mathbf{A}$. As I've said, $\mathbf{A}$ is an I x J matrix. The i-th row of $\mathbf{A}$ is the J-dimensional vector, $\mathbf{x}_j$, scaled by the i-th component of the I-dimensional vector $\delta$. That is, $\mathbf{A}_{i,j}$ is, in pseudo-code, delta[i] * \mathbf{x}_{i}[j]... Does that make more sense? – StevieP Jun 19 '15 at 6:41
• @StevieP: Then the vectors in $\mathbf{A}$ should be transposed. The usual convention in mathematics is that $n$-vectors are $n \times 1$ matrices. In any case, yes, you could premultiply by a diagonal matrix, but it would probably be better for you to use BLAS level 1 primitives or a loop. – Geoff Oxberry Jun 19 '15 at 6:51
• @FedericoPoloni: Ah, yes, right. – Geoff Oxberry Jun 19 '15 at 11:14