1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ possible duplicate of What is the best way to multiply a diagonal matrix (in fortran) $\endgroup$ Jun 18, 2015 at 12:10
  • 2
    $\begingroup$ 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. $\endgroup$ Jun 18, 2015 at 13:09
  • $\begingroup$ 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... $\endgroup$
    – StevieP
    Jun 18, 2015 at 14:08
  • $\begingroup$ @FedericoPoloni see edit. $\endgroup$
    – StevieP
    Jun 18, 2015 at 15:07

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – StevieP
    Jun 19, 2015 at 6:41
  • 1
    $\begingroup$ @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. $\endgroup$ Jun 19, 2015 at 6:51
  • $\begingroup$ @FedericoPoloni: Ah, yes, right. $\endgroup$ Jun 19, 2015 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.