Optimize table layout

Question

I have a set of cells arranged in a table. I need to minimize the table's height by adjusting column widths.

Each cell has an area such that its area is not encroached upon as its width and height are adjusted. In other words, given a final row height $h_i$, final column width $w_j$, and initial cell area $a_{ij} \in A$, this must hold true:

$$ \forall a_{ij} \in A : a_{ij} \leq h_i \cdot w_j $$

Since it's a table, each cell in a column has the same width and each cell in a row has the same height. Additionally, each row has the same width which is a chosen parameter $W$ (the width of the table). Thus:

$$ W = \sum_j w_j $$

...and the table will have this overall height:

$$ H = \sum_i h_i $$

So given $A$ (and knowing its dimensions), I need to compute all $w_j$ such that $H$ is minimized.

Minimum height of two-column table

Consider a two-column table with cell areas like the below. For simplicity, the table has a total width of 1. $p$ is the width of the first column; $1-p$ is the width of the second column; and column widths cannot be zero (so $0 < p < 1$):

  p  1-p
|<->|<->|

+---+---+
| a | b |
+---+---+
| c | d |
+---+---+
| e | f |
+---+---+
|...etc |

The height of the first row will be: $$ \cases{ p \leq \frac{a}{a+b} : \frac{a}{p} \\ p > \frac{a}{a+b} : \frac{b}{1 - p} } $$

...and of the second: $$ \cases{ p \leq \frac{c}{c+d} : \frac{c}{p} \\ p > \frac{c}{c+d} : \frac{d}{1 - p} } $$

...and so on. Notice how the left cell's area is considered (with one denominator) when $p$ is small enough; otherwise the right cell's area is used (with a different denominator).

Let's suppose that things are such that for a given $p$ these cells' areas are used: $( a, d, e, \ldots )$. This would be the table's height: $$ \frac{a}{p} + \frac{d}{1 - p} + \frac{e}{p} + \ldots $$

Let's take a moment to generalize this. Add up all the areas chosen from the left side and call that $l$, and $r$ for all areas from the right side. Thus: $$ H = \frac{l}{p} + \frac{r}{1 - p} $$

Now we want to minimize the table's height by finding the best $p$. So take the derivative and set it to zero: $$ 0 = \frac{d}{dp} H = \frac{r}{(1-p)^2} -\frac{l}{p^2} $$ $$ = r \cdot p^2 - l \cdot (1 - p)^2 $$ $$ = (r - l) \cdot p^2 + 2l \cdot p - l $$

Here are the solutions to this quadratic equation: $$ p = \cases{ l \neq r : \frac{-2l \pm \sqrt{4l^2 +4l(r-l)}}{2(r - l)} \\l = r : 0.5 }$$

Plug each of the solved $p$ back into $H$ to figure out which is best.

So now all you have to do is decide, for a given range of $p$, which cells contribute to $l$ and which cells contribute to $r$, and then use the above equations. The best table height from all ranges of $p$ is the global minimum table height.

I say "for a given range of $p$" because for every row we know the range of $p$ for which the left cell should be considered versus the right cell. For example, we know that cell $a$ should be added to $l$ when $p \leq \frac{a}{a + b}$. That means the first row contributes two possible ranges of $p$ that need to be checked (and $\frac{a}{a + b}$ is the boundary). The second row contributes another two possible ranges (with the boundary at $\frac{c}{c + d}$). And so on. In each range different cell areas are contributing to $l$ and the rest are contributing to $r$.

In other words, if there are $x$ table rows then there are up to $2x$ different equations for $H$ that you need to solve to find the minimum height of a two-column table.

But I do not know how to generalize this into more columns

False starts

#1

Here's an algorithm which at first glance might seem to do the trick. But it only works for certain table configurations. For example, this does not work when the diagonal cells begin as "king" cells.

1. Lay out the table so that all rows are tightly stacked (meaning no row exists in which all cells in that row have elbow room). At this point it doesn't matter how wide the table is. As a consequence of this some columns will be too wide
2. Select the first column
3. For every cell in the column, calculate the amount the column can be shrunk $\Delta w = w_y - a_i / h_x$ such that this cell will have no elbow room
4. Find the minimum $\Delta w > 0$ (if any) of the column
5. Shrink the column by that amount
6. Select the next column and repeat from #3 until all columns have been adjusted
7. Scale the table to the desired width, preserving relative column proportions
8. Recalculate row heights based on the final column widths

This comes from the intuition that when a table's rows all have minimum height then each row will have at least one "king" cell which has no elbow room and will only increase that row's height if its column is collapsed further. Therefore the table has to get taller if any "king" cell's column is shrunk. But that only covers columns in which a "king" cell is present. The goal of this algorithm is to get "king" cells in all columns.

Once there's a "king" cell in each row and in each column then one would think that no column can be shrunk without a net increase in table height. One would think that increasing a row's height cannot be compensated by a decrease in another row's height because one would think all other rows already have minimum height.

But that's an incorrect intuition. While it may be true that no column may be shrunk in isolation, there may still exist the possibility to alter the widths of several columns together in such a way that the total table height is reduced.

Regardless, I do believe it's the case that the optimal column widths are still optimal when scaled together. So I believe steps 7 and 8 are valid.

To illustrate why this algorithm does not work, consider this 2x2 table:

+---+---+
| a |   |
+---+---+
|   | b |
+---+---+

In this case, the table has two empty cells on one diagonal and two populated cells on the other diagonal. Thus these two cells are guaranteed to be king cells, and the algorithm will traverse the columns without altering anything. In other words, the original column arrangement (whatever that happens to be) is the final column arrangement. The algorithm does nothing but push the problem of optimizing table layout elsewhere.

In this specific case it's possible to demonstrate that the ideal ratio of first column width to second column width is $\sqrt{a} : \sqrt{b}$. Yet this isn't the ideal ratio for all tables. So the problem remains unsolved in general.

#2

Given that the optimal column distribution for a two-column table may be found in O(rows^2) time (see above), I was hopeful for an easy way to append columns. But this doesn't appear to be feasible.

To illustrate this, consider this optimal table (roughly at scale):

