Matrixwise subtraction yields vector - how to reshape to original logical matrix?

Hi,

I am working with a large dataset for which I need to do multiple substractions between vectors, then summing across logically indexed values. Currently, I index the vector subtraction in a for loop, but it takes forever as I have to run that loop ~ 1.5M, several times, so I would like to move to matrixwise calculations.

Below is a short version of the crux of my question. If I have a matrix A and a matrix B, and I index the values I am interested in using a logical matrix (cmpMatrix), then I subtract B from A and want to get the mean per column. However, the logical indexing always reshapes into a vector, not a matrix:

A = rand(3);
B = rand(3);
cmpMatrix = A>=0.5 %in my situation, there are always unequal numbers of 0s and 1s per column
cmpMatrix =

  3×3 logical array

   1   0   1
   0   1   1
   1   0   0
   
%Issue: ans for the following line is a vector.
A(cmpMatrix)

ans =

    0.9612
    0.7870
    0.9437
    0.9813
    0.5702

I'd like to take the mean of all the values that are indexed per column, so I'd need to be able to backtrack where the values in the vector came from. Ideally, I'd be able to do something like:

mean(A(cmpMatrix) - B(cmpMatrix),2)

And receive a 1 x 3 vector where B was subtracted from A at each logical 1, and then the mean was taken across the matrix columns, so as to get a mean per column.

How do I achieve that?

Thank you for your help.

Answer 1

[Stephen23]

A = rand(3);
B = rand(3);
X = A>=0.5;
[~,C] = find(X);
M = accumarray(C,A(X)-B(X),[],@mean).'

M = 1×3

0.6699 0.0500 0.5998

[Michaela Warnecke]

Huh.. OK. Thanks for this answer - love the oneliner.

Profiler shows that in a trial with 10k selections (instead of 1.5M) a for loop is actually fast (29.93s) than the matrix solution (49.61s).

But this is still good to know - I'd been struggling with the "matrix to vector" issues for other items as well. Thank you!

[Dyuman Joshi]

Most likely culprit for that would be find().

Try this -

A = rand(3)

A = 3×3

0.7416 0.4844 0.4685 0.6996 0.2079 0.9116 0.7800 0.6658 0.7368

B = rand(3)

B = 3×3

0.1584 0.4105 0.0889 0.7934 0.5021 0.2475 0.3965 0.1938 0.7773

Y = meshgrid(1:size(A,1))

Y = 3×3

1 2 3 1 2 3 1 2 3

X = A>=0.5

X = 3×3 logical array
   1   0   0
   1   0   1
   1   1   1

M = accumarray(Y(X),A(X)-B(X),[],@mean).'

M = 1×3

0.2909 0.4720 0.3119

[Michaela Warnecke]

