MATLAB: Preserving numerical symmetry in large nxn matrix

matrix multiplicationsparse matricessymmetricsymmetric matrices

Let W be a sparse, symmetric nxn matrix and D a sparse, diagonal nxn matrix. These matrices are very large (n~250,000) so both must be stored in sparse format to analyze them.
I would like to compute the following matrix: . I compute L in MALTAB in the following way:
Dinv = sparse(1:n,1:n, 1./diag(D)); % theoretical inverse of D, stored in a sparse matrix
L = Dinv*W*D;
Theoretically, this matrix should be symmetric. However, it isn't when I compute it numerically. Is this because I am somehow computing L wrong? Or does it have to do with sparsity? Thank you.

Best Answer

Here is a mex routine to do this calculation. It relies on inputting the diagonal matrix as a full vector 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:  spdimd.c                                                                 */
/* Compile:  mex spdimd.c                                                          */
/* Syntax  C = spdimd(M,D)                                                         */
/* Does  C = D^-1 * M * D                                                          */
/* where M = double real sparse NxN matrix                                         */
/*       D = 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, *Dpr, *Cpr;
    mwIndex *Mir, *Mjc, *Cir, *Cjc;
/* Argument checks */
    if( nlhs > 1 ) {
        mexErrMsgTxt("Too many outputs");
    }
    if( nrhs != 2 ) {
        mexErrMsgTxt("Need exactly two inputs");
    }
    if (!mxIsDouble(prhs[0]) || !mxIsSparse(prhs[0]) || mxIsComplex(prhs[0])) {
        mexErrMsgTxt("1st argument must be real double sparse matrix");
    }
    if( !mxIsDouble(prhs[1]) || mxIsSparse(prhs[1]) || mxIsComplex(prhs[1]) ||
	mxGetNumberOfDimensions(prhs[1]) != 2 || (mxGetM(prhs[1]) != 1 && mxGetN(prhs[1]) != 1)) {
        mexErrMsgTxt("2nd argument must be real double full vector");
    }
    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);
    if( m!=n || mxGetNumberOfElements(prhs[1])!=n ) {
        mexErrMsgTxt("Matrix dimensions must agree.");
    }
/* Sparse info */
    Mir = mxGetIr(prhs[0]);
    Mjc = mxGetJc(prhs[0]);
/* Create output */
    plhs[0] = mxCreateSparse( m, n, Mjc[n], mxREAL);
    
/* Get data pointers */
    Mpr = (double *) mxGetData(prhs[0]);
    Dpr = (double *) mxGetData(prhs[1]);
    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++ * (Dpr[j] / Dpr[*Cir++]);  
	}
    }
}
Related Question