Diagprod()

function out=Diagprod(a,B)

Input:

• a: $n \times 1$ column vector
• B: $n \times m$ matrix

Output:

• out: $\left({\mbox{Diag\,}}\vec a\right) B$

The expression for the gradient in Eq. (3) ultimately involves products of the form

$$C=\left({\mbox{Diag\,}}\vec a \right) B .$$ (11)

Here, ${\mbox{Diag\,}}$ takes a vector of length $\prod S_k$ and forms a very large, diagonal $\prod S_k \, \times \, \prod S_k$ matrix, which then multiples the matrix $B$. Again, this direct evaluation of the expression is extremely wasteful of both time and space, and so we shall provide an function which performs the operation in terms of a more efficient selection of BLAS routines.

In terms of components, Eq. (11) becomes

$$C_{i,n}=\sum_j a_j \delta_{j,i} B_{i,n} = a_i B_{i,n}.$$

Thus, each column of $C$ may be computed independently as the ``.*'' product of $\vec a$ with the corresponding column of $B$, a series of BLAS1 operations. Alternately, by using a little more memory, the same operation may be carried out with BLAS2 operations (which in this case run about 2 times faster) by first forming a matrix of the same size as $B$ with copies of $\vec a$ in each column and then taking the ``.*'' product, c=(a*ones(1,size(B,2))).*B. Using this strategy provide a function of the above prototype in the file ``Diagprod.m''.