how to find maximum value in a matrix?

Now, there is reason why I don't want to use a max() or similar built-in function to find the maximum value of a matrix. My question is following:

Imagine, we have a matrix, let's call it 'gauss' and looks like this:

 gauss =

         0    0.2000         0
    0.5000    1.0000    0.6000
         0    0.3000         0

What I want is to find the maximum, i.e. 1 or the value at gauss(2,2). I need to automize a step search from gauss(1,1) so that I end up finding gauss(2,2) in this entire matrix.

Would be glad if you could give some idea regarding the coding of the steps.

Thanks!

[Anton Semechko]

It is in general not possible to find (global) maximum value of the matrix when performing local steps, from one matrix element to its adjacent neighbours. If you do so, then for any given starting position, the algorithm will converge onto a local maximum, which may or may not be a global maximum.

[kat001]

That is correct! However, if you imagine a 2D map of a gaussian profile, where the peak value is in the center of map, then it should be possible to take x number of steps and then y number of steps from a random point on the map until it reaches the peak. Right?

Answer 1

[Anton Semechko]

Demo of local search based on 4-connected neighbourhood:

Demo of local search based on 8-connected neighbourhood:

Code used to generate these figures:

function Gauss_map_demo

% Sampling grid
N=30; % grid size
dx=2*rand(1)-1;
dy=2*rand(1)-1;
x=linspace(-2,2,N)+dx;
y=linspace(-2,2,N)+dy;
ds=x(2)-x(1);
[X,Y]=meshgrid(x,y);

% Random precision matrix (inverse of covariance)
t=rand(1)*pi;
R=[cos(t) -sin(t);sin(t) cos(t)];
D=sort((rand(1,2)+1)/2,'descend');
S=R'*diag(1./D)*R;

% Evaluate (zero-centred) Gaussian on specified grid
P=[X(:) Y(:)]';
G=exp(-0.5*sum(P.*(S*P),1));
G=reshape(G,size(X));

% Visualize 
figure('color','w')
imagesc([1 N],[1 N],G)
axis equal
set(gca,'XLim',[0 N+1],'YLim',[0 N+1])
hold on
colorbar;

% Identify corner with smallest value as a starting point; this is done for
% demonstration purposes. In practice,  you should select corner with 
% largest value, as it will reduce number of steps required to reach 
% absolute maximum.
ij=[1 1;N 1;N N;1 N]; % corner indices
id=sub2ind([N N],ij(:,1),ij(:,2));
[~,id_min]=min(G(id));
seed=ij(id_min,:);

% Search using on 4-connected neighbourhood
conn=false; % set conn=true to use 8 connected neighbourhood

[G_max,ij_max,trj]=GetMax(G,seed,conn);
fprintf('Maximum value found %.3f at G(%u,%u) in %u steps\n',G_max,ij_max(1),ij_max(2),size(trj,1))

plot(trj(:,2),trj(:,1),'-r','LineWidth',2);                      % ascent trajectory
plot(trj(1,2),trj(1,1),'ok','LineWidth',3,'MarkerSize',10);      % start
plot(trj(end,2),trj(end,1),'xk','LineWidth',3,'MarkerSize',10);  % end

function [G_max,ij_max,trj]=GetMax(G,seed,conn)
% Get maximum value of an image using local search
%
%   - G     : N-by-M image
%   - seed  : starting point
%   - conn  : set conn=true to use 8-connected neighbourhood {default}, 
%             and set conn=false to use 4-connected neighbourhood 

if nargin<2 || isempty(seed), seed=[1 1]; end
if nargin<3 || isempty(conn), conn=true; end

siz=size(G);

% Initial value
ij_max=seed;
G_max=G(ij_max(1),ij_max(2));
trj=ij_max;

% Neighbour offsets
ngb_x=repmat([-1 0 1],[3 1]);
ngb_y=ngb_x';
d_ij=[ngb_y(:) ngb_x(:)];
if conn
   d_ij(5,:)=[];
else
    d_ij([1 3 5 7 9],:)=[];
end

% Search
while true

      % Indices of neighbours around ij_max  
      i_ngb=ij_max(1)+d_ij(:,1);
      i_ngb=max(min(i_ngb,siz(1)),1);

      j_ngb=ij_max(2)+d_ij(:,2);
      j_ngb=max(min(j_ngb,siz(2)),1);

      id_ngb=sub2ind(siz,i_ngb,j_ngb);

      % Select neighbour with largest value; I am using 'max' function here, but you can loop through the neighbouring values if you want 
      [g_max,id_max]=max(G(id_ngb));
      if g_max>G_max
          G_max=g_max;
          ij_max_new=[i_ngb(id_max) j_ngb(id_max)];
          if nargout>2
              trj=cat(1,trj,ij_max_new);
          end            
          ij_max=ij_max_new;
      else
          break
      end

