Using duplicate linear indices for updating a single position in a matrix

1
votes

I am currently working on a vectorizing code that assembles a big sparse matrix where given linear indices and a matrix of increments, I would like all values in the increments matrix that coincide to the same linear index to be all accumulated together so that the position at this linear index in the matrix get updated.

In order to do that, I have came up with two structures.

  • GlobalPoint has is in reference to a big sparse zeros matrix A. GlobalPoint is a n x 4 matrix that contains linear indices to access A.
  • LM has the increments to be added in the big sparse matrix

For example:

GlobalPoint

1   730 10  739
83  2   82  1
83  812 92  821
165 84  164 83
165 894 174 903

LM

-0,531249999999875  0,531249999999875   0,531249999999875   -0,531249999999875
-0,531250000000125  0,531250000000125   0,531250000000126   -0,531250000000126
-0,343749999999875  0,343749999999875   0,343749999999875   -0,343749999999875
-0,249999999999812  0,249999999999812   0,250000000000332   -0,250000000000332
-0,249999999999909  0,249999999999909   0,249999999999909   -0,249999999999909

Initially A contains all zeroes. To perform what I would like, I have done the following update on A:

A(GlobalPoint) = A(GlobalPoint) + LM;

The linear index of a means location (1,1) in the sparse matrix A. Since the lear index of 1 shows up in the first and second row in the first column ((1,1) and (2,1), the expected result should be: (-0.531249) + (-0.53125) = -1.6025. However, this update shows up only with -0.53124999. I'm quite confident the references in GlobalPoints are right. I am also sure that LM is calculated correctly but I do not manage to see how this update is working.

So how do I update with multiple references at the same time? Are these interfering with my update?

1
matlabmatrixvectorization
A is indeed a big sparse zeroes matrix at the beginning but since I am using reference 1 twice that could be only -0,53125. This code should add -0,531249999999875 + -0,531250000000126. – Artur Castiel
ah I see. No, that's not how that works unfortunately. If you have multiple access values to a single coordinate, it will only remember the last update... so the latest index of 1 is seen at position (2,4) in your matrix which is 0.53125. The update only uses this value. Actually, the behaviour of specifying multiple values to reference a single location is undefined. If you want to accumulate these you need to use something else. – rayryeng
I'll write you an answer. Hold on. – rayryeng

1 Answers

2
votes

What's going on in your code is that if you specify duplicate linear indices into the same location in a matrix, you expect all of the values that map to the same index to accumulate in sum. Unfortunately that's not how it works in MATLAB. It will only remember the update corresponding to the last time it sees the linear index.

Here is a very small example. Suppose I have this small 2 x 2 matrix:

>> A = [1 2; 3 4]

A =

     1     2
     3     4

Let's also say that I have the following 2 x 2 matrix of linear indices B as well as a matrix C that I want to accumulate with:

>> B = [1 1; 2 2]

B =

     1     1
     2     2

>> C = [10 20; 30 40]

C =

    10    20
    30    40

If I am interpreting your question correctly, the result that you are expecting should be:

A = [1 + 10 + 20, 2;
     3 + 30 + 40, 4]

This is because the first two rows should map to linear index 1, which is the top left corner and we should accumulate those together. The same is said for linear index 2, which is the bottom left corner. However, this is what you get if you try to do what you coded above:

>> A(B) = A(B) + C

A =

    21     2
    43     4

As you can see, the top left corner only gets updated with 20 + 1 as the value 20 coincides with the last time we saw index 1 for the first row. Similarly, the bottom left corner only gets updated with 40 + 3, as 40 is the last time we saw index 2 for the second row.

If you want to perform this accumulation, one thing I can suggest is to use accumarray to bin all of the values with the same linear index, add them up, then write them back to where they're supposed to go.

Something like this should work:

[vals,~,ind] = unique(GlobalPoint);
sums = accumarray(ind, LM(:));
A(vals) = A(vals) + sums;

The above code takes a bit of imagination, but it does work. The function unique determines all unique linear indices in the matrix GlobalPoint and returns this as a 1D vector stored in vals. The third parameter ind reparameterizes all values in GlobalPoint so that they're changed into indices from 1 up to as many unique values there are in vals. We then use accumarray where ind specifies which bin each corresponding value in LM should go to. The result will be an array where each element corresponds to the sum of all values that corresponded to a single linear index. The last step is to take this sum array and update the right positions in A.

Here's a quick example run of the steps I just talked about. The output of the first line of code given your data gives me this:

>> vals.'

ans =

     1     2    10    82    83    84    92   164   165   174   730   739   812   821   894   903

>> ind.'

ans =

  Columns 1 through 17

     1     5     5     9     9    11     2    13     6    15     3     4     7     8    10    12     1

  Columns 18 through 20

    14     5    16

vals contains all of the unique values in GlobalPoint and ind contains an ID which tells you which value in GlobalPoint each value maps to. Therefore, the index of 1 maps to the value 1 in GlobalPoint, the value 3 maps to the value of 10 in GlobalPoint and so on. Take note that this is sorted and there are 16 unique points in total.

Once I run through accumarray, I now get this for my sums:

>> sums.'

ans =

  Columns 1 through 4

                   -1.0625         0.531250000000125         0.531249999999875         0.531250000000126

  Columns 5 through 8

         -1.12500000000033         0.249999999999812         0.343749999999875         0.250000000000332

  Columns 9 through 12

        -0.499999999999721         0.249999999999909         0.531249999999875        -0.531249999999875

  Columns 13 through 16

         0.343749999999875        -0.343749999999875         0.249999999999909        -0.249999999999909

We can see that each bin for a linear index is summed properly. The last step is to take our corresponding unique linear indices with this sum array and update the right positions. I don't have your data to do this to show you your final result, but you can trust me that it works. To be absolutely sure, here's what the linear index and corresponding sums are side by side:

>> [vals sums]

ans =

                         1                   -1.0625
                         2         0.531250000000125
                        10         0.531249999999875
                        82         0.531250000000126
                        83         -1.12500000000033
                        84         0.249999999999812
                        92         0.343749999999875
                       164         0.250000000000332
                       165        -0.499999999999721
                       174         0.249999999999909
                       730         0.531249999999875
                       739        -0.531249999999875
                       812         0.343749999999875
                       821        -0.343749999999875
                       894         0.249999999999909
                       903        -0.249999999999909

Taking a look at linear index 1 already gives me confidence. We can see that both values that mapped to index 1 are summed properly. If you look at linear index 83 for example, we see that there are three values that coincide to this position: (2,1), (3,1) and (4,4). The values for these are -0.53125, -0.34375 and -0.25. Accumulating these gives roughly -1.125 and that's what we see in the result. If you repeat the calculation for the rest of the values, we can verify that the sums are calculated correctly.