Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Upwind Gauss-Seidel algorithm

Subject: Upwind Gauss-Seidel algorithm

From: R

Date: 22 Aug, 2012 00:38:08

Message: 1 of 7

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.

Subject: Upwind Gauss-Seidel algorithm

From: R

Date: 22 Aug, 2012 01:07:07

Message: 2 of 7

"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.

Subject: Upwind Gauss-Seidel algorithm

From: Bruno Luong

Date: 22 Aug, 2012 12:14:08

Message: 3 of 7

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

Subject: Upwind Gauss-Seidel algorithm

From: Bruno Luong

Date: 22 Aug, 2012 13:46:09

Message: 4 of 7

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

Subject: Upwind Gauss-Seidel algorithm

From: Bruno Luong

Date: 22 Aug, 2012 16:01:07

Message: 5 of 7

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

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

Bruno

Subject: Upwind Gauss-Seidel algorithm

From: R

Date: 22 Aug, 2012 22:48:08

Message: 6 of 7

"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.

Subject: Upwind Gauss-Seidel algorithm

From: Bruno Luong

Date: 23 Aug, 2012 06:26:10

Message: 7 of 7

"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

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us