Disclaimer: I need guidance and help with where to start looking for solution of the problem I have described below. My background is in optimization and I am new to statistical methods, so there is a good chance that I am asking the wrong question or/and used wrong terminologies (please correct if that is the case).

Below I setup my problem:


  • the set of $n\times n$ matrices.

  • two functoins, $f: \mathbb{R}^{n^2} \rightarrow\mathbb{R}$ and ${\bf{g}}: \mathbb{R}^{n^2}\rightarrow Symm.(n\times n)$

  • Two constraints as follows:

$$1 - \epsilon_1<f({\bf{M}}) < 1 + \epsilon_1$$

$$-\epsilon_2< \lvert\lvert{{\bf{g}}}({\bf{M}})\rvert\rvert_{\infty} < \epsilon_2$$

for ${\bf{M}} \in \mathbb{R}^{n^2}$ and $\epsilon_1$ and $\epsilon_2$ both fixed small positive numbers.

Here are my questions:

  1. Can I find a model, or a distribution, which when I sample from it, it produces $n\times n$ matrices that satisfy the above two constraints (most of the time)? The sampled data needs to be close the real distribution in order to be representative.
  2. Is question (1) a well formulated question?
  3. If the answer to (2) is yes, what method(s) should I look into in order to work towards a solution?

For those who are interested in more concrete realizations of functions $f$ and $\bf{g}$, $f({\bf{M}})=\mathrm{det}({\bf{M}})$ and ${{\bf{g}}}({\bf{M}})=\frac{1}{2}({\bf{M}}^T{\bf{M}} - {\bf{I}})$.

I appreciate any hint or help with this problem. Tags mentioned below are speculative.

  • 2
    $\begingroup$ The answer for question (1) is yes : you can for example just generate random matrices and then filter out the ones that don't satisfy your constraints. This will clearly generate what you want, just not very efficiently. Are you asking for ways to achieve it more efficiently ? $\endgroup$
    – J. Delaney
    Feb 23, 2022 at 14:35
  • $\begingroup$ Your two functions generically define a region in $\mathbb{R}^{n^2}.$ Apart from that, you don't provide any information at all about what the distribution ought to be within that region. Developing probability models and generating random variables is a rich subject, with a very large number of different techniques (discussing just some univariate distributions the handbook requires two large, expensive volumes), Therefore, please tell us what distribution you need $M$ to have. $\endgroup$
    – whuber
    Feb 23, 2022 at 16:31
  • $\begingroup$ @J.Delaney Yes, I am looking for ways to generate these matrices more efficiently. Plus, having a model that satisfies these constraints, I can "query" the model for more matrices during an online training, and not rely on "random generation -> filtering" scheme. $\endgroup$
    – s.yadegari
    Feb 24, 2022 at 8:02
  • $\begingroup$ @whuber Thanks for your input. A prior distribution on $M$ could be of the form uniform $U[-\epsilon_2, \epsilon_2]$ for non-diagonal entries and $U[1-\epsilon_2, 1+\epsilon_2]$ for diagonal entries. $\endgroup$
    – s.yadegari
    Feb 24, 2022 at 8:19

1 Answer 1


Your question defines a set of matrices you want to sample from. However as pointed out in the comments, the notion of sampling depends on some distribution or probability measure over the set.

As a simple example, consider the $\delta$ distribution around some matrix $M_0$ that satisfies your constraints - this is a perfectly valid distribution and it answers your requirements because 'sampling' from it will always give you the matrix $M_0$, but this is probably not what you want.

Now having that in mind we can try to come up with of with some distribution that might be reasonable for your purpose (with the caveat that without knowing the purpose of this exercise it is hard to tell what is 'reasonable'): note that your specific set of constraints defines matrices that are close to unitary, since if $M$ is unitary ($M^TM=I$) then both $f(M)=1$ and $g(M)=0$. So, let

$$M = U(I + \epsilon X)$$

Where $U$ is unitary, $X$ is a random matrix (with some distribution) and $\epsilon$ is a small real number, therefore

$$ \det(M)=1 + \epsilon \text{tr}(X) + O(\epsilon^2) $$

$$M^TM - I = \epsilon(X+X^T) + O(\epsilon^2) $$

For a given $X$ this defines a possible range for $\epsilon$, since you require that

$$ -\epsilon_1 < \epsilon \text{tr}(X) < \epsilon_2 $$ $$ -\epsilon_1 < |\epsilon| \cdot \|\frac{X+X^T}{2}\|_{\infty} < \epsilon_2 $$

(For the second constraints only the upper bound is actually relevant since the norm can only be positive). So, generate a random matrix $X$ form some distribution, then generate a random $\epsilon$ from the range defined above (or alternatively you can use different distributions for the diagonal and off diagonal elements of $X$ in a way that will automatically satisfy the constraints): this will give you a matrix $M$ that satisfies the constraints with high probability (not surely because of the $O(\epsilon^2)$ terms).

A random unitary matrix $U$ can be generated in many ways, for example using Cayley transform: $U = (I - A)(I + A)^{-1}$ is unitary if $A$ is anti-symmetric. Therefore generate another random matrix $X$ and use its anti-symmetric part $X-X^T$ to get $U$.

You can generalize the above procedure to other functions $f$ and $g$ by expanding them in a similar way with a small parameter $\epsilon$ around $M_0$ that satisfies $f(M_0)-1=g(M_0)=0$ (because you assume $\epsilon_{1,2}$ are small). However if this condition defines a set of matrices that you want to sample from, you will have to find some parametrization that will allow you to do that (the equivalent of generating the unitary matrices in your case). There is probably no universal recipe for that so it will have to be worked out case by case.


Your Answer

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

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