Sorting to have continuous non negative values next to diagonal of square matrix

Hello

I have a square matrices of different sizes. I want to reorder the rows of the matrix in such a way that I get the maximum number of continuous nonnegative values at the positions occurring immediately to the next of the diagonal matrix. The diagonal of the matrix will always be -1.

For example in square matrix 5*5

    % 1  2  3  4  5

    [-1  1  2 -1  3    % 1
      0 -1 -1 -1  4    % 2 
      5 -1 -1  6  7    % 3 
     -1 -1  8 -1  9    % 4
     10 11 12 13 -1]   % 5

In this case I want maximum number of continuous non negative values on the positions (1,2),(2,3),(3,4),(4,5). Currently we can see that (1,2) has a positive value and since (2,3) has a negative value. But if we rearrange the rows then we can have more continuous values.

One more point to be noted is that the values in the rows and column are interlinked. What I mean by this is that suppose you swap rows 2 and 5 , then you need to swap columns 2 and 5 as well. This will ensure that diagonal matrix will remain -1.

If we do this we get

    % 1  2  3  4  5

    [-1  3  2 -1  1    % 1
     10 -1 12 13 11    % 2 
      5  7 -1  6  7    % 3 
     -1  9  8 -1 -1    % 4
      0  4 -1 -1 -1]   % 5

We swap rows 2 and 5 first and then followed by swapping columns 2 and 5, we have 3 continuous non negative values at (1,2),(2,3) and (3,4).

    % 1  2  3  4  5

    [-1  9  8 -1 -1    % 1
     13 -1 12 10 11    % 2 
      6  7 -1  5 -1    % 3 
     -1  3  2 -1  1    % 4
     -1  4 -1  0 -1]   % 5

And if we swap rows 1 and 4 and also columns 1 and 4 , we get continuous values at all 4 positions and by sorting in such a way we get maximum of continuous non negative values.

Note that I just want to get the maximum of continuous non negative values, may be in another example for a 5*5 matrix, it could just have 1 or 2 continuous values, but that should be the maximum.

I hope I have explained my problem and any help would be appreciated. This is just an example and I have square matrices of the size 30*30 on which I need to perform the operation.

[Manoj]

I would like to know if the problem could be solved at all in Matlab or maybe any toolboxes would be of any help since for a 30*30 square matrix , there are 30! possible ways the rows can be sorted to find out this maximum non negative values next to the diagonal

[Manoj]

Hello

Is my question unclear because I haven't figured out a way to solve this problem. Any input would be helpful at this time

[Guillaume]

I don't think your question is unclear at all. In fact, it's a lot clearer than many questions we get.

