NMFcomputations.mws

 >

This sheet performs some computations needed to verify various claims in the paper.

The non-math text in this sheet refers to objects discussed in the paper, e.g., r_1, ..., r_6, P, q_1^*, ...

The sheet is organized in sections that roughly correspond to sections of the paper.

The matrix factorizations claimed in the paper are verified in "Section 4.1 -- Geometry".

 > restart;

 > with(LinearAlgebra):

The following function maps a (nonnegative) matrix to a stochastic (i.e., column-stochastic) one, by normalizing the columns.

 > makeStoch := lA -> Matrix([seq(Normalize(lA[1..RowDimension(lA),i],1), i=1..ColumnDimension(lA))]):

Section 4.1 -- Geometry

The columns of the following matrix are the "inner points" r_1, ..., r_6.

 > InnerPoints := Transpose(Matrix([ [3/4,1/8,0], [3/4,1/2,0], [3/11,17/22,0], [2,0,1/2], [1/2,0,3/4], [1/6,0,7/12]]));

The columns of the following matrix contain the vertices of the polyhedron P.

 > OuterPoints := Transpose(Matrix([ [0,0,0], [1,0,0], [1,1/2,0], [0,1,0], [9/4,0,1/2], [0,0,8/7]]));

The columns of the following matrix contain the "intermediate points" q_1^*, ..., q_5^*.

 > MediatePoints := Transpose(Matrix([ [2-sqrt(2),0,0], [(3-sqrt(2))/7,(11+sqrt(2))/14,0], [1,(3+sqrt(2))/14,0], [0,0,(10+sqrt(2))/14], [(26+7*sqrt(2))/17,0,(12-2*sqrt(2))/17]])); evalf(%);

 > Hpre := Matrix([ [h11,0  ,h13,h14,0  ,h16], [0  ,h22,h23,0  ,0  ,0  ], [h31,h32,0  ,0  ,0  ,0  ], [0  ,0  ,0  ,0  ,h45,h46], [0  ,0  ,0  ,h54,h55,0  ]]);

The following matrix (denoted H below) is the matrix H' from the paper.

 > ({seq(seq(InnerPoints[i,j] = (MediatePoints . Hpre)[i,j], i=1..3), j=1..6)}): solve(%): H := subs(%, Hpre);

The following matrix is the matrix C.

 > C := Matrix([ [0, 10/11, 0], [0, 0    , 4/11], [-1/11, -2/11, 1/22], [-1/11, 0, 5/22], [4/11, 0, 0], [-2/11, -8/11, -7/11]]);

The following vector is the vector d.

 > d := Vector([0, 0, 2, 1, 0, 8]) / 11;

 > Vector[row]([1,1,1,1,1,1]) . C; Vector[row]([1,1,1,1,1,1]) . d;

The following matrix (denoted M below) is the matrix M' from the paper. We verify that it is stochastic.

 > M := C . InnerPoints + Matrix([d,d,d,d,d,d]); makeStoch(M) - M;

The following matrix is the matrix W. We verify that it is stochastic. We verify that M' = W * H'.

 > W := C . MediatePoints + Matrix([d,d,d,d,d]);   makeStoch(W)-W, simplify(M - W . H); evalf(W);

The columns of the following matrix contains the points q_1^epsilon, ..., q_5^epsilon.

 > EpsPoints := Matrix([ [99/169, 121/534, 9337/9338, 1/42216, 813/385], [0, 133/150, 64/203, 0, 0], [1/40560, 0, 0, 17209/21108, 997/1848]]);

 > Hepspre := Matrix([ [1-h41-h51,1-h22-h32,1-h23-h33,1-h44-h54,1-h45-h55], [0  ,h22,h23,0  ,0  ], [0  ,h32,h33,0  ,0  ], [h41,0  ,0  ,h44,h45], [h51,0  ,0  ,h54,h55] ]); Vector[row]([1,1,1,1,1]);

 > EpsPoints - MediatePoints . Hepspre: seq(seq(%[i,j], j=1..5), i=1..3): solHeps := solve({%}):

The following matrix is the matrix H_epsilon. We verify that it is stochastic.

 > Heps := subs(solHeps, Hepspre); makeStoch(Heps) - Heps; evalf(Heps);

The following matrix is the matrix W_epsilon. We verify that it is stochastic.

 > Weps := C . EpsPoints + Matrix([d,d,d,d,d]); makeStoch(Weps) - Weps; evalf(Weps);

