MATLAB data indexing issue. What is going on here?

Question

Let I be the identity, D an orthonormal projection, and p a vector.

I realized that several of my lines of code combined to be (I-(I-D))(p) and I could just simplify it to D(p). In replacing it, I computed the new method along-side the old to double check I was computing the same thing (Earlier in my code I had a line that was D = I - D. The D you see here is that D.) I wasn't getting the same answer, and traced it to an error in indexing D.

Here you can see I'm using the debugger and checking portions of D and getting the wrong data returned.

The values in the data explorer on the right are what I'd expect them to be. Sometimes I get what I'd expect from D(:,:,k,1), and sometimes I don't, even when I make the queries right after each other.

The vectors those red arrows are pointing to should be the same. Nothing else changed or was computed between those lines, and k = 2 when the first line was run. I've closed MATLAB and restarted it and get the same issue every time. (D depends on random input, but I'm not altering the seed, so I get the same thing every first run after newly opening MATLAB. The way D is computed, I do expect D(:,:,1,1) to be the identity matrix.)

What in the world is going on? Any help is appreciated.

I have wondered if MATLAB is messing with me on purpose. Sometimes when I open it, a pop-up dialog box says I need to update my student license. I click the update button, but nothing ever happens and the dialog box never closes, so I click cancel.

Edit:

K>> whos D P
  Name      Size             Bytes  Class     Attributes

  D         4-D               4608  double              
  P         4x1x6              192  double              

K>> size(D)
ans =
     4     4     6     6

I've been playing around with A and B a bit, and I get the same thing. Sometimes it computes correctly and sometimes it doesn't.

K>> B=permute(P,[1,3,2])
B =
       0.4155      0.27554      0.52338       0.6991     -0.11346      0.20999
      0.53573     -0.83781      0.53182    -0.022364      0.60291     -0.62601
     -0.49246     -0.46111     -0.39168      0.45919      0.42377      0.47074
      0.54574     0.097595      0.53835     -0.54763      0.66637      0.58516


K>> A=D
A(:,:,1,1) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
A(:,:,2,1) =
      0.99071    -0.091198    0.0020814    -0.029755
    -0.091198      0.10503     0.020426       -0.292
    0.0020814     0.020426      0.99953    0.0066643
    -0.029755       -0.292    0.0066643      0.90473
A(:,:,3,1) =
      0.46769     0.019281     -0.49725     0.036486
     0.019281       0.9993     0.018011   -0.0013215
     -0.49725     0.018011      0.53551     0.034083
     0.036486   -0.0013215     0.034083       0.9975
A(:,:,4,1) =
      0.96774     0.063488     -0.10826      0.12438
     0.063488      0.87506      0.21304     -0.24477
     -0.10826      0.21304      0.63673      0.41737
      0.12438     -0.24477      0.41737      0.52047
A(:,:,5,1) =
       0.7542     0.031217      0.42575     0.056052
     0.031217      0.99604    -0.054071   -0.0071187
      0.42575    -0.054071      0.26255    -0.097088
     0.056052   -0.0071187    -0.097088      0.98722
A(:,:,6,1) =
       0.9818     -0.10286     0.085279    0.0034902
     -0.10286      0.41855      0.48208      0.01973
     0.085279      0.48208      0.60031    -0.016358
    0.0034902      0.01973    -0.016358      0.99933
A(:,:,1,2) =
      0.99071    -0.091198    0.0020814    -0.029755
    -0.091198      0.10503     0.020426       -0.292
    0.0020814     0.020426      0.99953    0.0066643
    -0.029755       -0.292    0.0066643      0.90473
A(:,:,2,2) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
A(:,:,3,2) =
      0.97125     -0.15889   -0.0080537    -0.051131
     -0.15889      0.12194    -0.044507     -0.28256
   -0.0080537    -0.044507      0.99774    -0.014323
    -0.051131     -0.28256    -0.014323      0.90907
A(:,:,4,2) =
      0.91488     -0.16388     -0.18495      0.12967
     -0.16388       0.6845     -0.35607      0.24964
     -0.18495     -0.35607      0.59815      0.28174
      0.12967      0.24964      0.28174      0.80247
