Robert Israel

6522 Reputation

21 Badges

18 years, 182 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are replies submitted by Robert Israel

Just to expand on that.  Let A be the 50 x 50 Matrix with entries A[i,j] = 1 if your matrix has a nonzero entry in position [i,j], 0 otherwise.  If P is a 50 x 50 permutation matrix, you want to minimize
f = sum(sum(C[i,j]*P[j,i],i=1..50),j=1..50)  where C[i,j] = sum(A[k,j], k= i+1 .. 50).  The constraints on the entries
of P are sum(P[i,j], j=1..50) = 1 for each i, sum(P[i,j], i=1..50) = 1 for each j, and all P[i,j] >= 0.   A basic feasible solution will automatically have all P[i,j] in {0,1}, so we don't have to specify that.  However, LPSolve's solution might not be basic (it uses an interior-point method), so I'll check the result. 

For example:

with(LinearAlgebra): with(Optimization):  
# random sparse Matrix with entries 0 and 1
A:= RandomMatrix(50,50,density=0.1, generator=1);
U:= Matrix(50,50,(i,j) -> `if`(i<j,1,0));
C:= U.A;
P:= Matrix(50,50,symbol=p);
  # the objective


f:= Trace(C.P):

# the constraints
cons:= {seq(add(P[i,j],i=1..50)=1, j=1..50),
seq(add(P[i,j],j=1..50)=1, i=1..50)}: 

# solve the linear programming problem
S:=LPSolve(f, cons, assume=nonnegative,maximize=false);
PM:=map(round,eval(P,S[2]));
# this should be the permutation matrix
MatrixNorm(PM . PM^%T - IdentityMatrix(50));
# this should be 0 if PM is a permutation matrix

        0

AP := A . PM;     # the permuted Matrix

Just to expand on that.  Let A be the 50 x 50 Matrix with entries A[i,j] = 1 if your matrix has a nonzero entry in position [i,j], 0 otherwise.  If P is a 50 x 50 permutation matrix, you want to minimize
f = sum(sum(C[i,j]*P[j,i],i=1..50),j=1..50)  where C[i,j] = sum(A[k,j], k= i+1 .. 50).  The constraints on the entries
of P are sum(P[i,j], j=1..50) = 1 for each i, sum(P[i,j], i=1..50) = 1 for each j, and all P[i,j] >= 0.   A basic feasible solution will automatically have all P[i,j] in {0,1}, so we don't have to specify that.  However, LPSolve's solution might not be basic (it uses an interior-point method), so I'll check the result. 

For example:

with(LinearAlgebra): with(Optimization):  
# random sparse Matrix with entries 0 and 1
A:= RandomMatrix(50,50,density=0.1, generator=1);
U:= Matrix(50,50,(i,j) -> `if`(i<j,1,0));
C:= U.A;
P:= Matrix(50,50,symbol=p);
  # the objective


f:= Trace(C.P):

# the constraints
cons:= {seq(add(P[i,j],i=1..50)=1, j=1..50),
seq(add(P[i,j],j=1..50)=1, i=1..50)}: 

# solve the linear programming problem
S:=LPSolve(f, cons, assume=nonnegative,maximize=false);
PM:=map(round,eval(P,S[2]));
# this should be the permutation matrix
MatrixNorm(PM . PM^%T - IdentityMatrix(50));
# this should be 0 if PM is a permutation matrix

        0

AP := A . PM;     # the permuted Matrix

@alex_01 : since the candidates are in random order, the best candidate has the same probability (1/n) of being first, second, ..., n'th.  That's also where I get the (k-1)/(b-1): b-1 places where the best of the first b-1 candidates could be, all with the same probability, and k-1 of them are numbers 1 to k-1.

@alex_01 : since the candidates are in random order, the best candidate has the same probability (1/n) of being first, second, ..., n'th.  That's also where I get the (k-1)/(b-1): b-1 places where the best of the first b-1 candidates could be, all with the same probability, and k-1 of them are numbers 1 to k-1.

1) B is the random variable, b is a value of the random variable.

2) There are b-1 equally likely positions for the best of the first b-1 candidates.  If that candidate is in positions 1 through k-1, the strategy works (i.e. candidates number k to b-1 are not chosen since they are worse than that candidate, and then candidate b who is the overall best will be chosen).  If that candidate is in position k through b-1, he/she will be chosen instead of candidate b, and the strategy doesn't work.

3) You mean 1/n.  There are n possible positions for the best candidate, and they are all equally likely.

4) It's sometimes called the "total probability formula".  If A is any event, and B_1, ..., B_n are events forming a partition of the sample space (i.e. every possible outcome belongs to exactly one of B_1, ..., B_n, then
P(A) = sum(P(A|B_i)*P(B_i), i=1..n).  This is very basic probability theory: look in any probability text.

5) Because P(A_k | B = b) = 0 if b < k.

