Random Sampling Many-Dimensional Sets Arising in Control †

: Various Monte Carlo techniques for random point generation over sets of interest are widely used in many areas of computational mathematics, optimization, data processing, etc. Whereas for regularly shaped sets such sampling is immediate to arrange, for nontrivial, implicitly speciﬁed domains these techniques are not easy to implement. We consider the so-called Hit-and-Run algorithm, a representative of the class of Markov chain Monte Carlo methods, which became popular in recent years. To perform random sampling over a set, this method requires only the knowledge of the intersection of a line through a point inside the set with the boundary of this set. This component of the Hit-and-Run procedure, known as boundary oracle, has to be performed quickly when applied to economy point representation of many-dimensional sets within the randomized approach to data mining, image reconstruction, control, optimization, etc. In this paper, we consider several vector and matrix sets typically encountered in control and speciﬁed by linear matrix inequalities. Closed-form solutions are proposed for ﬁnding the respective points of intersection, leading to efﬁcient boundary oracles; they are generalized to robust formulations where the system matrices contain norm-bounded uncertainty.


Introduction
One of the first issues in data mining and pattern recognition is an economy representation of implicitly specified massive data arrays with the subsequent extraction of specific features and classification [1][2][3]. Every element of the data set can be associated with a vector of which the components are the numerical values of certain properties of the phenomenon of interest under one or another operating condition. As a rule, a complete description of the data is not available due to high dimensions and cardinality, and the very absence of an explicit "generating" mechanism. Instead, a representative enough sample from the data set may be considered, leading to a point representation of the respective continuum domain, and the desired characteristics can then be measured just at these sample points.
There exist many ways to organize such samples. The most obvious one is to arrange a deterministic rectangular grid over the respective many-dimensional domain. However, this approach suffers several drawbacks. First, for high dimensions, the density of the grid needed for a reasonable representation of the set is very large; second, it can be efficiently implemented only for box-shaped sets; finally, sets that differ from rectangles may be embedded into larger boxes, however, the so-called rejection rate (the amount of "idle" grid points lying outside the domain of interest) may grow very quickly with the growth of dimension. A more efficient technique relates to the generation of quasi-random points in the feasible domain; they are also referred to as LP τ sequences introduced in [4]; however, sets [23] and, most efficiently, in the probabilistic approach to numerical optimization. For the latter, see, for example, [24,25], where the authors proposed new versions of the cutting plane scheme for use in efficient methods of optimization. Namely, the HR points were generated with the aim of approximating the center of gravity of the support set of the optimized function. Other applications include random walks in graph spaces [26], machine learning [16,27], etc. In particular, for a specific determinantal point process problem, a version of the HR sampler based on the zonotopic structure of the problem was proposed in [16]. This approach, though being faster than the conventional HR, requires solving additional linear programs and is not easy to implement.
However, in spite of a wide range of applications of the Hit-and-Run sampler, to the best of our knowledge, almost no research was devoted to its use in control. Perhaps the only paper in this direction is [28], where this method was applied to the design of stabilizing controllers in the simplest form.
In this paper we do not deal with applications of the HR to specific control problems; this is the subject for future research. Nor do we analyze the convergence and optimization of the performance of the HR algorithm. Instead, we limit ourselves to the considerations related to its BO component, and propose several numerical schemes for the implementation of this component.
In principle, instead of having a BO, the availability of a weaker membership oracle might be assumed, and a straightforward line search can then be used. However, this search has to be performed many times, so that it is important to provide a "closed-form" solution to this operation in order to improve the performance of the algorithm. By closedform solution we also mean various numerical procedures such as finding the roots of a polynomial or the eigenvalues of a matrix, solving semidefinite programs (SDPs) of moderate size and the like. These operations are indeed very fast and numerically efficient in the standard MATLAB implementation.
In this paper, we propose efficient computations of several BOs for various vector and matrix sets defined by linear matrix inequalities (LMIs).
The sets that we analyze in the paper appear in control problems in a natural way, and we consider sets specified by LMIs in the canonical form, the classical matrix Lyapunov inequality, and algebraic matrix Riccati inequalities. We also pay attention to various robust versions of these inequalities, where the matrix coefficients contain additive uncertainty.
We focus on the sets defined by LMIs because of the generality and flexibility of this technique and it is wide spread in systems theory; see the classical book [29] and the recent survey paper [30]. Overall, LMI-sets describe a wast variety of convex sets. Notably, because of convexity, only two intersection points are to be found.
One of the few works in this research direction is the conference paper [31], where first attempts were made to the implementation of such boundary oracles. Since then, possible applications of the HR algorithm became broader, and some of them motivated us to get back to this subject and consider it in a different context. For instance, one of the novel applications of this technique might be processing huge data arrays obtained from ultrasound computed tomography in medical diagnostics; for example, such as in [32]. This may lead to fast image reconstruction algorithms and improve the quality of the image.
In the current paper several major changes are performed as compared to [31]. The introduction section is completely re-written to better substantiate the importance of the research; for the same purpose, the order of the exposition and the overall wording is changed, the bibliography list is considerably extended to account for the new available developments, some proofs are added and some of them are corrected and made more transparent, several mathematical inaccuracies in the exposition are also corrected, many new notations are introduced to facilitate the understanding of the material, and typos are removed.
The notation used in the paper is very standard. Namely, R n and R k×m denote the space of real column vectors of length n and real k × m matrices; is the transposition sign; ⊕ and ⊗ stand for the Kronecker sum and Kronecker product of two matrices; the signs ≺, , , denote the negative (positive) definiteness (semidefiniteness) of a matrix; · is the Euclidean vector norm or the spectral or the Frobenius matrix norm (this will be clear from the context); W 1/2 is the matrix square root of the positive-definite matrix W.

