<s>[INST] Create matrix containing values in [0,1] where sum of all diagonals and anti-diagonals is fixed The problem I am facing sounds at first glance pretty simple. However, as very often, it seems more complicated than I first assumed:
I want to calculate a matrix $P = (p_{j,k}) \in \mathbb{R}^{n \times n}$, $n\in\mathbb{N}$, satisfying the following constraints: 


*

*$p_{j,k} \in [0,1]$,

*$\sum_{j,k} \ p_{j,k} = 1$,

*the sums of all $2n-1$ diagonals are fixed by $b_{1-n},\ldots,b_{n-1}\in (0,1)$, e.g. $\sum_{j} p_{j,j} = b_0$ or $\sum_{j} p_{j,j+1} = b_{-1}$, and finally

*the sums of all $2n-1$ antidiagonals are fixed by $a_1,\ldots,a_{2n-1}\in (0,1)$.


Obviously, the problem is easy to solve for $n=2$, since the "corner elements" of the matrix are directly given from 3. and 4. However, for $n>2$, it becomes more difficult. 
A simple linear algebra approach for $n=3$ leads to the problem of "solving" $Mp=c$ with 
$$
M = \left(
\begin{array}{ccccccccc}
1 &  &  & 0  & 0  & 0  & 0  & 0  & 0 \\
& 1 & & 1 & & & 0 & 0 & 0 \\  
& & 1 & & 1 & & 1 & & \\
0 & 0 & 0 & & & 1 & & 1 & \\ 
0 & 0 & 0 & 0 & 0 & 0 & & & 1\\
0 & 0 & 0 & 0 & 0 & 0 & 1 & &\\
0 & 0 & 0 & 1 & & & & 1 & \\ 
1 & & & & 1 & & & & 1\\
& 1 & & & & 1 & 0 & 0 & 0 \\  
&  &  1 & 0  & 0  & 0  & 0  & 0  & 0 \\
1& 1 & 1& 1 & 1 & 1 & 1& 1&1 \\
\end{array}
\right), \quad p = 
\left( 
\begin{array}{c} 
p_{1,1} \\ 
p_{1,2} \\
p_{1,3} \\ 
p_{2,1} \\
p_{2,2} \\ 
p_{2,3} \\
p_{3,1} \\ 
p_{3,2} \\
p_{3,3}
\end{array} \right), \quad c = 
\left( 
\begin{array}{c} 
a_1 \\
a_2 \\
a_3 \\
a_4 \\
a_5 \\
b_{-2} \\
b_{-1} \\
b_{0} \\
b_{1} \\
b_{2} \\
1
\end{array} \right) 
$$
where $p \in \mathbb{R}^{n^2}$, $b \in \mathbb{R}^{4n-1}$ and $M$ is of size $(4n-1)\times n^2$ with entries in $\{0,1\}$. Note that $M$ is rank-deficient and the rows that specify the values for the "corner elements" are clearly visible. Dropping these rows leaves
$$
\tilde{M} = \left(
\begin{array}{ccccccccc}
& 1 & & 1 & & & 0 & 0 & 0 \\  
& & 1 & & 1 & & 1 & & \\
0 & 0 & 0 & & & 1 & & 1 & \\ 
0 & 0 & 0 & 1 & & & & 1 & \\ 
1 & & & & 1 & & & & 1\\
& 1 & & & & 1 & 0 & 0 & 0 \\  
1& 1 & 1&1 &1 & 1 & 1 & 1 & 1 \\  
\end{array}
\right), \quad \tilde{p} = 
\left( 
\begin{array}{c} 
p_{1,2} \\
p_{2,1} \\
p_{2,2} \\ 
p_{2,3} \\
p_{3,2} \\
\end{array} \right), \quad \tilde{c} = 
\left( 
\begin{array}{c} 
a_2 \\
a_3 \\
a_4 \\
b_{-1} \\
b_{0} \\
b_{1} \\
1-a_1-a_5-b_{-2}-b_2
\end{array} \right) .
$$
The question is, if there exists a solution (for $n \in \mathbb{N}$) for this problem satisfying all properties and how to calculate it (least-squares, SVD,...)? Or is there an approach different from the one I chose which is suitable to compute a solution iteratively? I have the intuition, that due to the somewhat special constraints, there is a close connection to probability ("values are positive, sum equals one") or graph (adjacency matrices) theory.
 [/INST] The given constraints are a system of linear inequalities, so you can find a feasible solution (or prove that none exists) by feeding these constraints to a linear program (LP) solver with some arbitrarily chosen objective function.
</s>