3
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Sep 20, 2016 at 18:18
  • 1
    $\begingroup$ 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}$.) $\endgroup$
    – whuber ♦
    Commented Sep 20, 2016 at 18:34
  • $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Sep 20, 2016 at 18:42
  • 1
    $\begingroup$ Sorry, I mixed up the parameters in my latest comment (but not in the preceding ones). I hope you get the idea. $\endgroup$
    – whuber ♦
    Commented Sep 20, 2016 at 18:53
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Sep 20, 2016 at 19:03

1 Answer 1

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ (+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. $\endgroup$
    – whuber ♦
    Commented Sep 23, 2016 at 15:27
  • $\begingroup$ @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. $\endgroup$
    – Xi'an
    Commented Sep 23, 2016 at 16:34
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Sep 23, 2016 at 16:40
  • 1
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Sep 26, 2016 at 23:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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