Iterate over diagonal elements of a Matrix in MatLab

Question

I need to get the indexes of all diagonals in a matrix. The matrix can be non square.

The diag function gives the values, but I need the coords. So for instance with this: [1 2 3; 4 5 6; 7 8 9], I want [1 1; 2 2; 3;3] PLUS [2 6] and 3 because they are the upper diagonals of the matrix, and same for below i.e. [4 8] and 7. So the full list of indexes is: [1 1; 2 2; 3 3], [1 2; 2 3], [1 3], [2 1], [3 2], [3 1].

And I need this in the other diagonal direction too...

And it needs to work for matrices of any shape i.e. not just square ones, including vectors (i.e. 1xn and nx1) - but I guess I can test for these last 2 cases and treat separately.

EDIT: An example of a non square matrix is this 3x2 one:

1   2
3   4
5   6

In this case, the indexes I need are:

[1 1; 2 2]
[2 1; 3 2]
[1 2]
[3 1]

And then the same for diagonals running the other way...

I am aware of the diag function and can write a bit of code to iterate over the matrices, but I can't then see how to easily turn those values into indexes.

As for why I need this: It's an online course I'm doing in MatLab. An introduction to programming, using MatLab. I'm a reasonably experienced programmer but I'm not familiar with MatLab and I want to learn it more.

The question asks us to find the biggest product of any n numbers in any direction, i.e. rows, cols, and diagonals (both ways). And the thing is we don't return the product, we return the indexes of the elements making up this product.

I hope this provides more information.

EDIT 2: Sorry for not getting back to you guys, I've been very busy doing other things. When I first posted here, I assumed I was missing something and this would be something that took a few lines. But it's trickier than that, so the code you guys are posting will take me a little more time to go through and understand etc. But I hope to do this today and accept one as an answer. Thanks for all the replies.

Answer 1