end

[kat001]

Thank you very much for the code.

However, I would like to understand the following parts of the code from the first function:

 % Random precision matrix (inverse of covariance)
 t=rand(1)*pi;
 R=[cos(t) -sin(t);sin(t) cos(t)];
 D=sort((rand(1,2)+1)/2,'descend');
 S=R'*diag(1./D)*R;

 % Evaluate (zero-centred) Gaussian on specified grid
 P=[X(:) Y(:)]';
 G=exp(-0.5*sum(P.*(S*P),1));
 G=reshape(G,size(X));

Why generating the R and D and how did you calculate G? I couldn't really understand that. Do you have any material regarding this which you can refer me to?

Many thanks again!

[Anton Semechko]

Gaussian is parameterized by mean and covariance. Precision is inverse of covariance. In the demo, I set mean to the origin, and just shifted the sampling grid with respect to the origin.

https://en.wikipedia.org/wiki/Multivariate_normal_distribution

[Anton Semechko]

Additional notes: (1) Covarince matrix, C, can be decomposed as C=R'*D*R, where R is a rotation matrix and D is a diagonal matrix. (2) Inverse of C = R'*inv(D)*R ==> precision

[kat001]

Thank you. As I am quite new to Matlab, I am learning a lot from this example.

Now, I am also going through the second function or tracking function I should say. Would you please provide a very short explanation of the following part of the code, please:

% Neighbour offsets
 ngb_x=repmat([-1 0 1],[3 1]);
 ngb_y=ngb_x';
 d_ij=[ngb_y(:) ngb_x(:)];
 if conn
    d_ij(5,:)=[];
 else
    d_ij([1 3 5 7 9],:)=[];
 end

In the while-loop, you have presented the neighbor indices of the starting index. Would be grateful if you could also explain inline of the following code:

% Indices of neighbours around ij_max  
      i_ngb=ij_max(1)+d_ij(:,1);
      i_ngb=max(min(i_ngb,siz(1)),1);
      j_ngb=ij_max(2)+d_ij(:,2);
      j_ngb=max(min(j_ngb,siz(2)),1);

Thanks!

[Anton Semechko]

Suppose that i and j are indices specifying location of an element in array G of size M-by-N. What are the indices of all elements in the immediate neighbourhood around G(i,j)?

[kat001]

So,

ngb_x =

      -1     0     1
      -1     0     1
      -1     0     1

ngb_y =

      -1    -1    -1
       0     0     0
       1     1     1

And then, for example, you are "erasing" linear index number 1, 3, 5, 7 and 9 from ngb_x and ngb_y and to me, it looks like I am getting two vectors, i.e.

d_ij = 

-1    0
 0   -1
 0    1
 1    0

Why doing that? Is there a way so that you could comment right next to the codes which I highlighted in my previous post? Thank you again for your patience with me!

[Anton Semechko]

Lets say you are at a "point" (i,j) in array G. In a 4-connected neighbourhood, neighbouring points around (i,j) are

(i-1,j+0)
(i+0,j-1)
(i+0,j+1)
(i+1,j+0)

So, the 1st index of all neighbours is i_ngb=i+d_ij(:,1), and corresponding 2nd index is j_ngb=j+d_ij(:,2). Makes sense?

Finally, recall that G is a N-by-M array. So if i=1 (resp. i=N) then regardless of j, point (i,j) doesn't have any neighbours above (resp. below) it. The statement

i_ngb=max(min(i_ngb,N),1);

will enforce this requirement.

Likewise, if j=1 (resp. j=M) then regardless of i, point (i,j) doesn't have any neighbours to the left (resp. right) of it. The statement

j_ngb=max(min(j_ngb,M),1);

will enforce this requirement.

Answer 2

[Walter Roberson]

gr = size(gauss,1);
gc = size(gauss,2);
rowpos = 1; colpos = 1;

improved = true;
while improved
  bestval = gauss(rowpos, colpos);
  improved = false;
  for rowtry = max(1, rowpos-1) : min(rowpos+1, gr)
    for coltry = max(1, colpos-1) : min(colpos+1, gc)
      newval = gauss(rowtry, coltry);
      if newval > bestval
        bestval = newval;
        rowpos = rowtry;
        colpos = coltry;
        improved = true;
        break;
      end
    end
    if improved; break; end
  end
end

If it did not get caught in a local minima, then afterwards rowpos and colpos will be the location of the peak, and bestval will be the value there.