Rearranging rows of a boolean matrix so that the diagonals are all non-zero if possible

Hello, just as a disclaimer I have a question about improving my code, not on how to get it to work.

I wrote a code that rearranges just the rows of a boolean matrix so that the diagonal values are all non-zero. This code gets pretty slow as the matrix being fed in gets pretty large. It starts lagging pretty bad with 25x25 matrices. I'm curious what I can do to speed my code up.

Here is the function that I have right now:

function [reAdmis, unmatchedTaskList] = rearrange(admis)
    [n,~] = size(admis);
    potAssignment = zeros(1,n);%potential assignments
    dList=1:n; %Decision list: elements are zeroed if chosen
    reAdmis = zeros(n,n); % reAdmis will be the answer
    AdmisTemp = admis;
    
    CC = sum(admis,1); %column count
    RC = sum(admis,2); %row count
    [~, colIndex] = sort(CC);
    [~, rowIndex] = sort(RC);
    [reAdmisTemp, dIndex] = reorganizeRow(AdmisTemp, rowIndex, n); %reorganizes the rows based on RC
    
    for i = 1:n
        for j = 1:n
            if reAdmisTemp(j,colIndex(i))
                reAdmis(colIndex(i),:)=admis(rowIndex(j),:);
                AdmisTemp(rowIndex(j),:)=zeros; %used row gets zeroed
                dList(dIndex(j))=0; % chosen element in dList is zeroed
                potAssignment(dIndex(j))=colIndex(i); %Potential assignment is updated
                AdmisTemp(:,colIndex(i))=zeros; %used column gets zeroed
                break
            end
        end
        
        CC = sum(AdmisTemp,1);
        RC = sum(AdmisTemp,2);
        [~, colIndex] = sort(CC);
        [~, rowIndex] = sort(RC);
        [reAdmisTemp, dIndex] = reorganizeRow(AdmisTemp, rowIndex, n);
    end
end

function [organized, dIndex] = reorganizeRow(A, dIndex, n)
organized = zeros(n,n);
for i = 1:n
    organized(i,:)=A(dIndex(i),:);
end
end

An example boolean matrix to use and its output could be:

admis = [1     0     1     0     0
         0     0     0     0     1
         1     1     0     0     1
         0     1     0     1     1
         1     1     1     1     1]
     
reAdmis = [1     1     0     0     1
           0     1     0     1     1
           1     0     1     0     0
           1     1     1     1     1
           0     0     0     0     1]
       

My hope is that this can be as fast as possible. I have no hard number for a time requirement, however my clocked speed using tic-toc was about 0.0003 seconds, which I'd like to be able to beat.

Thank you in advance for any help that you are able to provide

[Matt J]

This code gets pretty slow as the matrix being fed in gets pretty large.

Can you put a number to "pretty large"? Are we talking 10000x10000?

[2Lindsay2]

As I mentioned above, I start to notice it getting very slow when I use matrices above about 25x25 in size. The speed is somewhat dependant on how "nice" the matrix is so this is just a general size that I'm giving here.

[Matt J]

Yes, but how large are you hoping to make the matrix? Would it get much larger than 25x25?

[2Lindsay2]

I'd like to be able to solve 200x200 reasonably fast (as fast as possible)

Answer 1

[Bruno Luong]

Your problem belong to the so-called "Assignment problem", that has dedicated algorithm called Hungarian assigment.

Take a look on the formal definition, try to understand it then find some implementation on File Exchange.

[Matt J]

I don't think it applies, because here we are permuting only the rows. There is no re-ordering of the columns.

[Bruno Luong]

@Matt not sure if you where you read assigment problem seach for double permutation? Have you read what it solve?

% Generate random binary matrix B that suitable for existing solution
n = 200;
B = (sprand(n,n,0.1)+speye(n)) > 0;
B = full(B(randperm(n),:));

