Apply an operation/function to each 2D slice of a 3D matrix when the size of the slice varies?

I have a matrix X of dimensions 1000,75 representing 1000 observations of 75 variables. And I also have a matrix. P, of logical values with dimensions 128,75 representing 128 different combinations of the 75 variables in X (either they should be included, 1, or not, 0).

What I want to do is calculate the pseudo inverse for each combination of observations and variable combination using the pinv command.

First I did it by placing each possible X(:,P(i,:)) combination in a cell array of length i and use cellfun like so:

for i = 128 : -1 : 1
CellX{i} = X(:,P(i,:));
end
Inverse = cellfun(@(x) pinv(x),CellX,'UniformOutput',false);

but its slow so I would like to find something else, preferably by using the gpu. Making it a bit more complicated is the fact that the number of logical ones and zeroes in P (what I mean is that the sum of each row in P) is not the same in each row so the output size of pinv will vary a bit.

I was thinking that I could do a repmat of X so it gets dimensions 1000,75,128 and then use some command in combination with the logical indexing of P to do the calculation for each of the 128 2D matrices in the 3D matrix in one go, but I can’t quite figure out how to do it without looping though the third dimension.

Does anyone have an idea on how to do it? Arrayfun, pagefun etc.?

Thanks