# Sampling from an (almost) multivariate normal over matrices

Consider $n$ points in the euclidean plane, $p_i = (x_i,y_i)_{1\leq i \leq n}$. Now consider a $2 \times 2$ matrix $M = \left(\begin{array}{cc}a & b\\c& d\end{array}\right)$ a vector $r = \left(\begin{array}{c}e\\f\end{array}\right)$, and the function $\Phi$ representing the c.d.f of a standard normal.

Let $$f:v \rightarrow \Phi(Mv+r)$$ be a mapping between the plane and the unit square. $\Phi$ is taken component by component. We are interested in the product of the jacobian of $f$ across all points $p$.

The higher this product, the more spread out the mapping of the points over the unit square. Intuitively, the points will be most spread out when $r$ and $M$ respectively demean and perform a PCA decomposition of the set of points $p$. I think I may have proved it at some point but I forgot, it definitely sounds truthy.

Beyond the optimal value of $M$ and $r$, one can interpret the product of the jacobians as a probability density, and we would like to sample values of $M$ and $r$ with respect to this density.

Let $$A = \sum_{i=1}^n \left(\begin{array}{ccc}x_i^2& x_i y_i & x_i\\ x_i y_i & y_i^2 & y_i\\ x_i & y_i & 1\end{array}\right)$$

I think that the log product of the jacobians can be written as

$$\sum_{i=1}^n \log |J(f)(p_i)| = C(n) -\frac{1}{2}(a~b~e~c~d~f)\left(\begin{array}{cc}A & 0\\0 & A\end{array}\right)(a~b~e~c~d~f)^{\top}+n \log |a d - b c|$$

where $C(n)$ is a constant. The first term represents the contribution of $\Phi$, and it tries to make the coefficients small to stay in the region where $\Phi$ is most expansive. The second represents the contribution of the multiplication by $M$, and it tries to make the determinant large to spread out the points.

Here's where I'm stuck...

As it is, I have the product of a multivariate gaussian density in $(a~b~e~c~d~f)$ (which is great and easy), and a function of the determinant of $M$, (which isn't).

1) Is there any trick to sample from this density, or to get rid of this determinant?

2) Am I even correct that this is maximized with PCA?

• Since M is 2*2, you could just analytically find its eigenvalues (it will be some simple quadratic equation) and then use the fact that the determinant is the product of the eigenvalues – Sid May 13 '15 at 18:51
• how would that help me exactly ? – Arthur B. May 13 '15 at 19:14

$f(x) = e^{-\frac{n}{2} x^2} x^n$
For $n>10$ or so, $f$ is very closely approximated by a gaussian distribution with mean $$\sqrt{\frac{n}{2}} \frac{\Gamma\left(\frac{n}{2}\right)}{\Gamma \left(\frac{n+1}{2}\right)} ~= 1+\frac{1}{4n} + O\left(\frac{1}{n^2}\right)$$ and variance $$1+\frac{1}{n}-\frac{\Gamma\left(1+\frac{n}{2}\right)\Gamma\left(\frac{n}{2}\right)}{\Gamma \left(\frac{n+1}{2}\right)^2}~= \frac{1}{2n} + O\left(\frac{1}{n^2}\right)$$