% Hungarian assignnment, or Munkres algorithm
% https://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
tic
p = assignmentoptimal(1-B.');
toc % Elapsed time is 0.010675 seconds.

% Check
ax=subplot(1,2,1);
imagesc(ax,B);
title(ax,'Before permutation')
ax=subplot(1,2,2);
imagesc(ax,B(p,:));
title(ax,'After permutation')
colormap gray

Results

[Matt J]

not sure if you where you read assigment problem seach for double permutation?

I read it here:

https://en.wikipedia.org/wiki/Hungarian_algorithm#Matrix_formulation

[Bruno Luong]

I think I see what disturbs you. The equation of the perturbed matrix

L*C*R 

seems to indicate to allow both row and column permutations.

BUT actually it does only matter if C is not square.

However when C is square, then indeed

Trace(L*C*R)= = Trace((R*L)*C) 

So actually the algorithm finds also the row-only permutation (R*L) that minimizes the trace.

In other word, we don't have more effective freedom to have double permutation than single permutation, as the intuition apparently guides.

[Matt J]

Very good! Except I think you really mean,

Trace(L*C*R)= = Trace((R*L)*C) 

[Bruno Luong]

Correct, I'll edit my post. Thanks

Answer 2

[Matt J]

When your version works, it often works faster than the version below, however your version doesn't always get it right. I've attached a .mat file with an example matrix admis for which your version does not manage to get all non-zero diagonal elements.

load admis700 admis

N=size(admis,1);     
     
P=optimvar('P',[N,N],'LowerBound',0,'UpperBound',1,'Type','integer');

prob=optimproblem('ObjectiveSense','max');
prob.Constraints.row=sum(P,1)==ones(1,N);
prob.Constraints.col=sum(P,2)==ones(N,1);

E=speye(N);

opts=optimoptions('intlinprog','Display','none');

tic;
prob.Objective=sum(sum((P*admis).*E));

sol=solve(prob,'options',opts);
reAdmis=round(sol.P)*admis;
toc;

[2Lindsay2]

I'm interested in exploring this more. I'm unable to get solve() to run, I get an error saying to use linprog or intlinprog instead, so I haven't been able to verify what you've said.

I havn't been generating any matrices that cause it to fail, that's not to say that it doesn't. How are you generating the matrix to verify that there is a solution? If there is a matrix that you've found that is smaller in size (below 10x10 so I can go through it by hand easier) that causes it to fail I'd appreciate if you could share that.

[Matt J]

I'm interested in exploring this more. I'm unable to get solve() to run, I get an error saying to use linprog or intlinprog instead, so I haven't been able to verify what you've said.

I can't comment on that without seeing a copy/paste of the error message. What version of Matlab are you using? Presumably you have the Optimization Toolbox if the error message would suggest intlinprog.

How are you generating the matrix to verify that there is a solution?

I came across it by accident as I was running repeated comparisons on randomly generated matrices,

admis=double(sprand(N,N,7/N)>0);

If there is a matrix that you've found that is smaller in size (below 10x10 so I can go through it by hand easier) that causes it to fail I'd appreciate if you could share that.

Attached. In my tests, your method achieves 9 positive diagonal elements whereas my proposal finds 10. For N=10, I had to increase the density of ones in my random search,

N=10;
admis=double(sprand(N,N,3/N)>0);

[2Lindsay2]

The exact error is as follows:

==============================================

Error using optim.problemdef.OptimizationProblem/solve

options is not a valid solver. Use 'linprog' or 'intlinprog' instead.

Error in choiceMaker (line 13)

sol=solve(prob,'options',opts);

==============================================

I am using Matlab 2017b

[Matt J]

Here's a version requiring only intlinprog

N=size(admis,1);     
     
E=speye(N);
f=kron(admis.',-E);
 f=sum(f(1:N+1:end,:),1);
intcon=1:N^2;
Aeq=[kron(E,ones(1,N));kron(ones(1,N), E)];
beq=ones(size(Aeq,1),1);
lb=zeros(N);
ub=ones(N);
opts=optimoptions('intlinprog','Display','none');

tic;
P=intlinprog(f,intcon,[],[],Aeq,beq,lb,ub,[],opts);
   P=reshape(P,N,N);

reAdmis=round(P)*admis;
toc;

[2Lindsay2]

Thank you! I notice though that there is no guarantee that the sparse matrix that gets produced actually has a solution where all the diagonal elements are one. When there is not a solution where all the diagonal elements are one my code stops inserting the rows into the empty matrix "reAdmis", which is why the results would be different. When you tested this did you verify there was a solution where all the diagonal elements were one?

[Matt J]

The randomization procedure that I used to explore different choices of admis did not ensure this, however, the example failure cases admis10.mat and admis700.mat that I provided for you did have solutions where all the diagonal elements are one, as you will see when you run either my code or Bruno's.

In any case, I do not think there is any longer any reason for you to focus on either my solution or your own. Bruno's solution works and, for large N, is the fastest among all 3 methods. To demonstrate, I attach a 3rd data example, for which your solution does happen to work:

 load admis200

 tic;
  R1=Lindsay(admis);
 toc;
 
 tic
  R2=MattJ(admis);
 toc
 
 tic
  R3=Bruno(admis);
 toc
 

 
Tr1=trace(R1)
Tr2=trace(R2)
Tr3=trace(R3)

results in

Elapsed time is 0.030986 seconds.
Elapsed time is 0.135317 seconds.
Elapsed time is 0.015373 seconds.

Tr1 =

   200


Tr2 =

   200


Tr3 =

   200