I think you didn't get any reply because it's not a trivial problem to solve. (or maybe it's just bad luck and your question was missed). It's certainly not a problem domain I'm familiar with.

[Matt J]

What I mean by this is that suppose you swap rows 2 and 5 , then you need to swap columns 2 and 5 as well.

Suppose I break this rule and use some arbitrary permutation of rows and columns, but still manage to keep -1's on the diagonal. Would you still reject that as a solution?

[Manoj]

Yes the solution would not be what I want. The rule has to be maintained because I get those particular values because of its positions and hence the swapping of both the rows and cols is mandatory.

[Manoj]

And I also need to know the initial and final positions of the rows. As in the original rows have ended up where after we are done with all the swapping.

Answer 1

[Matt J]

The code below does the job, I believe. Or at least, I find that it works with all of your sample inputs, and with a worst case time of 8.5 seconds on my machine (for Input_six).

Note that the sorting problem really turns out to be a Longest Path Problem which is known to be NP-hard. Therefore, sparsity of the non-negative values is key to keeping the computation time low.

    function [Msorted,rowperm] = mySort(M)
    %Sorts input M to achieve a maximum length string of non-negative values on
    %diag(Msorted,1). Optionally, also, returns "rowperm", an index vector
    %satisfying
    %
    %  Msorted=M(rowperm,rowperm)

            n=size(M,1);

            p=LongestPath(M>=0);

            rowperm=[p,setdiff(1:n,p)];

            Msorted=M(rowperm,rowperm);

    function out=LongestPath(G,options)
    %Finds a longest path through adjacency matrix, G

         if nargin<2
           options=1:size(G,1);  
         end

          m=length(options);

          c=cell(m,1);

          for i=1:m

             j=options(i);
             tmp=G;
             tmp(:,j)=0;
             suboptions=find(tmp(j,:));
             if isempty(suboptions)%nowhere to go
              c{j}=j;   
             else    
              c{j}=[j,LongestPath(tmp,suboptions)]; 
             end
          end

          L=cellfun('length',c);

          [~,idx]=max(L);

          out=c{idx};

[Manoj]

Hello Matt

I ran the code and it seems to be working perfectly for all my input matrices so far. Thanks a ton for all your time and effort to solve this problem, really appreciate it.

Best Regards

Manoj

Answer 2

[Guillaume]

I'm really not familiar with this problem domain (optimisation). The only way I can think of doing what you want is to try all the permutations and see which one is the best. This is only reasonable for small matrices otherwise computation time will be prohibitive.

You would use perms to generate all possible combinations of row/column indices and diag to extract the diagonal for testing:

function bestm = findbest()
  m =[-1  1  2 -1  3    % 1
       0 -1 -1 -1  4    % 2 
       5 -1 -1  6  7    % 3 
      -1 -1  8 -1  9    % 4
     10 11 12 13 -1];   % 5 
  diagcount = countcontinuous(diag(m, 1));
  bestm = m;
  allcombs = perms(1:size(m, 1));
  for comb = allcombs'
     mswappedrow=m(comb, :);
     mswapped = mswappedrow(:, comb);
     count = countcontinuous(diag(mswapped, 1));
     if count > diagcount
        diagcount = count;
        bestm = mswapped;
     end
     if count == size(m, 1)-1
        break; %it's the maximum possible anyway
  end
end

function count = countcontinuous(diagonal)
   positive = diagonal' > 0;
   transitions = diff([0 positive 0]); % a 1 is a transition from 0 to 1, a -1 from 1 to 0
   %because we've put a 0 before and after we're sure to have the same number of transitions
   %from 0 to 1 and from one to 0
   runlength = find(transitions == -1) - find(transitions == 1);
   count = max(runlength);
end

[Manoj]

Hello , thank you for your reply and providing me with a solution . Yes , the problem domain would be quite big. That is why I was unsure of going ahead with permutation of all the rows especially for 30*30 matrix which cannot be solved with respect to computational time and also space requirements.

I will try to see if I can work with what you have posted and try to adapt it to my needs. Thanks once again

[Guillaume]

Yes, perms is certainly not feasible for 30 rows. As an alternative you could use randperm to try different combinations at random and stop after a given number of iterations.

Note that my algorithm above will stop immediately if it manages to have only positive number on the 2nd diagonal. Therefore, if your 30*30 matrix is mostly positive numbers, with randperm it should hopefully stop quickly.

If your matrix is mostly negative numbers, you could change the stop condition

if count == size(m, 1)-1

to

if count == maxpossiblecount

with

maxpossiblecount = min([size(m, 1)-1 sum(m(:) >= 0)]);

[Manoj]

Yes, my matrices will generally not have all non negative numbers on the 2nd diagonal. There will be mostly few negative numbers on the second diagonal for matrices having size 20 and above.

[Guillaume]

My solution does not really care how the matrix starts, since it tries random (with randperm) or ordered (with perms) permutations of all the rows and columns.

What I was saying is that if your matrix is mostly positive numbers anywhere, then there would be many permutations that result in all numbers positive on the 2nd diagonal and the algorithm can stop quite quickly.

For example, modifying my code to use randperm

function [bestm, diagcount, step] = findbest(m, maxstep)
  diagcount = countcontinuous(diag(m, 1));
  bestm = m;
  %allcombs = perms(1:size(m, 1));
  %for comb = allcombs'
  step = 0;
  while step < maxstep
     step = step+1;
     comb = randperm(size(m, 1)); 
     mswappedrow=m(comb, :);
     mswapped = mswappedrow(:, comb);
     count = countcontinuous(diag(mswapped, 1));
     if count > diagcount
        diagcount = count;
        bestm = mswapped;
     end
     if count == size(m, 1)-1
         break
     end
  end
  if step == maxstep
      warning('%d steps reached. calculation stopped', maxstep);
  end
end

function count = countcontinuous(diagonal)
   positive = diagonal' > 0;
   transitions = diff([0 positive 0]); % a 1 is a transition from 0 to 1, a -1 from 1 to 0
   %because we've put a 0 before and after we're sure to have the same number of transitions
   %from 0 to 1 and from one to 0
   runlength = find(transitions == -1) - find(transitions == 1);
   count = max(runlength);
end

And using it on a 30*30 matrix of mostly positive numbers

m = randi(20, 30)-6; %uniformly distributed between -5 and 14
[bestm, runcount, step] = findbest(m, 100000);

returns quite quickly.

When the number of possible answer diminishes, the algorithm becomes too expensive. At least, it'll still give you a better combination if not the best.

[Manoj]

Oh I understood what you were saying about positive numbers now. Unfortunately , my problem has very few solutions since there are quite a lot of negative numbers. So the solution which you have provided with randperm may not be ideal in this case. That is why I do not have any shortcut methods to reduce the search area for the solution.

[Manoj]

I was thinking of other ways to solve this problem and reduce the search area. What if I set a criteria for the bigger square matrices(>20*20) where the criteria would be a number for these continuous values.

Eg: I set criteria==15, then I just need a sorted matrix which will give me 15 continuous values next to the diagonal.

Would this make it easier or would I be complicating it further?

Thanks for all your help

Answer 3

[Matt J]

I'm not sure how robust an approach this is, but it did give me one of the non-unique solutions to your posted example

    result =

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

Basically, the idea is to use fmincon from the Optimization Toolbox to solve this as a nonlinear minimization problem. The code below searches for the permutation matrix P such that if A is the original matrix, then the permuted off-diagonal

     b=diag(P*(A>=0)*P.',1)

will have a large sum and be smooth. I would have liked to have added the nonlinear constraint

     P*P.'-eye(size(P))=0

to enforce the fact that P should be orthogonal, but when I do that, I find the optimization gets trapped very easily at local minima or non-solutions. Not so surprising, since this is a highly non-convex constraint. Without the non-linear constraint, though, the whole thing is really a quadratic program, and quadprog() would probably perform even better. I wouldn't take the trouble to recast the input data into quadprog's syntax until you're sure that it's working, though.

    A=   [-1  1  2 -1  3    % 1
          0 -1 -1 -1  4    % 2 
          5 -1 -1  6  7    % 3 
         -1 -1  8 -1  9    % 4
         10 11 12 13 -1]   % 5 ]

     n=size(A,1);

     map=(A>=0);

     Aeq=  [kron(speye(n), ones(1,n));...  %column sum
            kron(ones(1,n), speye(n))];    %row sum
     beq=ones(2*n,1);   

     lb=zeros(size(A));
     ub=ones(size(A));

      [P0,fval] = fmincon(@(P) cost(P,map),speye(n),[],[],Aeq,beq,lb,ub);

     P=double(bsxfun(@ge,P0,max(P0,[],2))); %round to nearest permutation matrix

     if any(Aeq*P(:)-beq)  %check that it is a permutation
        error 'Approximate permutation not found' 
     end

     result=P*A*P.',

    function  f=cost(P,map) 

      b=diag(P*map*P.',1);

      f=sum(diff(b).^2)-sum(b);

[Manoj]

First of all thank you for your time and efforts to provide me with a solution to my problem. The code with fmincon seems to work for the example I have given here but unfortunately it does not seem to work for other square matrices which are larger (> 10*10).

I get the error ' Approximate permutation not found ' and no result. But I know that there is a solution for the problem when I go about doing it manually.

Could you please let me know why this is happening or suggestions to go about to get the solution.

Thanks once again

[Matt J]

Manoj,

The version below might be more robust, but I say again, I am not sure how confident we can be in this approach. I'd be curious to have one of your large matrix examples to play with, however, if you're willing to attach one or several of them in a .mat file and tell us what the optimal number of continuous non-negative elements are on the diagonal. Part of the problem might be that non-default options may need to be passed to fmincon.

   n=size(A,1);

   map=(A>=0);

   Aeq=  [kron(speye(n), ones(1,n));...  %column sum
          kron(ones(1,n), speye(n))];    %row sum
   beq=ones(2*n,1);   

   lb=zeros(size(A));
   ub=ones(size(A));

   P0=(speye(n));
   nlc=[];

   for i=1:2

    M = fmincon(@(P) cost(P,map),P0,[],[],Aeq,beq,lb,ub,nlc,optimset('Display','none'));

    P=double(bsxfun(@ge,M,max(M,[],2))); %round to nearest permutation matrix

       %second pass
       P0=P; 
       nlc=@nonlcon;

   end

    if any(Aeq*P(:)-beq)  %check that it is a permutation
      warning 'Approximate permutation not found' 
    end

   result=P*A*P.',

function  f=cost(P,map) 

    b=diag(P*map*P.',1);

      m=length(b);

    f=sum(diff(b).^2)-sum(b);

function [c,ceq]=nonlcon(P)

       ceq=[];
       c=P*P.'-eye(size(P));

[Manoj]

Hello

Thanks once again , I will attach some of these matrices in .mat file as early as possible once I get the optimal solutions for each of them manually.

[Manoj]

Hello Matt

I am attaching few of the large matrices in a mat file and also an excel file for better viewing. My aim apart from getting these diagonal to be non negative would also be to know the initial and final positions of the original rows.

I am not sure if this is of help but the way those numbers works is that suppose a position for ex. (1,2) has a -1 then the position (2,1) also will be a -1.

Similarly for ex. (3,4) is a non negative number say 1, then the position (4,3) will also be a non negative number (0,1,2..)

Thanks once again