A(:,:,5,2) =
      0.95461      0.16812      0.10326     0.066372
      0.16812      0.37733     -0.38244     -0.24582
      0.10326     -0.38244      0.76511     -0.15098
     0.066372     -0.24582     -0.15098      0.90295
A(:,:,6,2) =
      0.99628     0.012018     0.052874     0.027665
     0.012018      0.96117     -0.17085    -0.089393
     0.052874     -0.17085      0.24833     -0.39329
     0.027665    -0.089393     -0.39329      0.79422
A(:,:,1,3) =
      0.46769     0.019281     -0.49725     0.036486
     0.019281       0.9993     0.018011   -0.0013215
     -0.49725     0.018011      0.53551     0.034083
     0.036486   -0.0013215     0.034083       0.9975
A(:,:,2,3) =
      0.97125     -0.15889   -0.0080537    -0.051131
     -0.15889      0.12194    -0.044507     -0.28256
   -0.0080537    -0.044507      0.99774    -0.014323
    -0.051131     -0.28256    -0.014323      0.90907
A(:,:,3,3) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
A(:,:,4,3) =
      0.98622     0.043449    -0.066709     0.085142
     0.043449      0.86297      0.21038     -0.26852
    -0.066709      0.21038      0.67698      0.41227
     0.085142     -0.26852      0.41227      0.47382
A(:,:,5,3) =
      0.62859     0.041458      0.47558     0.074661
     0.041458      0.99537    -0.053085   -0.0083339
      0.47558    -0.053085      0.39105      -0.0956
     0.074661   -0.0083339      -0.0956      0.98499
A(:,:,6,3) =
      0.95505     -0.16608      0.12371    0.0067153
     -0.16608      0.38639      0.45705      0.02481
      0.12371      0.45705      0.65956     -0.01848
    0.0067153      0.02481     -0.01848        0.999
A(:,:,1,4) =
      0.96774     0.063488     -0.10826      0.12438
     0.063488      0.87506      0.21304     -0.24477
     -0.10826      0.21304      0.63673      0.41737
      0.12438     -0.24477      0.41737      0.52047
A(:,:,2,4) =
      0.91488     -0.16388     -0.18495      0.12967
     -0.16388       0.6845     -0.35607      0.24964
     -0.18495     -0.35607      0.59815      0.28174
      0.12967      0.24964      0.28174      0.80247
A(:,:,3,4) =
      0.98622     0.043449    -0.066709     0.085142
     0.043449      0.86297      0.21038     -0.26852
    -0.066709      0.21038      0.67698      0.41227
     0.085142     -0.26852      0.41227      0.47382
A(:,:,4,4) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
A(:,:,5,4) =
      0.73864      0.20112    -0.011394      0.39048
      0.20112      0.84524    0.0087678     -0.30047
    -0.011394    0.0087678       0.9995     0.017023
      0.39048     -0.30047     0.017023      0.41662
A(:,:,6,4) =
      0.87322     -0.15647    0.0029936      0.29363
     -0.15647      0.80689    0.0036946      0.36238
    0.0029936    0.0036946      0.99993   -0.0069332
      0.29363      0.36238   -0.0069332      0.31996
A(:,:,1,5) =
       0.7542     0.031217      0.42575     0.056052
     0.031217      0.99604    -0.054071   -0.0071187
      0.42575    -0.054071      0.26255    -0.097088
     0.056052   -0.0071187    -0.097088      0.98722
A(:,:,2,5) =
      0.95461      0.16812      0.10326     0.066372
      0.16812      0.37733     -0.38244     -0.24582
      0.10326     -0.38244      0.76511     -0.15098
     0.066372     -0.24582     -0.15098      0.90295
A(:,:,3,5) =
      0.62859     0.041458      0.47558     0.074661
     0.041458      0.99537    -0.053085   -0.0083339
      0.47558    -0.053085      0.39105      -0.0956
     0.074661   -0.0083339      -0.0956      0.98499
A(:,:,4,5) =
      0.73864      0.20112    -0.011394      0.39048
      0.20112      0.84524    0.0087678     -0.30047
    -0.011394    0.0087678       0.9995     0.017023
      0.39048     -0.30047     0.017023      0.41662