6) To make it a function of k.

7) I'm getting asymptotic formulas as n -> infinity for P(k) (which, when n and k are integers, is the probability that the strategy succeeds), for the cases k = n/e and k = n/e + 1 - 1/(2*e), and their difference.

8) Say k_1 = n/e and k_2 = n/e + 1 - 1/(2*e).  Then P(k_2) - P(k_1) is approximately  0.9051258449/n^2 as n -> infinity, verifying that A_{k_2} is a slightly better strategy than A_{k_1}.

1) B is the random variable, b is a value of the random variable.

2) There are b-1 equally likely positions for the best of the first b-1 candidates.  If that candidate is in positions 1 through k-1, the strategy works (i.e. candidates number k to b-1 are not chosen since they are worse than that candidate, and then candidate b who is the overall best will be chosen).  If that candidate is in position k through b-1, he/she will be chosen instead of candidate b, and the strategy doesn't work.

3) You mean 1/n.  There are n possible positions for the best candidate, and they are all equally likely.

4) It's sometimes called the "total probability formula".  If A is any event, and B_1, ..., B_n are events forming a partition of the sample space (i.e. every possible outcome belongs to exactly one of B_1, ..., B_n, then
P(A) = sum(P(A|B_i)*P(B_i), i=1..n).  This is very basic probability theory: look in any probability text.

5) Because P(A_k | B = b) = 0 if b < k.

6) To make it a function of k.

7) I'm getting asymptotic formulas as n -> infinity for P(k) (which, when n and k are integers, is the probability that the strategy succeeds), for the cases k = n/e and k = n/e + 1 - 1/(2*e), and their difference.

8) Say k_1 = n/e and k_2 = n/e + 1 - 1/(2*e).  Then P(k_2) - P(k_1) is approximately  0.9051258449/n^2 as n -> infinity, verifying that A_{k_2} is a slightly better strategy than A_{k_1}.

@kmadkins : did you install the Maple Advisor Database?  It sounds to me like you didn't.  Do other Maple Advisor Database commands work for you?

@kmadkins : did you install the Maple Advisor Database?  It sounds to me like you didn't.  Do other Maple Advisor Database commands work for you?

@Kamel : the problem with Classic is that it doesn't fill non-convex (actually non-starshaped) polygons correctly.

@Kamel : the problem with Classic is that it doesn't fill non-convex (actually non-starshaped) polygons correctly.

@hermitian: if you don't have the Student:-NumericalAnalysis package (which was introduced, I think, in Maple 13), you can use ApproximateInt in the Student:-Calculus1 package. So replace the first and last lines with

with(Student:-Calculus1);

and

ApproximateInt(f(x), x = 0 .. 9, partition = 9, method=trapezoid);

The other lines are just to produce a function where f(0), ..., f(10) have been defined, but f(x) for symbolic x has not.  You could do the same by, somewhat tediously, assigning the values one-by-one:

f(0) := 0;
f(1) := 1;
...
f(9) := 81;
f(10) := 100;

@hermitian: if you don't have the Student:-NumericalAnalysis package (which was introduced, I think, in Maple 13), you can use ApproximateInt in the Student:-Calculus1 package. So replace the first and last lines with

with(Student:-Calculus1);

and

ApproximateInt(f(x), x = 0 .. 9, partition = 9, method=trapezoid);

The other lines are just to produce a function where f(0), ..., f(10) have been defined, but f(x) for symbolic x has not.  You could do the same by, somewhat tediously, assigning the values one-by-one:

f(0) := 0;
f(1) := 1;
...
f(9) := 81;
f(10) := 100;

The solutions of the equation t = f(t) are fixed points of the recurrence t[n] = f(t[n-1]), but in order for the recurrence, started sufficiently close to a solution s, to converge to that solution, you want |f'(s)| < 1.  You'll get really fast (quadratic rather than linear) convergence if f'(s) = 0.  Try it with B = -2*C*s and y = 1+s*A-C*s^3,
where the desired solution is s.

The solutions of the equation t = f(t) are fixed points of the recurrence t[n] = f(t[n-1]), but in order for the recurrence, started sufficiently close to a solution s, to converge to that solution, you want |f'(s)| < 1.  You'll get really fast (quadratic rather than linear) convergence if f'(s) = 0.  Try it with B = -2*C*s and y = 1+s*A-C*s^3,
where the desired solution is s.

Sorry, that should be Surfplot.  And the correct link (from the posting by jakubi in
<http://www.mapleprimes.com/questions/38398-surface-Plots>) is
<http://www.mapleprimes.com/view.aspx?sf=96307/280321/surfplot.zip>.  That is a zip archive containing files Surfplot.ind, Surfplot.lib and 8933_data2.mw.  Extract these to a directory, open the worksheet 8933_data2.mw, and follow the directions there.  Here's the plot I get:

First 19 20 21 22 23 24 25 Last Page 21 of 187