+-+-------------+
|1|             |
+-+-------------+
| |             |
| |             |
| |     169     |
| |             |
| |             |
+-+-------------+

Since it's optimal, the first row height is $\sqrt{1} / \sqrt{169} = 7.7\%$ of the table height.

What happens when we append the following column to it?

+-----+
| 1e6 |
+-----+
|  0  |
+-----+

169 is peanuts compared to 1e6. And what are we going to do—place it in a row that's only 7.7% of the total table height while the other 92.3% goes to the cell with 169? Of course not! We'll want to give the second column proportionally more space so that it gets shorter and the 1e6 can grow taller/skinnier. And as the 1e6 grows taller we can give the first column proportionally less space (so that the height of the 1 is equal to the height of the 1e6).

In other words, appending a column requires laying out the entire table again. Meaning to lay out a three-column table you need to know how to lay out a three-column table. That doesn't really get us anywhere. And for the general case I think that would work out to O(rows^2 * columns!) time complexity.

Answer 1

I tried to implement Rahul's suggestion to view it as a convex optimization problem. The results are mixed. I can easily do small tables like 30 by 30, but 300 by 300 can be done with only about 1% precision if you are willing to wait 1 minute and getting down from there will take eternity. That is primarily due to the inefficiency of the solution finder I'm using (which is more or less just cycling over variables and optimizing over certain subsets of them; I wish I could find a better way or, at least, accelerate convergence somewhat). Nevertheless it is a good exercise in convex programming, so I'll post the details here. The algorithm can be modified to take into account "natural" restrictions of the kind $w_j\ge W_j$ or $h_i\ge H_i$ (width/height should not be too small) and the modification has pretty much the same rate of performance as far as I can tell from simulations, but I'll restrict myself to the original question here.

Let $w_j$ be the unknown widths and $a_{ij}$ be the known areas. We want to minimize $\sum_i\max_j \frac{a_{ij}}{w_j}$. It is useful to consider the dual problem. I will spare you from the general theory of duality and will just note that $$ \max_j \frac{a_{ij}}{w_j}=\max\left\{\sum_j b_{ij}\frac{a_{ij}}{w_j}:b_{ij}\ge 0, \sum_j b_{ij}=1\right\} $$ so if we consider all admissible vectors $w=(w_1,\dots,w_n)$ (non-negative entries, total sum $1$) and all admissible matrices $b=(b_{ij})$ (non-negative entries, all row sums equal to $1$), we can write our problem as that of finding $$ \min_w\max_b \sum_{i,j} b_{ij}\frac{a_{ij}}{w_j}\,. $$ The dual problem to that is finding $$ \max_b \min_w\sum_{i,j} b_{ij}\frac{a_{ij}}{w_j}\,. $$ The inner $\min_w$ is here easy to find: if we denote $S_j=\sum_i b_{ij}a_{ij}$, then it is just $(\sum_j \sqrt{S_j})^2$ with unique optimal $w_j$ proportional to $\sqrt{S_j}$.

There are two things one should understand about duality. The first one is that every admissible matrix $b$ (computed or just taken from the ceiling) can serve as the certificate of the impossibility to do better than a certain number in the original problem, i.e., the minimax is never less than the maximin. This is pretty trivial: just use the given $b$ to estimate the minimax from below. The second one is that the true value of minimax is actually the same as the true value of maximin (under some mild assumptions that certainly hold in our case). This is a somewhat non-trivial statement.

Together these two observations allow one to use the following strategy. We shall try to solve the dual problem. For every approximation $b$ to the solution, we will look at the easily computable lower bound $(\sum_j\sqrt{S_j})^2$ it produces and at the corresponding minimizer $w$. For that $w$ we can easily compute the original function $\sum_j\max_i\frac{a_{i,j}}{w_j}$. If its value is reasonably close to the lower bound, we know that we should look no further.

Now, of course, the question is how to maximize $\sum_j\sqrt S_j$ under our constraints on $b$. It doesn't look like an attractive problem because the number of unknowns increased from $n$ to $mn$. Still, one can notice that if we fix all rows of $b$ except, say, the $i$'th one, then the optimization of the $i$'th row is rather straightforward. Indeed, the corresponding problem is of the following kind:

**Find $\max\sum_j\sqrt{a_jb_j+c_j}$ where $a_j,c_j\ge 0$ are given and $b_j\ge 0$ are the unknowns subject to the constraint $\sum_j b_j=1$. Using the standard Lagrange multiplier mumbo-jumbo, we conclude that the optimal $b_j$ must satisfy the equations $\frac{a_{j}}{\sqrt{a_jb_j+c_j}}=\lambda$ whenever $b_j>0$ and the inequalities $\frac{a_{j}}{\sqrt{a_jb_j+c_j}}\le \lambda$ whenever $b_j=0$. Thus, the optimizer is just a vector $b_j=\max(\Lambda a_{j}-\frac{c_j}{a_j},0)$ with an unknown $\Lambda=\frac 1{\lambda^2}>0$ that should be found from the constraint $\sum_j b_j=1$. This is a one-variable equation for the root of a monotone function, so it can be easily solved in various ways.

Thus, we can optimize each row of $b$ with other rows fixed rather quickly. The natural idea is then just to cycle over rows optimizing each one in turn. It takes about 20 full cycles to get the lower bound and the value of the function within 1% range from each other on a random matrix (structured matrices seem to be even better) up to the size of 300 by 300 at least.

This is the description. The code (in Asymptote) is below.

srand(seconds());

int m=50, n=55;

real[][] a, b;
for(int i=0;i<m;++i)
{
    a[i]=new real[]; b[i]=new real[];
    for(int j=0; j<n; ++j)
    {
        a[i][j]=unitrand();
        b[i][j]=1/n;
    }
    //a[i][rand()%n]=2;
    a[i]*=unitrand();
}

real[] p, S;