The Hit-and-Run Procedure
We briefly formulate the main components of the HR procedure as applied to convex sets and uniform distributions.
Denote by x 0 ∈ R n an initial point in the interior of a closed convex bounded domain D ⊂ R n , and let x j denote the point generated at the jth step of the procedure. The next point is generated as follows. First, a random unit-length vector y ∈ R n is generated uniformly on the sphere in R n and the line x j + λy through this point in the direction y is considered. Denote by x j and x j the points of intersection of this line and the boundary of D; they are assumed to be obtained with the use of BO. Then the next point x j+1 is generated uniformly randomly on the segment [x j , x j ], a new random direction is generated, etc. The scheme of the HR algorithm is presented in Figure 1 for n = 2. In the HR-related literature it is shown that, under mild conditions on the set D and the choice of initial x 0 , the resulting sequence of random points tend to the uniform distribution over D as the number of points grow.
It is interesting to note that the HR method and particularly, its BO component can be used not only for generating points inside sets, but also for approximate point description of the boundary of implicitly specified sets. This description is obtained at no extra cost by simply memorizing the "side-product" endpoints x j and x j . We note that there exist specialized methods that generate points exactly on the boundary of the sets of interest [33].

Illustration: Static Output Feedback Stabilization
Let us consider a linear time invariant system in the state-space descriptioṅ Here the system matrices A ∈ R n×n , B ∈ R n×m , C ∈ R k×m are fixed and known. Assume that there exists a static output control u = Ky that makes the closed-loop matrix A + BKC Hurwitz stable, and let the set D ⊂ R k×m of all such stabilizing gain matrices K be bounded. A typical control problem is to optimize certain "engineering" performance indices over the set of these matrix gains, such as overshoot, oscillation, damping time.
First, it is well known that the very existence of a stabilizing output feedback is hard to establish [34], not talking about optimizing such a control. Second, the engineering performance indices mentioned above are not easy to formalize and optimize even in a much simpler situation where state feedback is considered.
As an alternative to the existing techniques we can instead arrange random sampling over the domain D, compute the value of the performance index for each sample, and pick the best one. This approach to optimal controller design suggests the availability of a representative enough sample set in D, and it can be obtained by using the HR-algorithm. In its simplest form this approach was considered in [35,36], where only low-order controllers were considered and the domain D was two-dimensional. No random sampling was used, but the authors suggested to pick a point in D with a mouse on the computer screen and compute the respective performance index. Such a "manual" semi-heuristic approach, being very simple, often leads to controllers that outperform those obtained by the methods available in the literature.
To arrange the HR algorithm in the general many-dimensional setting, assume that we possess a feasible K ∈ D, denote by L ∈ R k×m a matrix direction, take K + λL instead of K, and consider the λ-dependent n × n matrix A + B(K + λL)C, where λ ∈ R. Then the goal of the BO is to find the minimum λ and maximum λ values of λ ∈ R that guarantee the stability of A + B(K + λL)C over the whole segment [λ, λ].
To further simplify notation, denote A 0 : = A + BKC and A 1 : = BLC; then we have A + B(K + λL)C = A 0 + λA 1 , where A 0 is stable. Hence, the boundary oracle reduces to finding the so-called unidirectional stability segment of maximum length [37].
For the clarity of the overall exposition, we present the corresponding result.

