Preserving numerical symmetry in large nxn matrix
3 次查看(过去 30 天)
显示 更早的评论
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/529294/image.png)
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.
0 个评论
采纳的回答
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++]);
}
}
}
0 个评论
更多回答(1 个)
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 个评论
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!