for(int k=0;k<101;++k)
{
    for(int j=0;j<n;++j)
    {
        real s=0;
        for(int i=0;i<m;++i)
            s+=a[i][j]*b[i][j];
        S[j]=s;
        p[j]=sqrt(S[j]);
    }
    if(k%10==0)
    {
        write("*** Iteration "+string(k)+" ***");
        write(sum(map(sqrt,S))^2);
    }

    p/=sum(p);

    real s=0; 
    for(int i=0;i<m;++i)
    {
        real M=0; 
        for(int j=0;j<n;++j)
        {
            real h=a[i][j]/p[j];
            if(h>M)
                M=h;
        }
        s+=M;
    }
    if(k%10==0)
        write(s);
    //pause();

    for(int i=0;i<m;++i)
    {
        real[] A,V,C,B;
        for(int j=0;j<n;++j)
        {
            A[j]=a[i][j];
            V[j]=S[j]-a[i][j]*b[i][j];
            C[j]=V[j]/a[i][j];
        }
        real aa=(sum(C)+1)/sum(A);
        real da=1;
        while(da>1/10^10)
        {
            for(int j=0;j<n;++j)
            {
                B[j]=aa*A[j]-C[j];
                if(B[j]<0)
                {
                    A[j]=0;
                    B[j]=0;
                }
            }
            da=sum(B)-1; aa-=da/sum(A); 
        }
        for(int j=0;j<n;++j)
        {
            b[i][j]=B[j];
            S[j]=V[j]+a[i][j]*B[j];
        }
    }
}

write("************");

pause();

Answer 2

This problem can be solved readily using a convex programming library like CVX or CVXPY, after applying the transformation $a_{ij}\le h_i w_j \iff \log a_{ij} \le \log h_i + \log w_j$ to convert it to a convex problem. Here is CVXPY code for fedja's problem:

import cvxpy as cp
import numpy as np
from math import *

# Problem data.
m = 50
n = 37
# np.random.seed(0)
A = np.ones((m,n))
for i in range(m):
    for j in range(n):
        A[i,j] = 1 + 0.0001*sin(i + j*j)
wmax = 1

# Construct the problem.
h = cp.Variable((m,1))
w = cp.Variable((1,n))
objective = cp.Minimize(cp.sum(h))
H = cp.hstack([h for _ in range(n)])
W = cp.vstack([w for _ in range(m)])
constraints = [cp.log(A) <= cp.log(H) + cp.log(W), cp.sum(w) <= wmax]
problem = cp.Problem(objective, constraints)

problem.solve(verbose=True)

print("Optimal value", problem.value)
# print("Optimal widths", w.value)
# print("Optimal heights")
# print(h.value)

[...solver output, remove 'verbose=True' to hide...]
Maximum number of iterations reached, recovering best iterate (98) and stopping.

Close to OPTIMAL (within feastol=2.5e-07, reltol=5.8e-12, abstol=1.1e-08).
Runtime: 0.491104 seconds.

Optimal value 1850.1460524691356

(Note that this is not exactly a feasible solution: some constraints are violated by ${\sim}10^{-5}$. A feasible solution could be recovered by increasing the row heights slightly.)

Upper and lower bounds on $h_i$ and $w_j$ can easily be added as well.

Answer 3

Here is just the code that is (presumably) finding the exact answer for not too large matrices. All explanations will come later. Again it is in Asymptote. Why not in C? I know C++ but having a slow computer with an interpreted rather than compiled language allows you to watch the program as it goes (if you bother to output the work protocol to the screen, of course) and see many nuances that otherwise would be easily missed. I found at least 6 logical bugs in the original version this way (I hope the remaining number is less). The progress is not linear because the program tries to locate what Matt calls "king cells" in the matrix and the "vert=" line is a better mark of progress than the difference between the upper and the lower bounds. The final array of widths is $P$ and the program terminates when the relative error is $10^{-15}$ (be careful here: due to the rounding errors it may continue to run beyond that point, but once you see "no conflict" and "nothing to do" repeating again and again, the job is actually done and you are watching the geometric convergence game end (which can be also played differently, but who cares?).

srand(seconds());

int m=30, n=17, K=100001, k, SK=10, Count=0, proccnt=0, Failtime=1000000, I=0,J=0, cycletime=0; 
real M=0, Mdel=0, TTT=0, B;
int time=0, failtime=0, successtime=0; 
int tt=1, ttt=1, blcount=0, ulttol=3;


int start=seconds();

int[][] blacklist;
for(int i=0;i<m;++i) blacklist[i]=array(n,1);

real[][] a, aaa,  b , bb, bbb, db, dbb, bres;
real[] AA;

/// Initiating the matrix ////////////

real delta=0.0001;

for(int i=0;i<m;++i)
{
real u=unitrand(), v=unitrand();
a[i]=new real[]; b[i]=new real[];
for(int j=0; j<n; ++j) 
{
a[i][j]=1+delta*sin((i+j^2)); 
b[i][j]=1/n;
}
//a[rand()%(i+1)][rand()%n]=2;
//a[i]*=unitrand();
}

////////////////////////////////////////////


aaa=copy(a); bres=copy(b);
real kill=1/1000;



int [][] temp; bool[] conf=array(n,true); 
bool fast=true, notreset=true, confl=true;

for(int i=0;i<m;++i) temp[i]=array(n,0);

int[] Ind; for(int i=0;i<m;++i) Ind[i]=i;

real Norm(real[][] b)
{
real[] S;
for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=a[i][j]*b[i][j]; S[j]=s;
}
return sum(map(sqrt,S))^2;
}


void shuffle()
{
for(int kk=0;kk<m;++kk) {int a=rand()%m, b=rand()%m; int II=Ind[a]; Ind[a]=Ind[b]; Ind[b]=II;}
}

