Preserving numerical symmetry in large nxn matrix

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.

Answer 1

[James Tursa]

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++]);  
	}
    }
}

Answer 2

[James Tursa]

Why do you think L should be symmetric? E.g.,

(1) L = D^-1 * W * D

(2) L^T = (D^-1 * W * D)^T = D^T * W^T * (D^-1)^T = D * W * D^-1

In (1) above, the D factors are applied to the columns of W and the D^-1 factors are applied to the rows of W.

In (2) above, the D factors are applied to the rows of W and the D^-1 factors are applied to the columns of W.

E.g., the L(1,2) element is (1/D(1)) * W(1,2) * D(2)

And the L(2,1) element is (1/D(2)) * W(2,1) * D(1) = (1/D(2)) * W(1,2) * D(1)

I don't see why you would think L should be symmetric.

P.S. You might be able to speed up that triple matrix multiply with a simple mex routine that expoits the diagonal D. Is D originally a full vector? It would actually simplify things for the mex routine if D was passed in as a full vector. Indexing for the individual element multiplies would be trivial. It would still work if D were passed in as a 2D sparse matrix but would only be as fast if the assumption of D being diagonal is made without checking.

[Samuel L. Polk]

Sorry, I was trying to make this problem as general as possible so that others can use it. But I see what you mean that my generality was making it a ill-posed problem. The matrices W and D that I am working with are such that is indeed symmetric theoretically. I think the issue I'm running into is numerical. Despite L being theoretically symmetric, it is not when I compute it, probably due to roundoff error.

To your answer post-script question, that's a really good point! D is indeed a full vector. I will implement this computation using a .mex file, as you suggested. Thank you for your thoughtful comment.

[James Tursa]

From what I have shown above, it seems the symmetry of L will only be guaranteed if D(i)=D(j) for all non-zero W(i,j). Is that what you have?

For the mex routine, have you ever coded with MATLAB format sparse matrices before? It can be a bit tricky. Let me know if you need help.

[Samuel L. Polk]

According to Wikipedia, ifA and B are symmetric, then is symmetric too iff Aand B commute. I.e., . This is the condition I have.

I haven't used sparse matrices much, so I definitely would be open to some help. One way I was thinking I could compute L without three matrix multiplications is the following:

d = diag(D);  % diagonal of the matrix D.
DinvW = W./d; % Equivalent to dividing the ith row of W by D(i,i). Yields D^{-1}W.
DinvWD = (d').*DinvW; % Equivalent to multiplying the ith row of DinvW by the ith column of D. Yields D^{-1}WD.

Do you have any thoughts on the above? I'm least confident that the last step will work. Thank you for your help!

[James Tursa]

I don't know how your commuting AB = BA applies in this case. As I pointed out, the only way the result can be symmetric is if a certain condition is applied to specific D elements for the non-zero W elements. Maybe you have this condition.

Regardless, I think your approach is a good one. It replaces the matrix multiplies with element-wise multiplies, but does seem to involve an extra temporary sparse matrix in the process. It will be interesting to see how this compares to the mex routine.