Lemma 1 ([37]
). Let A 0 , A 1 ∈ R n×n , and let A 0 be Hurwitz stable. The minimal and the maximal values of the parameter λ ∈ R that guarantee the stability of A 0 + λA 1 are defined by where λ i are real eigenvalues of the matrix −(A 0 ⊕ A 0 ) −1 (A 1 ⊕ A 1 ), and A ⊕ A denotes the Kronecker sum.
We recall that the Kronecker sum of the two square matrices A and B is defined as where I a , I b are identity matrices of appropriate dimensions and ⊗ stands for the Kronecker product of matrices.
Back to our problem, we see that the matrix segment [K + λL; K + λL] obtained with the use of the lemma above belongs to the set D of output stabilizing gain matrices. Therefore, the next matrix (i.e., the next HR point) is then generated randomly on this segment.
Unfortunately, in Lemma 1 one has to compute eigenvalues of n 2 × n 2 matrices, which makes computations slow for high dimensions. For instance, already for n = 25, the computations would take about a second on a standard laptop, which is considered rather slow, keeping in mind multiple executions of this procedure as required in the HR-algorithm. On the other hand, (i) typically, problems arising in control are of lower dimensions, and (ii) the availability of boundary oracle can be relaxed to the assumption of having a membership oracle, so that a boundary oracle can always be constructed by means of a linear search. However, for a wide range of problems, a boundary oracle can be easily formulated in closed form, and below we present particular boundary oracles for matrix sets specified by linear matrix inequalities.

LMI Sets in Canonical Form
We turn to formation of boundary oracles related to sets defined by linear matrix inequalities and consider LMIs in the canonical form.

