Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Fast matrix copying and array appending

Asked by nedo nodo on 30 Jan 2013

Hi,

I would like to know if there exists a fast way to copy a subpart of a matrix into a new matrix and to append an array to the same matrix. For instance, something like this: if true A=[A(2:end,2:end),B;B',b]; end where A is nxn matrix, B is nx1 array and b is a scalar.

Thank you Best Regards Salvo

0 Comments

nedo nodo

Products

No products are associated with this question.

2 Answers

Answer by Jan Simon on 30 Jan 2013
Edited by Jan Simon on 30 Jan 2013

I'm convinced that this is a good idea already:

A = [A(2:end, 2:end), B; B', b];

You save some micro-seconds by avoiding end, because the explicit definition if the size is faster:

[s1, s2] = size(A);
A = [A(2:s1, 2:s2), B; B', b];

In my measurements this is remarkably slower:

A(1:s1-1, 1:s2-1) = A(2:s1, 2:s2);
A(s1,     1:s2-1) = B';
A(1:s1-1, s2)     = B;
A(s1, s2)         = b;

This might save some temporary memory, but the indexing needs more time.

 

[EDITED]

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
  double *A, *B, *Bp, b, *R;
  mwSize i, n, n_m1;
    % Get inputs: 
    A = mxGetPr(prhs[0]);
    B = mxGetPr(prhs[1]);
    b = mxGetScalar(prhs[2]);
    n = mxGetM(prhs[0]);
    % Create outputs:
    plhs[0] = mxCreateDoubleMatrix(n, n, mxReal);
    R = mxGetPr(plhs[0]);
    Bp = B;
    n_m1 = n - 1;
    for (i = 0; i < n_m1; i++) {
      memcpy(R, A+1, n_m1 * sizeof(double));
      R[n_m1] = *Bp++;
      R += n;
      A += n;
    }
    memcpy(R, B, n_m1 * sizeof(double));
    R[n_m1] = b;
  }

Use this with great care and after exhaustive testing only! This is written in the forum's editor and neither compiled nor tested. Do you know how to compile a C-Mex file?

4 Comments

José-Luis on 30 Jan 2013

It depends whether you want your data in a single matrix. Depending on the application, the size of your matrix, your access patterns, etc., it might be more efficient to use a cell array. It is very difficult to know what will be faster without actually testing it. And even then, the results will only be true for your computer.

Jan Simon on 30 Jan 2013

@nedo nodo: I do not believe, that this assignment needs the most time in your program. Optimizing a not relevant part of the code does not lead to a significant acceleration of the whole program. So please use the profiler to find the bottlenecks, if you did not do this already.

You could create a C-MEX function for the creation of the matrix, but it is not sure if this is faster at all, because most of Matlab's indexing functions are very efficient already. If you find out, that this operation needs 90% of the overall processing time, I could offer to create a small C-code.

nedo nodo on 30 Jan 2013

@Jan Simone I have already done the profiling of my code and for this reason I have said that it is fundamental to improve this row. Your solution without the esplicit use of end is sligthly faster. Thank you

Jan Simon
Answer by Matt J on 30 Jan 2013
Edited by Matt J on 30 Jan 2013

If this really is the performance critical part of the code, you should consider redesigning so that A(2:end,2:end) gets prepended instead of appended. It leads to orders of magnitude speed difference, as shown below. The redesign should be trivial if you are doing likely operations like A*x. Just do A*[x(end);x(1:end-1)] instead.

    N=3000;
    A=rand(N);
    B=rand(N-1,1);
    b=1;

#########

    tic %Appending
     A=[A(2:end,2:end),B;B',b];
    toc
    %Elapsed time is 0.086694 seconds.

############

    tic%Prepending
      A(1)=b;
      A(2:end,1)=B;
      A(1,2:end)=B';
    toc
    %Elapsed time is 0.000218 seconds.

1 Comment

Matt J on 30 Jan 2013

You could also use my ProdCascade class so that even though you're prepending A, it will automatically pre- and post-permute things you want to multiply with A to give the same I/O behavior as if you had appended to A. The class introduces some overhead, but it's still close to three times as fast as your original approach when doing A*x operations:

tic %Appending
 A=[A(2:end,2:end),B;B',b];
 z1=A*x;
toc

Elapsed time is 0.079278 seconds.

P=circshift(speye(N),[1,0]);
AA=ProdCascade({P.',A,P});
tic;
    A(1)=b;
    A(2:end,1)=B;
    A(1,2:end)=B';
    AA{2}=A;
    z2=AA*x;
  toc;

Elapsed time is 0.027279 seconds.

Matt J

Contact us