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.

 A: 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


*

*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}$;

*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).