bool[] conflict(real[][] b)
{
bool[] conf=array(n,false);

int count=0; 

for(int i=0;i<m;++i) 
{
if(min(b[i])<0) {write("karaul"); pause();}
b[i]=max(b[i],array(n,0));
count+=sum(map(sgn,b[i]));
}
int[] pres=array(m,1);
int[][] sb;
for(int i=0;i<m;++i) {sb[i]=map(sgn,b[i]); sb[i][n]=1;}


for(int I=1;I<m;++I)
for(int i=0; i<I; ++i)
{
if(pres[i]>0 && sum(sb[i]*sb[I])>sb[i][n]*sb[I][n]) {pres[i]=0; sb[I]=sb[i]+sb[I];}
}

int vert,edgecnt,Vert=0,Edgecnt=0; int comp=sum(map(sgn,pres));
for(int i=0;i<m;++i) 
{
if(pres[i]>0) 
{
vert=sum(sb[i])-sb[i][n];
Vert+=vert;
edgecnt=-sb[i][n];
for(int j=0;j<n;++j) edgecnt+=max(2*sb[i][j]-1,0); 
Edgecnt+=edgecnt;
if(edgecnt>vert-1) for(int j=0;j<n;++j) {if(sb[i][j]>0) conf[j]=true;}
}
}
int alive=0; for(int i=0;i<m;++i) for(int j=0;j<n;++j)
if(conf[j] && b[i][j]>0 && blacklist[i][j]<=ttt) ++alive;
write("vert="+string(Vert)+"("+string(alive)+") edgecnt="+string(Edgecnt)+" comp="+ string(comp));
return conf;
}





real[] p, P, S;

for(k=0;k<K;++k)
{

void procedure()
{
for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=aaa[i][j]*b[i][j]; S[j]=s;
}
for(int i:Ind)
{
real aa;
real[] A,V,C,B;
for(int j=0;j<n;++j) {A[j]=aaa[i][j]; V[j]=S[j]-aaa[i][j]*b[i][j]; C[j]=V[j]/aaa[i][j];}
real aa=(k==0?(sum(C)+1)/sum(A):AA[i]);

int countbound=40;

for(int j=0;j<n;++j) B[j]=max(aa*A[j]-C[j],0);
if(sum(B)>1/2)
{
if(sum(B)<1)
{
real sl=0;
for(int j=0;j<n;++j) sl+=A[j]*sgn(B[j]);
aa+=1.0001*((1-sum(B))/sl); countbound=4;
}
}
else aa=(sum(C)+1)/sum(A);

real da=1;
int count=0;

while(da>0 && count<countbound)
{
++count; 
//write(da,count); //pause();
for(int j=0;j<n;++j) {B[j]=aa*A[j]-C[j]; if(B[j]<0) {B[j]=0; A[j]=0; C[j]=0;}}
if(sum(A)>0) {da=sum(B)-1; aa-=da/sum(A);}
else {write("alert"); pause(); for(int j=0;j<n;++j) {if(b[i][j]>0) A[j]=aaa[i][j];} aa=(sum(C)+1)/sum(A); } 
//write(sum(B),aa,da,sum(A),sum(C));
}
for(int j=0;j<n;++j) {b[i][j]=B[j]; S[j]=V[j]+aaa[i][j]*B[j];}
Count+=count; 

if(abs(sum(b[i])-1)>0.1) {write("rough!"); pause();}
AA[i]=aa; b[i]/=sum(b[i]);
}
++proccnt;
}

bool check()
{
bool check=false;
for(int i=0;i<m && !check;++i) for(int j=0;j<n;++j) check=check || (bres[i][j]>0 && b[i][j]==0);
return check;
}




void fix()
{
for(int i=0;i<m;++i) for(int j=0;j<n;++j) 
{
if(b[i][j]==0 && conf[j]) aaa[i][j]=a[i][j]*kill;
//if(b[i][j]==0) blacklist[i][j]=1;
}
}


void advance(bool adv=true)
{
for(int kk=0;kk<(adv?ttt:tt)*SK;++kk) procedure(); bres=copy(b); if(adv) {write("advancing with speed "+string(TTT)); fix();}
}


void reset(bool hard=true)
{
if(!confl) write("nothing to do"); else write("RESETTING "+(hard?"HARD":"SOFT")); 
fast=true; if(hard) blcount=0;   
//shuffle();
aaa=copy(a); for(int kk=0;kk<(confl && hard?ttt:1)*SK;++kk) procedure(); 
if(confl && hard) ttt*=2;  
fix(); 
}

real minb=1, shift=0;

TTT=1;

while (TTT>1/3) 
{ 
TTT=0;
//bbb=copy(b); advance(false); 
bb=copy(b); advance(false); bres=copy(b);

for(int i=0;i<m;++i) 
{
db[i]=b[i]-bb[i]; 
//dbb[i]=bb[i]-bbb[i]; 
shift=max(shift,max(map(abs,db[i]))); temp[i]=array(n,0);
}

for(int i=0;i<m;++i) for(int j=0;j<n;++j)
{
if(b[i][j]>0 && db[i][j]<0 && bb[i][j]>0) 
{
real u=-db[i][j]/b[i][j];
//v=-dbb[i][j]/bb[i][j]; 
if(u>TTT && u>0 && aaa[i][j]>a[i][j]/2 && blacklist[i][j]<=ttt) {TTT=u; I=i; J=j; minb=min(minb,b[i][j]);}
}
}
tt=(confl?blacklist[I][J]:1);
if(TTT>1/3) advance(); 
else if(TTT==0 || blcount>ulttol) reset();
else write('\n \naccelerating from speed '+string(TTT)+
"; position=("+string(I)+","+string(J)+"); cycle count= "+string(2*tt*SK)); 
}

time=seconds()-start; if(time>Failtime) {write('\n\nTotal failure'); pause(); Failtime*=2;} 

write("time= "+string(time)+", cycling "+string(cycletime)+
" seconds, failures =  "+string(failtime)+ ", successes= "+string(successtime));

write("count="+string(Count/m/proccnt)); 

conf=conflict(b);

for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=aaa[i][j]*b[i][j]; S[j]=s; p[j]=sqrt(s);  
}

p/=sum(p); 
if(k==0) P=copy(p); 
write(Mdel); 

{
real s=0, sss=0; 
for(int i=0;i<m;++i)
{
real M=0; 
for(int j=0;j<n;++j) {real h=a[i][j]/p[j]; if(h>M) M=h;}
sss+=M;
}


for(int i=0;i<m;++i)
{
real M=0; 
for(int j=0;j<n;++j) {real h=a[i][j]/P[j]; if(h>M) M=h;}
s+=M;
}

if(sss<s) P=copy(p); 
write(s,s-Mdel); 
if(s-Mdel<1/10^15*s) {write("******it took "+string(seconds()-start)+" seconds******");
pause();}
}

