Bug when using min() with sparse matrices?

Hi,

I translated one of my Matlab routines to C and now I've encountered a really strange behaviour which seems to be a bug. M is a sparse matrix and I printed M, M', min(M, M') and min(full(M), full(M')), but all of them as full matrices to see the problem better:

http://dl.dropbox.com/u/12356013/weirdMatlab.png

The last matrix is what I want and what I'm supposed to already get by using min(M, M'), but for some reason min(M, M') shows some really weird behaviour. Any gueses what causes this?

Edit: Here's the mex file:

#include "mex.h"

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray *prhs[])
{
  mxArray *tempArray, *outputDist[1], *inputDist[2], *outputSort[2];
  double *tempVal, *pr, *sr, *sorted, *order;
  mwSize n, nzmax;
  mwIndex *irs, *jcs, i, j, k;

  // distance function call
  tempArray   = mxCreateDoubleMatrix(1, 1, mxREAL);
  tempVal     = mxGetPr(tempArray);
  inputDist[0] = prhs[0];
  inputDist[1] = tempArray;

  // read input matrix size
  n  = mxGetN(prhs[0]);

  // read real part of input matrix
  pr = mxGetPr(prhs[0]);

  // read parameter for k
  k = (mwIndex)*mxGetPr(prhs[1]);

  // create output sparse matrix
  nzmax   = (mwSize)(n*k);
  plhs[0] = mxCreateSparse(n, n, nzmax, mxREAL);
  sr      = mxGetPr(plhs[0]);
  irs     = mxGetIr(plhs[0]);
  jcs     = mxGetJc(plhs[0]);

  // go through data points
  for (i = 0; i < n; i++)
  {
    // calculate distance vector
    tempVal[0] = (double)(i+1);
    mexCallMATLAB(1, outputDist, 2, inputDist, "distEuclidean2");

    // sort vector
    mexCallMATLAB(2, outputSort, 1, &outputDist[0], "sort");
    sorted = mxGetPr(outputSort[0]);
    order  = mxGetPr(outputSort[1]);

    // assign values to sparse matrix
    jcs[i] = i*k;
    for (j = 0; j < k; j++)
    {
      irs[i*k+j] = (mwIndex)(order[j]-1);
      sr[i*k+j]  = sorted[j];
    }
  }
  jcs[n] = n*k;
}

And here is the called distEuclidean2:

function [ d ] = distEuclidean2( M, ii )
  d = sqrt(sum((repmat(M(:, ii), 1, size(M, 2)) - M) .^ 2, 1));
emd

[Ingo]

I should add, that creating M with the Matlab version of my routine causes absolutely no such problem.

[Andrew Newell]

How did you translate it to C?

[Ingo]

I used mex files and created the sparse matrix there. I can check the matrix created with the mex file in Matlab and it seems fine. But using min(M, M') on it causes this.

Of course I could just paste the complete mex source code, if you want me to? It's not that much, like 55-60 lines.

[Andrew Newell]

That might help - also the commands you used to create the mex files.

[Walter Roberson]

Suggest you edit the source in to your Question, as comments don't support formatting (or indentation even)

[Ingo]

I added the mex source code now. The function is called with a d-by-n matrix of n data points and a number k as the second parameter. It then calculates the n-by-n distance matrix, but only keeps the k smallest distances, all other values are set to zero (thus making it sparse).

I added distEuclidean2.m as well. I compile the mex file with "mex filename.c -largeArrayDims" and I'm running a 64-bit Matlab. One can then run it with

A = (1:5);

M = whateveryoucalledit(A, 3);

for example.

Answer 1

[James Tursa]

You might try using this to check the matrix:

href=""<http://www.mathworks.com/matlabcentral/fileexchange/20186-spok-checks-if-a-matlab-sparse-matrix-is-ok</a>>

Not sure if it works for 64-bit.

*** EDIT ***

Note that spok has been removed from the FEX. You can use my replacement function sarek instead:

https://www.mathworks.com/matlabcentral/fileexchange/158161-sarek-sparse-analyzer-real-et-komplex?s_tid=srchtitle

[Ingo]

Thanks. I could use it by modifying the installation to 64-bit. I knew about having explicit zero entries (didn't they they'd be a problem), but I also have row index problems. After using M = 1*M and M = M'' the min(M, M') command works. Now I know where the problem is!