The Uncertainty-Free Setup
Introduce the following matrix-valued function of a vector variable: Consider the convex domain defined as where the sign denotes the negative semidefiniteness of a matrix, and we assume that D as nonempty. Let x, y ∈ R n be, respectively, a point in the interior of D and a direction. To formulate a boundary oracle for this domain, we need to compute the minimal and the maximal values λ, λ of the scalar parameter λ ∈ R that guarantee the sign-definiteness of F(x + λy) over the interval [λ, λ]. We have where we denoted A 0 : = F(x) and A 1 : = ∑ n i=1 y i F i . It is seen that A 0 ≺ 0 and A 1 = A 1 ; hence, similarly to the above, the goal is to compute the minimum and the maximum values of λ ∈ R that guarantee the signdefiniteness of the affine matrix function A 0 + λA 1 . The lemma below is a slight modification of the result in [7] it shows that this can be done by computing the generalized eigenvalues of the matrix pair (A 0 , −A 1 ). Lemma 2. Let A 0 ≺ 0 and A 1 = A 1 , then the critical values of the scalar λ that guarantee negative semidefiniteness of the matrix function A 0 + λA 1 are obtained as: where λ i are the generalized eigenvalues of the matrix pair (A 0 , −A 1 ); i.e., A 0 e i = −λ i A 1 e i .
Proof. The proof is nearly trivial: The critical value of λ that makes the matrix A 0 + λA 1 sign indefinite is also responsible for the loss of nonsingularity; by definition, this means that there exists a nonzero e ∈ R m such that (A 0 + λA 1 )e = 0.
This result of the lemma above applies to symmetric matrices; i.e., it is a special case of Lemma 1, since a symmetric matrix is stable if and only if it is negative definite.
As far as the initial point x 0 ∈ intD is considered, it can be found as a solution {x,η} of the following SDP: min η subject to F(x) η I.

A Robust Formulation
We now consider the situation where the matrix coefficients F i in (4) contain uncertainties. Specifically, they have the form F i + ∆ i , where ∆ i ∈ R m×m are bounded in the spectral norm ∆ i ≤ ε i , i = 0, . . . , n, and mutually independent, and ε i ≥ 0 are given levels of uncertainty. To retain the symmetric structure of the matrices (F i + ∆ i ), the ∆ i 's are also assumed to be symmetric. We call such uncertainties admissible.
As a result, we obtain the uncertain linear matrix-valued function of the vector variable x ∈ R n and the respective (convex) robustly feasible domain D rob = {x ∈ R n : F(x, ∆) 0 for all admissible ∆ i }.
Likewise the uncertainty-free setup, having an x ∈ intD rob and a directional vector y ∈ R n , we aim to compute the minimal λ rob and maximal λ rob values of the parameter λ that guarantee that the uncertain LMI F(x + λy, ∆) 0 holds for all admissible uncertainties ∆ i over the segment [λ rob , λ rob ]. This procedure will be referred to as robust boundary oracle.
The construction of the feasible and the associated robustly feasible domains is depicted in Figure 2 for x ∈ R 2 , the same symmetric matrices F 0 , F 1 , F 2 ∈ R 3×3 that were used in the example of Figure 1, and certain numerical values of ε i .
Here, λ and λ denote the critical values of the parameter λ obtained via formulae (6) and (7) for the corresponding uncertainty-free formulation in Lemma 2.
From the lemma below, Reference [7] presents the robust BO for the uncertain LMI F(x, ∆) 0; it amounts to solving certain nonlinear equations in one scalar variable. Lemma 3. Let x ∈ intD rob and y ∈ R n be given. Consider the two functions in the variable λ ∈ R: The critical values λ rob and λ rob of the parameter λ, which guarantee that the matrix F(x + λy, ∆) is negative definite for all admissible ∆, are found as the two solutions of the equation φ(λ) = ε(λ) on the segment [λ, λ] obtained in (6) and (7). Clearly, the robust BO is more laborious than its uncertainty-free analog; indeed, not only a specific nonlinear equation is to be solved numerically, but the nonrobust bounds λ and λ are to be computed as well.
The proof of Lemma 3 uses the same basic reasonings: The loss of robust signdefiniteness of the matrix function F(x + λy, ∆) coincides with the loss of nonsingularity for some admissible ∆ i s. The proof heavily exploits the independence of the ∆ i s and uses the symmetric analog of Theorem 3 in [38] on the radius of matrix nonsingularity.
In Figure 3, the plots of the functions φ(λ) and ε(λ) are depicted for an uncertain LMI with x ∈ R 4 , and the desired segments for λ are shown.

Quadratic Stabilization
From this section onward we switch our attention from "abstract" LMIs in vector variables to those encountered in control problems; namely, to specific LMIs in matrix variables.

