# Randomly constructing a bounded ellipsoid

In a project, I am working with constraints of the following type

$$\frac{1}{2}{x}^\top Q x + q^\top x + q_0 \leq 0$$

where I randomly generate the data by (randn stands for random standard normal sampling)

q = randn(n)
Q = randn(n,n); Q = Q'*Q;
q0 = - 1.0 * n + randn();


However, most of the time

$$\left\{x \in \mathbb{R}^n \ : \ \frac{1}{2}{x}^\top Q x + q^\top x + q_0 \leq 0 \right\}$$

is not bounded. Hence, my optimization problems under this constraint are producing errors. Is there a way of randomly generating such data while ensuring that this feasible set stays bounded?

• Do you want the ellipsoids to be centered at the origin? Commented Jun 15, 2023 at 8:04
• @RodrigodeAzevedo would be good to have that, but I can also parametrize that center and randomize it. Commented Jun 15, 2023 at 14:01
• I'm not sure it is possible (meaning probability > 0) for your feasible region to be unbounded. How are you detecting unboundedness? (If via a solver, what solver?)
– prubin
Commented Jun 15, 2023 at 21:48
• $Q$ being all zeros has probability 0. In fact, $Q$ having rank $< n$ has probability 0, which is key to my skepticism at it being unbounded. Also, I've tried three examples ($n=5$, different random seeds) and all were bounded.
– prubin
Commented Jun 16, 2023 at 2:49
• I think I have a proof that the ellipsoid is bounded with probability 1. At the same time, I did some experiments with two optimization models (using CPLEX). Given random problem instances ($n=5$), the first model sometimes said they were bounded but frequently said they were unbounded, while the second always said they were bounded. So if you are getting unboundedness from a solver, perhaps it is a solver issue?
– prubin
Commented Jun 16, 2023 at 21:35

As noted in a comment to the original question, I have (I believe) a proof that the ellipsoid is bounded with probability 1 (assuming the pseudorandom number generator is, well, random). I've also tested a few different optimization models, attempting to maximize either the $$L_1$$ or $$L_\infty$$ norm of $$x$$ over the ellipsoid, with mixed results. I published the proof (maybe a bit too long to include here) in one blog post and some computational results in another post.

• oh wow, such a great post, many thanks! I will investigate this and get back to you :) Commented Jun 19, 2023 at 15:45

# Mathematical boundedness

The model is bounded if all eigenvalues of $$Q$$ are positive. As @prubin already answered this is the case with high probability.

Mathematically, the distribution of eigenvalues for this matrix has been analyzed already (the Wishart distribution). This distribution is continuous, so the probability of an eigenvalue being exactly $$0$$ is nil.

# Solver troubles

The solver may have trouble with it because the generated matrix has really bad conditioning. The ratio between larger and smaller eigenvalues is typically around $$10^5$$ for a $$100 \times 100$$ matrix, and $$10^6$$ for a $$1000 \times 1000$$ matrix, but is often as high as $$10^8$$ and $$10^{10}$$.

In Python:

H = np.random.randn(100, 100)
Q = H @ H.transpose()
eig = np.sort(np.linalg.eigvalsh(Q))
print(f"{eig[-1] / eig[0]}")


This is called the condition number. A high condition number means that the smallest eigenvalue is very close to $$0$$ compared to the largest. Depending on the solver, this can cause numerical issues, and this is probably why the model is declared unbounded.

A plot of the condition number across matrices:

# Avoiding the issue

If this is indeed the issue for your solver, you could generate matrices with better eigenvalue distributions. A very simple way to improve the conditioning (flatten the eigenvalue distribution) is to make your initial matrix $$H$$ non-square ($$n \times m$$ with $$m > n$$).

I plotted the condition number for varying $$m$$, averaged over 100 matrices. We see that it improves the condition number very quickly (note the log scale!).

• Interesting. If H is n by m with m < n, $Q = HH^T$ would be at most rank m, which is less than full rank.. With m > n, Q is "super full rank", which is a term I just made up, and denotes Q being full rank, but with extra margin of being so. Commented Jun 19, 2023 at 18:04
• Indeed, it seems that there is still a log going on after the $n = m$ threshold. I was surprised by how steep the change is Commented Jun 19, 2023 at 18:25
• I ran a handful of simulations in R, with $H$ $100\times 100$, and the ratio of maximum (absolute) eigenvalue to minimum (absolute) eigenvalue came out between 60,000 and 90,000. That's not great, but I would not expect a good quality solver to have major problems with it. Clearly, though, something is throwing off the solver.
– prubin
Commented Jun 19, 2023 at 21:42
• I just saw that your tests are for $n=5$, and it seems very early to have problems due to condition numbers too, so I plotted the conditioning of many random 5x5 matrices. It can still get bad quickly if you are unlucky, but for most of them I wouldn't expect CPLEX to fail either Commented Jun 20, 2023 at 5:49
• Beautiful analysis @Ggouvine! Commented Jun 20, 2023 at 11:07

MOSEK offers a whitepaper on quadratic constraints in which it is concluded that if $$Q$$ has a suitable factor model, then it is often beneficial to represent your model in a way that exploits this factor model. In your case you have $$Q=H^T H$$ and although $$H$$ has no particular properties of interest (besides full rank with high probability), I suggest you try to reformulate $$\frac{1}{2} x^T Q x + q^T x + q_0 \leq 0$$

as $$\frac{1}{2} y^T y + q^T x + q_0 \leq 0,\quad y = Hx,$$

or perhaps even better as $$\frac{1}{2} y^T y + p^T y + q_0 \leq 0,\quad y = Hx,$$

where $$p$$ is the solution to $$H^T p = q$$, such that $$p^T y = p^T H x = q^T x$$ by definition. In your case, the potential benefit comes from having the quadratic coefficient matrix (the only nonlinear part) becoming an identity matrix which is numerically stable. There is no eigenvalue distribution to worry about.

• That's a good thing to try. I'm not convinced it will work in practice, as any numerical issues would be due to the interaction of the constraint with branch-and-cut, and it remains ill-conditioned (large variations in $x$ for small variations of constraint value) Commented Jun 21, 2023 at 10:19
• Definitely an overstatement, @Ggouvine, as the question doesn't even mention MIP, but let's go with it. If you solve nonlinear MIP by outer approximation, and the OA is unbounded, you try to generate cuts to reject this conclusion. If these cuts are rejected as numerical garbage, unboundedness is accepted. This is likely in x-space, but not in y-space where all OA-cuts are reasonable. Hence you obtain a relaxation bounded in y-space, where the conclusion of unboundedness is much harder to reach. Commented Jun 23, 2023 at 7:33
• Maybe @prubin has time to test this change as a follow-up to his blog posts? Commented Jun 23, 2023 at 7:34