# Multinomial transformation for matrices

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:

$r_i=\log(p_i/p_n)$

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 http://en.wikipedia.org/wiki/Multinomial_logit>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)},$

$p_i=\exp(r_i)p_n$.

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.

-
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: mathoverflow.net/questions/156983/… (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: github.com/mpiktas/retacoro. 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

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.

-
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}$.

-
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