Subtract corresponding nonzero element from each element in a vector

Ok, I want to create an M x M matrix A with the following.

Vec1=[0 1 1 0] Vec2=[1 3 5 7]

Then, matrix A has size M = size(Vec1) = size(Vec2) and all elements start at some constant value c.

I want to update the matrix A such that, for each index where Vec1 is not zero, create a matrix row where the row values are obtained by subtracting the value of Vec2 at that same index from each other value and getting the absolute values.

So in the above example, if I initialize matrix A as:

c c c c
c c c c
c c c c
c c c c

then, the first row stays the same (Vec 1 is zero) the second row becomes [abs(3-1), abs(3-3), abs(3-5), abs(3-7)], Third row becomes [abs(5-1), abs(5-3), abs(5-5), abs(5-7)]. Fourth row stays the same, so we end up with:

c c c c
2 0 2 4
4 2 0 2
c c c c

I would like to accomplish this without loops or possibly without even repmats and the like, since the matrices and vectors in question are huge.

Thanks!

Answer 1

[Sean de Wolski]

%EDIT again
Vec1=[0 1 1 0]; 
Vec1(~Vec1) = nan;
Vec2=[1 3 5 7];
C2 = bsxfun(@times,Vec1.',abs(bsxfun(@minus,Vec2.',Vec2)));
C2(isnan(C2)) = A(isnan(C2));

[Walter Roberson]

That doesn't take in to account the existing values in A, which are required to stay the same when the vec1 location is non-zero.

[CP]

This is close. It doesn't leave the elements of C unchanged when Vec1 is 0 though. It sets them to 0 instead.

[CP]

Sorry, mabye I wasn't clear. it isn't as simple as just setting everything to a constant. I want the values of C where Vec 1 is 0 to retain their previous values. It's true I initialize to a constant c, but the matrix is updated again and again, each time with a different Vec1.

[Sean de Wolski]

Okay, how about this one? Each time you repeat the update you save C2 as A? If not; then I'm still not clear...

[CP]

This seems to work. Thanks!

[CP]

Hmm I tried profiling it and I was surprised to find that the line C2(isnan(C2)) = A(isnan(C2)) is taking up 40% of the entire calculation time of my whole code (which includes a lot of numeric integration and the like), while the bsxfun line is taking 23% of the time.

If you can think of any more efficient ways of getting it done it would be appreciated, but apparently moving 1024x1024 values around from matrix to matrix is more of a performance hit than I thought.

[Sean de Wolski]

You could remove one call to isnan by defining an index matrix:

idx = isnan(C2);

C2(idx) = A(idx);

Answer 2

[Walter Roberson]

In order to avoid having at least one temporary matrix the same size as A, you would need to update A "in-place", which would require using a loop.

Possibly you might be able to do the whole calculation with a single bsxfun() call; I would have to think more about how to handle the vectorization efficiently in the face of the fact that bsxfun() expects user-provided functions to work on column vectors.

[CP]

I think a temporary matrix would be preferable to a loop in this case in terms of speed, so if you have a way to do it with a temporary matrix that could work.

[Sean de Wolski]

I wasn't able to figure out how do it all with one call to BSXFUN. I'd like to see it since it seems easily possible.

Answer 3

[Matt Fig]

Try this for speed. On my machine it is much faster than the double BSXFUN call. I assume random V1 (your Vect1). .

.

.

EDIT Putting the call to ABS inside the loop is faster(?).

A = 3*ones(length(V1));  % This is your starting matrix... (Example)
TMP = 1:length(V1);
TMP = TMP(logical(V1));

for ii = TMP
     A(ii,:) = abs(V2-V2(ii));
end

.

.

EDIT2

What was I thinking with that inner FOR loop? I was thinking I need to brew some more coffee...

[CP]

Well I tried it with the inner for loop and I could not tell the difference between having the code or not having the code (e.g. compared to the rest of the code it's a trivial fraction). Pretty great given that it was like 75% of the overall code with the earlier method. I'll update to this version. Thanks a lot :)

[Matt Fig]

Never underestimate the power of the FOR loop! You're welcome. ;-)