We verify that W_epsilon = W * H_epsilon.

 > simplify(Weps - W . Heps);

Section 4.2 -- Type 1

In the following we perform determinant computations pertaining to the polygon P_0 (type-1 case).

 > Matrix([[u,0,1],[1,v/2,1],[3/4,1/8,1]]); det1 := Determinant(%); Matrix([[1, v/2, 1], [1-w, 1/2 + w/2, 1], [3/4, 1/2, 1]]); det2 := Determinant(%); Matrix([[1-w, 1/2 + w/2, 1], [u, 0, 1], [3/11, 17/22, 1]]); det3 := Determinant(%); eliminate({det1, det2}, {v, w}); f := simplify(subs([%[1]][1], det3)); solve(%);

 > plot(f, u=0..0.6);

In the following we perform determinant computations pertaining to the polygon P_1 (type-1 case).

 > Matrix([[u,0,1], [(9-9*v)/4, (7+9*v)/14, 1], [2, 1/2, 1]]); det1 := Determinant(%); Matrix([[(9-9*v)/4, (7+9*v)/14, 1], [0, (8-8*w)/7, 1], [1/2, 3/4, 1]]); det2 := Determinant(%); Matrix([[(0,8-8*w)/7, 1], [u, 0, 1], [1/6, 7/12, 1]]); det3 := Determinant(%); eliminate({det1, det2}, {v,w}); f := simplify(subs([%[1]][1], det3)); solve(%);

 > plot(f, u=0.5..1);

Section 4.5 -- Type 4

The following functions up and lo round up and down to n decimal digits.

 > up := (x,n) -> evalf(ceil(x*10^n)/10^n): lo := (x,n) -> evalf(floor(x*10^n)/10^n):

 > up(0.123456,3); lo(0.123456,3);

In the following we verify computations of bounds for the type-4 case.

 > W12lo := lo(W[1,2],2);

 > W13lo := lo(W[1,3],3);

 > W13up := up(W[1,3],3);

 > W24lo := lo(W[2,4],2);

 > W25lo := lo(W[2,5],3);

 > W33lo := lo(W[3,3],4);

 > W34lo := lo(W[3,4],2);

 > W35up := up(W[3,5],3);

 > W42lo := lo(W[4,2],3);

 > W44lo := lo(W[4,4],2);

 > W45up := up(W[4,5],3);

 > W55lo := lo(W[5,5],3);

 > W61lo := lo(W[6,1],2);

 > W63up := up(W[6,3],2);

 > W64up := up(W[6,4],2);

 > ex := 5; eps := 10.^(-ex);

 > R24lo := W24lo;

 > R25lo := W25lo;

 > L32up := up(W35up/R25lo,3);

 > L42up := up(W45up/R25lo,2);

 > L62up := up(eps/R25lo,ex);

 > L52up := up(eps/R24lo,ex);

 > L22lo := lo(1 - L32up - L42up - L52up - L62up,2);

 > R21up := up(eps/L22lo,ex);

 > L63lo := lo(W61lo - L62up*R21up,2);

 > L54lo := lo((W55lo - L52up)/(1-R25lo),4);

 > L34up := 1 - L54lo;

 > 2*L32up, 2-3*L63lo, 2*L34up;

 > L35lo := lo( (2*W34lo - W64up - max(2*L32up, 2-3*L63lo, 2*L34up))/2, 3);

 > L44up := 1 - L54lo;

 > 2*L42up, 2 - 3*L63lo, 2*L44up;

 > L45lo := lo((2*W44lo - W64up - max(2*L42up, 2 - 3*L63lo, 2*L44up))/2,3);

 > L11lo := W12lo;

 > R12lo := W12lo;

 > R13lo := W13lo;

 > L31up := up(eps/R12lo,ex);

 > R13up := up(W13up/L11lo,2);

 > R33up := up(W63up/L63lo,2);

 > R53up := up(eps/L45lo,ex);

 > R43lo := lo(1 - R13up - R33up - R53up,2);

 > L44up := up(eps/R43lo,ex);

 > L41up := up(eps/R13lo,ex);

 > R52up := up(eps/L35lo,ex);

 > L43lo := lo( (W42lo - L41up - L44up - R52up) / (1-R12lo), 3);

 > maxlo := lo( (W33lo - L31up - R53up) / (1-R13lo), 4);

 > is(maxlo + L43lo + L63lo > 1);

 > is(maxlo + L54lo > 1);

 >