confl=false; for(int j=0;j<n;++j) confl=confl || conf[j]; 
if(!confl) {write("no conflict"); reset();} else fix();

if(fast)
{
for(int i=0;i<m;++i) for(int j=0;j<n;++j)
{
if(conf[j] && b[i][j]>0 && bb[i][j]>0) 
{
real u=-db[i][j]/b[i][j]; 
//v=-dbb[i][j]/bb[i][j]; 
if(u>TTT/10 && aaa[i][j]>a[i][j]/2 && blacklist[i][j]<=ttt) temp[i][j]=1;
}
}
}

if(confl) temp[I][J]=1;

void loop()
{
bres=copy(b); Mdel=Norm(b); B=b[I][J]; if(B==0) B=1;

int cyclestart=seconds();

for(int i=0; i<m;++i) for(int j=0; j<n; ++j) if(temp[i][j]>0) aaa[i][j]=a[i][j]*kill;

for(int kk=0;kk<tt*SK;++kk) procedure(); 

if(b[I][J]>0 && confl) {write("Too weak killing!"); pause(); kill/=10;}

for(int i=0; i<m ;++i) for(int j=0; j<n; ++j) if(temp[i][j]>0) aaa[i][j]=a[i][j];

for(int kk=0;kk<tt*SK;++kk) procedure();

cycletime+=seconds()-cyclestart+1;

M=Norm(b); 
}

loop(); real rat=b[I][J]/B;

while (rat>0 && rat<0.9 && M>Mdel) {write("Repeating..."); loop(); rat=b[I][J]/B;}

if(confl && rat>0 && M>Mdel) {write("BLACKLISTING!"); blacklist[I][J]=2*ttt; ++blcount; if(blcount>0) reset((blcount>4?true:false));} 


int bl=0; for (int i=0;i<m;++i) 
bl+=sum(map(sgn,max(blacklist[i]-array(n,ttt),array(n,0)))); 
write(string(bl)+"  vertices blacklisted");


if(M>Mdel) 
{
if(rat==0) {fast=true; blcount=0;}
if(confl) write("Success!"+(b[I][J]==0?" Vertex is gone": "Vertex stays with ratio "+string(b[I][J]/B)+
" and abs value "+string(b[I][J])));
if(!check()) tt*=2; 
Mdel=M; successtime+=2*tt*SK; notreset=true;} 
else 
{
b=copy(bres); fast=false; failtime+=2*tt*SK;
blacklist[I][J]=2*tt;
if(confl) write("Failure! "+string(Mdel-M)+" short...");   
if (tt<ttt) tt*=2; else 
if (TTT>0 && confl) 
{
write("BLACKLISTING!"); blacklist[I][J]=2*ttt; ++blcount; if(blcount>0) reset((blcount>ulttol?true:false));
//write(tt,ttt); pause(); 
} 
else reset(); 
//else {tt*=2;}
}


}

Answer 4

I know that adding a second answer to the same thread is somewhat frowned upon but I felt like a couple of things here merit a special discussion. To avoid any issues with undeserved reputation points, etc., I'll make it a community wiki. Also I apologize in advance that I do not have a continuous stretch of time to type the full thing in one go, so I'll type it by parts, which will, probably, bump it to the front page more often than necessary.

Before I go into mathematics, let me just say that Rahul's answer is both excellent and terrible. It is excellent because it allows one to draw from readily existing sources and just avoid any hard thinking and it is terrible for the same very reason. The code he offers doesn't solve the problem. It merely restates it in the language understandable to the machine, after which the problem is delegated to a black box that spits out an uncheckable answer (even apparently bogus sometimes, as our discussion with Rahul shows, though I still believe that it might be an issue with human programming rather than with the solver itself) and you are left with no better understanding of the matters than you had in the first place. Of course, most of the available solvers are far superior to anything you or I can write ourselves when we have a whole bunch of complicated problems with some crazy constraints and objective functions and need one solver that works for all of them. However I'm really curious what is the price one has to pay for creating a Universal Monster instead of a small application that is aimed at a specific question (and what is the price one has to pay for delegating tasks to such a monster instead of trying to find one's own approach if not to all, then at least to some questions). That's why I wanted to see what is the best precision one can obtain using the standard software on a particular matrix for which I can find an exact solution using a few tricks.

So, the questions I'm going to address now are adding natural additional constraints and the speed of convergence. Note that I can easily add only lower bounds $w_j\ge W_j$ and $h_i\ge H_i$ but not the upper ones. You'll see why in a minute.

Adding the height restrictions is easy. The duality is ultimately just a statement that you need to consider all "trivial lower bounds" and switch from minimax to maximin (the devil is, of course, in the details starting with the exact definition of "trivial lower bounds"). The objective function now is $\sum_i\max(H_i,\max_j\frac {a_{ij}}{w_j})$ and we can use the same trick to estimate it from below by $\sum_{i}[c_iH_i+\sum_jb_{ij}\frac {a_{ij}}{w_j}]$ where now $c_i,b_{ij}\ge 0$ and $c_i+\sum_j b_{ij}=1$. If we had no width restrictions, the story would be almost exactly as before: we would just add terms with the same $j$, use the relation between $c$ and $b$ and get $$ \sum_i H_i+\sum_j \frac{S_j}{w_j}-\sum_i H_i\sum_j b_{ij} $$ with $S_j=\sum_i a_{ij}b_{ij}$ as before. The minimum over $w_j$ is again attained when they are proportional to $\sqrt{S_j}$, so the functional to maximize is $$ \left[\sum_j\sqrt{S_j}\right]^2-\sum_i H_i\sum_{j}b_{ij}=\sigma^2-\sum_i H_i\sum_{j}b_{ij} $$ We can consider one row and take the derivatives, as before, and see that two cases are possible: either we have $\sum_{j}b_{ij}<1$, in which case we have the equations $\frac \sigma{\sqrt{S_j}}=H_i$ for all $j$ with $b_{ij}>0$ and the corresponding inequalities for $b_{ij}=0$, or we have the inequalities everywhere but the constraint $\sum_j b_{ij}=1$ instead. Both cases result in a one-parametric family of vectors to consider and we just should check which constraint is stronger. Note also that we do not need to get an exact maximizer in the row at each step. It is enough to move in the direction of the maximizer and not to overshoot. So, in effect we can view $\sigma$ as a constant when recalculating $b_{ij}$ (the non-overshooting property requires a proof, of course). That's what I'm using in the code though, of course, it is still a story about finding the root of a monotone function of one variable. Since we'll not get a final answer in one step, we'd better not to waste two much time on finding that root with high precision.

