# Simulation under Wishart-like constraint in $\mathbb{R}^{k\times p}$

Given a $(p,p)$ symmetric positive semi-definite matrix $\mathbf{H}$ of rank $k\le p$, I am looking for a (possibly efficient) way of generating a set of $k$ vectors $\alpha_i\in\mathbb{R}^p$ uniformly distributed under the constraint $$\sum_{i=1}^k \alpha_i \alpha_i^\text{T} = \mathbf{H}$$

If needed I can further assume that there exists a set of $k$ vectors $\beta_i\in\mathbb{R}^p$ such that $\mathbf{H}$ is constructed as $$\sum_{i=1}^k \beta_i \beta_i^\text{T} = \mathbf{H}$$

Note that, despite the title, this is unrelated to simulating a Wishart in that $\mathbf{H}$ is fixed and the $\alpha_i$ are not $\text{N}_p(0,\mathbf{I}_p)$ variates.

• Ultimately, this is a question of generating points from a uniform distribution on an ellipsoid in $\mathbb{R}^k$. What can you tell us about $p$ and the relative sizes of the eigenvalues of $H$? That information would be useful guidance concerning which approaches are worth considering. – whuber Sep 20 '16 at 18:18
• I believe it is an ellipsoid because it is readily parameterized as such. Adapting a basis of $\mathbb{R}^p$ to $H$ we may as well assume it is of full rank. Letting $$H=UD^2U^\prime$$ be the SVD, the matrices $A=(\alpha_1,\ldots,\alpha_k)$ are parameterized by the orthogonal matrices $Q$; the parameterization is $$Q\to UDQ=A.$$ This works because $$(UDQ)(UDQ)^\prime=UD^2U^\prime=H$$ and all such decompositions arise this way. Thus the set of points is a linear transformation of the set of orthogonal matrices $O(k)$: that's the "ellipsoid." (It lives in $\mathbb{R}^{k^2}$.) – whuber Sep 20 '16 at 18:34
• Choose an orthogonal basis for $\mathbb{R}^k$ in which the first $p$ elements span the image of $H$. All the action occurs within that $p$-dimensional subspace. The SVD will give you such a basis. – whuber Sep 20 '16 at 18:42
• Sorry, I mixed up the parameters in my latest comment (but not in the preceding ones). I hope you get the idea. – whuber Sep 20 '16 at 18:53
• I'm sorry, you lost me: your link goes right back to this thread and I don't see any related questions in your user profile. I'm not saying you're wrong: I have read so many questions and answered so many others that sometimes I completely forget what I have written. – whuber Sep 20 '16 at 19:03

To sum up W. Huber's comments, if I may, a way to simulate $k$ vectors $a_i$'s under the constraint $$\sum_{i=1}^k \alpha_i \alpha_i^\text{T} = \mathbf{H}$$ is to

1. obtain a singular value decomposition (SVD) of $\mathbf{H}$ with eigenvectors represented as the $p\times k$ matrix $\mathbf{U}$ (each column being one such eigenvector) and eigenvalues $\lambda_1,\ldots,\lambda_k$, stored in the $k\times k$ diagonal matrix $\mathbf{D}$;
2. generate an orthogonal matrix $\mathbf{Q}$ in $\mathcal{O}(k)$ and compute the $p\times k$ matrix $\mathbf{A}=\mathbf{U}\mathbf{D}^{1/2}\mathbf{Q}$, with each column of $\mathbf{A}$ being an $\alpha_i$.

Note that, to generate an orthogonal matrix, the following applies (quoting from Wikipedia):

To generate an (n + 1) × (n + 1) orthogonal matrix, take an n × n one and a uniformly distributed unit vector of dimension n + 1. Construct a Householder reflection from the vector, then apply it to the smaller matrix (embedded in the larger size with a 1 at the bottom right corner).

• (+1) (1) Note that this does not produce vectors that are "uniformly" distributed unless all the eigenvalues of $H$ are equal. (That's why I left my reply in a comment). (2) A quick and dirty way to construct orthogonal matrices is to generate an $m\times n$ array of iid standard Normal variates with $m \ge n$ (preferably a little larger to make sure the array is not rank deficient) and diagonalize its covariance matrix. The diagonalization yields an orthogonal transformation of the basis represented as an orthogonal $n\times n$ matrix. – whuber Sep 23 '16 at 15:27
• @whuber: (1) Thank you for pointing this out: I am interested in one parameterisation with a manageable prior, so picking a uniform prior on $\mathcal{O}(k)$ is good enough. (2) I was thinking of simulating a matrix of $k^2$ iid standard normal and of brute-force orthonormalising the columns. I may be missing the uniform distribution though. – Xi'an Sep 23 '16 at 16:34
• That approach in (2) should work--it ought to be the equivalent of my quick-and-dirty method. There's a tiny chance the columns will appear to be collinear (up to floating point error). That's why I recommend using one or a few extra rows: the collinearity becomes that much less likely. The proof of uniformity is inductive (on the number of columns). The base case reduces to the fact that the standard Normal distribution in $k$ dimensions is spherically symmetric. The induction step relies on the independence of the columns and the retention of symmetry upon projection to a subspace. – whuber Sep 23 '16 at 16:40
• See stats.stackexchange.com/questions/2746 concerning generating random orthogonal matrices. Many of the techniques described there for generating random correlation matrices can be adapted. – whuber Sep 26 '16 at 23:34