The Lyapynov Boundary Oracle
First of all we consider the Lyapunov inequality where G = G and A are given and X = X is the n × n matrix variable.
This inequality naturally appears in the design of stabilizing state feedback controllers for continuous-time systems. Specifically, consider the systemẋ = Ax + Bu (1) with controllable pair (A, B); then the gain matrix K for the stabilizing control u = Kx may be found as K = −B X −1 , where X 0 is a solution of the LMI (9) with G = −2BB , and V(x) = x X −1 x is a quadratic Lyapunov function for the closed-loop system [29].
Therefore, the convex closed set of matrices specified by the LMI (9) is of obvious interest in control design, and we formulate its boundary oracle. Such an oracle is seen to be computed with the use of Lemma 2. Given an X ∈ intD and a matrix direction Y = Y , we obtain A(X + λY) + (X + λY)A + G 0.
Denoting A 0 : = AX + XA + G and A 1 : = AY + YA and keeping in mind that A 0 ≺ 0 and A 1 = A 1 , we find ourselves in the conditions of Lemma 2, which gives us the interval Λ A of sign-definiteness of the λ-dependent function A 0 + λA 1 . Since the matrix X + λY is to be positive definite (as the matrix of a Lyapunov function), we again use Lemma 2 to calculate the corresponding range Λ X of λ retaining the sign-definiteness of X + λY. Then adopt Λ = Λ A ∩ Λ X as the desired maximal interval of λ for which the matrix X + λY belongs to D (10).

A Robust Formulation
The Lyapunov BO of the previous subsection can be generalized to the situation where the matrix A in (9) contains uncertainty. We consider the model of structured uncertainty, which often appears in control applications: Here the matrix perturbation ∆ ∈ R p×q is only known to be bounded in some norm ∆ ≤ 1, and the matrices M ∈ R n×p , N ∈ R q×n are given.
In such a setup, the uncertain Lyapunov inequality is used for the construction of a common quadratic Lyapunov function with matrix X for the uncertain system and yields robustly stabilizing controllers. This gives rise to the robust Lyapunov set D rob = X 0 : inequality (12) holds for all ∆ ≤ 1 , (13) and the goal is to devise a boundary oracles for it. The main technical trick for operating with uncertainty (11) is the so-called Petersen lemma proposed in [39]; it provides necessary and sufficient conditions for an uncertain matrix to be sign-definite.

Lemma 4.
Assume that R = R ∈ R n×n , P ∈ R n×p , Q ∈ R q×n . Then, for the inequality R + P∆Q + (P∆Q) 0 (14) to hold for all ∆ ∈ R p×q , ∆ ≤ 1, it is necessary and sufficient that there exist ε > 0 such that We stress that, with the lemma above, checking the robust sign-definiteness of a continuum matrix family reduces to a convex problem in one scalar variable. To see this, using the Schur complement, represent inequality (15) in the equivalent form which is an LMI in the scalar variable ε. Moreover, with this machinery, the maximal admissible span for the uncertainty ∆ can be found: γ max = sup{γ : inequality (14) holds for all ∆ ≤ γ}. (17) It is sometimes called the radius of robust sign-definiteness of the uncertain matrix R + P∆Q + (P∆Q) ; computing this quantity reduces to the solution of a semidefinite program in two scalar variables and a "scaled" LMI constraint of the form (16), see [40].
With the Petersen lemma in mind, we are ready to formulate the boundary oracle for the set (13); it will be referred to as the robust Lyapunov BO.

