Calculate the mean in a single operation with matrices

Given 2 matrices M and P: I want to calculate the weighted geometric mean. For a single column vector p i would do like these:

M = rand(3,4);
P = (perms(1:3))';
p_single = P(:,1);                             % And iterate with a for loop through P matrix columns
result = prod(M.^p_single).^(1/sum(p_single)); % This gives me a single line

I want to get Result matrix where each row represents a Mean vector with the corresponding column from P matrix.

How can I do this without iteration in a single or few lines of code ?

EDIT (rewrote the complete loop):

M = rand(3,4);
P = (perms(1:3))';
[~,sz_C] = size(P);
result = [];
for k = 1: sz_C
  p_single = P(:,k);
  mean_single = prod(M.^p_single).^(1/sum(p_single));
  result = [result; mean_single];
end

This code appends to the 'result' matrix every calculated geometric mean combination for P columns.

Is it possible to make it one liner or to eliminate the loop for speed (in case of big data set) ?

[Rik]

Except that your code doesn't work. If you show how you would do it in a complete working loop, we might be able to help you do it without a loop.

[Marco Salvo Cottone]

Just an error typing here the code, corrected.

[dpb]

Preallocate and populate instead of augment...since the sizes of M and P aren't commensurate in that size(result) is

size(P,2) x size(M,2)

where P,M aren't directly commensurate since size(P,2)-->factorial(M,1). This means if M gets very big you're going to need a whole bunch of memory real soon now...just the output of perms for N is

>> factorial(10)*8/1024/1024
ans =
            27.685546875
>> factorial(11)*8/1024/1024
ans =
           304.541015625
>>

[Marco Salvo Cottone]

M is static data it will stay small or at least not that big, what will really grow here (when I refer to big data set) is W matrix here named P and the Result matrix will be as you wrote size(P,2) x size(M,2).

I agree I can preallocate:

result = zeros(size(P,2), size(M,2));

But still the for loop problem remains.

Answer 1

[dpb]

mngeo=geomean(M.*W);

where M, W are the matrix and associated weights matrix

[Marco Salvo Cottone]

Not correct, I updated the code with the whole loop to explain the process, your code gives error, which it shouldn't be there:

Matrix dimensions must agree.

[dpb]

That code isn't really a geometric mean; it's got some similarity to one but that's about all...

[Marco Salvo Cottone]

Every row in M must be considered a single value cause it represents a single data in my case, that's why I apply same weight to each element in a row. It may look different but it's actually the same.

[dpb]

Same as what? I can see randomizing weights as you're apparently doing but fail to see how having a value for every element of M qualifies as somehow averaging. Looks to me like for this case there should be 1xlength(P) end results.

I can't (yet, anyway) wrap my head around what it is that you think you're calculating; perhaps if you more fully explained what is the underlying problem it would become more clear and that could also lead to better coding solution.

[Marco Salvo Cottone]

Each row of M represents a Color expressed as reflectance value at a distance of 10nm each in the visible spectrum of light. I wrote a simple example to make it more easy to follow.

Real M matrix is like 6x36 and Weights/Parts matrix alot bigger, to simulate my requirements think about a W matrix like permutations with repetition with 10^6 columns. (In my example I used different permutation method, but it still applies to the theory, not the function).

To find a subractive color mixing recipe for a color I must create the Weighted Geometric Mean of the 6 rows of M (6 colors at a time) for 10^6 times with the columns of P (each column represent the proportion of parts of those 6 colors), the resulting matrix will give me the mixes, those mixes will be compared to a reference color (the recipe color) after being converting to Lab color space.

You can read it more at Subtractive Color Mixture Computation.

Currently the calculation is too slow because of the loop.

[dpb]

"Currently the calculation is too slow because of the loop."

I don't think it's the loop itself that is the bottleneck nearly as much as that you're augmenting the output array every pass instead of preallocating and filling. If it really is calculating what you want and you're absolutely sure of that, fix the above and then see if it isn't adequate--loops are much more well-optimized now in ML than in years gone by but the reallocation within the loop is a real killer.

I'll have to come back and read up a little more on the actual computation; sorry I don't have the time at the moment to do more than just make a quick comment.

[Marco Salvo Cottone]

No I tested, it's pretty slow, even after preallocation.

Preallocation:

result = zeros(size(P,2), size(M,2));

Copying the data:

result(k,:) = mean_single;

It's the huge number of loops the problem.

Thanks for your efforts, are much appreciated.