diagouter()

function out=diagouter(A,B)

Input:

• A,B: $n \times m$ matrices

Output:

• out: ${\mbox{diag\,}}\left(A B^\dagger\right)$

The expression for the density in Eq. (1) is in the form of taking as a column vector the diagonal elements of the ``outer product'' of two matrices,

$$\vec c ={\mbox{diag\,}}\left( A B^\dagger \right).$$ (10)

Because $A$ and $B$ are of dimension $\prod S_k\,\times\,N_s$, it is extremely wasteful (of time and memory) to form directly the matrix $\prod S_k \times \prod S_k$ matrix $A B^\dagger$ only to then take its diagonal elements. In this case, it is better to provide our own function, diagouter(), to perform the using other BLAS operations. In terms of components Eq. (10) is

$$c_i = \sum_n A_{i,n} B^*_{i,n}.$$

Thus, we can perform this operation by first forming the matrix elements $C_{i,n} \equiv A_{i,n} B^*_{i,n}$ using octave's ``.*'' and ``conj()'' operators, and then summing along the rows (second index). In octave notation, this becomes simply

c=sum(A.*conj(B),2);

Use this approach to produce a function of the above prototype in the file ``diagouter.m''.