Theorem 1.
Assume that X ∈ D rob (11)-(13) and Y = Y ∈ R n×n are given. Denote Λ A = [λ, λ], where λ and λ are found as the solutions of the two semidefinite programs: in the scalar variables λ, ε. The maximum interval of λ ∈ R that guarantees X + λY ∈ D rob is given by Λ = Λ A ∩ Λ X , where Λ X ⊂ R denotes the maximal interval of sign-definiteness of the pencil X + λY.
Proof. Having a feasible X ∈ D rob defined by (11)- (13), and a matrix direction Y, we are aimed at finding the minimum and the maximum values of λ for which the inequality is satisfied for all ∆ ≤ 1.
By Petersen lemma, this inequality holds for all admissible ∆ ≤ 1 and for some λ if and only if there exist ε > 0 such that Since the functions R(λ) and Q(λ) are affine in λ, the last inequality is seen to be an LMI in the two scalar variables λ and ε. Therefore, finding the critical values of λ reduces to the minimization and maximization of the variable λ subject to the LMI constraint above. Finally, the maximum interval of the sign-definiteness of the matrix X + λY is found using Lemma 2.
We make two comments at this point. First, an initial X ∈ intD rob can be obtained from the solution {X,η,ε} of the following SDP in the matrix variable X and two scalar variables η, ε: In complete analogy with (8), ifη < 0, adopt X =X; otherwise the set D rob (13) is empty.
This observation leads to the second comment. If the LMIs (19) are infeasible, we decrease the level of uncertainty and re-solve the SDP to test the nonemptiness of D rob . On the other hand, by solving a similar semidefinite program, the maximal allowed level γ rob = max{γ : ∃X 0 : inequality (12) holds for all ∆ ≤ γ}, the radius of robust quadratic stabilizability can be calculated in advance; cf. (17) and see [30,40] for the details.

Optimal Control
Assume now that, on top of stabilizing system (1) by means of linear control u = Kx, we want to optimize a certain performance index.
In one of the classical optimal control problems, the linear quadratic regulation (LQR) problem, the quadratic cost functional has the form where W x 0 and W u 0 are weight matrices. The minimization of this cost leads to the algebraic Riccati equation in the matrix variable X; under certain conditions on A, B its unique positive-definite solution defines the optimal gain K = −W −1 u B X −1 .

The Riccati Boundary Oracle
In a slightly more general formulation consider the set where G 0 is a given real n × n matrix. This set accumulates all possible solutions X 0 of the algebraic Riccati inequality, and the boundary oracle for this set is formulated in much the same way as in Section 4.1, with Lemma 2 as the main tool. Indeed, let us first take the Schur complement and represent the quadratic inequality in the LMI form: Next, considering X + λY instead of X, we arrive at Hence we are again in the conditions of Lemma 2, which gives us the associated segment Λ A .
Next, accounting for the sign definiteness of the matrix X + λY and using Lemma 2 once again, we obtain the corresponding feasible segment Λ X and adopt Λ A ∩ Λ X as the final answer.
Overall, the computation of the boundary oracle for the set defined by the quadratic inequality reduces to the Lyapunov BO described in Section 4.1; the only essential difference is that now the size of the constraint matrices is twice as large as before.

A Robust Formulation
Let now the matrix A be subjected to the uncertainty of the same structured form (11). The corresponding robust counterpart of the LQR problem can be formulated and solved in several ways; for example, see [41]. However, following the logic of the current paper we consider the set defined by the uncertain Riccati inequality and propose a boundary oracle for it. From the exposition of the two previous subsections it is seen that the robust version of the Riccati boundary oracle can be formulated by exploiting the Schur complement (to transform quadratic inequalities into linear ones) and Petersen lemma (to handle the uncertainty). We omit the obvious and not quite insightful derivations and formulate the result in the theorem below. Theorem 2. Let X ∈ intD rob (21), (11) and let Y = Y ∈ R n×n . The maximal and the minimal values of λ, which guarantee that the matrix segment X + λY belongs to the set D rob (21) are obtained from the solutions of the two semidefinite programs min / max λ subject to and the optimization is performed in the two scalar variables λ and ε.

