Thread Subject:

Upwind Gauss-Seidel algorithm

I wonder if someone can help with an algorithm I'm trying to write. This is in the context of a dynamic programming problem, in which I'm trying to implement an Upwind Gauss-Seidel (UGS) algorithm. I understand that there are alternative methods like Simulated Upwind Gauss-Seidel (SUGS) that avoid the re-ordering problem, but right now I want to try this algorithm.

PROBLEM:
I have a square matrix Q with n rows. Q is a matrix of transition probabilities, so that all elements are between 0 and 1, and the elements of each row sum to 1. I would like to re-order the rows of Q such that Q_ij >= Q_ji for all i ~= j. Any suggestions on an efficient algorithm?

In the general case I'm thinking about, I don't know the size of Q, but it could easily be large. It's not obvious to me that a solution will exist for all Q, so the sooner I can recognize this the better. If I can get 'close' to a solution, that would also be okay, though I have to define 'close': if there are n rows then there are a total n-1 comparisons, so I might define 'close' as the number of comparisons that do satisfy the constraint divided by the total number of comparisons.

It might also be worth noting that to solve a given problem I'll need to re-order the matrix many times. The first time, there's no telling much re-ordering will be needed. But I believe that in later iterations, the amount of re-ordering needed will be less, so an algorithm that knows how to start by looking at 'minor' re-orderings before major re-orderings could be useful.

ALGORITHM 1: n small
Use 'perms' to generate all possible orderings of the n rows. Loop through each ordering, and keep the 'closest' ordering found so far, or break after the first ordering that satisfies the constraint. Obviously, this is constrained by the use of perms.

ALGORITHM 2: n large
While the constraint is unsatisfied, generate a random ordering of the rows. If the ordering satisfies the constraints, end the loop. Else remember that ordering if it 'closer' than any previous ordering and continue the loop. But this algorithm also seems crummy.

CODE (Algorithm 1 for different size n):
% Functions to check how close we are to a solution
ij_fun = @(n)(n+1)*[1:n-1]';
ji_fun = @(n)2+(n+1)*[0:n-2]';
close_fun = @(Qij,Qji)sum(Qij >= Qji) / length(Qij);
tests = 100;

% Dimensions
for n = 4:9
  tic;
  % Simulate problem
  ij = ij_fun(n); ji = ji_fun(n);
  P = perms([1:n])'; numP = size(P,2);
  cnt = 0; satisfied = 0;

  % Generate a bunch of random probability matrices
  for ii = 1:tests
    Q = rand(n); Q = bsxfun(@times,Q,1./sum(Q,2));

    % Algorithm 1: Check all permutations
    close = 0; order = []; jj = 0;
    for p = P
      jj = jj+1;
      Qp = Q(p,:); Qp_ij = Qp(ij); Qp_ji = Qp(ji);

      closep = close_fun(Qp_ij,Qp_ji);
      if (closep > close), close = closep; order = p; end
      if (close == 1), break; end
    end

    cnt = cnt + jj;
    satisfied = satisfied + (close == 1);
  end
  avg = cnt / tests;
  disp(['n = ' num2str(n)]);
  disp(['On average, searched ' num2str(avg) ' out of ' num2str(numP) ' permutations.']);
  disp([num2str(satisfied) ' out of ' num2str(tests) ' random matrices satisfied.']);
  toc;
  disp('=====================================================');
end

Any suggestions welcome. Thanks.

"R" wrote in message <k119lg$9m2$1@newscl01ah.mathworks.com>...
> PROBLEM:
> I have a square matrix Q with n rows. Q is a matrix of transition probabilities, so that all elements are between 0 and 1, and the elements of each row sum to 1. I would like to re-order the rows of Q such that Q_ij >= Q_ji for all i ~= j. Any suggestions on an efficient algorithm?
>

Quick clarification: That is, I need Q_ij >= Q_ji for all i,j such that j = i+1. So I need Q_12 > Q_21, Q_23 > Q_32, etc.

I tried to find an algorithm of this problem, and so far I haven't find anything better than brute-force method. Here is a vectorized version, which is faster by permutation analysis, but does not save much time (since it can't break when a solution is found). I still provide it bellow:

tests = 100;