The tricky part is to incorporate the width restrictions. Of course, I can formally write $\min_w$ with the restricted domain but then I'll not be able to compute it easily and all my nice formulae and the speech about admissible $b_{ij}$ forming a one-parameter family will go down the drain. So we need to be a bit inventive here. Note that we can add any sum $\sum_j\tau_j(\frac{W_j}{w_j}-1)$ with non-negative $\tau_j$ to our lower bound because this quantity is never positive for $w_j\ge W_j$. This will allow us to bring $\tau$'s and $b$'s together and to redefine $S_j$ as $\tau_jW_j+\sum_{i}a_{ij}b_{ij}$, so that we would have the expression $$ \left[\sum_j\sqrt{S_j}\right]^2-\sum_i H_i\sum_{j}b_{ij}-\sum_j\tau_j $$ to maximize. Again, it is quite a story about why the minimax is the same as the maximin here, but it is at least clear that any such expression can serve as a lower bound for the original problem. Note that $\tau$ enters it in exactly the same way as each row of $b$ and the only difference is that we do not have the restriction that the sum of $\tau_j$ is bounded by $1$ (in fact, those numbers can be as large as they wish), so updating $\tau$'s can be done in pretty much the same way as updating $b$'s.

There is one important catch in this new setup however. Note that we may have the situation when all $b$'s and $\tau$'s are $0$, in which case $w_j$ cannot be determined as "proportional to $\sqrt{S_j}$" because anything is proportional to a string of zeroes. It really happens if (and only if) the constant height restrictions are the strongest constraint, so all weight goes on them. In this case we have no real competition between $w_j$, just the restriction that they shouldn't force the height of any cell to be above the corresponding $H_i$, so we can just initially put $w_j=\max_i \frac{a_{ij}}{H_i}$. The sum will be automatically not greater than $1$ and we can then just scale it to $1$ by enlarging each $w_j$.

The code is below (again in Asymptote, and again not combed, but, apparently, working). Feel free to edit and rewrite it in C#, etc. if you are still interested in how it all works :-). The next question to discuss is the convergence rate. With this simple iteration scheme, it is not good at all (something like $1$ over the number of iterations). I was curious for a while whether one could invent something that would facilitate finding the "exact" (technically machine precision) solutions for reasonable size matrices and after experimenting with few ideas I found something that works at least up to size 50 by 50 on my old laptop though, to be honest, I do not quite understand why exactly it works (however, as before, it outputs both the answer and the certificate of optimality, so technically it doesn't matter how it finds them; the result is completely verifiable when it is achieved).

srand(seconds());

int m=50, n=75, K=201, cc=20;

real[] H,P;
for(int i=0;i<m;++i) H[i]=n*unitrand();
for(int j=0;j<n;++j) P[j]=unitrand();
P*=unitrand()/sum(P);

real[][] a, b;
for(int i=0;i<m;++i)
{
a[i]=new real[]; b[i]=new real[];
if(i<m) {for(int j=0; j<n; ++j) {a[i][j]=unitrand(); b[i][j]=1/n;}}
//a[i][rand()%n]=2;
a[i]*=unitrand();
}

real[] p,t,S;
for(int j=0;j<n;++j) t[j]=0;

for(int k=0;k<K;++k)
{
for(int j=0;j<n;++j)
{
real s=P[j]*t[j]; for(int i=0;i<m;++i) s+=a[i][j]*b[i][j]; S[j]=s;
}


for(int j=0;j<n;++j)
{
p[j]=P[j]; for(int i=0;i<m;++i) p[j]=max(p[j],a[i][j]/(H[i]+1/10^10));
}
if(sum(p)<1) p/=sum(p);
else {p=map(sqrt,S); p/=sum(p);}





if(k%cc==0)
{
write("*** Iteration "+string(k)+" ***");
real s=sum(map(sqrt,S))^2-sum(t)+sum(H);
for(int i=0;i<m;++i) s-=H[i]*sum(b[i]);
write(s);
}

for(int kk=0;kk<20;++kk)
{
p=max(p,P);
p/=sum(p);
}
real s=0; 
for(int i=0;i<m;++i)
{
real M=H[i]; 
for(int j=0;j<n;++j) {real h=a[i][j]/p[j]; if(h>M) M=h;}
s+=M;
}
if(k%cc==0) write(s);
//pause();

real SRS=sum(map(sqrt,S));
for(int kk=0;kk<5;++kk)
{
real[] V,T;
for(int j=0;j<n;++j) {V[j]=S[j]-t[j]*P[j]; t[j]=(P[j]>0?max(SRS^2*P[j]-V[j]/P[j],0):0); S[j]=V[j]+t[j]*P[j];}
SRS=sum(map(sqrt,S));
}

for(int i=0;i<m;++i)
{
real[] A,V,C,B;
for(int j=0;j<n;++j) {A[j]=a[i][j]; V[j]=S[j]-a[i][j]*b[i][j]; C[j]=V[j]/a[i][j];}
if(H[i]>0) 
{
for(int j=0;j<n;++j) B[j]=max(SRS^2/H[i]^2*A[j]-C[j],0);
}
if(H[i]==0 || sum(B)>1)
{
real aa=(sum(C)+1)/sum(A);
real da=1;
while(da>1/10^10)
{
for(int j=0;j<n;++j) {B[j]=aa*A[j]-C[j]; if(B[j]<0) {A[j]=0;B[j]=0;}}
da=sum(B)-1; aa-=da/sum(A); 
}
}
for(int j=0;j<n;++j) {b[i][j]=B[j]; S[j]=V[j]+a[i][j]*B[j];}
SRS=sum(map(sqrt,S));
}


}

write("************");

write(t,P,p);

pause();

