Most Efficient Way of Using mexCallMATLAB in Converting Double* to mxArray*

Question

I am writing a MEX code in which I need to use pinv function. I am trying to find a way to pass the array of type double to pinv using mexCallMATLAB in the most efficient way. Let's for the sake of example say the array is named G and its size is 100.

double *G = (double*) mxMalloc( 100 * sizeof(double) );

where

G[0] = G11; G[1] = G12;
G[2] = G21; G[3] = G22;

Which means every four consecutive elements of G is a 2×2 matrix. G stores 25 different values of this 2×2 matrix.

I should note that these 2×2 matrices are not well-conditioned and they may contain all zero in their element. How can I use pinv function to calculate the pseudoinverse in the elements of G? For example, how can I pass the array to mexCallMATLAB in order to calculate the pseudoinverse of the first 2×2 matrix in G?

I thought of the following approach:

mxArray *G_PINV_input     = mxCreateDoubleMatrix(2, 2, mxREAL);
mxArray *G_PINV_output    = mxCreateDoubleMatrix(2, 2, mxREAL);
double  *G_PINV_input_ptr = mxGetPr(G_PINV_input);

memcpy( G_PINV_input_ptr, &G[0], 4 * sizeof(double));
mexCallMATLAB(1, G_PINV_output, 1, G_PINV_input, "pinv");

I am not sure how good this approach is. Copying the values is not economical at all because the total number of elements in G in my actual application is large. Is there anyway to skip this copying?

Answer 1

Here is my implementation of the MEX-function:

my_pinv.cpp

#include "mex.h"

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
    // validate arguments
    if (nrhs!=1 || nlhs>1)
        mexErrMsgIdAndTxt("mex:error", "Wrong number of arguments");
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]))
        mexErrMsgIdAndTxt("mex:error", "Input isnt real dense double array");
    if (mxGetNumberOfElements(prhs[0]) != 100)
        mexErrMsgIdAndTxt("mex:error", "numel() != 100");

    // create necessary arrays
    mxArray *rhs[1], *lhs[1];
    plhs[0] = mxCreateDoubleMatrix(100, 1, mxREAL);
    rhs[0] = mxCreateDoubleMatrix(2, 2, mxREAL);
    double *in = mxGetPr(prhs[0]);
    double *out = mxGetPr(plhs[0]);
    double *x = mxGetPr(rhs[0]), *y;

    // for each 2x2 matrix
    for (mwIndex i=0; i<100; i+=4) {
        // copy 2x2 matrix into rhs
        x[0] = in[i+0];
        x[2] = in[i+1];
        x[1] = in[i+2];
        x[3] = in[i+3];

        // lhs = pinv(rhs)
        mexCallMATLAB(1, lhs, 1, rhs, "pinv");

        // copy 2x2 matrix from lhs
        y = mxGetPr(lhs[0]);
        out[i+0] = y[0];
        out[i+1] = y[1];
        out[i+2] = y[2];
        out[i+3] = y[3];

        // free array
        mxDestroyArray(lhs[0]);
    }

    // cleanup
    mxDestroyArray(rhs[0]);
}

Here is a baseline implementation in MATLAB so that we can verify the results are correct:

my_pinv0.m

function y = my_pinv0(x)
    y = zeros(size(x));
    for i=1:4:numel(x)
        y(i:i+3) = pinv(x([0 1; 2 3]+i));
    end
end

Now we test the MEX-function:

% some vector
x = randn(100,1);

% MEX vs. MATLAB function
y = my_pinv0(x);
yy = my_pinv(x);

% compare
assert(isequal(y,yy))

---

EDIT:

Here is an another implementation:

my_pinv2.cpp

#include "mex.h"

inline void call_pinv(const double &a, const double &b, const double &c,
                      const double &d, double *out)
{
    mxArray *rhs[1], *lhs[1];

    // create input matrix [a b; c d]
    rhs[0] = mxCreateDoubleMatrix(2, 2, mxREAL);
    double *x = mxGetPr(rhs[0]);
    x[0] = a;
    x[1] = c;
    x[2] = b;
    x[3] = d;

    // lhs = pinv(rhs)
    mexCallMATLAB(1, lhs, 1, rhs, "pinv");

    // get values from output matrix
    const double *y = mxGetPr(lhs[0]);
    out[0] = y[0];
    out[1] = y[1];
    out[2] = y[2];
    out[3] = y[3];

    // cleanup
    mxDestroyArray(lhs[0]);
    mxDestroyArray(rhs[0]);
}

void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
{
    // validate arguments
    if (nrhs!=1 || nlhs>1)
        mexErrMsgIdAndTxt("mex:error", "Wrong number of arguments");
    if (!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) || mxIsSparse(prhs[0]))
        mexErrMsgIdAndTxt("mex:error", "Input isnt real dense double array");
    if (mxGetNumberOfElements(prhs[0]) != 100)
        mexErrMsgIdAndTxt("mex:error", "numel() != 100");

    // allocate output
    plhs[0] = mxCreateDoubleMatrix(100, 1, mxREAL);
    double *out = mxGetPr(plhs[0]);
    const double *in = mxGetPr(prhs[0]);

    // for each 2x2 matrix
    for (mwIndex i=0; i<100; i+=4) {
        // 2x2 input matrix [a b; c d], and its determinant
        const double a = in[i+0];
        const double b = in[i+1];
        const double c = in[i+2];
        const double d = in[i+3];
        const double det = (a*d - b*c);

        if (det != 0) {
            // inverse of 2x2 matrix [d -b; -c a]/det
            out[i+0] =  d/det;
            out[i+1] = -c/det;
            out[i+2] = -b/det;
            out[i+3] =  a/det;
        }
        else {
            // singular matrix, fallback to pseudo-inverse
            call_pinv(a, b, c, d, &out[i]);
        }
    }
}

This time we compute the determinant of the 2x2 matrix, if is non-zero, we calculate the inverse ourselves according to:

Otherwise we fallback to invoking PINV from MATLAB for the pseudo-inverse.

Here is quick benchmark:

% 100x1 vector
x = randn(100,1);           % average case, with normal 2x2 matrices

% running time
funcs = {@my_pinv0, @my_pinv1, @my_pinv2};
t = cellfun(@(f) timeit(@() f(x)), funcs, 'Uniform',true);

% compare results
y = cellfun(@(f) f(x), funcs, 'Uniform',false);
assert(isequal(y{1},y{2}))

I get the following timings:

>> fprintf('%.6f\n', t);
0.002111   % MATLAB function
0.001498   % first MEX-file with mexCallMATLAB
0.000010   % second MEX-file with "unrolled" matrix inverse (+ PINV as fallback)

The error is acceptable and within machine precision:

>> norm(y{1}-y{3})
ans =
   2.1198e-14

You could also test the worst case, when many of the 2x2 matrices are singular:

x = randi([0 1], [100 1]);

Answer 2

You don't need to allocate the output. Just make the pointer and let pinv create the mxArray automatically.

mxArray *lhs;

Then just use & like,

mexCallMATLAB(1, &lhs, 1, &rhs, "pinv");