A(:,:,5,5) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1
A(:,:,6,5) =
      0.93556      0.24481   -0.0093576     0.016177
      0.24481     0.069855     0.035553    -0.061461
   -0.0093576     0.035553      0.99864    0.0023492
     0.016177    -0.061461    0.0023492      0.99594
A(:,:,1,6) =
       0.9818     -0.10286     0.085279    0.0034902
     -0.10286      0.41855      0.48208      0.01973
     0.085279      0.48208      0.60031    -0.016358
    0.0034902      0.01973    -0.016358      0.99933
A(:,:,2,6) =
      0.99628     0.012018     0.052874     0.027665
     0.012018      0.96117     -0.17085    -0.089393
     0.052874     -0.17085      0.24833     -0.39329
     0.027665    -0.089393     -0.39329      0.79422
A(:,:,3,6) =
      0.95505     -0.16608      0.12371    0.0067153
     -0.16608      0.38639      0.45705      0.02481
      0.12371      0.45705      0.65956     -0.01848
    0.0067153      0.02481     -0.01848        0.999
A(:,:,4,6) =
      0.87322     -0.15647    0.0029936      0.29363
     -0.15647      0.80689    0.0036946      0.36238
    0.0029936    0.0036946      0.99993   -0.0069332
      0.29363      0.36238   -0.0069332      0.31996
A(:,:,5,6) =
      0.93556      0.24481   -0.0093576     0.016177
      0.24481     0.069855     0.035553    -0.061461
   -0.0093576     0.035553      0.99864    0.0023492
     0.016177    -0.061461    0.0023492      0.99594
A(:,:,6,6) =
     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

---

Edit 2: Added relevant code. I've been pausing the code and getting the errors inside the for loops at the end. (I believe it's also giving errors in S, but I've been focusing on D trying to figure it out.)

mtimesx is from here.

n = 4;
M = 6;

P = Normalize(2*rand(n,1,M)-1);

%differences between p_i and p_j
%sum of p_i and p_j
d = Normalize(repmat(permute(P,[1,3,2]),[1,1,M]) - repmat(P,[1,M,1]));
s = Normalize(repmat(permute(P,[1,3,2]),[1,1,M]) + repmat(P,[1,M,1]));

d(isnan(d)) = 0;

%orthogonal projection onto d(:,i,j), i.e. outer product of differences
%orthogonal projection onto s(:,i,j), i.e. outer product of sums
D = mtimesx(permute(d,[1,4,2,3]), permute(d,[4,1,2,3]));
S = mtimesx(permute(s,[1,4,2,3]), permute(s,[4,1,2,3]));
D2 = D;
S2 = S;

%projection onto the complement of d(:,i,j)
%projection onto the complement of s(:,i,j)
D = repmat(eye(n),[1,1,M,M]) - D;
S = repmat(eye(n),[1,1,M,M]) - S;

%total distance to the nearest subspace
PDist = zeros([1,M]);
PDist2 = PDist;
for j = 1:M
    for k = 1:M-1
        for l = k:M
            if j~=k && j~=l
                PDist(j) = PDist(j) + min(norm(P(:,1,j) - mtimes(D(:,:,k,l),P(:,1,j))), norm(P(:,1,j) - mtimes(S(:,:,k,l),P(:,1,j))));
                PDist2(j) = PDist2(j) + min(norm(D2(:,:,k,1)*P(:,1,j)),norm(S2(:,:,k,1)*P(:,1,j)));
            end
        end
    end
end
PDist-PDist2

---

Normalize.m

%Normalize
%Accepts an array (of column vectors) and normalizes the columns

function B = Normalize(A)
    B = A./repmat(sqrt(sum(A.*A)),size(A,1),1);
end

Answer 1

The problem is that you indexed the matrices using the constant 1 instead of the variable l (lowercase L), both in the first example and in the code for computing PDist2.

In general it is good to avoid using variable names that look similar to each other and/or similar to numbers.

This can be avoided by using an editor that highlights uses different colors for variables and constants (I don't know if this is possible in MATLAB). In fact, this is how I found the error in your code. As you can see, when indexing D2 for the computation of PDist2 the number 1 is colored red.