Fast matrix multiplication with diagonal matrices

Question

Samuel L. Polk 2021-2-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/755294-fast-matrix-multiplication-with-diagonal-matrices

评论： z cy 2022-7-28

采纳的回答： James Tursa

在 MATLAB Online 中打开

Let Wbe a large, sparse matrix. Let

and

be diagonal matrices of the same size. I would like to calculate

. However, these matrices are large enough that matrix multiplication is very expensive. I would like to speed up the calculation of L.

I know that computing L can be sped up by utilizing the fact that

and

are diagonal. For example, I know that I can compute

as follows.

diagD1 = diag(D1);  % diagonal of the matrix D1.
D1W = W.*diadD1; % Equivalent to multiplying the ith row of W by D(i,i). Yields D1*W.

My question is whether there is a similar exploitation of the diagonality of

that will allow me to avoid matrix multiplication to compute L.

Thank you.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Matt J 2021-2-26

编辑：Matt J 2021-2-26

在 MATLAB Online 中打开

It doesn't appear to be beneficial to avoid matrix multiplication, though I am surprised the latter is so slow (in R2020b).

N=1e5;
W=sprand(N,N,100/N);
d=rand(N,1);
D=spdiags(d,0,N,N);
tic
 L1=D*W*D;
toc
Elapsed time is 0.572682 seconds.
tic;
 L2=(d.*W).*d.';
toc
Elapsed time is 13.663614 seconds.

saskia leary 2022-3-20

(diag(D))'.*A works for right mulipltication - i.e.

A*D=(diag(D)).*A

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

James Tursa 2021-2-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/755294-fast-matrix-multiplication-with-diagonal-matrices#answer_632589

编辑：James Tursa 2021-2-25

在 MATLAB Online 中打开

Here is a mex routine to do this calculation. It relies on inputting the diagonal matrices as full vectors of the diagonal elements. It does not check for underflow to 0 for the calculations. A robust production version of this code would check for this and clean the sparse result of 0 entries, but I did not include that code here. It also does not check for inf or NaN entries. This could be made faster with parallel code such as OpenMP, but I didn't do that either.

/* File:  spdmd.c                                                                   */
/* Compile:  mex spdmd.c                                                            */
/* Syntax  C = spdmd(D1,M,D2)                                                       */
/* Does  C = D1 * M * D2                                                            */
/* where M  = double real sparse NxN matrix                                         */
/*       D1 = double real N element full vector representing diagonal NxN matrix    */
/*       D2 = double real N element full vector representing diagonal NxN matrix    */
/*       C  = double real sparse NxN matrix                                         */
/* Programmer:  James Tursa                                                         */
/* Date: 2/24/2021                                                                  */
/* Includes ----------------------------------------------------------- */
#include "mex.h"
#include <string.h>
/* Gateway ------------------------------------------------------------ */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize m, n, j, nrow;
    double *Mpr, *D1pr, *D2pr, *Cpr;
    mwIndex *Mir, *Mjc, *Cir, *Cjc;
/* Argument checks */
    if( nlhs > 1 ) {
        mexErrMsgTxt("Too many outputs");
    }
    if( nrhs != 3 ) {
        mexErrMsgTxt("Need exactly three inputs");
    }
    if (!mxIsDouble(prhs[1]) || !mxIsSparse(prhs[1]) || mxIsComplex(prhs[1])) {
        mexErrMsgTxt("2nd argument must be real double sparse matrix");
    }
    if( !mxIsDouble(prhs[0]) || mxIsSparse(prhs[0]) || mxIsComplex(prhs[0]) ||
	mxGetNumberOfDimensions(prhs[0]) != 2 || (mxGetM(prhs[0]) != 1 && mxGetN(prhs[0]) != 1)) {
        mexErrMsgTxt("1st argument must be real double full vector");
    }
    if (!mxIsDouble(prhs[2]) || mxIsSparse(prhs[2]) || mxIsComplex(prhs[2]) ||
	mxGetNumberOfDimensions(prhs[2]) != 2 || (mxGetM(prhs[2]) != 1 && mxGetN(prhs[2]) != 1)) {
	mexErrMsgTxt("3rd argument must be real double full vector");
    }
    m = mxGetM(prhs[1]);
    n = mxGetN(prhs[1]);
    if (m != n || mxGetNumberOfElements(prhs[0]) != n || mxGetNumberOfElements(prhs[2]) != n) {
        mexErrMsgTxt("Matrix dimensions must agree.");
    }
/* Sparse info */
    Mir = mxGetIr(prhs[1]);
    Mjc = mxGetJc(prhs[1]);
/* Create output */
    plhs[0] = mxCreateSparse( m, n, Mjc[n], mxREAL);
    
/* Get data pointers */
    Mpr = (double *) mxGetData(prhs[1]);
    D1pr = (double *) mxGetData(prhs[0]);
    D2pr = (double *) mxGetData(prhs[2]);
    Cpr = (double *) mxGetData(plhs[0]);
    Cir = mxGetIr(plhs[0]);
    Cjc = mxGetJc(plhs[0]);
/* Fill in sparse indexing */
    memcpy(Cjc, Mjc, (n+1) * sizeof(mwIndex));
    memcpy(Cir, Mir, Cjc[n] * sizeof(mwIndex));
/* Calculate result */
    for( j=0; j<n; j++ ) {
        nrow = Mjc[j+1] - Mjc[j];  /* Number of row elements for this column */
        while( nrow-- ) {
            *Cpr++ = *Mpr++ * (D2pr[j] * D1pr[*Cir++]);  
	}
    }
}