Implementation Issues
We briefly discuss some technical implementation issues and report on the processor time required for the computation of the boundary oracles proposed above.
First, to generate a vector uniformly distributed on the surface of the unit ball in R n , we use the MATLAB routine randn to obtain a random vector y with the standard Gaussian distribution and then normalize it to obtain y = y/ y , which is distributed uniformly on the unit sphere. As far as the generation of matrix directions Y in Sections 4 and 5 is considered, we generate them randomly uniformly on the surface of the unit ball in the Frobenius norm (i.e., on the surface of the ball in R n 2 ) with the subsequent symmetrization.
The main tool for computing the Lyapunov and Riccati BOs, Lemma 2 requires finding generalized eigenvalues of a pair of matrices, and it can be efficiently applied to matrices of rather high dimensions.
The robust boundary oracles in Sections 4.2 and 5.2 require solving semidefinite programs in two scalar variables with LMI constraints having dimension 2n × 2n (Theorem 1) or 4n × 4n (Theorem 2). To solve these problems, we used the MATLAB-compatible software cvx [42].
The following illustrative experiment was performed. We generated randomly N = 1000 samples of the data required in Lemma 2 and Theorems 1 and 2. For each of these samples, we measured the execution time for the four boundary oracles proposed in the paper, and averaged over N samples. The results obtained on a standard laptop are presented in Table 1. From the first two rows of the table we see that nonrobust oracles are indeed very fast, since they exploit just the eigenvalue computation, and this operation is "perfectly" implemented in various available toolboxes. Instead, robust boundary oracles are based on solving a specific optimization problem, which is more time-consuming. The last row of the table clearly shows that up to n = 30, the robust Riccati oracle is fast enough, whereas for n = 50 the required execution time is approximately one and a half second, which may be considered slow. Note however, that in the latter case, the matrices in the LMI constraints are of quite large size 200 × 200.
In the experiment, we considered "full" matrices and used the standard freeware cvx for solving semidefinite programs. To speed up the computations in higher dimensions, one may exploit a sparse structure of the system matrices, which is typically the case for control problems, or try to use alternative optimization methods and software. On the other hand, practical control problems have lower dimensions, usually not exceeding n = 15; see [43], a popular collection of control-related test problems. In other words, in practice, all the proposed BOs seem to be fast enough to be used in sampling the respective sets.
Another observation is that the execution time grows nonlinearly with increase of problem dimension n. This is because the coefficient matrices in the respective LMIs contain n 2 elements and because of a specific representation of the LMI constraints in the internal language of cvx solvers. It is also seen that, for a given n, the execution time of the robust Riccati oracle is not exactly the same as that of the robust Lyapunov oracle for the doubled dimension (recall the relative sizes of the respective LMI constraints). This is because the structure of these constraints in the robust Riccati oracle is slightly more complicated.

Concluding Remarks and Directions for Future Research
In this paper, we focused on boundary oracles, one of the key components of Hitand-Run, a popular approach to random generation of points inside bounded sets. We proposed simple computational schemes that implement boundary oracles for several typical control-related matrix sets defined by linear matrix inequalities. At the nucleus of these oracles is the computation of the generalized eigenvalues of a pair of matrices or solving low-dimensional semidefinite programs in scalar variables. These operations are fast as implemented in the MATLAB environment, which makes the proposed BOs efficient for use in various random walk methods.
Future research will be targeted at the implementation of boundary oracles for broader classes of sets (e.g., the sets of Schur stable polynomials, positive polynomials, uncertain affine matrix families), which account for the presence of exogenous unknown-but-bounded disturbances, and speeding up the numerical implementation of the HR procedure in these specific problems for higher dimensions. We also note that the results in Sections 4 and 5 can be easily extended to discrete-time systems by using the discrete-time Lyapunov and Riccati inequalities.
Of the most interest is the application of the machinery presented here to problems of practical interest, in particular, to the design of stabilizing controllers subject to several engineering specifications as mentioned at the beginning of Section 2.2.