# 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 ... > 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" wrote in message ... > 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 ... > > 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

### Everyone's Tags:

Separated by commas
Ex.: root locus, bode

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