Thank you ! I tested it, but it's still 64s using your solution. Looks like the culprit is @mean, which I'm not sure I can do much about: In a call of 10,000 selections (loops), mean gets called 1,040,104 times (my variable "X" is 104 x 8193, so it's a calculation of 10,000 x 104 steps = 1,040104 calls to mean (across 8193 logical indeces)).

@mean is also the major time suck in the for-loop solution, despite more calls:

Matrix-wise solution (64s):

For-loop solution (27s):

[Voss]

@Michaela Warnecke: It would be interesting to see timings for all 5 approaches suggested so far, running in your actual code.

Here they are for random A and B of size 104x8193:

A = rand(104,8193);
B = rand(104,8193);

timeit(@()Stephen(A,B))

ans = 0.0411

timeit(@()Dyuman(A,B))

ans = 0.0397

timeit(@()Walter1(A,B))

ans = 0.0464

timeit(@()Walter2(A,B))

ans = 0.0436

timeit(@()Voss(A,B))

ans = 0.0149

function M = Stephen(A,B)
X = A>=0.5;
[~,C] = find(X);
M = accumarray(C,A(X)-B(X),[],@mean).';
end

function M = Dyuman(A,B)
[m,n] = size(A);
Y = repmat(1:n,m,1); % modification for non-square A
X = A>=0.5;
M = accumarray(Y(X),A(X)-B(X),[],@mean).';
end

function result = Walter1(A,B)
cmpMatrix = A>=0.5;
result = arrayfun(@(IDX) mean(A(cmpMatrix(:,IDX),IDX)-B(cmpMatrix(:,IDX),IDX)), 1:size(A,2));
end

function result = Walter2(A,B)
cmpMatrix = A>=0.5;
result = cellfun(@mean, mat2cell(A(cmpMatrix) - B(cmpMatrix), sum(cmpMatrix,1), 1)).'; % removing temp variable, and transposing for consistency with other methods
end

function result = Voss(A,B)
A_temp = A;
B_temp = B;
idx = A >= 0.5;
A_temp(~idx) = NaN;
B_temp(~idx) = NaN;
result = mean(A_temp-B_temp,1,'omitnan');
end

[Dyuman Joshi]

Comparing using isequal does not work for Voss's solution, but does for the rest.

Weird?

A = rand(104,8193);
B = rand(104,8193);

timeit(@()Stephen(A,B))

ans = 0.0426

timeit(@()Dyuman(A,B))

ans = 0.0399

timeit(@()Walter1(A,B))

ans = 0.0452

timeit(@()Walter2(A,B))

ans = 0.0435

timeit(@()Voss(A,B))

ans = 0.0129

S = Stephen(A,B);
D = Dyuman(A,B);
W1 = Walter1(A,B);
W2 = Walter2(A,B);
V = Voss(A,B);

isequal(S,D,W1,W2,V)

ans = logical
   0

isequal(S,D,W1,W2)

ans = logical
   1

function M = Stephen(A,B)
X = A>=0.5;
[~,C] = find(X);
M = accumarray(C,A(X)-B(X),[],@mean).';
end

function M = Dyuman(A,B)
[m,n] = size(A);
Y = repmat(1:n,m,1); % modification for non-square A
X = A>=0.5;
M = accumarray(Y(X),A(X)-B(X),[],@mean).';
end

function result = Walter1(A,B)
cmpMatrix = A>=0.5;
result = arrayfun(@(IDX) mean(A(cmpMatrix(:,IDX),IDX)-B(cmpMatrix(:,IDX),IDX)), 1:size(A,2));
end

function result = Walter2(A,B)
cmpMatrix = A>=0.5;
result = cellfun(@mean, mat2cell(A(cmpMatrix) - B(cmpMatrix), sum(cmpMatrix,1), 1)).'; % removing temp variable, and transposing for consistency with other methods
end

function result = Voss(A,B)
A_temp = A;
B_temp = B;
idx = A >= 0.5;
A_temp(~idx) = NaN;
B_temp(~idx) = NaN;
result = mean(A_temp-B_temp,1,'omitnan');
end

[Voss]

They're equal up to round-off error:

A = rand(104,8193);
B = rand(104,8193);

max(abs(Stephen(A,B)-Voss(A,B)),[],'all')

ans = 1.6653e-16

function M = Stephen(A,B)
X = A>=0.5;
[~,C] = find(X);
M = accumarray(C,A(X)-B(X),[],@mean).';
end

function M = Dyuman(A,B)
[m,n] = size(A);
Y = repmat(1:n,m,1); % modification for non-square A
X = A>=0.5;
M = accumarray(Y(X),A(X)-B(X),[],@mean).';
end

function result = Walter1(A,B)
cmpMatrix = A>=0.5;
result = arrayfun(@(IDX) mean(A(cmpMatrix(:,IDX),IDX)-B(cmpMatrix(:,IDX),IDX)), 1:size(A,2));
end

function result = Walter2(A,B)
cmpMatrix = A>=0.5;
result = cellfun(@mean, mat2cell(A(cmpMatrix) - B(cmpMatrix), sum(cmpMatrix,1), 1)).'; % removing temp variable, and transposing for consistency with other methods
end

function result = Voss(A,B)
A_temp = A;
B_temp = B;
idx = A >= 0.5;
A_temp(~idx) = NaN;
B_temp(~idx) = NaN;
result = mean(A_temp-B_temp,1,'omitnan');
end

[Dyuman Joshi]

I know and understand that they are round-off error; my point is why does it only occur in a single method and not in others, when all the methods use the same function (essentially the same operation)?

Answer 2

[Voss]

A_temp = A;
B_temp = B;
idx = A >= 0.5;
A_temp(~idx) = NaN;
B_temp(~idx) = NaN;
result = mean(A_temp-B_temp,1,'omitnan');

Answer 3

[Walter Roberson]

One way:

result = arrayfun(@(IDX) mean(A(cmpMatrix(:,IDX),IDX)-B(cmpMatrix(:,IDX),IDX)), 1:size(A,2))

Chances are that this will be slower than looping.

Another way:

temp = A(cmpMatrix) - B(cmpMatrix);
result = cellfun(@mean, mat2cell(temp, sum(cmpMatrix,1), 1));

It would not surprise me if this is slower than looping.