Preserving numerical symmetry in large nxn matrix

1 次查看（过去 30 天）
Samuel L. Polk
Samuel L. Polk 2021-2-23
编辑： James Tursa 2021-2-25
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.

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

采纳的回答

James Tursa
James Tursa 2021-2-25
编辑：James Tursa 2021-2-25
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++]);
}
}
}

更多回答（1 个）

James Tursa
James Tursa 2021-2-23
编辑：James Tursa 2021-2-23
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.
  4 个评论
显示 2更早的评论
Samuel L. Polk
Samuel L. Polk 2021-2-24
编辑：James Tursa 2021-2-25
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
James Tursa 2021-2-25
编辑：James Tursa 2021-2-25
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by