Here is just the code that is (presumably) finding the exact answer for not too large matrices. All explanations will come later. Again it is in Asymptote. Why not in C? I know C++ but having a slow computer with an interpreted rather than compiled language allows you to watch the program as it goes (if you bother to output the work protocol to the screen, of course) and see many nuances that otherwise would be easily missed. I found at least 6 logical bugs in the original version this way (I hope the remaining number is less). The progress is not linear because the program tries to locate what Matt calls "king cells" in the matrix and the "vert=" line is a better mark of progress than the difference between the upper and the lower bounds. The final array of widths is $P$ and the program terminates when the relative error is $10^{-15}$ (be careful here: due to the rounding errors it may continue to run beyond that point, but once you see "no conflict" and "nothing to do" repeating again and again, the job is actually done and you are watching the geometric convergence game end (which can be also played differently, but who cares?).

srand(seconds());

int m=30, n=17, K=100001, k, SK=10, Count=0, proccnt=0, Failtime=1000000, I=0,J=0, cycletime=0; 
real M=0, Mdel=0, TTT=0, B;
int time=0, failtime=0, successtime=0; 
int tt=1, ttt=1, blcount=0, ulttol=3;


int start=seconds();

int[][] blacklist;
for(int i=0;i<m;++i) blacklist[i]=array(n,1);

real[][] a, aaa,  b , bb, bbb, db, dbb, bres;
real[] AA;

/// Initiating the matrix ////////////

real delta=0.0001;

for(int i=0;i<m;++i)
{
real u=unitrand(), v=unitrand();
a[i]=new real[]; b[i]=new real[];
for(int j=0; j<n; ++j) 
{
a[i][j]=1+delta*sin((i+j^2)); 
b[i][j]=1/n;
}
//a[rand()%(i+1)][rand()%n]=2;
//a[i]*=unitrand();
}

////////////////////////////////////////////


aaa=copy(a); bres=copy(b);
real kill=1/1000;



int [][] temp; bool[] conf=array(n,true); 
bool fast=true, notreset=true, confl=true;

for(int i=0;i<m;++i) temp[i]=array(n,0);

int[] Ind; for(int i=0;i<m;++i) Ind[i]=i;

real Norm(real[][] b)
{
real[] S;
for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=a[i][j]*b[i][j]; S[j]=s;
}
return sum(map(sqrt,S))^2;
}


void shuffle()
{
for(int kk=0;kk<m;++kk) {int a=rand()%m, b=rand()%m; int II=Ind[a]; Ind[a]=Ind[b]; Ind[b]=II;}
}

bool[] conflict(real[][] b)
{
bool[] conf=array(n,false);

int count=0; 

for(int i=0;i<m;++i) 
{
if(min(b[i])<0) {write("karaul"); pause();}
b[i]=max(b[i],array(n,0));
count+=sum(map(sgn,b[i]));
}
int[] pres=array(m,1);
int[][] sb;
for(int i=0;i<m;++i) {sb[i]=map(sgn,b[i]); sb[i][n]=1;}


for(int I=1;I<m;++I)
for(int i=0; i<I; ++i)
{
if(pres[i]>0 && sum(sb[i]*sb[I])>sb[i][n]*sb[I][n]) {pres[i]=0; sb[I]=sb[i]+sb[I];}
}

int vert,edgecnt,Vert=0,Edgecnt=0; int comp=sum(map(sgn,pres));
for(int i=0;i<m;++i) 
{
if(pres[i]>0) 
{
vert=sum(sb[i])-sb[i][n];
Vert+=vert;
edgecnt=-sb[i][n];
for(int j=0;j<n;++j) edgecnt+=max(2*sb[i][j]-1,0); 
Edgecnt+=edgecnt;
if(edgecnt>vert-1) for(int j=0;j<n;++j) {if(sb[i][j]>0) conf[j]=true;}
}
}
int alive=0; for(int i=0;i<m;++i) for(int j=0;j<n;++j)
if(conf[j] && b[i][j]>0 && blacklist[i][j]<=ttt) ++alive;
write("vert="+string(Vert)+"("+string(alive)+") edgecnt="+string(Edgecnt)+" comp="+ string(comp));
return conf;
}





real[] p, P, S;

for(k=0;k<K;++k)
{

void procedure()
{
for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=aaa[i][j]*b[i][j]; S[j]=s;
}
for(int i:Ind)
{
real aa;
real[] A,V,C,B;
for(int j=0;j<n;++j) {A[j]=aaa[i][j]; V[j]=S[j]-aaa[i][j]*b[i][j]; C[j]=V[j]/aaa[i][j];}
real aa=(k==0?(sum(C)+1)/sum(A):AA[i]);

int countbound=40;

for(int j=0;j<n;++j) B[j]=max(aa*A[j]-C[j],0);
if(sum(B)>1/2)
{
if(sum(B)<1)
{
real sl=0;
for(int j=0;j<n;++j) sl+=A[j]*sgn(B[j]);
aa+=1.0001*((1-sum(B))/sl); countbound=4;
}
}
else aa=(sum(C)+1)/sum(A);

real da=1;
int count=0;

while(da>0 && count<countbound)
{
++count; 
//write(da,count); //pause();
for(int j=0;j<n;++j) {B[j]=aa*A[j]-C[j]; if(B[j]<0) {B[j]=0; A[j]=0; C[j]=0;}}
if(sum(A)>0) {da=sum(B)-1; aa-=da/sum(A);}
else {write("alert"); pause(); for(int j=0;j<n;++j) {if(b[i][j]>0) A[j]=aaa[i][j];} aa=(sum(C)+1)/sum(A); } 
//write(sum(B),aa,da,sum(A),sum(C));
}
for(int j=0;j<n;++j) {b[i][j]=B[j]; S[j]=V[j]+aaa[i][j]*B[j];}
Count+=count; 

if(abs(sum(b[i])-1)>0.1) {write("rough!"); pause();}
AA[i]=aa; b[i]/=sum(b[i]);
}
++proccnt;
}

bool check()
{
bool check=false;
for(int i=0;i<m && !check;++i) for(int j=0;j<n;++j) check=check || (bres[i][j]>0 && b[i][j]==0);
return check;
}




