MathOverflow is a question and answer site for professional mathematicians. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Suppose we have a vector of probabilities $\mathbf{p}=(p_1,...,p_n)$, where $p_i>0$ for $i=1,...n$ and $\sum p_i=1$. Define new vector $\mathbf{r}=(r_1,...,r_{n-1})$ in a following way:


This defines the transformation $T:(0,1)^n\to\mathbb{R}^{n-1}$, $\mathbf{r}=T\mathbf{p}$. This transformation can be called multinomial transformation (or to be more precise inverse multinomial transformation), since similar formula is used in>multinomial logit model.

This transformation is useful for modelling, since resulting $r_i$ can be any real number, and there is an easy way to transform $r_i$ back to probabilities:

$p_n=\dfrac{1}{1+\sum \exp(r_i)},$


My question is whether there exists a similar transformation for matrices. Suppose we have two probability vectors $\mathbf{p}=(p_1,...,p_n)$, $\mathbf{q}=(q_1,...q_m)$ and $n\times m$ matrix $P=(p_{ij})$, satisfying

$\sum p_i=1$, $\sum q_i=1$

$\sum_{j=1}^m p_{ij}=p_i$, for each $i=1,...,n$, (1)

$\sum_{i=1}^n p_{ij}=q_j$, for each $j=1,...,m$, (2)

(what we actualy have is a bivariate discrete probability distribution with given marginal distributions).

Now what I am looking for is a transformation which transforms $p_{ij}$ to unbounded real numbers, but such that the inverse would satisfy constraints (1) and (2). In effect I am looking for the bijection from subset of $(0,1)^{nm}$ to $R^{k}$, where $k$ should be $(n-1)(m-1)$.

I suspect that maybe copulas can be involved here, or some properties of stochastic matrices. If somebody could give me any pointers I would be very grateful.

share|cite|improve this question
Constraints 1 and 2 are not enough to describe your bijection. Marginals can be represented as two vectors of log-odds, n+m-2 parameters total, then you can get the joint satisfying the constraints by multiplying the marginals – Yaroslav Bulatov Aug 20 '10 at 16:36
For an alternate solution, see:… (and take the log of the entries of $T$.) – Bill Bradley Feb 10 '14 at 21:25
Judging from the link it seems that the transformation $r_{ij}=\log(p_{ij})+\log(p_{i-1,j-1})-\log(p_{i-1,j})-\log(p_{j,i-1})$ is the desired transformation. Intuition says that this is a correct formula, I'll check how to recover $p_{ij}$, given $r_{ij}$ and the marginal distributions. – mpiktas Feb 11 '14 at 8:39
@BillBradley Your solution works. Set $r_{ij}=\log(p_{i1})+\log(p_{1j})-\log(p_{11})-\log(p_{ij})$ and given row and column sums it is possible to recover all the matrix. I did not managed to get a closed form solution, but it is possible to solve for the answer numerically. I've made this into R package: If you write this as an answer, I will accept it. Thanks for your help! – mpiktas Nov 23 '15 at 12:14
@mpiktas I'm happy to write up a description of the linked method, but just to check, my answer below is (also) correct, right? I'm happy to describe the linked method because it seems aesthetically better; is that why you prefer it? – Bill Bradley Nov 25 '15 at 22:48
up vote 2 down vote accepted

At the suggestion of the original poster, I am summarizing an alternate answer that has a few strengths relative to my original answer. It is related to the question and answer at this MathOverflow Question.

For any matrix $A$, define $$ L_A(i,j)=\log\left(\frac{A_{i,j}A_{i+1,j+1}}{A_{i+1,j}A_{i,j+1}} \right)$$

Let $S$ be the set of $n\times m$ matrices with fixed row and columns sums $p_i$ and $q_j$ and positive entries. Let $Q:S\rightarrow R^{(n-1)\times (m-1)}$ be the map where the $(i,j)$-th entry of $Q(A)$ is $L_A(i,j)$.

Then $Q$ is a bijection between $S$ and $R^{(n-1)\times (m-1)}$. Moreover, the inverse is efficiently computable. Mpiktas wrote an R package called retacoro that appears to recover the inverse through gradient descent. Alternately, the inversion of this projection can be expressed as a geometric program, and geometric programs can be solved in polynomial time (and efficiently in practice).

Establishing the bijectivity of $Q$ and defining the geometric program precisely seem a bit long for a MathOverflow post. I'll post a more detailed description to ArXiv and put a link here when I'm done.

share|cite|improve this answer
I look forward to see how the inversion can be expressed as a geometric problem. – mpiktas Nov 26 '15 at 16:20

(Edited to fix a bug.)

I think the following bijection will do what you want.

For $1\leq i,j\leq n-1$, define $$r_{ij}=\log(p_{ij}/p_i)$$

Given the $r_{ij}$ and the marginals, we can recover the $p_{ij}$ as follows:

$p_{ij}=p_i \exp(r_{ij})$ for $1\leq i,j \leq n-1$.

$p_{in}=p_i - \sum_{j=1}^{n-1}p_{ij}$ for $1 \leq i \leq n-1$.

$p_{nj}=q_j - \sum_{i=1}^{n-1}p_{ij}$ for $1 \leq j \leq n$.

Note that we do not use the fact that $\sum p_i=\sum q_j=1$, only that the individual $p_i$ and $q_j$ are known and non-zero. For example, we could apply this transformation to doubly stochastic matrices with non-zero entries. (In that case, $p_i=q_j=1$ for all $i$ and $j$.)

Also, if we replaced the matrix by an order $k$ tensor with the same type of constraints, an analogous argument gives a bijection to $R^{(n-1)^k}$.

share|cite|improve this answer
Shouldn't $d$ depend on $j$ given the transformation? Also what about constraints $p_i$ and $q_j$? Note that if I do not care about the constraints then I can simply convert bivariate discrete distribution to univariate and use the defined multinomial transformation. – mpiktas Feb 3 '14 at 18:36
Ah, I somehow misread your question and answered it for doubly stochastic matrices. Sorry-- I'm not quite sure how I did that. I'll see if I can fix the post. – Bill Bradley Feb 3 '14 at 22:21
I fixed the post above to address the correct question. I tried to remove the asymmetry between the rows and columns, but didn't have any luck. – Bill Bradley Feb 4 '14 at 15:26
If we perturb initial matrix (by adding random matrix where elements have small variance for example) and use newly calculated $r_{ij}$ but $p_i$ and $q_j$ from the initial matrix, we can get negative entries for the last column and rows. So it seems that $r_{ij}$ have certain restrictions, which is not the case for the univariate case. – mpiktas Nov 26 '15 at 8:37
@mpiktas You mean if you use (slightly) incorrect values for the $p_i$ and $q_j$, the recovery fails? That doesn't seem too surprising itself. Are you concerned that it points to a numerical stability issue? – Bill Bradley Nov 26 '15 at 15:35

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.