All combinations of differences

Hey, I have the following problem:

I have two vectors that each are made up of a series of 3d sub-vectors, so something like: [v1; v2; v3] and [w1; w2], where vi and wi are column vectors. In my case, there is a very large number of subvectors vi and wi.

What I would like is to have some array of all combinations of differences between wi and vj, ie something like: [v1-w1, v1-w2, v2-w1, v2-w2, v3-w1, v3-w2].

Any ideas would be very much appreciated.

PS: In fact my actual problem is somewhat different. I have a program that is running very slowly, and I would like to try to vectorize it with the hopes of speeding it up (although maybe it is just very slow because it involves a lot of computation). The program calculates the energy between two chains of atoms. To do that it takes in two 3d arrays of different sizes of the following format:

Rot1_DevK_Atom1_x / Rot1_DevK_Atom2_x / ... / Rot1_DevK_AtomM_x

... / ... / ... / ...

Rot1_Dev3_Atom1_x / Rot1_Dev3_Atom2_x / ... / Rot1_Dev3_AtomM_x

Rot1_Dev2_Atom1_x / Rot1_Dev2_Atom2_x / ... / Rot1_Dev2_AtomM_x

Rot1_Dev1_Atom1_x | Rot1_Dev1_Atom2_x | ... | Rot1_Dev1_AtomM_x

Rot1_Dev1_Atom1_y | Rot1_Dev1_Atom2_y | ... | Rot1_Dev1_AtomM_y

Rot1_Dev1_Atom1_z | Rot1_Dev1_Atom2_z | ... | Rot1_Dev1_AtomM_z

Rot2_Dev1_Atom1_x | Rot2_Dev1_Atom2_x | ... | Rot2_Dev1_AtomM_x

Rot2_Dev1_Atom1_y | Rot2_Dev1_Atom2_y | ... | Rot2_Dev1_AtomM_y

Rot2_Dev1_Atom1_z | Rot2_Dev1_Atom2_z | ... | Rot2_Dev1_AtomM_z

... | ... | ... | ...

RotN_Dev1_Atom1_x | RotN_Dev1_Atom2_x | ... | RotN_Dev1_AtomM_x

RotN_Dev1_Atom1_y | RotN_Dev1_Atom2_y | ... | RotN_Dev1_AtomM_y

RotN_Dev1_Atom1_z | RotN_Dev1_Atom2_z | ... | RotN_Dev1_AtomM_z

N runs over the number of 'ROTamers' (that are different basic conformations of the sidechains), K runs over different DEViants (those are for each rotamer a set of slightly modified structures) and M runs over the number of atoms in each chain. N, K and M are different for the two sidechains.

Now to get the energy between two sidechains I need to calculate the energy between all combinations of atoms. I would like to do that over all combinations of rotamers and deviants. If I do that over a large number of randomly generated rotamers and deviants, involving Boltzmann distribution, that should give me an estimate of the average energy between the two sidechains, given that they are left to assume any of many random conformations.

My approach to vectorizing this problem was to reshape the 3d arrays into vectors and apply the above operation which I have yet to find, but if you have a different idea of how to speed this problem up other than doing 3 for loops, please let me know! My current program looks something like this:

    for n = 1:numbCombRotamers
        
        nIndex1 = ((combinationsRotamers(n,1)-1)*3+1):(combinationsRotamers(n,1)*3);
        nIndex2 = ((combinationsRotamers(n,2)-1)*3+1):(combinationsRotamers(n,2)*3);
        
        energyContributionDeviant = zeros(1,numbCombDeviants);
        for m = 1:numbCombDeviants
            
            energyContributionAtomPair = zeros(1,numbCombAtoms);
            for k = 1:numbCombAtoms
                
                R1 = DeviantsPositions1(nIndex1, combinationsAtoms(k,1), combinationsDeviants(m,1)); %VERY slow
                R2 = DeviantsPositions2(nIndex2, combinationsAtoms(k,2), combinationsDeviants(m,2)); %VERY slow
                
                energyContributionAtomPair(k) = input_EnergyLJ*(input_R0^12/(sum((R1-R2).^2))^6 ...
                    - 2*input_R0^6/(sum((R1-R2).^2))^3);
                                
            end
            
            energyContributionDeviant(m) = ...
                sum(energyContributionAtomPair);
            
        end
        
        energyContributionRotamer(n) = sum(energyContributionDeviant.*...
            exp(-energyContributionDeviant/gasConst/input_temperature))/...
            sum(exp(-energyContributionDeviant/gasConst/input_temperature));
               
    end

Here, numbCombDeviants is all the combinations of Deviants between the two sidechains, numbCombAtoms the combination of all atoms, combinationsThing are all combinations of things (rotamers, atoms or deviants) saved as a n*2 matrix. DeviantsPositionsN are the input arrays as above. I have timed the program and by far, far the very much worst part in terms of how much time it takes (by at least 2 orders of magnitude more than anything else) are the two rows defining R1 and R2.

Sorry for the long question and I hope my problem is clear!

Answer 1

[Steven Lord]

In your code it's not clear whether the identifiers you're using refer to data or functions. But to answer the question, if you're using release R2016b or later use implicit expansion.

% Sample data
A = randi(20, 10, 3);
B = randi(20, 10, 2);

% Preallocate
nA = size(A, 2);
nB = size(B, 2);
C = zeros(size(A, 1), nA*nB);
beforeStartOfNext = 0;
for whichColumnOfB = 1:nB
    C(:, beforeStartOfNext+(1:nA)) = A + B(:, whichColumnOfB);
    beforeStartOfNext = beforeStartOfNext + nA;
end

Alternately if you're okay with storing your results as a 3-D array, it's even easier. The sum of columns of A with the second column of B (for example) is in the second page of D. Remember that numeric arrays in MATLAB have implicitly stored singleton dimensions past ndims, so even though B is 2-dimensional we can permute it to make it 3-dimensional.

B3 = permute(B, [1 3 2]); % Makes it 10-by-1-by-2
D = A + B3; % The result is 10-by-3-by-2

isequal(C(:, 4:6), D(:, :, 2)) % true