Don't forget that Matlab understands linear indexing into multi-dimensional matrices. For OP's example of a 3*3 (let's call it a) matrix the expression a(3) identifies the same element as a(3,1) (Matlab uses column-major layout). So, for a k*k matrix, the expression

[1, k+2, 2k+3, 3k+4, ...]

gives the indices of the elements on the main diagonal. There is a range of ways to generate that vector. Then, a judicious application of ind2sub will give the 2-d indices into the same matrix.

It shouldn't be too hard to generalise this approach to the other diagonals and to non-square matrices.

Answer 2

Here is a solution for the "northwest to southeast" diagonals.

A = [1,2,3,4,5;6,7,8,9,10;11,12,13,14,15];  
rows = size (A,1);   cols = size (A,2);   padding = rows - 1; 
padded = [ones(rows, padding), A, ones(rows, padding)];

% Create a Helper array containing column/row pairs and a product
Helper = zeros (cols + padding, rows * 2 + 1);
for i = 1 : cols + padding;  % i moves along columns of the padded array
  i0 = i - padding; % "actual" index w.r.t A, where leftmost index is 1-padding
  c  = 1;       % a counter used for filling current "Helper" row appropriately
  for j = 1 : rows;       % j moves along rows
    Helper(i, c)   = i0;  % store current "actual" column
    Helper(i, c+1) = j;   % store current row
    c  = c  + 2;          % move to next 2 slots in "Helper"
    i0 = i0 + 1;          % move to next "actual" column for processing
  end 

  % Process product at appropriate indices
  Helper(i, end) = prod (padded (sub2ind (       ...
                     size (padded),              ...
                     Helper(i,2:2:end-1),        ...
                     Helper(i,1:2:end-1)+padding ...
                   )));
end

% Isolate the correct indices for which the product is largest
[~, biggest] = max (Helper(:, end)); % column index of largest product in Helper
Pairs = reshape (Helper(biggest, 1:end-1), 2, []); % collect indices
Pairs = fliplr (Pairs.');                          % express as rows/cols

% Remove pairs where column occurs in the "padded" zone
for i = 1 : size (Pairs, 1)
  if Pairs(i, 1) < 1 || Pairs(i, 1) > cols; Pairs(i, :) = []; end
end

display(Pairs)

If you understand the logic behind this, the "NorthEast to SouthWest" diagonals should be a piece of cake to obtain too. I'll leave that for you to figure out :)

Answer 3

The easy way to do this is to construct an array the same shape as your input where the value of each element is its index. Then you can get the values (indices) of each diagonal.

function diag_cell = alldiags(A)    
   [m,n] = size(A);

   M = reshape(1:m*n, m, n);   % construct array of indices
   % there are (m-1) diagonals below the main (k parameter negative)
   %   and (n-1) diagonals above, plus 0 for the main diagonal
   for d = 1-m:n-1
      C1{d+m} = diag(M, d);   % store in C1{1}..C1{m+n-1}
   end

   % to get the antidiagonals, rotate the index matrix `M` and repeat
   %   (remember to swap `m` and `n`):
   M = rot90(M);
   for d = 1-n:m-1
      C2{d+n} = diag(M, d);
   end

   diag_cell = {C1{:}, C2{:}};   % concatenate the cell arrays and return
end

Here's a test run with the 3x3 sample matrix:

A =

   1   2   3
   4   5   6
   7   8   9

>> alldiags(A)
ans =
{
  [1,1] =  3
  [1,2] =

     2
     6

  [1,3] =

     1
     5
     9

  [1,4] =

     4
     8

  [1,5] =  7
  [1,6] =  1
  [1,7] =

     4
     2

  [1,8] =

     7
     5
     3

  [1,9] =

     8
     6

  [1,10] =  9
}

---

If you really have to, you can convert all of the indices to subscripts as you go, but if you're accessing the original matrix it's probably easier to do it with indices. Here's an example converting the main diagonal to subscripts:

>> C{3}
ans =

   1
   5
   9

>> [r,c] = ind2sub(size(A), C{3});
>> subs = [r c]
subs =

   1   1
   2   2
   3   3

---

Here's an example to illustrate why this works. Let's take a random 6x3 matrix A:

A =

   0.634825   0.560609   0.926716
   0.049504   0.226152   0.748278
   0.746754   0.493896   0.773868
   0.993245   0.457522   0.412092
   0.430003   0.695042   0.641917
   0.935298   0.960166   0.022872

The first thing the function does is get the size of the input matrix A and create a new matrix of the same size such that the value of each element is equal to its linear index. Linear indices start at 1 at subscript [1,1], increase down the column to m at [m,1], then go to m+1 for subscript [1,2], ending at index m*n at subscript [m,n].

>> [m,n] = size(A)
m =  6
n =  3
>> M = reshape(1:m*n, m, n)
M =

    1    7   13
    2    8   14
    3    9   15
    4   10   16
    5   11   17
    6   12   18

The output diag_cell for this matrix looks like this (heavily reformatted so it's easier to see all at once):

diag_cell =
{
  [1,1]  =   6
  [1,2]  =   5    12
  [1,3]  =   4    11    18
  [1,4]  =   3    10    17
  [1,5]  =   2     9    16
  [1,6]  =   1     8    15
  [1,7]  =   7    14
  [1,8]  =  13
  [1,9]  =   1
  [1,10] =   7     2
  [1,11] =  13     8     3
  [1,12] =  14     9     4
  [1,13] =  15    10     5
  [1,14] =  16    11     6
  [1,15] =  17    12
  [1,16] =  18
}

The first half of the entries are the normal diagonals (top-left to bottom-right), and the second half are the antidiagonals (top-right to bottom-left).

As an example, take the diagonal diag_cell{4} = [3; 10; 17]. You can look at the matrix M to see which elements these indices refer to.

A =

   0.634825   0.560609   0.926716
   0.049504   0.226152   0.748278
  *0.746754*  0.493896   0.773868
   0.993245  *0.457522*  0.412092
   0.430003   0.695042  *0.641917*
   0.935298   0.960166   0.022872

This diagonal starts at [3,1] and ends at [5,3]. If we recover the values from A, we get:

>> A(diag_cell{4})
ans =

   0.74675
   0.45752
   0.64192

Answer 4

This solution goes all the way to your final task - finding the indices of the maximum product of any vector in the matrix:

function out = maxVecProd(A)
% validate the input is a vector\matrix:
assert(ismatrix(A),'Input must be a vector or matrix')
[cmax,cind] = max(prod(A,1)); % find where the maximum prod in the columns
[rmax,rind] = max(prod(A,2)); % find where the maximum prod in the rows
if all(size(A))>1 % if its a matrix:
    all_d = -(size(A,1)-1):(size(A,1)-1); % get all the valid diagonals
    pd = zeros(2*numel(all_d),2); % the solution of all products
    c = 1;
    for d = all_d
        pd(c,1) = prod(diag(A,d)); % get the diagonal product
        pd(c,2) = prod(diag(fliplr(A),d)); % get the anti-diagonal product
        c = c+1;
    end
    [dmax,dind] = max(pd(:)); % find where the maximum prod in the diagonals
    % get the orientation of the maximum prod (column, row, diagonal)
    [~,dirmax] = max([rmax,cmax,dmax]);
else
    [~,dirmax] = max([rmax,cmax]);
end
switch dirmax
    case 1
        % create the indices for the maximum row:
        out = [repelem(rind,size(A,2)).' (1:size(A,2)).'];
    case 2
        % create the indices for the maximum column:
        out = [(1:size(A,1)).' repelem(cind,size(A,1)).'];
    case 3
        [r,dir] = ind2sub(size(pd),dind); % convert to [row,column]
        d_ind = all_d(r); % find which value of d gave the maximum product
        [R,C] = ndgrid(1:size(A,1),1:size(A,2)); % create a grid of all indices
        if dir==1
            % find the indices for a diagonal:
            out = [diag(R,d_ind),diag(C,d_ind)];
        else
            % find the indices for an anti-diagonal:
            out = [diag(fliplr(R),d_ind),diag(fliplr(C),d_ind)];
        end
end
end

and you get for maxVecProd(A):

out =
     1     2
     2     2
     3     2

which means the largest prod in is of the elements in (1,2), (2,2) and (3,2).

Answer 5

Ok, so I coded this up myself last night. It's very long and I'm sure it can be done in better ways. I don't know MatLab well, so my solution is pretty naive. But it works.

function indicies = maxproduct(A, n)


function diagIndicies = getAllReverseDiagonalIndicies()
    % This function returns lists of indicies in the reverse diagonal direction (top-right to bottom-left)
    if rows == 1 
        numCells = 0;
        for j=1:length(A)
            temp = [1, j];
            numCells = numCells + 1;
            diagIndicies{numCells} = temp;
        end
        return
    end

    % This loop adds all diagonals down and to the right, to the end of the matrix.
    %fprintf('Diagonals down to main one are:\n');
    numCells = 0;
    for x=1:rows
        rowNum = x;
        colNum = 1;
        temp = [];
        while rowNum >= 1 && rowNum <= rows && colNum <= cols
            temp = [temp; [rowNum colNum]];
            rowNum = rowNum - 1;
            colNum = colNum + 1;
        end
        numCells = numCells + 1;
        temp = flipud(temp); % Need row major order for assignment
        %disp(temp);
        diagIndicies{numCells} = temp;
    end

    % Now go along bottom row
    %fprintf('Diagonals along bottom are:\n');
    for y=2:cols
        rowNum = rows;
        colNum = y;
        temp = [];
        while rowNum >= 1 && colNum <= cols
            temp = [temp; [rowNum colNum]];
            rowNum = rowNum - 1;
            colNum = colNum + 1;
        end
        numCells = numCells + 1;
        temp = flipud(temp); % Need row major order for assignment
        %disp(temp);
        diagIndicies{numCells} = temp;
    end                

end

function diagIndicies = getAllDiagonalIndicies()
    % This function returns lists of indicies in the main diagonal direction (top-left to bottom-right)
    if rows == 1 
        numCells = 0;
        for j=1:length(A)
            temp = [1, j];
            numCells = numCells + 1;
            diagIndicies{numCells} = temp;
        end
        return
    end

     % This loop adds all diagonals down and to the left, to the lhs of the matrix.
    %fprintf('Diagonals down to main one are:\n');
    numCells = 0;
    for x=1:rows
        rowNum = x;
        colNum = cols;
        temp = [];
        while rowNum >= 1 && rowNum <= rows && colNum >= 1
            temp = [temp; [rowNum colNum]];
            rowNum = rowNum - 1;
            colNum = colNum - 1;
        end
        numCells = numCells + 1;
        temp = flipud(temp); % Need row major order for assignment
        %disp(temp);
        diagIndicies{numCells} = temp;
    end

    % Now go along bottom row...
    %fprintf('Diagonals along bottom are:\n');
    for y=cols-1:-1:1
        rowNum = rows;
        colNum = y;
        temp = [];
        while rowNum >= 1 && colNum >= 1
            temp = [temp; [rowNum colNum]];
            rowNum = rowNum - 1;
            colNum = colNum - 1;
        end
        numCells = numCells + 1;
        temp = flipud(temp); % Need row major order for assignment
        %disp(temp);
        diagIndicies{numCells} = temp;
    end

 end


 % Main function starts here.

[rows, cols] = size(A);

theMaxProduct = -10000000;
indicies = [];

if isscalar(A)
    if A == 1 && n == 1
        indicies = [1,1];
        return;
    else
        return;
    end
end



% Find max product in each row of A
for i=1:rows
    for j=1:cols-n+1
        theProduct = 1;
        for k=j:j+n-1
            theProduct = theProduct * A(i, k);
        end
        if theProduct > theMaxProduct
            theMaxProduct = theProduct;
            indicies = [];
            for k=j:j+n-1
                indicies = [indicies; [i, k]];
            end

        end
    end
end
fprintf('theMaxProduct after examining rows = %15.15f\n', theMaxProduct);

% Find max product in each column of A
for i=1:cols
    for j=1:rows-n+1
        theProduct = 1;
        for k=j:j+n-1
            theProduct = theProduct * A(k, i);
        end
        if theProduct > theMaxProduct
            theMaxProduct = theProduct;
            indicies = [];
            for k=j:j+n-1
                indicies = [indicies; [k, i]];
            end

        end
    end
end
fprintf('theMaxProduct after examining cols = %15.15f\n', theMaxProduct);



% Find max product along reverse diagonals of A
diagIndicies = getAllReverseDiagonalIndicies();
%disp(diagIndicies);
for i=1: length(diagIndicies)
    [numIndicies, ~] = size(diagIndicies{i});
    for j=1:numIndicies-n+1
        theProduct = 1;
        for k=j:j+n-1
            theProduct = theProduct * A(diagIndicies{i}(k, 1), diagIndicies{i}(k, 2));
        end
        if theProduct > theMaxProduct
            theMaxProduct = theProduct;
            indicies = [];
            for k=j:j+n-1
                indicies = [indicies; [diagIndicies{i}(k, 1), diagIndicies{i}(k, 2)]];
            end

        end
    end
end
fprintf('theMaxProduct after examining reverse diag = %15.15f\n', theMaxProduct);


% Find max product along diagonals of A
diagIndicies = getAllDiagonalIndicies();
%disp(diagIndicies);
for i=1: length(diagIndicies)
    [numIndicies, ~] = size(diagIndicies{i});
    for j=1:numIndicies-n+1
        theProduct = 1;
        for k=j:j+n-1
            theProduct = theProduct * A(diagIndicies{i}(k, 1), diagIndicies{i}(k, 2));
        end
        if theProduct > theMaxProduct
            theMaxProduct = theProduct;
            indicies = [];
            for k=j:j+n-1
                indicies = [indicies; [diagIndicies{i}(k, 1), diagIndicies{i}(k, 2)]];
            end

        end
    end
end
fprintf('theMaxProduct after examining main diag = %15.15f\n', theMaxProduct);

end