void fix()
{
for(int i=0;i<m;++i) for(int j=0;j<n;++j) 
{
if(b[i][j]==0 && conf[j]) aaa[i][j]=a[i][j]*kill;
//if(b[i][j]==0) blacklist[i][j]=1;
}
}


void advance(bool adv=true)
{
for(int kk=0;kk<(adv?ttt:tt)*SK;++kk) procedure(); bres=copy(b); if(adv) {write("advancing with speed "+string(TTT)); fix();}
}


void reset(bool hard=true)
{
if(!confl) write("nothing to do"); else write("RESETTING "+(hard?"HARD":"SOFT")); 
fast=true; if(hard) blcount=0;   
//shuffle();
aaa=copy(a); for(int kk=0;kk<(confl && hard?ttt:1)*SK;++kk) procedure(); 
if(confl && hard) ttt*=2;  
fix(); 
}

real minb=1, shift=0;

TTT=1;

while (TTT>1/3) 
{ 
TTT=0;
//bbb=copy(b); advance(false); 
bb=copy(b); advance(false); bres=copy(b);

for(int i=0;i<m;++i) 
{
db[i]=b[i]-bb[i]; 
//dbb[i]=bb[i]-bbb[i]; 
shift=max(shift,max(map(abs,db[i]))); temp[i]=array(n,0);
}

for(int i=0;i<m;++i) for(int j=0;j<n;++j)
{
if(b[i][j]>0 && db[i][j]<0 && bb[i][j]>0) 
{
real u=-db[i][j]/b[i][j];
//v=-dbb[i][j]/bb[i][j]; 
if(u>TTT && u>0 && aaa[i][j]>a[i][j]/2 && blacklist[i][j]<=ttt) {TTT=u; I=i; J=j; minb=min(minb,b[i][j]);}
}
}
tt=(confl?blacklist[I][J]:1);
if(TTT>1/3) advance(); 
else if(TTT==0 || blcount>ulttol) reset();
else write('\n \naccelerating from speed '+string(TTT)+
"; position=("+string(I)+","+string(J)+"); cycle count= "+string(2*tt*SK)); 
}

time=seconds()-start; if(time>Failtime) {write('\n\nTotal failure'); pause(); Failtime*=2;} 

write("time= "+string(time)+", cycling "+string(cycletime)+
" seconds, failures =  "+string(failtime)+ ", successes= "+string(successtime));

write("count="+string(Count/m/proccnt)); 

conf=conflict(b);

for(int j=0;j<n;++j)
{
real s=0; for(int i=0;i<m;++i) s+=aaa[i][j]*b[i][j]; S[j]=s; p[j]=sqrt(s);  
}

p/=sum(p); 
if(k==0) P=copy(p); 
write(Mdel); 

{
real s=0, sss=0; 
for(int i=0;i<m;++i)
{
real M=0; 
for(int j=0;j<n;++j) {real h=a[i][j]/p[j]; if(h>M) M=h;}
sss+=M;
}


for(int i=0;i<m;++i)
{
real M=0; 
for(int j=0;j<n;++j) {real h=a[i][j]/P[j]; if(h>M) M=h;}
s+=M;
}

if(sss<s) P=copy(p); 
write(s,s-Mdel); 
if(s-Mdel<1/10^15*s) {write("******it took "+string(seconds()-start)+" seconds******");
pause();}
}

confl=false; for(int j=0;j<n;++j) confl=confl || conf[j]; 
if(!confl) {write("no conflict"); reset();} else fix();

if(fast)
{
for(int i=0;i<m;++i) for(int j=0;j<n;++j)
{
if(conf[j] && b[i][j]>0 && bb[i][j]>0) 
{
real u=-db[i][j]/b[i][j]; 
//v=-dbb[i][j]/bb[i][j]; 
if(u>TTT/10 && aaa[i][j]>a[i][j]/2 && blacklist[i][j]<=ttt) temp[i][j]=1;
}
}
}

if(confl) temp[I][J]=1;

void loop()
{
bres=copy(b); Mdel=Norm(b); B=b[I][J]; if(B==0) B=1;

int cyclestart=seconds();

for(int i=0; i<m;++i) for(int j=0; j<n; ++j) if(temp[i][j]>0) aaa[i][j]=a[i][j]*kill;

for(int kk=0;kk<tt*SK;++kk) procedure(); 

if(b[I][J]>0 && confl) {write("Too weak killing!"); pause(); kill/=10;}

for(int i=0; i<m ;++i) for(int j=0; j<n; ++j) if(temp[i][j]>0) aaa[i][j]=a[i][j];

for(int kk=0;kk<tt*SK;++kk) procedure();

cycletime+=seconds()-cyclestart+1;

M=Norm(b); 
}

loop(); real rat=b[I][J]/B;

while (rat>0 && rat<0.9 && M>Mdel) {write("Repeating..."); loop(); rat=b[I][J]/B;}

if(confl && rat>0 && M>Mdel) {write("BLACKLISTING!"); blacklist[I][J]=2*ttt; ++blcount; if(blcount>0) reset((blcount>4?true:false));} 


int bl=0; for (int i=0;i<m;++i) 
bl+=sum(map(sgn,max(blacklist[i]-array(n,ttt),array(n,0)))); 
write(string(bl)+"  vertices blacklisted");


if(M>Mdel) 
{
if(rat==0) {fast=true; blcount=0;}
if(confl) write("Success!"+(b[I][J]==0?" Vertex is gone": "Vertex stays with ratio "+string(b[I][J]/B)+
" and abs value "+string(b[I][J])));
if(!check()) tt*=2; 
Mdel=M; successtime+=2*tt*SK; notreset=true;} 
else 
{
b=copy(bres); fast=false; failtime+=2*tt*SK;
blacklist[I][J]=2*tt;
if(confl) write("Failure! "+string(Mdel-M)+" short...");   
if (tt<ttt) tt*=2; else 
if (TTT>0 && confl) 
{
write("BLACKLISTING!"); blacklist[I][J]=2*ttt; ++blcount; if(blcount>0) reset((blcount>ulttol?true:false));
//write(tt,ttt); pause(); 
} 
else reset(); 
//else {tt*=2;}
}


}