% Dimensions
for n = 4:9
  tic;
  % Simulate problem
  P = perms(uint32(1:n))'; numP = size(P,2);
  satisfied = 0;

  % Generate a bunch of random probability matrices
  for ii = 1:tests
    Q = rand(n);
    Q = bsxfun(@times,Q,1./sum(Q,2));

    % Algorithm 1: Check all permutations
    ij = bsxfun(@plus,uint32(n*(1:n-1)'),P(1:end-1,:));
    Qij = Q(ij);
    ji = bsxfun(@plus,uint32(n*(0:n-2)'),P(2:end,:));
    Qji = Q(ji);
    closep = sum(Qij>=Qji, 1) / (n-1);
    [close jbest] = max(closep);
    order = P(:,jbest);
    satisfied = satisfied + (close == 1);
  end
  disp(['n = ' num2str(n)]);
  disp([num2str(satisfied) ' out of ' num2str(tests) ' random matrices satisfied.']);
  toc;
  disp('=====================================================');
end

% Bruno

Here is a better method, significantly faster in average.

tests = 100;

% Dimensions
for n = 4:9
  tic;
  % Simulate problem
  P = perms(uint32(1:n))'; numP = size(P,2);
  satisfied = 0;

  % Generate a bunch of random probability matrices
  for ii = 1:tests
    Q = rand(n);
    Q = bsxfun(@times,Q,1./sum(Q,2));

    % Algorithm 1: Check all permutations
    [found order] = rowperm(1, Q);% mfile bellow
    satisfied = satisfied + found;
  end
  disp(['n = ' num2str(n)]);
  disp([num2str(satisfied) ' out of ' num2str(tests) ' random matrices satisfied.']);
  toc;
  disp('=====================================================');
end

%% Save this in an m-file
function [found p] = rowperm(r, Q)

persistent p_
persistent Q_
persistent n_

if r==1
    % initialization persistent variables
    n_ = size(Q,1);
    Q_ = Q;
    p_ = zeros(1,n_);
    keep = 1:n_;
elseif r > n_
    p = p_;
    found = true;
    return
else
    keep = Q_(:,r-1) <= Q_(p_(r-1),r);
    keep(p_(1:r-1)) = false;
    keep = find(keep)';
end

for k = keep
    p_(r) = k;
    % call recursive
    [found p] = rowperm(r+1);
    if found
        return
    end
end

% fail
p = p_;
found = false;

end % rowperm

% Bruno

To avoid noise, this statement should be removed in the last algorithm.

> P = perms(uint32(1:n))'; numP = size(P,2);

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <k12vo3$egf$1@newscl01ah.mathworks.com>...
> To avoid noise, this statement should be removed in the last algorithm.
>
> > P = perms(uint32(1:n))'; numP = size(P,2);
>
> Bruno

This is nicely done. In the brute force algorithm for an 8x8 matrix, I might check all 5040 permutations in which row 1 was the first row of the solution matrix before ruling out row 1 as a candidate. If I understand Bruno's algorithm correctly, I might be able to rule out row 1 after making as few as 7 comparisons. The recursive program then extends this logic for all following rows of the solution matrix.

Also, Bruno's algorithm uses the initial matrix as a starting point, so if we have reason to think the initial matrix is already 'close' to a good solution, the search should be short.

This is great.

I might point out that a trade-off of the smarter design is that for any permutation that falls short of a full-solution, you only have a lower bound estimate for how 'close' it is to a complete solution. Which is too bad if it's useful or necessary to cut the search short for large n while tracking the best candidate solution up to that point (though it's a very reasonable trade-off).

One follow-up question. When I tried Bruno's code for a medium size Q matrix (n=500), I ran up against the recursion limit. How worried should I be about just increasing the recursion limit? I think the concern is memory, but since Q_ is persistent and that's probably the only thing that requires memory, can I think of this as a negligible problem? Or is there a lot of memory overhead from the recursion process anyway? In which case for large Q it would be necessary to re-write rowperms() as iterative rather than recursive?

In any case, thanks to Bruno, and to anyone else who took a look.

"R" wrote in message <k13nj8$g9t$1@newscl01ah.mathworks.com>...
>
> One follow-up question. When I tried Bruno's code for a medium size Q matrix (n=500), I ran up against the recursion limit. How worried should I be about just increasing the recursion limit? I think the concern is memory, but since Q_ is persistent and that's probably the only thing that requires memory, can I think of this as a negligible problem? Or is there a lot of memory overhead from the recursion process anyway? In which case for large Q it would be necessary to re-write rowperms() as iterative rather than recursive?

Actually the goal of using persistent is that there is only one copy of Q regardless the number of levels of recursion. So it shouldn't have concern abut Q. The only parameters that are stacked are r, k, and keep.

I have tested with n=500 without any problem on my 2Gbytes-RAM laptop.

Bruno